VBAでの拡散計算

エクセルのVBAを使って添付画像のようなグラフを作成しようと考えています。

以下の計算で作成できそうなのですが、⊿tを10-5より小さく設定し、1000sec～の濃度変化が知りたいため表計算ではなくVBAを使ってみました。

表計算では、

t=0のときx0=C0(飽和濃度)、x>x0でC=0とし(初期条件)

x0では、

C(j+1,i) = C(j,i) +D * ( C(j,i-1) - C(j,i) ) / dx / dx * dt

x >x0では、

C(j+1,i) = C(j,i) +D * ( C(j,i-1) - 2 * C(j,i) + C(j,i+1) ) / dx / dx * dt

の計算を行い、セル表記が以下のようになりました。どの時間も物質量は一定です。(たぶん・・・)

　　t0、t1、t2、t3t、・・・
x1,60、58.8、57.6・・・
x2, 0、 1.2、2.38・・・
x3, 0、 0、0.02・・・

上の計算をVBAで以下のように書いてみました。

Sub diffusion_dry_up2()
Dim n As Integer, nt As Integer
Dim i As Integer, j As Integer
Dim b As Double, te As Double, dt As Double
Dim c0 As Double, cs As Double, d As Double
Dim x As Double, dx As Double, t As Double
Dim a As Double, cjp As Double, cj0 As Double
Dim cjm As Double
b = InputBox("配管の長さb(m)")
n = InputBox("配管の長さの方向分割数n")
te = InputBox("計算する時間長t(sec)")
dt = InputBox("時間増分dt(m)")
c0 = InputBox("配管底部の濃度c0(vol.%)")
cs = InputBox("時刻t=0の時の配管内の濃度cs(vol.%)")
d = InputBox("拡散係数d(m^2/sec)")
Sheet1.Cells(1, 2) = "配管の長さb(m)"
Sheet1.Cells(2, 2) = "配管の長さの方向分割数n"
Sheet1.Cells(3, 2) = "計算する時間長t(sec)"
Sheet1.Cells(4, 2) = "時間増分dt(sec)"
Sheet1.Cells(5, 2) = "配管底部の濃度c0(vol.%)"
Sheet1.Cells(6, 2) = "時刻t=0の時の配管内の濃度cs(vol.%)"
Sheet1.Cells(7, 2) = "拡散係数(m^2/sec)"
Sheet1.Cells(1, 3) = b
Sheet1.Cells(2, 3) = n
Sheet1.Cells(3, 3) = t
Sheet1.Cells(4, 3) = dt
Sheet1.Cells(5, 3) = c0
Sheet1.Cells(6, 3) = cs
Sheet1.Cells(7, 3) = d
nt = te / dt
dx = b / n
Sheet1.Cells(1, 1) = nt
t = 0
a = d * dt / dx ^ 2
Sheet1.Cells(1, 5) = t
Sheet1.Cells(1 + 1, 4) = -dx
Sheet1.Cells(1 + 2, 5) = c0
Sheet1.Cells(1 + n + 2, 4) = b + dx
Sheet1.Cells(1 + n + 3, 5) = cs
Sheet1.Cells(3, 5 + i) = a * (-cj0 + cjp) + cj0
For i = 1 To nt
t = t + dt
Sheet1.Cells(1, 5 + i) = t
Sheet1.Cells(1 + n + 3, 5 + i) = a * (cjp - 2 * cjo + cjm) + cj0
For j = 2 To n + 2
cjp = Sheet1.Cells(1 + j + 1, 5 + i - 1)
cj0 = Sheet1.Cells(1 + j + 0, 5 + i - 1)
cjm = Sheet1.Cells(1 + j - 1, 5 + i - 1)
Next j
Next i
linegraph 2, 4, 4 + n, nt + 2
End Sub
Sub linegraph(sr As Integer, sc As Integer, lr As Integer, lc As Integer)
ActiveSheet.ChartObjects.Add(200, 10, 240, 200).Select
ActiveChart.ChartWizard _
Source:=Range(Cells(sr, sc), Cells(lr, lc)), _
gallery:=xlXYScatter, _
Format:=3, _
PlotBy:=xlColumns, _
categoryLabels:=1, _
SeriesLabels:=0, _
HasLegend:="false", _
Title:="ex2", _
categoryTitle:="t", _
ValueTitle:="y", _
ExtraTitle:=""
End Sub

しかしまったく表計算のようになりませんでした。
a = d * dt / dx^2以降の書き込みが変だと思うのですが、どのようにすればよいのでしょうか。
また、上のような表記ではtを大きくdtを小さくするとエラーになってしまいます。

