[LibreOffice/OpenOffice Basic] 微分方程式の数値的解法 (2)

本ページは「Excel で学ぶ Visual Basic for Applications (VBA)」内のページ「微分方程式の数値的解法 (2) VBA + Excel」を
LibreOffice/OpenOffice Basic に対応させたものである。

前ページに引続き、微分方程式の数値解法を学ぶ。
本ページでは LibreOffice/OpenOffice Basic のプログラミング機能を用いて実現する。

Excel / OpenOffice で学ぶフーリエ変換入門 表紙 本章で解説する微分方程式の数値的解法を用いた 「Excel / OpenOffice で学ぶフーリエ変換入門」では
マクロをソース閲覧可能な形で用いています。

Euler 法

Euler 法で数値積分を実現するのが以下の LibreOffice/OpenOffice Basic プログラムである。

ポイントは 記述して Main 関数を実行してみよう。
(Basic 用 Editor の起動方法は「第一回：プログラムを書き始めるまでの準備」参照) 

Option Explicit

Sub Main()
Dim t As Double, tmax As Double
Dim dt As Double

Dim x1 As Double, x2 As Double    ' Euler 法用の変数
Dim count As Long, cellindex As Long

  t = 0 : tmax = 30 : dt = 0.01
  count = 0 : cellindex = 0
  x1 = 0 : x2 = 0 ' Euler 法用変数の初期化
 
  ThisComponent.addActionLock() 
 	
  While t < tmax
    If count Mod 10 = 0 Then  ' 10 回に 1 回だけ描画
      ThisComponent.Sheets(0).getCellByPosition(0, cellindex).Value = t
      ThisComponent.Sheets(0).getCellByPosition(1, cellindex).Value = x1
      ThisComponent.Sheets(0).getCellByPosition(2, cellindex).Value = 0.5 * t * Sin(2 * t)
    	
      cellindex = cellindex + 1
    End If
 	
    Euler x1, x2, t, dt
 		
    count = count + 1
    t = t + dt 	
  Wend
 	
  ThisComponent.removeActionLock() 
End Sub

Sub Euler(x1 As Double, x2 As Double, t As Double, dt As Double)
Dim f1 As Double, f2 As Double

  f1 = x2
  f2 = -4 * x1 + 2 * Cos(2 * t)
	
  x1 = x1 + dt * f1
  x2 = x2 + dt * f2
End Sub


実際にプログラムを書いてみて、以下のことに気づいただろうか。

f2(x1, x2, t) = -4 x1 + 2 cos(2t) という式を実現するため、
前ページの シート のみの手法では

-4*B2+2*cos(2*A2)

と記述したのに対し、今回の VBA では、
-4*x1+2*cos(2*t)

と記述できる。

つまり、シート では変数名として B2 や A2 といったセル番号を用いねばならないのに対し、
Basic では t や x1 などという、元の式に近い変数を (宣言すれば) 使える。

シート のみの時のように、B2 や A2 という抽象的な変数名で式を記述することは (実際にやってみるとすぐにわかることだが)、
非常にわずらわしく、間違いが混入しやすい。

一方、Basic のようにある程度分かりやすい変数名でプログラミングできれば、 間違いが起こりにくく、また計算の内容も非常にわかりやすい。
(もちろん、Basic というプログラミング言語にある程度慣れていなければならないという制約はある)

そのようなわけで、ある程度複雑な処理をしようとすると、シート のみで頑張るよりも、
Basic を利用した方が逆に簡単になる場合がある、ということに気づいてもらえると嬉しい。

さて、プログラムを記述して実行してみよう。
Basic 用 Editor 上では何も変化が起こらないが、
通常のシートに戻ってみると、以下のように 3 列のデータがプログラムにより 書き込まれている。

これらは、左から順に「時刻」、「Euler 法による数値解」、「解析解」である。



このグラフを描くために、3 列を選択しよう。
「A」から「C」をマウスをドラッグするように選択するか、
「A」を一回クリックした後、「Ctrl キーを押しながら、『B』、『C』を一回ずつクリック」すると良いのだった。



そしてメニューから「挿入」→「グラフ」を選択し、散布図を選択する。
得られるグラフは以下のようになる。(下図は若干の調整を済ませてある)

シートのみの時と同様、Euler 法による数値解と解析解とには大きな誤差があることがわかる。

なお、凡例の記述は図の用に１行目に名前を記述する行を挿入することで行った。
Basic による凡例の記述は「実用例(2) グラフを自動的に描く」で紹介する。



なお、LibreOffice/OpenOffice では、グラフを表示したままグラフのデータを更新しようとすると、
処理に非常に時間がかかるようになってしまう。

その回避方法は「実用例(2) グラフを自動的に描く」で紹介するが、
ここではひとまず、「データ更新の前にグラフは削除する」ことを心がけて欲しい。

4 次の Runge-Kutta 法

誤差の大きい Euler 法に代わって広く使われているのが4 次の Runge-Kutta 法である。

