Scilabで熱拡散方程式 その3 (陰解法)

Scilabで熱拡散方程式 その1 (陽解法)Scilabで熱拡散方程式 その2 (無次元化)では陽解法で偏微分方程式を解く方法を書くと同時に、時間と空間の刻み幅を自由に選べないという問題点を指摘しました。今回はその問題点を解決した陰解法(Crank-Nicolson法)を紹介します。


001_20140215031531793.png
Fig.1: 一次元熱伝導問題の解



偏微分方程式の陰解法


Scilabで熱拡散方程式 その1 (陽解法)では、Octaveの精義を参考にして、一次元の熱伝導を題材に偏微分方程式を陽解法で解くということを行いました。しかしながらScilabで熱拡散方程式 その2 (無次元化)で紹介したとおり、陽解法では λ ≦ 1/2 という安定性の条件から空間の分解能と時間の分解能を独立に決めることが出来ないという欠点があることが確認されました。

Octaveの精義では、この問題を陰解法のひとつであるCrank-Nicolson法を用いて解決する方法を紹介しています。そこで今回は、その1の計算をCrank-Nicolson法で計算します。

計算


計算結果が冒頭に示したFig.1です。
そのスクリプトはCN_sce.txtです。

clear;

// *** 時間と空間の分割 ***
n = 129; // 空間の分割数
dx = 1 / (n - 1); // 空間の刻み幅

r = 4

dt = r * dx ^ 2; // 時間の刻み幅
m = 1001; // 時間の分割数
dm = 100; // プロットする m の間隔

// *** 初期条件と境界条件 ***
// 初期条件
u = zeros(n,m);
// 境界条件
u(1,1) = 0;
u(n,1) = 1;

// *** 行列C ***
C = sparse(toeplitz([-2, 1, zeros(1,n - 2)]));
A = sparse(2 * eye(n,n) - r * C);
B = sparse(2 * eye(n,n) + r * C);
b = zeros(n,1);
b($) = 2 * r;

// *** 偏微分方程式の計算 ***
for j = 1:(m - 1)
// 境界以外の計算
u(:, j + 1) = A \ (B * u(:,j) + b);
// 境界条件
u(1, j + 1) = 0;
u(n, j + 1) = 1;
end

// *** グラフのプロット ***
plot([0:dx:1],u(:,1:dm:m));


陰解法の利点


Scilabで熱拡散方程式 その1 (陽解法)Scilabで熱拡散方程式 その2 (無次元化)では、陽解法の問題点として時間と空間の刻み幅を独立に決定することが出来ないという点を挙げました。

Crank-Nicolson法はこの問題点を改善したものです。


002_201402150315301ea.png
Fig.2: 1次元熱伝導方程式の陽解法(FTCS法)による数値解で、Octaveの精義のコードをScilabへ移植したもの。全てh2=1/322に固定しr = 1/6 (左上), 1/2 (右上), 65/128 (左下)と変化させた結果。r = 64/128 > 1/2 の場合には解が不安定になっていることが読み取れる。

003_20140215031530c4c.png
Fig.3: 1次元熱伝導方程式の陰解法(Crank-Nicolson法)による数値解で、Octaveの精義のコードをScilabへ移植したもの。全てh2=1/1282に固定しr=1 (左上),2 (右上), 4(左下)と変化させた結果。全てのrに対して安定な解が得られていることが分かる。



関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器






フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 偏微分方程式 熱伝導 陰解法 Crank-Nicolson法 

Scilabで熱拡散方程式 その2 (無次元化)

Scilabで熱拡散方程式 その1 (陽解法)ではOctaveの精義を参考にして以下に示す一次元の熱拡散方程式を無次元化してScilabで計算しました。

\frac{\partial T'}{\partial t'} = D^2 \frac{\partial^2 T'}{\partial {x'}^2}

ただし 0 ≦ x' ≦ λ, t' ≧ 0

