「みその計算物理学」では,一次元の拡散方程式(Diffusion equation)を,
CプログラムやJavaプログラムでの例が紹介されている。
http://www.geocities.jp/supermisosan/kakusan.html

本稿では,Igorで同じ問題を解く(間違っている可能性もあるが)。

一定の温度に熱せられた1次元の系が、両端が温度0にふれているという条件で解く。

Diffusion equation procedure window

拡散方程式(Diffusion equation)のプロシージャ

Function Diffusion_1D_main()
Variable i,j
for (i=0;i<1000;i+=1)
Make /O/N=100 $("wave"+num2str(i))=100 //初期値として100℃にする。
Bound_routine($("wave"+num2str(i))) //境界条件を適用する関数
endfor
for(i=0;i<999;i+=1) //
for(j=1;j<99;j+=1) // 拡散方程式を解く関数Diffuを呼び出す
Diffu($("wave"+num2str(i+1)),$("wave"+num2str(i)),j)
endfor
if (mod(i,100)==0) //iが100でと割り切れる場合に,Tableに追加
AppendToTable $("wave"+num2str(i))
endif
endfor
AppendToTable $("wave"+num2str(i)) //テーブルにwaveを追加する。
Displaywave0,wave100,wave200,wave300,wave400,wave500,wave600,
wave700,wave800,wave900,wave999 //グラフ化する。
End

Function Bound_routine(w) //境界条件を適用する関数
wave w
w[0]=0 //両端の温度を0℃にする。
w[99]=0
end

Function Diffu(w0,w1,jj) //拡散方程式を解く関数
Variable jj
wave w0,w1
w0[jj]=w1[jj]+0.5*(w1[jj+1]+w1[jj-1]-2*w1[jj])
end

図のように,
プロシージャーウィンドウで,
Diffusion_1D_main()
を実行する。

テーブルに計算結果であるwaveが追加され,
グラフ表示される。

Diffusion equation graph and table

グラフとテーブル

広告