拡散方程式のプログラミング

質問

拡散方程式を種々の初期条件、境界条件の下にプログラミング化(できればfortranがいいですが、その他でも大歓迎です)
まず、前回の質問たくさんの方に返事していただきました。
本当にありがとうございます！！！前回の質問を整理しました。

厚さがq(μm)、底面が気体が透過しないように設定されたペイントの表面の酸素濃度C(t)を時間幅0<=t<=t0(最大時間t0まで)周期的に変えます。
ここで、拡散方程式
σC/σt=Dm*(σC^2)/(σZ^2)
C：気体濃度、Z：厚さ方向、t：時間、Dm(μm2/s)：拡散係数
というのがあり、これを中心差分近似したものが、
(C(i,t+1)-C(i,t))/Δt=(C(i+1,t)-2C(i,t)+C(i-1,t))/(ΔZ^2)
となります。この拡散方程式に加えて、初期条件、境界条件を以下のように加えます。
・初期条件
(初期状態は厚さ方向に一様の真空状態)
→t=0;C=0(この式は正しいかわかりません↓↓)
・境界条件
ペイントの底面では厚さ方向にて濃度が同じ。
→Z=h;σC/σZ=0
この条件を離散化したものが、
→(C(1,t)-C(0,t))/ΔZとなります。
また、境界条件といたしまして、ペイント表面の濃度変化を正弦関数C(t)=Ksin2πft
(K：定数、f：周波数、t：時間)
として与えます。
ここで、補足をいたしますと、
ペイント表面の気体濃度を変化させると、その表面の気体はペイントの拡散係数に応じてペイント内部に拡散します(ペイント底面は気体は透過しない)。
ここで求めたい問題は、ペイント内のある時間での気体の厚さ方向での濃度分布、かつ、その厚さ方向の濃度を合計したペイント全体の気体濃度、さらに、一番求めたいのはペイント全体の気体濃度の時間変化(0<=t<=t0まで)プログラミング化するのですが。。。
様々な方のアドバイスを聞きたくてたくさん質問しています。。。申し訳ありません。

回答No.1

FT56F001

プログラムの足しになりませんが，この問題には理論解がありますね。
C(z,t)=Im[ C1*exp{z(1+i)/δ+iωt}+C2*exp{-z(1+i)/δ+iωt}]
C1=K/{1+exp[2h(1+i)/δ]}
C2=K/{1+exp[-2h(1+i)/δ]}
δ=√{2Dm/ω}(浸透深さ)，h:ペイント厚み
z=0:ペイント表面(C=Ksinωt，z=h:ペイント底面(∂C/∂z=0)
i=√(-1)，Imは虚数部(Im(x+iy)=y)，ω=2πf

できたプログラムの答えと理論解が合うか，確認してみてください。
/*　大きな計算プログラムを作ると，それっぽい答えが出てきても，
　　ホントに正しいのか確認するのは，いつも苦労します。*/

お礼

ありがとうございます！！
プログラムが完成した際の参考にぜひさせていただきます！！