(T': 温度, t': 時間, D2: 熱拡散率, x': 位置)

今回はもう少し具体的に、現実に即した問題を無次元化し、更にもう一度具体的な次元を持った値に直すという手順をやってみます。


問題設定


Scilabで熱拡散方程式 その1 (陽解法)では一次元の熱伝導の問題を解くために熱拡散方程式を無次元化して陽解法でで計算を行いました。
無次元化を行ったということは、計算結果は実際の長さや時間とは違う目盛りで書かれているということなので、無次元化とは逆の処理をして実際のものに直す必要があります。

今回は以下のような具体的な条件でScilabで熱拡散方程式 その1 (陽解法)と全く同じ計算を行います。

長さ1.5mの銅でできた棒の左端を0℃の氷に、右端を100℃のお湯に接しさせたとき、1時間後の温度分布はどのようになるか。ただし、最初の温度分布は全て0℃であったとする。

銅の熱拡散率D2 = 116 × 10-6 m2/sは以下の計算式とパラメータから計算しました。

D^{2} = \frac{k}{\rho C_p}

(D2: 熱拡散率, ρ: 密度, CP: 定圧比熱)

パラメータ
熱伝導度 k398 W/m/K
密度 ρ8.92×103 kg/m3
定圧モル比熱 CP,m24.44 J/mol/K
原子量 M6.3546×10-2 kg/mol
table.1: 熱拡散率D2=116×10-6 m2/s を計算するためのパラメータ。定圧比熱は密度と単位をそろえるためにCP=CP,m/Mとする。各パラメータの出展はWikipedia。


分子のおもちゃ箱 熱拡散長には色々な金属の熱拡散率が表にまとめられています。

熱拡散方程式


Scilabで熱拡散方程式 その1 (陽解法)と全く同じ内容をもう一度書いておきます。今回解くべき偏微分方程式は以下のものになります。

\frac{\partial T'}{\partial t'} = D^2 \frac{\partial^2 T'}{\partial {x'}^2}

ただし 0 ≦ x' ≦ λ, t' ≧ 0