4 次の Runge-Kutta 法の細かな式はここでは重要ではないが、興味のある人は このリンク先を見て欲しい。

さて、プログラム全体を記述すると以下のようになる。 

Option Explicit

Sub Main()
Dim t As Double, tmax As Double
Dim dt As Double

Dim x1 As Double, x2 As Double    ' Euler 法用の変数
Dim y1 As Double, y2 As Double    ' Runge-Kutta 法用の変数
Dim count As Long, cellindex As Long

  t = 0 : tmax = 30 : dt = 0.01
  count = 0 : cellindex = 0
  x1 = 0 : x2 = 0 ' Euler 法用変数の初期化
  y1 = 0 : y2 = 0 ' Runge-Kutta 法用変数の初期化
 
  ThisComponent.addActionLock() 
 	
  While t < tmax
    If count Mod 10 = 0 Then  ' 10 回に 1 回だけ描画
      ThisComponent.Sheets(0).getCellByPosition(0, cellindex).Value = t
      ThisComponent.Sheets(0).getCellByPosition(1, cellindex).Value = x1
      ThisComponent.Sheets(0).getCellByPosition(2, cellindex).Value = 0.5 * t * Sin(2 * t)
      ThisComponent.Sheets(0).getCellByPosition(3, cellindex).Value = y1
    	
      cellindex = cellindex + 1
    End If
 	
    Euler x1, x2, t, dt
    RK_4 y1, y2, t, dt
 		
    count = count + 1
    t = t + dt 	
  Wend
 	
  ThisComponent.removeActionLock() 
End Sub

Sub Euler(x1 As Double, x2 As Double, t As Double, dt As Double)
Dim f1 As Double, f2 As Double

  f1 = x2
  f2 = -4 * x1 + 2 * Cos(2 * t)
	
  x1 = x1 + dt * f1
  x2 = x2 + dt * f2
End Sub	

Sub RK_4(x1 As Double, x2 As Double, t As Double, dt As Double)
Dim k11 As Double, k12 As Double
Dim k21 As Double, k22 As Double
Dim k31 As Double, k32 As Double
Dim k41 As Double, k42 As Double

    k11 = dt * func1(x1, x2, t)
    k12 = dt * func2(x1, x2, t)
    
    k21 = dt * func1(x1 + k11 / 2, x2 + k12 / 2, t + dt / 2)
    k22 = dt * func2(x1 + k11 / 2, x2 + k12 / 2, t + dt / 2)
    
    k31 = dt * func1(x1 + k21 / 2, x2 + k22 / 2, t + dt / 2)
    k32 = dt * func2(x1 + k21 / 2, x2 + k22 / 2, t + dt / 2)
    
    k41 = dt * func1(x1 + k31, x2 + k32, t + dt)
    k42 = dt * func2(x1 + k31, x2 + k32, t + dt)

    x1 = x1 + (k11 / 6 + k21 / 3 + k31 / 3 + k41 / 6)
    x2 = x2 + (k12 / 6 + k22 / 3 + k32 / 3 + k42 / 6)

End Sub

Function func1(x1 As Double, x2 As Double, t As Double)
    func1 = x2
End Function
 
Function func2(x1 As Double, x2 As Double, t As Double)
    func2 = -4 * x1 + 2 * Cos(2 * t)
End Function


このプロシージャ RK_4 の意味を理解する必要は本講義ではないが、以下のことは注意しておいて欲しい。

RK_4 は見たところ行数も多く複雑に思えるが、式と比較しながら注意深く観察すれば、
理解できないほどの複雑さではない。変数名と式にある程度対応がつけられるからである。

一方、これを シートのみで実現しようとすると、不可能ではないが非常に大変である。
(Euler 法の計算のために 「t」、「x1」、「x2」、「f1」、「f2」の５列を利用したが、
4 次の Runge-Kutta 法には何列のデータが必要だろうか？)

このように、ある程度簡単な計算ならシートのみの機能で十分であるが、
複雑な計算になると、Basic を用いた方が簡単な場合がある。

さて、Main 関数を実行してシートに戻ると、以下のように 4 列のデータが完成している (文字は出ないので注意)。
左から順に 「時刻」、「Euler 法による数値解」、「解析解」、「4 次の Runge-Kutta 法による数値解」である。

グラフを描くために 4 列を同時選択して、メニューから「挿入」→「グラフ」を実行して散布図を選択する。





すると、以下のようなグラフが得られる (下図は若干の調整済み)。

「解析解」と「4 次の Runge-Kutta 法による数値解」はほぼ同じ値を取って重なっているため、
解析解しか見えなくなっている。

すなわち、4 次の Runge-Kutta 法は Euler 法に比べると精度が非常に高いことがわかる。






[LibreOffice/OpenOffice Basic] 微分方程式の数値的解法 (3) 汎用性を高める→

Excel で学ぶ Visual Basic for Applications (VBA)に戻る