12. 常微分方程式

常微分方程式の初期値問題のXLPackによる解き方の例を示します。

例題

次の初期値問題を解く。

    dy/dt = 1/(2 - t)2   (y(0) = 0.5)

この例題は解析解がわかっており、次のようになります。

    y = 1/(2 - t)

この関数の形は次のとおりです。

12-01

VBAプログラムを使用した解き方 (1)

VBAサブルーチンDerkfを使ったプログラム例を示します。

Dim Neval As Long

Sub F(N As Long, T As Double, Y() As Double, Yp() As Double)
    Yp(0) = 1 / (2 - T) ^ 2
    Neval = Neval + 1
End Sub

Sub Start()
    Const NMax = 10, RowR = 4, ColR = 1
    Dim N As Long, Info As Long
    Dim Y(NMax) As Double, T As Double, Tout As Double
    Dim RTol(0) As Double, ATol(0) As Double, I As Integer
    N = 1
    RTol(0) = 0.000000000001 '** 1.0e-12
    ATol(0) = RTol(0)
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Info = 0
    For I = 1 To 19
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        Call Derkf(N, AddressOf F, T, Y(), Tout, RTol(), ATol(), Info)
        If Info <> 1 And Info <> 3 Then
            Cells(RowR + I, ColR + 2) = "Error: " + Str(Info)
            Exit For
        End If
        Cells(RowR + I, ColR + 1) = Y(0)
        Cells(RowR + I, ColR + 2) = Neval
    Next
End Sub

Derkfは4(5)次ルンゲ・クッタ・フェールベルグ法の適応ルーチンで、目標精度を与えるときざみ幅を自動的に調節して解を求めます。微分方程式を外部関数として用意し、初期値と要求精度を指定して呼び出します。なお、Nevalは単に関数評価回数をカウントするための変数です。

このプログラムを実行すると、次の結果が得られました。

12-02

11~12桁の精度で正しい結果が得られています。なお、tが大きくなり変化が急になってくると、精度を保つために関数評価回数が増えているのがわかります。

VBAプログラムを使用した解き方 (2)

リバースコミュニケーション版のVBAサブルーチンDerkf_rを使ったプログラム例を示します。

Dim Neval As Long

Sub F(N As Long, T As Double, Y() As Double, Yp() As Double)
    Yp(0) = 1 / (2 - T) ^ 2
    Neval = Neval + 1
End Sub

Sub Start()
    Const NMax = 10, RowR = 4, ColR = 1
    Dim N As Long, Info As Long
    Dim Y(NMax) As Double, T As Double, Tout As Double
    Dim RTol(0) As Double, ATol(0) As Double, I As Integer
    Dim TT As Double, YY(NMax - 1) As Double, YYp(NMax - 1) As Double, IRev As Long
    N = 1
    RTol(0) = 0.000000000001 '** 1.0e-12
    ATol(0) = RTol(0)
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Info = 0
    For I = 1 To 19
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        IRev = 0
        Do
            Call Derkf_r(N, T, Y(), Tout, RTol(), ATol(), Info, TT, YY(), YYp(), IRev)
            If IRev = 1 Then Call F(N, TT, YY(), YYp())
        Loop While IRev <> 0
        If Info <> 1 And Info <> 3 Then
            Cells(RowR + I, ColR + 2) = "Error: " + Str(Info)
            Exit For
        End If
        Cells(RowR + I, ColR + 1) = Y(0)
        Cells(RowR + I, ColR + 2) = Neval
    Next
End Sub

目的関数を外部関数として与えるのではなく、IRev = 1のときにTTおよびYYの値を使って関数値を計算しYYp()に入れて再度Derkf_rを呼び出します。リバースコミュニケーション版の詳細についてはこちらを参照してください。

このプログラムを実行すると上と同じ結果が得られます。

Top