(T': 温度, t': 時間, D2:熱拡散率, x': 位置)

なお文字の上についている'(アポストロフィー)は(後で行う)無次元化をしていないことを意味します(微分を意味するものではありません)。

方程式の無次元化


以下の様に方程式の無次元化することを考えます。

x' = λx
t' = τt
T'(x',t') = T0u(x,t)

ただし

\frac{D^2 \tau}{\lambda^2}=1

の関係に注意が必要です。

いま銅の棒の長さは1.5mなので λ=1.5 m です。
これに熱拡散率 D2 = 116 × 10-6 m2/s をあわせてτが必然的に決まります。

\tau = \frac{\lambda^2}{D^2}

計算結果は τ = 19394.442 s となりました。

最後に温度の無次元化ですが、これは単純に T0 = 100 ℃ と置けばよいでしょう。このように置けば 0 ≦ T' ≦ 100 に対して 0 ≦ u ≦ 1 となります。

結局、無次元化した熱拡散方程式は下記の様になります。

\frac{\partial u(x,t)}{\partial t} = \frac{\partial^2 u(x,t)}{\partial x^2}

0 ≦ x ≦ 1, t ≧ 0

方程式の差分化


無次元化された熱拡散方程式を差分化するに際して、時間と空間の刻み幅をまとめて

r = \frac{\Delta t}{(\Delta x)^2}

と置いたわけですが、これが r < 1/2 では計算が失敗し r = 1/6 のときだけ誤差が小さくなるという特徴があると前回書きました。

ここからは r = 1/6 を選ぶとして話を進めます。rにこれ以外の数字を選ぶにしても、必ず1/2 未満にならないようにチェックをする必要があります。

いま一次元空間 0 ≦ x ≦ 1 の範囲でn点ほど計算するとするならば、Δxの大きさは Δx = 1/(n-1) となります。
rとΔxが決まるとΔtも必然的に決まり Δt = r×(Δx)2 です。

時間はτを使って無次元化されているため、無次元化された時間の刻み幅Δtもτを使って実際の時間の刻み幅に戻すことが出来ます。

\tau \Delta t = \tau r (\Delta x)^2 = \frac{\tau r}{(n-1)^2}

仮に空間の計算点数を n=50 とすると、τΔt= 1.3462753 s となります。
今回の計算では1時間後(60×60=3600秒後)の温度分布を計算するので m ≒ 3600 / (τΔt) ≒ 2674 回の繰り返し計算をします。

計算結果


以下に示すのが 約1時間後(3599.9401秒後)の銅の棒の温度分布です。


001_20140123132748b05.png
Fig.1: 全体が0℃であった銅の棒の左端を0℃に固定、右端を100℃に固定した状態で約1時間(3599.9401秒)保持したときの内部の温度分布。


計算の手順はこのエントリに書いたとおりです。

clear;

// *** 銅の物性値 ***
k = 398; // 熱伝導度 398 W/m/K
rho = 8.92E3; // 密度 8.92×10^3 kg/m^3

cpm = 24.44; // 定積モル比熱 24.44 J/mol/K
m = 63.546; // 原子量 g/mol
mkg = m * 1E-3; // 原子量 kg/mol
cp = cpm / mkg; // 定積質量比熱 J/kg/K

d2 = k / (rho * cp) // 熱拡散率 m^2/s

// *** 無次元化前の時間と空間 ***
lambda = 1.5; // 棒の長さ λ (m)
tdmax = 60 * 60; // t'max (s) まで計算

// *** 無次元化 ***
tau = lambda ^ 2 / d2; // 時間の無次元化係数
tmax = tdmax / tau; // 無次元化時間で tmax まで計算
temp0 = 100; // u=1 に対してT'=100℃

// *** 時間と空間の分割 ***
n = 50; // 空間の分割数
dx = 1 / (n - 1); // 空間の刻み幅

r = 1 / 6;

dt = r * dx ^ 2; // 時間の刻み幅
m = round(tmax / dt); // 時間の分割数

// *** 初期条件と境界条件 ***
// 初期条件
U = zeros(n,m);
// 境界条件
U(1,1) = 0;
U(n,1) = 1;

// *** 行列P ***
s = 1 - 2 * r;
P = toeplitz([s, r, zeros(1, n - 2)]);

// *** 偏微分方程式の計算 ***
for j = 1:(m - 1)
// 境界以外の計算
U(:, j + 1) = P * U(:,j);
// 境界条件
U(1, j + 1) = 0;
U(n, j + 1) = 1;
end

// *** グラフのプロット ***
X = linspace(0,lambda,n);
plot(X,temp0 * U(:,m),'o-b');
zoom_rect([0,0,lambda,temp0]);

xlabel("Position (m)");
ylabel("Temperature (degree C)");


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器







フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 偏微分方程式 熱拡散方程式 熱伝導 

Scilabで熱拡散方程式 その1 (陽解法)

Scilabには常微分方程式ソルバや非線形方程式ソルバは存在しますが、偏微分方程式ソルバというものは存在しません。
そこで今回はOctaveの精義をを参考にして、偏微分方程式の最も簡単な例の一つである、一次元の熱伝導の問題を計算しました。

\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}


偏微分方程式


Scilabを利用すると色々な数値計算を行うことが出来ます。
この中で常微分方程式非線形方程式を解くのは比較的簡単でScilabがデフォルトで持っている命令をワンパターンで使うだけです。

これに対して偏微分方程式をワンパターンで解けるような命令はScilabには用意されていません。今回はOctaveの精義を参考に一次元の熱伝導を陽解法(FTCS法:Forward Time Space methods)で解いてみます。

熱拡散方程式


物質の拡散や熱伝導は熱拡散方程式と呼ばれる偏微分方程式で表されます。
今回扱う一次元の物体内部の温度分布は以下のような形になります。

\frac{\partial T'}{\partial t'} = D^2 \frac{\partial^2 T'}{\partial {x'}^2}

ただし 0 ≦ x' ≦ λ, t' ≧ 0