質問項目が多いですが、よろしくお願いします。

No.4 ベストアンサー 回答者： queuerev2 回答日時： 2012/10/09 12:40

No.2&3です。

質問者様のプログラムを修正して動作するようにしました。
ただし修正は最小限でありまだ改良の余地はかなりあります。
とはいえ、このプログラムの改良だけではNo.1様ご指摘の問題はあると思います。
そのため、質問者様の構想を実現するためには、プログラムをC言語などの高速な言語で書き直し、データは見たい部分だけを保存するようにする、などが考えられます。
なお、計算方法そのものに関することですが、tが大きくnが小さい場合、計算値の合計（物質量）が時間とともに少しずつ減っていくようです。原因は、xが最大のところで存在しないはずのより大きいxのデータを使用しているところにあるようです。今回のプログラムではこの点は修正しておりませんのでご了承ください。

Sub diffusion_dry_up2()
Dim n As Integer, nt As Integer
Dim i As Integer, j As Integer
Dim b As Double, te As Double, dt As Double
Dim c0 As Double, cs As Double, d As Double
Dim x As Double, dx As Double, t As Double
Dim a As Double, cjp As Double, cj0 As Double
Dim cjm As Double
b = InputBox("配管の長さb(m)")
n = InputBox("配管の長さの方向分割数n")
te = InputBox("計算する時間長t(sec)")
dt = InputBox("時間増分dt(m)")
c0 = InputBox("配管底部の濃度c0(vol.%)")
cs = InputBox("時刻t=0の時の配管内の濃度cs(vol.%)")
d = InputBox("拡散係数d(m^2/sec)")
Sheet1.Cells(1, 2) = "配管の長さb(m)"
Sheet1.Cells(2, 2) = "配管の長さの方向分割数n"
Sheet1.Cells(3, 2) = "計算する時間長t(sec)"
Sheet1.Cells(4, 2) = "時間増分dt(sec)"
Sheet1.Cells(5, 2) = "配管底部の濃度c0(vol.%)"
Sheet1.Cells(6, 2) = "時刻t=0の時の配管内の濃度cs(vol.%)"
Sheet1.Cells(7, 2) = "拡散係数(m^2/sec)"
Sheet1.Cells(1, 3) = b
Sheet1.Cells(2, 3) = n
Sheet1.Cells(3, 3) = te '修正
Sheet1.Cells(4, 3) = dt
Sheet1.Cells(5, 3) = c0
Sheet1.Cells(6, 3) = cs
Sheet1.Cells(7, 3) = d
nt = te / dt
dx = b / n
Sheet1.Cells(1, 1) = nt
t = 0
a = d * dt / dx ^ 2
Sheet1.Cells(2, 5) = t '修正
Sheet1.Cells(1 + 1, 4) = -dx
Sheet1.Cells(1 + 2, 5) = c0
Sheet1.Cells(1 + n + 2, 4) = b + dx
For j = 3 To n + 2 '追加
Sheet1.Cells(1 + j, 5) = cs '修正
Next j '追加
For i = 1 To nt
t = t + dt '移動
cjp = Sheet1.Cells(1 + 2 + 1, 5 + i - 1) '追加
cj0 = Sheet1.Cells(1 + 2 + 0, 5 + i - 1) '追加
Sheet1.Cells(3, 5 + i) = a * (-cj0 + cjp) + cj0 '移動修正
Sheet1.Cells(2, 5 + i) = t '修正
For j = 3 To n + 2 '修正
cjp = Sheet1.Cells(1 + j + 1, 5 + i - 1)
cj0 = Sheet1.Cells(1 + j + 0, 5 + i - 1)
cjm = Sheet1.Cells(1 + j - 1, 5 + i - 1)
Sheet1.Cells(1 + j + 0, 5 + i) = a * (cjp - 2 * cj0 + cjm) + cj0 '移動修正
Next j
Next i
linegraph 2, 5, 3 + n, nt + 5 '修正
End Sub
Sub linegraph(sr As Integer, sc As Integer, lr As Integer, lc As Integer)
ActiveSheet.ChartObjects.Add(200, 10, 240, 200).Select
ActiveChart.ChartWizard _
Source:=Range(Cells(sr, sc), Cells(lr, lc)), _
gallery:=xlXYScatter, _
Format:=3, _
PlotBy:=xlRows, _
categoryLabels:=1, _
SeriesLabels:=0, _
HasLegend:="false", _
Title:="ex2", _
categoryTitle:="t", _
ValueTitle:="y", _
ExtraTitle:="" '修正 (PlotBy)
End Sub

