1次元拡散方程式(陰解法)をExcel VBAで計算。
さて、昨日Numerical Recipesの三重対角行列を計算するルーチンを作ったのは、この計算をするため。
∂u(t,x) /∂t = D ∂2u(t,x) / ∂x2
を陰的に差分化すると、
uik+1 = uik + s * (ui-1k+1 +ui+1k+1 -2*uik+1 )
となる。なので三重対角行列が計算できればこれも計算できる。
こんなルーチンでどうでしょう。
Option Explicit
Private Sub CommandButton1_Click()
Dim n As Long
Dim a() As Double, b() As Double, c() As Double
Dim r() As Double, u() As Double
Dim i As Integer, k As Integer
Dim s As Double
Application.ScreenUpdating = False
s = 1#
n = 128
ReDim a(n), b(n), c(n), r(n), u(n)
For i = 1 To n
b(i) = 1# + 2# * s
a(i) = -s#
c(i) = -s#
Next i
c(1) = -2# * s
a(n) = -2# * s
For i = 1 To n
If i > 0.5 * n - 5 And i < 0.5 * n + 5 Then
r(i) = 1#
Else
r(i) = 0#
End If
Worksheets("Sheet1").Cells(i + 1, 2) = r(i)
Next i
For k = 1 To 500
Call tridag(a, b, c, r, u, n)
For i = 1 To n
r(i) = u(i)
Worksheets("Sheet1").Cells(i + 1, k + 2) = r(i)
Next i
Next k
Application.ScreenUpdating = True
End Sub
で計算した結果がこちら。s=1でも計算出来てるようだ。
で、まだこれで終わりじゃなくてまだまだ続きが、、、
« Numerical Recipes in Cの三重対角行列の方程式を解く関数をExcel VBAに移植 | トップページ | SMAPのソフトバンク プラチナバンドスタートのCMを見た。民放各局放映はチキンラーメンのCM以来かな? »
「学問・資格」カテゴリの記事
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2021年2月16日 IEEE Microwave Magazineの特集はオールデジタルのRFID、Microwave JournalはEバンド ミリ波通信に衛星や気球を使う話、アメリカの半導体企業がバイデンに投資を迫る、(2021.02.17)
- カオスを生じる電気回路、Chua’s circuitをLTspiceで回路シミュレーションしてみる。(2021.02.19)
- Labyrinth Chaos(迷宮カオス)を生むThomas-Rössler方程式のパラメータbを色々変えて、Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使って計算してGIFアニメ(2021.02.16)
- フィッツヒュー・南雲 (FitzHugh-Nagumo) 方程式をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算。(2021.02.23)
- 「水晶振動子の等価回路計算」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUP! インピーダンスの大きさと位相がグラフ化できる。(2021.02.12)
« Numerical Recipes in Cの三重対角行列の方程式を解く関数をExcel VBAに移植 | トップページ | SMAPのソフトバンク プラチナバンドスタートのCMを見た。民放各局放映はチキンラーメンのCM以来かな? »
コメント