(T': 温度, t': 時間, D2:熱拡散率, x': 位置)

なお文字の上についている'(アポストロフィー)は(後で行う)無次元化をしていないことを意味します(微分を意味するものではありません)。

方程式の無次元化


次にこれを数値計算が行いやすいように、以下の様に無次元化します。

x' = λx
t' = τt

するとxの範囲が 0 ≦ x ≦ 1 となります。

またτは

\frac{D^2 \tau}{\lambda^2}=1

すなわち

\tau = \frac{\lambda^2}{D^2}

です。

温度も計算結果の値が大きすぎたり小さすぎたりしないように適切なT0を選んで

T'(x',t') = T0u(x,t)

の様に無次元化します。

すると偏微分方程式は以下のような形になります。

\frac{\partial u(x,t)}{\partial t} = \frac{\partial^2 u(x,t)}{\partial x^2}

0 ≦ x ≦ 1, t ≧ 0

無次元化によって

T' → u
x' → x
t' → t

と3つの変数が変更されました。
以降ではこの無次元化された偏微分方程式について数値計算を行います。

方程式の差分化


更に前述の無次元化された偏微分方程式の微分を差分に変換します。

\frac{\partial u(x,t)}{\partial t} = \frac{\partial^2 u(x,t)}{\partial x^2}

を差分化すると以下の様になります。

\frac{\partial u(x,t)}{\partial t} \sim \frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}

\frac{\partial^2 u(x,t)}{\partial x^2} \sim \frac{u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x, t)}{(\Delta x)^2}

式が長くなるので、以下の様に ui,j を定義します。

u_{i,j} \equiv u((i-1)\Delta x, (j-1)\Delta t)

結局

\frac{u_{i,j+1}-u_{i,j}}{\Delta t} = \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}

ここで

r = \frac{\Delta t}{(\Delta x)^2}

とおくと

u_{i,j+1} = (1-2r)u_{i,j} + r(u_{i+1,j}+u_{i-1,j})

これには式の見た目が簡単になるということ以上に重要な意味があります。
rが小さいほど計算が高速になるのですが r < 1/2 では計算に失敗し正しい結果が得られないことが知られています。また r = 1/6 のときだけ打切り誤差が特別に小さくなるという性質があることも知られています。従って実際の計算では r = 1/2 または r = 1/6 のどちらかが採用されます。

この事は裏を返せば空間分解能を上げたいときは、同時に時間分解能を上げなければならない(逆も然り)ということでもあります。

Scilabでガウス型波束の散乱でハミルトニアンを行列で表したのと同様に

7861cb402820e6f58fddc8d3abb3e335_90_black.png

となる様な行列Pを考えると

f99827543db336a6148d03e0d830d8ab_90_black.png

ただし

r = \frac{\Delta t}{(\Delta x)^2}

s = 1 - 2r

です。

初期条件と境界条件


常微分方程式を解くために初期条件が必要であったのと同様に、偏微分方程式を解くためには初期条件と境界条件が必要になります。
問題設定によって色々な境界条件が考えられ、その境界条件に応じて異なったプログラムを書かなければならない事が偏微分方程式の難しさの元になっているのだと思うのですが、今回は入門書でよく見る以下のような簡単なものを考えます。

初期条件(t=0)
u(x,0) = f(x)

境界条件(x=0, x=1)
u(0,t) = 0
u(1,t) = 1

この初期条件と境界条件は、次のような状態を表していると考えることが出来ます。

まず、金属の棒全体を0℃の状態にしておく(初期条件)。
次に左端を0℃の氷に接しさせ、右端を100℃のお湯に接しさせる(境界条件)。

これを時間発展させれば、金属の棒の内部の温度分布がどうなるかを調べられます。
非常に長い時間を置いた後の最終的な温度分布は、直線的になるであろうことは、計算しなくても直感的に予想が付きます。

計算結果


以下に計算結果を示します。
各線は 200×Δt の時間ごとに結果をプロットしたもので一番最後のものが 1600×Δt の時間のものです。
この段階では温度分布がほとんど直線的になっており、直感的な予想が正しかったことが分かります。


001_20140119171908a1e.png
Fig.1: 計算結果


clear;

// *** 時間と空間の分割 ***
n = 33; // 空間の分割数
dx = 1 / (n - 1); // 空間の刻み幅

r = 1 / 6;

dt = r * dx ^ 2; // 時間の刻み幅
m = 1601; // 時間の分割数
dm = 200; // プロットする m の間隔

// *** 初期条件と境界条件 ***
// 初期条件
U = zeros(n,m);
// 境界条件
U(1,1) = 0;
U(n,1) = 1;

// *** 行列P ***
s = 1 - 2 * r;
P = toeplitz([s, r, zeros(1, n - 2)]);