この回答へのお礼

回答ありがとうございます。

上端の境界条件などについては自分でも検討してみます。

ありがとうございました。

お礼日時：2012/10/09 20:14

No.3 回答者： queuerev2 回答日時： 2012/10/05 21:57

No.2です。

回答し終わってみたらすでに先に回答を付けられた方がいらっしゃいました。
もしNo.1様の回答で十分であれば、表計算内容の補足説明をいただかなくとも結構です。

No.2 回答者： queuerev2 回答日時： 2012/10/05 21:48

もうちょっと解説いただけないでしょうか。

特に表計算で、どのセルにどんな数式を書けばよいのでしょうか。

質問者様は以下のような結果を得ておられます。
> 　　t0、t1、t2、t3t、・・・
> x1,60、58.8、57.6・・・
> x2, 0、 1.2、2.38・・・
> x3, 0、 0、0.02・・・
そして、(x=?)x0の場合とx>x0の場合について式を書かれています。
しかしこの式をどのセルにどのように書けばいいのかがわかりません。
（t0のときの60, 0, 0は定数でしょうか？　t1, t2の時の6個の数値、つまり58.8から0.02までですが、x1～x3ということは、x>x0の式を使用するのでしょうか？　C(j+1,i)等の２つの引数のうちどちらがxでどちらがtなのでしょうか？　などです）

表計算の式がわかれば、私のような拡散計算が全くわからない者でも同じ結果を得るVBAプログラムを作ることができると思うので、解説をお願いする次第です。

この回答への補足

回答ありがとうございます。

濃度計算についてですが、

　　t0、 t1、 t2、 t3t、・・・
x1,60、58.8、57.6 ・・・
x2, 0、 1.2、2.38 ・・・
x3, 0、 0、0.02 ・・・

t=0のときx1=C0(飽和濃度)、x>x1(x2,x3・・・)でC=0とし(初期条件)

x1(x1の右側の行のセル)では、

C(j+1,i) = C(j,i) +D * ( C(j,i-1) - C(j,i) ) / dx / dx * dt　　・・・(1)

j:時間tの引数、i:位置xの引数です。

行x0のセルは式(1)を用いて各時間の濃度を算出しています。(t1で58.8、t2で57.6の濃度は式(1)で算出しました。)


x >x1(x2,x3・・・)のセルでは、

C(j+1,i) = C(j,i) +D * ( C(j,i-1) - 2 * C(j,i) + C(j,i+1) ) / dx / dx * dt　・・・(2)

式(2)はx0以外の行(x2,x3・・・)のセルで用いました。(x2で1.2、2.38…、x3で0、0.02・・・は式(2)で算出しています。)

t=0における濃度は、式(1)、式(2)は用いておらず、自分で設定した値なので60, 0, 0は定数と考えて頂いても大丈夫です。

以上のように1行目は式(1)、2行目以降は式(2)で計算しました。

補足日時：2012/10/06 00:44

No.1 回答者： ki073 回答日時： 2012/10/05 21:32

まず確認ですが、

＞⊿tを10-5より小さく設定し、1000sec～の濃度
とありますが、これって10の-5乗ですよね。もしそうだとすると、
dt=10^-5でte=1000だと
nt = te / dt
からntが10^8になります。ntはInteger型なのでover flowしてしまいます。それと
For i = 1 To nt
からiも
http://www.geocities.jp/cbc_vbnet/kisuhen/hensuu …
まあLongにすれば計算はできますが、

もっと大きな問題ですが、10^8回ループを回る訳ですが、セルにどんどん書き込んでいるようですが、足りるのでしょうか。
結果を1秒ごとに書き出すとかにしないと足りなくなりますが。

もうひとつこちらの方が当たっていると深刻だと思いますが、
計算が実用的な時間で終わるのでしょうか？
VBAは結構遅いのでかなり難しいように思います。適当な回数で時間を測定してみてはいかがでしょうか。
Sheet1.Cells()のようにセルを直接アクセスするのをやめて、配列を使うと少しは速くなるようにも思いますが、それも限界が。

いろいろ手はあるのですが、

この回答への補足

回答ありがとうございます。

表計算で行うと1秒間のデータを作成しただけでもファイル容量が大きくなったため、別の方法(VBA)で試してみました。

ネットや拡散計算の文献では添付画像のようなグラフを見かけるため、出来ないことはないとは思うのですが・・・

試行錯誤しているところです。

補足日時：2012/10/06 00:40