1次元拡散方程式(陰解法)をExcel VBAで計算。

さて、昨日Numerical Recipesの三重対角行列を計算するルーチンを作ったのは、この計算をするため。

∂u(t,x) /∂t = D ∂²u(t,x) / ∂x²

を陰的に差分化すると、

u_i^k+1 = u_i^k + s * (u_i-1^k+1 +u_i+1^k+1 -2*u_i^k+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でも計算出来てるようだ。

で、まだこれで終わりじゃなくてまだまだ続きが、、、