// *** 偏微分方程式の計算 ***
for j = 1:(m - 1)
// 境界以外の計算
U(:, j + 1) = P * U(:,j);
// 境界条件
U(1, j + 1) = 0;
U(n, j + 1) = 1;
end

// *** グラフのプロット ***
plot([0:dx:1],U(:,1:dm:m));


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器






フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 偏微分方程式 熱伝導 

Ubuntu 12.04 64bitでHiLAPW

CMDワークショップでいただいて帰ったHiLAPW64bit版のUbuntu 12.04にインストールしてみました。

HiLAPWは計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)の中でも紹介されている第一原理計算パッケージのひとつで、その名前の通りLinearized Augmented Planewave(LAPW)法を使ったコードです。

私はCMDワークショップのアドバンストコースでHiLAPWの実習を受講させていただいた折に、バージョン1.13のコードをいただきました。HiLAPWのコードが使いたい方はCMDワークショップに参加するか、あるいはメールか何かの手段で小口多美夫先生に相談されるのが良いかと思います。

HiLAPWのインストール


さて、こういったコードのインストールにはmakefile(と場合によってはソースコード)の編集が必要になりますが、makefileの設定を見てみるに64bit版LinuxでIntel fortranコンパイラ(ifort)を使うように書かれていたので、Ubuntu 12.04の64bit版と非商用Intel Fortran Compilerの組み合わせで挑んだところ、何の編集もせずに動作させることが出来ました。手順は以下の通りです。

  1. Ubuntu 12.04 64bit版のダウンロード
  2. Ubuntu 12.04 64bitのインストール
  3. ifortのインストール
  4. LAPACKとBLASのインストール
  5. HiLAPWのインストール


Ubuntuは主にVNCのことを考えてちょっと古いバージョン12.04としました。最新版以外のダウンロード先は探しにくいかもしれませんが、12.04はhttp://releases.ubuntu.com/12.04/からダウンロードできます。この際64bit版を選んでください。
インストール方法自体は、バージョンや32bit/64bitの違いに関わらずほとんど同じなので悩むことはないと思います。

ifortのインストールはUbuntu12.04 IntelR Parallel Studio XE 2011 for Linuxのインストールの通りです。64bit用のオプションを選ぶことを忘れずに。

LAPACKとBLASのインストールはLinux で LAPACK バージョン 3.4.0 のダウンロードとビルドとインストール (Ubuntu や Fedora を使用)のシェルスクリプトで簡単に出来ます。
シェルスクリプトをlapackinst.shという名前で保存し、ホームディレクトリに置き、実行権限をつけた後、管理者権限で実行します。

chmod +x lapackinst.sh
sudo ./lapackinst.sh


最後にHiLAPWをインストールします。マニュアルPDFのはじめのほうにもインストール方法が書いてありますが、34ページ以降の新しいバージョンに向けて書かれたインストール方法を読むのが良いと思います。

テスト計算


差し当たりマニュアルPDFの例は全て出来ました。


ek.png
Fig.1: bcc鉄のバンド構造


参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: HiLAPW 

Scilabで二次元プロット

Scilabは数値計算をするために適したプログラミング言語です。
より具体的に言うと、Scilab「数値計算をするのに適した命令」と「グラフを描画するのに適した命令」の両方をあらかじめ持っているため、私たち人間がコードを書くために考えなければならない事が少なくてすみます。

「数値計算をするのに適した命令」の解説は、ねがてぃぶろぐのこれまでのエントリの中で色々と紹介してきました。例えば以下のようなものです。
これらの計算は、その都度適した方法でグラフ化をしてきました。
そこで今回は、グラフの描画方法に着目してエントリをまとめなおします。


目次



二次元プロットの基本


Scilabの二次元プロットはExcelで言うところの散布図が基本です。
例えば 0 < x < 1 の範囲で y = sin(2πx) をプロットするのは以下のようなプログラムになります。
X = linspace(0,1);
Y = sin(2 * %pi * X);

plot(X,Y);

001_201312281856542e7.png
Fig.1: 最も基本的な二次元プロット。


これがScilabのグラフ描画の基本です。plotでは青の実線がデフォルトになります。
Scilabには二次元グラフを描く命令としてplotのほかにplot2dというものもありますが、ねがてぃぶろぐではplotをメインに使用しています。


シンボルと線種


plotはオプションの引数を与えることでグラフの色やシンボルを指定することが出来ます。
// *** 正弦波の計算 ***
X = linspace(0,1);
Y1 = sin(2 * %pi * X);
Y2 = sin(3 * %pi * X);
Y3 = sin(4 * %pi * X);

// *** グラフのプロット ***
plot(X,Y1,'.r');
plot(X,Y2,'-*g');
plot(X,Y3,':b');

002_201312281936478ef.png
Fig.2: 色々なシンボル/線種の例


与えるオプションの識別子と描画の関係はScilab 入門: 6 グラフィックスに詳しく書かれています。以下はその引用です。

識別子識別子線種
r赤(Red)-実線
g緑(Green)-- 破線
b青(Blue)-.点破線
cシアン(Cyan):二点破線
m赤紫(Magenta)
y黄色(Yellow)
k黒(Black)
w白(White)
table.1: グラフの色と線種


識別子形状説明
++十字
x×クロス
**アスタリスク
o中抜きの円
.塗りつぶし円
^上向き三角
v下向き三角
>右向き三角
<左向き三角
s四角(Square)
dダイヤ(Diamond)
p星(Pentagram)
table.2: グラフのシンボル


Scilabで繰り返し計算: ロジスティック写像ではグラフを描いた時点でのScilabのバージョンが低かったため上手く行っていませんが、markersizeでマーカーの大きさを変更することが出来ます。
003_201312282216375a3.png
Fig.3: ロジスティック写像の例でマーカーの大きさを小さくしたもの

// *** ロジスティック写像の計算 ***
nmax = 200;
nout = 100;

A = [0:0.001:4];
X = 0.2 * ones(nmax,length(A));

for n = 1: nmax - 1 do
X(n+1,:) = A .* X(n,:) .* (1.0 - X(n,:));
end

// *** グラフのプロット ***
plot(A,X(nout:nmax,:),'.r','markersize',1);
zoom_rect([3,0,4,1]);


更に上記のサンプルではzoom_rect([xmin,ymin,xmax,ymax])でグラフの指定範囲のズームも行っています。


グラフの装飾


グラフの装飾として、グラフのタイトル、軸ラベル、凡例、グリッドの追加をします。以下に示すのがその例です。


004_20131228231755a1e.pngFig.4: グラフの装飾例

// *** 正弦/余弦の計算 ***
X = linspace(0,1);
Y1 = sin(2 * %pi * X);
Y2 = cos(2 * %pi * X);

// *** グラフのプロット ***
plot(X,Y1,'-r');
plot(X,Y2,'-b');

// *** グラフの装飾 ***
// タイトル
title("sin and cos");
// 軸ラベル
xlabel("x axis");
ylabel("y axis");
// 凡例
legend(["sin(x)","cos(x)"],3);
// グリッド
xgrid(color(128,128,128));


グラフのタイトルはtitleでx軸のラベルはxlabelで指定することが出来ます。y軸やz軸のラベルも同様にylabel,zlabelで指定できます。更にこれらのタイトルとラベルを一気に指定するのがxtitleです(が、私は別個に指定するほうが好きなのであまり使いません)。

グラフの凡例はlegendで描画することが出来ます。このlegendは必ずplotより後に書かなければなりません。
サンプルプログラムの最後のパラーメータである 3 は、凡例を配置する位置を指定するパラメータです。右上(1),左上(2),左下(3),右下(4)とマウスでクリックした場所に配置(5)が選べる他、マイナスの数値を指定するとグラフの外に配置することが出来ます。

xgridを使うとグラフ上の数値を読み取りやすくするためにグリッドを表示することが出来ます。引数なしでxgird()とすると黒の破線でグリッドが描かれます。
黒だと少し濃くて読みづらいときは別の色を割り当てます。これにはcolorを使いましょう。引数の3つの数値は色のRGB整数値です。上記のサンプルでは灰色にしてあります。


分割と差し込みグラフ


ひとつの図の中に複数のグラフを描くことが出来ます。
一番簡単な方法はsubplotを使う方法です。

subplot(m,n,p)とすることで、上下にm分割、左右にn分割した後、p番目のパネルにプロットを行います。


005_20131229022716c9a.png
Fig.5: subplotを使ったグラフの分割。下のパネルのx軸ラベルにも"x axis"と書いてあるはずだが切れてしまっている。

// *** 正弦/余弦の計算 ***
X = linspace(0,1);
Y1 = sin(2 * %pi * X);
Y2 = cos(2 * %pi * X);

// *** グラフのプロット ***
// 上部パネル
subplot(2,1,1);
plot(X,Y1,'-r');
// 軸ラベル
xlabel("x axis");
ylabel("y axis");
// 下部パネル
subplot(2,1,2);
plot(X,Y2,'-b');
// 軸ラベル
xlabel("x axis");
ylabel("y axis");


しかしながらsubplotでグラフを分割すると、Fig.5の様に軸ラベルが描画範囲からはみ出してしまうことがあるようです。
そういう場合は少し面倒ですがxsetech(wrect)を使います。
引数は wrect=[x,y,w,h]で x,y はグラフの左上の点の位置で w,h は幅と高さです。


006_20131229025350382.png
Fig.6: xsetechを用いたグラフの分割。subplotよりも細かい設定が可能。

// *** 正弦/余弦の計算 ***
X = linspace(0,1);
Y1 = sin(2 * %pi * X);
Y2 = cos(2 * %pi * X);

// *** グラフのプロット ***
// 上部パネル
//subplot(2,1,1);
xsetech([0,0,1,0.45]);
plot(X,Y1,'-r');
// 軸ラベル
xlabel("x axis");
ylabel("y axis");
// 下部パネル
//subplot(2,1,2);
xsetech([0,0.5,1,0.45]);
plot(X,Y2,'-b');
// 軸ラベル
xlabel("x axis");
ylabel("y axis");


この方法を応用すればScilabで差し込みグラフ:金属の比熱で描いたような差込グラフを作ることも出来ます。


グラフの縦横比


squareまたはisoviewを利用すれば、グラフの縦横比を1:1に固定することが出来ます。


007_20131230052535f40.png
Fig.7: isoviewを用いて縦横比を1:1としたプロット

// *** 円の計算 ***
T = linspace(0,1);
X = sin(2 * %pi * T);
Y = cos(2 * %pi * T);
// 描画範囲の計算
xmin = min(X);
xmax = max(X);
ymin = min(Y);
ymax = max(Y);

// *** グラフのプロット ***
// 円のプロット
plot(X,Y);
// 縦横比を1:1に設定する
//square(xmin,ymin,xmax,ymax);
isoview(xmin,xmax,ymin,ymax);


上記のサンプルはisoviewを用いてプロットしたものです。squareのコメントアウトをはずしてisoviewをコメントアウトすればsquareでの描画のサンプルにもなります。

これら二つの違いはウインドウサイズを変更するかしないかです。
なのですが、ややこしい事にx軸方向の最小値/最大値とy軸方向の最小値/最大値の設定する順番が異なっています。


複数枚のグラフ


subplotxsetechを利用して一枚の図を分割する他に複数枚の図にグラフをプロットする方法もあります。これにはscfを利用します。


008_20131230055603cfa.png
009_20131230055602c14.png
Fig.8-9: scfを使うと複数枚のグラフが描画できる

// *** 正弦/余弦の計算 ***
X = linspace(0,1);
Y1 = sin(2 * %pi * X);
Y2 = cos(2 * %pi * X);

// *** グラフのプロット ***
// 1枚目の図
scf(0);
plot(X,Y1,'-r');
// 軸ラベル
xlabel("x axis");
ylabel("y axis");
// 2枚目の図
scf(1);
plot(X,Y2,'-b');
// 軸ラベル
xlabel("x axis");
ylabel("y axis");



対数グラフなど


plot2dを利用すれば簡単に対数グラフを描くことが出来ます。しかしながらシンボルや線種の指定などの方法がplotと異なってしまいます。
そこでplotで対数グラフを描きたいわけですが、これはここまで紹介してきたほかの方法よりも少しだけややこしいです。

対数グラフに限った話ではないのですがgcaで軸の設定をしている行列を呼び出しaxes_propertiesに従って新しい設定を書き込めば、軸の設定を変更できます。以下は片対数グラフを描くサンプルです。


010_20131230071240a53.png
Fig.10: 片対数グラフ

// *** 常用対数の計算 ***
X = linspace(0.1,10);
Y = log10(X);

// *** グラフのプロット ***
subplot(2,1,1);
plot(X,Y);
xgrid(color(128,128,128));

// *** 片対数グラフのプロット ***
subplot(2,1,2);
a = gca();
a.data_bounds = [0.1,-1;10,1];
a.log_flags = "ln";
plot(X,Y);
xgrid(color(128,128,128));


このテクニックをマスターすれば、かなり色々なことに応用できるはずです。


関連エントリ



参考URL



付録


今回のエントリのサンプルは数が多い反面、ひとつひとつは短いので独立のファイルとしてアップロードすることはしません。各パートからコピー&ペーストしてください。


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab シンボル 線種 凡例 軸ラベル グラフの分割 差し込みグラフ 片対数グラフ 両対数グラフ 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性OPアンプPICCPA常微分方程式モンテカルロ解析ecaljodeトランジスタ状態密度インターフェースDOS定電流スイッチング回路PDS5022半導体シェルスクリプトレベルシフト乱数HP6632AR6452AI2C可変抵抗分散関係トランジスタ技術ブレッドボード温度解析反強磁性確率論バンドギャップセミナー数値積分熱設計非線形方程式ソルババンド構造絶縁偏微分方程式ISO-I2CLM358フォトカプラ三端子レギュレータカオスLEDシュミットトリガGW近似A/Dコンバータ発振回路PC817C直流動作点解析USBマフィンティン半径数値微分アナログスイッチTL43174HC4053カレントミラーサーボ量子力学単振り子チョッパアンプ補間2ちゃんねる開発環境bzqltyFFT電子負荷LDAイジング模型BSch基本並進ベクトルブラべ格子パラメトリック解析標準ロジックアセンブラ繰り返し六方最密充填構造SMPコバルトewidthFET仮想結晶近似QSGW不規則合金VCAMaximaGGA熱伝導cygwinスレーターポーリング曲線キュリー温度スイッチト・キャパシタ失敗談ランダムウォークgfortran抵抗相対論位相図スピン軌道相互作用VESTA状態方程式TLP621ラプラス方程式TLP552条件分岐NE555LM555TLP521マントル詰め回路MCUテスタFXA-7020ZR三角波過渡解析ガイガー管自動計測QNAPUPSWriter509ダイヤモンドデータロガー格子比熱熱力学awkブラウン運動起電力スーパーセル差し込みグラフ第一原理計算フェルミ面fsolveCIFxcrysden最大値最小値ubuntu最適化平均場近似OpenMPシュレディンガー方程式固有値問題井戸型ポテンシャル2SC1815TeX結晶磁気異方性OPA2277非線型方程式ソルバフラクタルFSM固定スピンモーメントc/agnuplotPGA全エネルギーfccマンデルブロ集合縮退正規分布キーボード初期値interp1multiplotフィルタ面心立方構造ウィグナーザイツ胞L10構造半金属二相共存ZnOウルツ鉱構造BaOSIC重積分磁気モーメント電荷密度化学反応クーロン散乱岩塩構造CapSenseノコギリ波デバイ模型ハーフメタルRealforceフォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt合金等高線凡例軸ラベル線種シンボルトラックボールグラフの分割MAS830LPIC16F785トランス入出力CK1026PC直流解析パラメータ・モデル等価回路モデル不規則局所モーメント関数フィッティング日本語ヒストグラムTS-112ExcelGimp円周率TS-110LMC662片対数グラフ三次元specx.fifortUbuntu文字列疎行列不純物問題ジバニャン方程式ヒストグラム確率論マテリアルデザインP-10境界条件連立一次方程式AACircuit熱拡散方程式HiLAPW両対数グラフ陰解法MBEナイキスト線図負帰還安定性Crank-Nicolson法EAGLE最小二乗法

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