y y = y(x) dy/dx = f(x, y) (2.1) 4 y = y(x) (x, y)=(x 0, y 0 ) x h (x 0, y 0 ) x 1 =x 0 +h y y 1 x 1 x 2 =x 1 + h y y 2 y 3 y 2 y 1 y 0 x 0 h

2

常微分方程式の数値解法

2.1

１階常微分方程式のオイラー法

x0 h x1 x2 x3 x y 厳密解曲線 オイラー法に による近似解曲線 y0 y = y(x) h h y1 y2 y3 図

1:

オイラー法

1

階常微分方程式

dy/dx = f (x, y)

(2.1)

を考えます4。この微分方程式の解を

y = y(x)

とします。初期 条件を

(x, y) = (x

0

, y

0

)

とし，

x

軸上の刻み幅を

h

とします。  オイラー法は初期条件

(x

0

, y

0

)

から出発して隣接点

x

1

= x

0

+h

での

y

の近似値

y

1を求め，次々に

x

1の隣接点

x

2

= x

1

+ h

での

y

の近似値

y

2を求め・・・と，次々に同じこと繰り返して得られ る

y

i値を結んで微分方程式の近似解を得ようというものです。  オイラー法は解析的にはテイラー展開の

1

次の項をとったも のとなります。

x

i近傍でテイラー展開し

1

次の項だけをとると

y

i+1

= y(x

i

+ h) = y(x

i

) +

dy

dx

¯

¯

¯

¯

x=xi

h +

1

2

d

2

_y

dx

2

¯

¯

¯

¯

x=xi

h

2

+ · · ·

(2.2)

; y

i

+ f (x

i

, y

i

)h

(2.3)

 次に幾何学的に見ていきます。初期値

(x

0

, y

0

)

での

y = y(x)

の接線の方程式は，傾き

f (x

0

, y

0

)

をもつ直線の方 程式なので

y = y

0

+ f (x

0

, y

0

)(x − x

0

)

となります。

x

1

= x

0

+ h

での

y

の値を

y

1とすると

y

1

= y

0

+ f (x

0

, y

0

)(x

1

− x

0

) = y

0

+ f (x

0

, y

0

)h

次に，点

(x

1

, y

1

)

における接線の方程式を求めると，接線の傾きは

f (x

1

, y

1

)

で点

(x

2

, y

2

)

を通るとすると

y

₂

= y

₁

+ f (x

₁

, y

₁

)h

以下，同様にして点

x

iの隣接点

x

i+1

(= x

i

+ h)

での

y

値を

y

i+1とすると

y

i+1

= y

i

+ f (x

i

, y

i

)h (i = 0, 1, 2, · · · , n − 1)

(2.4)

このようにして得られた点

(x

0

, y

0

), (x

1

, y

1

), · · · , (x

i

, y

i

), · · · , (x

n

, y

n

)

を順次結んでいって近似解曲線を得るわけ です。オイラー法は簡単でわかりやすいですが，刻み幅

h

を小さく取らないと誤差が大きくなるし，また刻み幅 を小さくとっても

x

の増加とともに誤差が累積していくという問題があります。そこで，次節ではこの点を改良 した修正オイラー法を紹介します。 ●オイラー法

VBA

プログラム 微分方程式

dy

dx

=

1

2

x(1 − y)

（初期値：

x = 0, y = 4

）

(2.5)

の近似解を求めます。ちなみに，この方程式の厳密解（解析解）は

y = e

−x2/4

(3 + e

x2

/4)

 刻み幅は

h = 0.01

としました。計算結果の一部を載せておきます。「

Calc

ボタン」で計算結果と解曲線を表示 し，「

Erase

ボタン」で消去します。厳密解と近似解を見比べてください。 4 幾何学的にはxy平面上に，点(x, y)における解曲線y = y(x)の接線の傾きを与えているということになります。

＜ソースコード＞ Option Explicit ’ モジュール内で有効な変数の定義 Private i As Integer Private cnt As Integer Private FL As Integer ’ グラフ消去のフラグ ’ 関数 f(x) ：dy/dx=f(x) Function Df(x, y) Df = x * (1 - y) / 2 End Function Function Sf(x)  ’ 厳密解 Sf = Exp(-x ^ 2 / 4) * (3 + Exp(x ^ 2 / 4)) End Function Sub オイラー法 ()

’---Dim x As Single, y As Single, h As Single, V As Single ’ 初期値 x = Range("C4").Value:y = Range("C5").Value FL = 1 h = Range("E4").Value ’ 刻み幅 cnt = Range("E5").Value ’ 繰り返し回数 ’ 逐次計算 For i = 0 To cnt Cells(i + 8, 1) = i Cells(i + 8, 2) = x Cells(i + 8, 3) = y Cells(i + 8, 4) = Sf(x) ’ 厳密解の値 y = y + h * Df(x, y) Cells(i + 8, 5) = Sf(x) - y ’ 厳密解との差 x = x + h Next i ’ 近似解曲線を散布図で描く V = 20 ’ 作図範囲 (厳密解との差がよく分かる範囲。V=cnt とすると全範囲） With ActiveSheet.Shapes.AddChart.Chart .ChartType = xlXYScatterSmoothNoMarkers ’ 平滑線付き散布図データマーカーなし .SetSourceData Range(Cells(7, 2), Cells(7 + V, 4)) ’ 作図範囲

End With End Sub

’ 計算結果の消去 Sub Erase_Click()

Range(Cells(8, 1), Cells(8 + cnt, 5)).Clear ’ セルに格納したデータを消去 If FL = 1 Then ActiveSheet.ChartObjects("1").Activate ’ グラフの消去 Selection.Delete FL = 0 End If End Sub

2.1.1

修正オイラー法 傾き同じ xi xi+1 yi h k1 k2 x y = dy dx ¯ ¯ ¯ ¯ x=xi 1 2(k1+ k2) 最終的な増分 図

2:

修正オイラー法 微分方程式の解曲線が比較的穏やかな曲線の場合は いいのですが，変化の激しい解曲線となると隣接点の 情報を取り込んでおかないと大きなズレが生じて精度 が落ちてしまいます。修正オイラー法はこの欠点を補 うもので，隣接点の情報を取り込んでいます。  図

2

をご覧ください。この方法は

3

ステップに分か れていて，

1.

まず点

(x

i

, y

i

)

での微係数

f (x

i

, y

i

)

から隣接点

x

i+1

(= x

i

+ h)

での

y

値の増分

k

1を求めます。

2.

次に増分から得られた座標点

(x

i

+ h, y

i

+ k

1

)

に 注目し，その点における微係数

f (x

i

+ h, y

i

+ k

1

)

をもつ接線を点

(x

i

, y

i

)

を通るように引き直して

y

_iからの増分

k

2を求めます。

3.

最後に

2

つの増分

k

1

, k

2の算術平均をとってこれ を最終的な増分とするものです。























k

1

= f (x

i

, y

i

)h

k

2

= f (x

i

+ h, y

i

+ k

1

)h

y

i+1

= y

i

+

1

₂

(k

1

+ k

2

)

(2.6)

●修正オイラー法

VBA

プログラム オイラー法と比較するために先程の微分方程式を取り上げます。

dy

dx

=

1

2

x(1 − y)

（初期値：

x = 0, y = 4

）

(2.7)

結果は次の通りでオイラー法と比べて厳密解との一致精度が高いことがわかります。 ＜ソースコード＞ Option Explicit ’ モジュール内で有効な変数の定義 Private i As Integer Private cnt As Integer Private FL As Integer ’ グラフ消去のフラグ ’ 関数

---Function Df(x, y) Df = x * (1 - y) / 2 End Function ’ 厳密解 ---Function Sf(x) Sf = Exp(-x ^ 2 / 4) * (3 + Exp(x ^ 2 / 4)) End Function ’**************************** Sub 修正オイラー法 () ’****************************

Dim x As Single, y As Single, V As Integer Dim k1 As Single, k2 As Single

Dim h As Single ’ 初期値 x = Range("C4").Value y = Range("C5").Value FL = 1 h = Range("E4").Value ’ 刻み幅 cnt = Range("E5").Value ’ 繰り返し回数 ’ 逐次計算 For i = 0 To cnt Cells(i + 8, 1) = i Cells(i + 8, 2) = x Cells(i + 8, 3) = y Cells(i + 8, 4) = Sf(x) ’ 厳密解の値 Cells(i + 8, 5) = Sf(x) - y ’ 厳密解との差 k1 = Df(x, y) * h k2 = Df(x + h, y + k1) * h y = y + (1 / 2) * (k1 + k2) x = x + h Next i ’ 近似解曲線を散布図で描く V = 20 ’ 作図範囲 (厳密解との差がよく分かる範囲。V=cnt とすると全範囲） With ActiveSheet.Shapes.AddChart.Chart .ChartType = xlXYScatterSmoothNoMarkers ’ 平滑線付き散布図データマーカーなし .SetSourceData Range(Cells(7, 2), Cells(7 + V, 4)) ’ 作図範囲

End With End Sub

’ 計算結果の消去 Sub Erase_Click()

Range(Cells(8, 1), Cells(8 + cnt, 5)).Clear ’ セルに格納したデータを消去 If FL = 1 Then ActiveSheet.ChartObjects("1").Activate ’ グラフの消去 Selection.Delete FL = 0 End If End Sub 【補足】ホイン法 同じ傾き x0 x1 y0 h 2 y1(O) x y f (x0, y0) h 2 傾き f (x1, y1O) 修正オイラー法の

1

つで，ドイツの数学者ホイン5が開発した方法で す。説明は省略しますが，

Excel

のワークシートをご覧ください。 オイラー法の基本式は

y

_i+1(O)

= y

_i

+ f (x

_i

, y

_i

)h

(2.8)

これを使って

y

i+1

= y

i

+

h

₂

f (x

i+1

, y

i+1(O)

)

(2.9)

とします。

5_{Karl Heun}

2.2

2 階常微分方程式のオイラー法

2

階常微分方程式の場合は

1

階連立常微分方程式に変換したあと各微分方程式にオイラー法を適用します。

d

2

_y

dx

2

= f

µ

x, y,

dy

dx

¶

,

初期値：

y(a) = y

0

, y

0

(a) = y

00

(2.10)

z = dy/dx

とおくと次の

1

階連立常微分方程式の変換されます。















dy

dx

= z = g(x, y, z),

y(a) = y

0

dz

dx

= f (x, y, z),

z(a) = y

0 0

(2.11)

したがって逐次近似式は











y

i+1

= y

i

+ g(x

i

, y

i

, z

i

)h

z

i+1

= z

i

+ f (x

i

.y

i

, z

i

)h

x

_i+1

= x

_i

+ h

(i = 0, 1, 2, · · · , n − 1)

(2.12)

それでは次の

2

階常微分方程式にチャレンジしましょう。

(1 − x

2

)

d

2

y

dx

2

− 2x

dy

dx

+ 6y = 0

  （初期条件：

y(0) = 1, y

0

_{(0) = 1)}

_(2.13)

この微分方程式の厳密解は

y = 1 − 3x

2

(2.14)

です。さて，

_













dy

dx

= z

dz

dx

=

2xz − 6y

1 − x

2

(2.15)

とおいて

1

階の連立微分方程式に変換すると逐次近似式は次のようになります6。















y

i+1

= y

i

+ g(x

i

, y

i

, z

i

),

g(x, y.z) = z

z

i+1

= z

i

+ f (x

i

.y

i

, z

i

)h,

f (x, y, z) =

2xz − 6y

_{1 − x}

₂

x

_i+1

= x

_i

+ h

(2.16)

  (※) エクセルの＜おまけ＞シートに次節で述べる 2 階常微分方程式の修正オイラー法を適用したものをのせておきました。 参考にしてください。なお，都合により変数名は変えてありますので留意のほどを（＾＾）；。 6_{f (x, y, z)} の分母1 − x2 に注意。h ∗ cntがx = 1を跨ぐと大変なことになります。xの範囲に注意！。

●

2

階常微分方程式のオイラー法・

VBA

ログラム   ＜ソースコード＞ Option Explicit ’ モジュール内で有効な変数の定義 Private cnt As Integer Private FL As Integer ’ グラフ消去のフラグ Function f(x, y, z) f = (2 * x * z - 6 * y) / (1 - x ^ 2) End Function Function g(x, y, z) g = z End Function ’ 厳密解 ---Function Sf(x) Sf = 1 - 3 * x ^ 2 End Function ’** 2 階常微分方程式 *********** Sub オイラー法 () ’******************************

Dim x0 As Single, y0 As Single, z0 As Single Dim x1 As Single, y1 As Single, z1 As Single Dim i As Integer, h As Single

’ 初期値 FL = 1

x0 = Range("C4").Value: y0 = Range("C5").Value: z0 = Range("C6").Value h = Range("E4").Value cnt = Range("E5").Value Cells(9, 1) = 0 Cells(9, 2) = x0 Cells(9, 3) = y0 Cells(9, 4) = Sf(x0) ’ 逐次計算 For i = 1 To cnt y1 = y0 + h * g(x0, y0, z0) z1 = z0 + h * f(x0, y0, z0) y0 = y1 z0 = z1 x0 = x0 + h Cells(i + 9, 1) = i Cells(i + 9, 2) = x0 Cells(i + 9, 3) = y0 Cells(i + 9, 4) = Sf(x0) ’ 厳密解 Cells(i + 9, 5) = y0 - Sf(x0) ’ 近似解と厳密解の差 Next i ’ 近似解曲線を散布図で描く With ActiveSheet.Shapes.AddChart.Chart .ChartType = xlXYScatterSmoothNoMarkers ’ 平滑線付き散布図データマーカーなし .SetSourceData Range(Cells(8, 2), Cells(8 + cnt, 4)) ’ 作図範囲

End With End Sub

’ 計算結果の消去 Sub Erase_Click()

Range(Cells(9, 1), Cells(9 + cnt, 5)).Clear ’ セルに格納したデータを消去 If FL = 1 Then ActiveSheet.ChartObjects("1").Activate ’ グラフの消去 Selection.Delete FL = 0 End If End Sub

2.2.1

２階常微分方程式・修正オイラー法 次の

2

階常微分方程式をオイラー法と修正オイラー法で計算し，厳密解と比較します。

d

2

_x

dt

2

= −x,

(x(0) = 1, x

0

(0) = 0)

／厳密解：

x = cos t

(2.17)

 さて，例によって

1

階常微分連立方程式に変換します。















dx

dt

= v

dv

dt

= −x

(2.18)

刻み幅を

h = 0.05

とした場合の計算結果を示します。修正オイラー法と厳密解の曲線はほとんど一致しますが， オイラー法の曲線はズレが次第に大きくなっていきます。オイラー法でも

h

を細かくとることでズレは抑えられ ますので確認してください。 ●

2

階常微分方程式の修正オイラー法・

VBA

ログラム   ＜ソースコード＞ Option Explicit ’ モジュール内で有効な変数の定義 Private cnt As Integer Private FL As Integer ’ グラフ消去のフラグ Function f(x) f = -x End Function ’**☆ 2 階オイラー法 *********** Sub オイラー法 ()

Dim t0 As Single, x0 As Single, v0 As Single Dim x As Single, v As Single, t As Single Dim x1 As Single, v1 As Single

Dim f0 As Single, f1 As Single Dim k1 As Single, k2 As Single Dim h As Single, i As Integer ’ 初期値

t0 = Range("C4").Value x0 = Range("C5").Value v0 = Range("C6").Value h = Range("E4").Value cnt = Range("E5").Value Cells(9, 1) = 0 Cells(9, 2) = t0 Cells(9, 3) = x0 Cells(9, 4) = x0 ’ ←修正オイラー用 ’ オイラー法逐次計算---For i = 1 To cnt x1 = x0 + h * v0 v1 = v0 + h * f(x0) t = t + h Cells(i + 9, 1) = i Cells(i + 9, 2) = t Cells(i + 9, 3) = x1 x0 = x1 v0 = v1 Next i ’**** 修正オイラー法逐次計算 **** t0 = Range("C4").Value x0 = Range("C5").Value v0 = Range("C6").Value x = x0 v = v0 t = t0 For i = 1 To cnt v0 = v f0 = f(x) x1 = x + h * v0 v1 = v + h * f0 f1 = f(x1) k1 = (v0 + v1) * h k2 = (f0 + f1) * h t = t + h x = x + k1 / 2 v = v + k2 / 2 Cells(i + 9, 4) = x Cells(i + 9, 5) = Cos(t) ’ 厳密解 Next i ’ 近似解曲線を散布図で描く With ActiveSheet.Shapes.AddChart.Chart .ChartType = xlXYScatterSmoothNoMarkers ’ 平滑線付き散布図データマーカーなし .SetSourceData Range(Cells(8, 2), Cells(8 + cnt, 5)) ’ 作図範囲

End With End Sub

’ 計算結果の消去 Sub Erase_Click()

Range(Cells(9, 1), Cells(9 + cnt, 5)).Clear ’ セルに格納したデータを消去 If FL = 1 Then ActiveSheet.ChartObjects("1").Activate ’ グラフの消去 Selection.Delete FL = 0 End If End Sub

2.3

ルンゲ・クッタ法

xi+1 x yi h/2 k1

➀

➁

_k2

➃

k3 h

➁

➂

➂

k4

➃

y xi k1/2 図

3: 4

次ルンゲ・クッタ法   近似精度が高く通常よく使われる

4

次のルンゲ・クッタ7法 を使って

1

階，

2

階，

3

階の常微分方程式を解いていきます。 なお，ルンゲ・クッタ公式の解析的な導出は相良紘「技術 者のための数値計算入門」（日刊工業新聞社

,2007

）に詳し く載っていますので興味がある方は一読ください。

2.3.1

1

階常微分方程式 次の

1

階常微分方程式を考えます。

dy

dx

= f (x, y)

(2.19)

 初期条件

(x

0

, y

0

)

からスタートして解の近似曲線の座標

(x

_i+1

, y

_i+1

) (i = 0, 1, 2, · · · , n)

を順次計算して求めていく わけですが，ルンゲ・クッタ法では近似精度が高くなる見 返りとして手順のステップは増えます。  刻み幅を

h

として具体的な手順は次の

5

ステップとなります。

1.

オイラー法と同じように，点

(x

i

, y

i

)

のおける接線を引いて隣接点

x

i+1での

y

値の増分

k

1を求めます。

k

1

= f (x

i

, y

i

)h

(2.20)

2. k

1を使って

x

iから

h/2

進んだ点

(x

i

+ h/2, y

i

+ k

1

/2)

での傾き

f (x

i

+ h/2, y

i

+ k

1

/2)

をもつ接線を点

(x

i

, y

i

)

を通るように引き直し，

x

iの隣接点

x

i+1での

y

値の増分を

k

2とします。

k

2

= f

µ

x

i

+

h

₂

, y

i

+

k

₂

1

¶

h

(2.21)

3.

ステップ

2

と同様に，

k

2を使って

x

iと

x

i+1の中点での傾き

f (x

i

+ h/2, y

i

+ k

2

/2)

をもつ接線を点

(x

i

, y

i

)

を通るように引き直し，

x

iの隣接点

x

i+1での

y

値の増分を

k

3とします。

k

3

= f

µ

x

i

+

h

₂

, y

i

+

k

₂

2

¶

h

(2.22)

4. k

3を使って

x

iから

h

進んだ点

(x

i+1

, y

i

+ k

3

)

での傾き

f (x

i+1

, y

i

+ k

3

)

を持つ接線を点

(x

i

, y

i

)

を通るよう に引き直し，

x

iの隣接点

x

i+1での

y

値の増分を

k

4とします。

k

4

= f (x

i

+ h, y

i

+ k

3

) h

(2.23)

5.

各ステップで増分

k

1

, k

2

, k

3

, k

4にぞれぞれ

1, 2, 2, 1

の重み8を付けて加重平均した値

1

6

(k

1

+ 2k

2

+ 2k

3

+ k

4

)

(2.24)

を

x

iと

x

i+1間の

y

値の最終的な増分とすると，

x

iの隣接点

x

i+1

(= x

i

+ h)

での

y

i+1値は

y

i+1

= y

i

+

1

₆

(k

1

+ 2k

2

+ 2k

3

+ k

4

) (i = 0, 1, 2, · · · , n − 1)

(2.25)

 それでは次の

1

階常微分方程式の近似解をルンゲ・クッタ法で求めてみましょう。

dy

dx

=

1

2

x(1 − y)

（初期値：

y(0) = 4)

(2.26)

この微分方程式の厳密解は

y = 1 + 3e

−x2/4です。刻み幅を

h = 0.1

とした結果を以下に示します。

7_{Carl David Tolme Runge}

：1856-1927.ドイツの数学者,物理学者,分光学者。Martin Wilhelm Kutta：1867-1944.ドイツの数学者。

8

●ルンゲ・クッタ法

1

階常微分方程式の

VBA

ログラム   ＜ソースコード＞ Option Explicit Private FL As Integer ’ グラフ消去のフラグ Private cnt As Integer Function f(x, y) f = x * (1 - y) / 2 End Function ’**ルンゲ・クッタ法 *********** Sub ルンゲクッタ法 ()

Dim x0 As Single, y0 As Single, y As Single Dim i As Integer, h As Single

Dim k1 As Single, k2 As Single, k3 As Single, k4 As Single ’ 初期値 FL = 1 x0 = Range("C4").Value y0 = Range("C5").Value h = Range("E4").Value cnt = Range("E5").Value ’ ルンゲクッタ法逐次計算 For i = 0 To cnt k1 = h * f(x0, y0) k2 = h * f(x0 + h / 2, y0 + k1 / 2) k3 = h * f(x0 + h / 2, y0 + k2 / 2) k4 = h * f(x0 + h, y0 + k3) y = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6 Cells(i + 8, 1) = i Cells(i + 8, 2) = x0 Cells(i + 8, 3) = y0 Cells(i + 8, 4) = 1 + 3 * Exp(-x0 ^ 2 / 4) ’ 厳密解 x0 = x0 + h ’x 値を更新 y0 = y ’y 値を更新 Next i ’ 近似解曲線を散布図で描く With ActiveSheet.Shapes.AddChart.Chart .ChartType = xlXYScatterSmoothNoMarkers ’ 平滑線付き散布図データマーカーなし .SetSourceData Range(Cells(8, 2), Cells(8 + cnt, 4)) ’ 作図範囲

End With End Sub

’ 計算結果の消去 Sub Erase_Click()

Range(Cells(8, 1), Cells(8 + cnt, 4)).Clear ’ セルに格納したデータを消去 If FL = 1 Then ActiveSheet.ChartObjects("1").Activate ’ グラフの消去 Selection.Delete FL = 0 End If End Sub

2.3.2

2

階常微分方程式

2

階常微分方程式

d

2

_y

dx

2

= f

µ

x, y,

dy

dx

¶

,

(

初期値：

y(x

0

) = y

0

, y

0

(x

0

) = y

00）

(2.27)

の近似解をルンゲクッタ法で求めます。例によって

dy/dx = z

とおいて

1

階連立常微分方程式に変換します。















dy

dx

= z = f (x, y, z)

dz

dx

= g(x, y, z)

（初期条件：

x = x

0のとき

y = y

0

, z = z

0）

(2.28)

上の

2

式にそれぞれルンゲクッタ法を適用すると各ステップの増分は

k

1

= f (x

i

, y

i

, z

i

)h,

j

1

= g(x

i

, y

i

, z

i

)h

k

2

= f

µ

x

i

+

h

₂

, y

i

+

k

₂

1

, z

i

+

j

₂

1

¶

h,

j

2

= g

µ

x

i

+

h

₂

, y

i

+

k

₂

1

, z

i

+

j

₂

1

¶

k

3

= f

µ

x

i

+

h

2

, y

i

+

k

₂

2

, z

i

+

j

₂

2

¶

h,

j

3

= g

µ

x

i

+

h

2

, y

i

+

k

₂

2

, z

i

+

j

₂

2

¶

k

4

= f (x

i

+ h, y

i

+ k

3

, z

i

+ j

3

)h,

j

4

= g (x

i

+ h, y

i

+ k

3

, z

i

+ j

3

) h

(2.29)

となり，初期値

(x = x

0

, y = y

0

, z = z

0

)

から出発して

x

の刻み幅

h

に対する

y

および

z

の増分は

y

i+1

= y

i

+

1

₆

(k

1

+ 2k

2

+ 2k

3

+ k

4

)

z

i+1

= z

i

+

1

₆

(j

1

+ 2j

2

+ 2j

3

+ j

4

)

(i = 0, 1, 2, · · · , n − 1)

(2.30)

から計算されます。  それでは次の

2

階常微分方程式の近似解をルンゲ・クッタ法で求めましょう。

d

2

_y

dx

2

− 2

dy

dx

+ y = 0 (

初期値：

x = 0

のとき

y = 1, dy/dx = 0)

(2.31)

この方程式の厳密解は

y = (1 − x)e

xです。

dy/dx = z

とおいて次の

1

階連立微分方程式に変換します。















dy

dx

= z = f (x, y, z)

dz

dx

= 2z − y = g(x, y, z)

(x

0

= 0, y

0

= 1, z

0

= 0)

(2.32)

x

の刻み幅を

h

とすると，各ステップの増分と逐次近似式は

(2.29), (2.30)

で与えられます。数値計算の結果（一 部）は次のとおりです。

●ルンゲ・クッタ法

2

階常微分方程式の

VBA

ログラム   ＜ソースコード＞ Option Explicit Private FL As Integer ’ グラフ消去のフラグ Private cnt As Integer Function f(x, y, z) f = z End Function Function g(x, y, z) g = 2 * z - y End Function ’**ルンゲ・クッタ法（2 階常微分方程式）*********** Sub ルンゲクッタ法 ()

Dim x0 As Single, y0 As Single, z0 As Single Dim y As Single, z As Single

Dim i As Integer, h As Single

Dim k1 As Single, k2 As Single, k3 As Single, k4 As Single Dim j1 As Single, j2 As Single, j3 As Single, j4 As Single ’ 初期値 FL = 1 x0 = Range("C4").Value y0 = Range("C5").Value z0 = Range("C6").Value h = Range("E4").Value cnt = Range("E5").Value ’ ルンゲクッタ法逐次計算 For i = 0 To cnt k1 = h * f(x0, y0, z0) j1 = h * g(x0, y0, z0) k2 = h * f(x0 + h / 2, y0 + k1 / 2, z0 + j1 / 2) j2 = h * g(x0 + h / 2, y0 + k1 / 2, z0 + j1 / 2) k3 = h * f(x0 + h / 2, y0 + k2 / 2, z0 + j2 / 2) j3 = h * g(x0 + h / 2, y0 + k2 / 2, z0 + j2 / 2) k4 = h * f(x0 + h, y0 + k3, z0 + j3) j4 = h * g(x0 + h, y0 + k3, z0 + j3) y = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6 z = z0 + (j1 + 2 * j2 + 2 * j3 + j4) / 6 Cells(i + 9, 1) = i Cells(i + 9, 2) = x0 Cells(i + 9, 3) = y0 Cells(i + 9, 4) = (1 - x0) * Exp(x0) ’ 厳密解 x0 = x0 + h ’x 値の更新 y0 = y ’y 値の更新 z0 = z ’z 値の更新 Next i ’ 近似解曲線を散布図で描く With ActiveSheet.Shapes.AddChart.Chart .ChartType = xlXYScatterSmoothNoMarkers ’ 平滑線付き散布図データマーカーなし .SetSourceData Range(Cells(8, 2), Cells(8 + cnt, 4)) ’ 作図範囲

End With End Sub

’ 計算結果の消去 Sub Erase_Click()

Range(Cells(9, 1), Cells(9 + cnt, 4)).Clear ’ セルに格納したデータを消去 If FL = 1 Then ActiveSheet.ChartObjects("1").Activate ’ グラフの消去 Selection.Delete FL = 0 End If End Sub

ファン・デル・ポールの方程式 非線形の減衰を受けた非保存系の振動子をファン・デル・ポール振動子といい， この力学系を記述するのがファン・デル・ポール9の方程式です。

d

2

x

dt

2

+ µ(x

2

− 1)

dx

dt

+ x = 0

(2.33)

左辺第

2

項の

µ

は非線形の減衰の強さを表すパラメータです。

4

次のルンゲ・クッタ法で解くために

dx/dt = y

と おいて

1

階の連立微分方程式に変換すると















dx

dt

= y = f (t, x, y)

dy

dt

= µ(1 − x

2

_{) − x = g(t, x, y)}

 （初期条件：

x(0) = 0.5, dx/dt|

_t=0

= y(0) = 0

）

(2.34)

あとは処方箋

(2.29), (2.30)

に従って計算するだけです。 ●

VanderPole

方程式の

VBA

ログラム   ＜ソースコード＞ Option Explicit Private FL As Integer ’ グラフ消去のフラグ Private cnt As Integer ’ カウンター Private m As Single ’ 減衰係数 Function f(t, x, y) f = y End Function Function g(t, x, y) g = m * y * (1 - x ^ 2) - x End Function ’**ルンゲ・クッタ法（ファン・デル・ポール方程式）*********** Sub ルンゲクッタ法 ()

Dim t0 As Single, x0 As Single, y0 As Single Dim x As Single, y As Single

Dim i As Integer, h As Single

Dim k1 As Single, k2 As Single, k3 As Single, k4 As Single Dim j1 As Single, j2 As Single, j3 As Single, j4 As Single ’ 初期値 m = Range("C6").Value ’ 減衰係数 FL = 1 x0 = Range("C4").Value y0 = Range("C5").Value t0 = 0 h = Range("E4").Value

9_{Balthasar van der Pol :1889-1959}

cnt = Range("E5").Value ’ ルンゲクッタ法逐次計算 For i = 0 To cnt k1 = h * f(t0, x0, y0) j1 = h * g(t0, x0, y0) k2 = h * f(t0 + h / 2, x0 + k1 / 2, y0 + j1 / 2) j2 = h * g(t0 + h / 2, x0 + k1 / 2, y0 + j1 / 2) k3 = h * f(t0 + h / 2, x0 + k2 / 2, y0 + j2 / 2) j3 = h * g(t0 + h / 2, x0 + k2 / 2, y0 + j2 / 2) k4 = h * f(t0 + h, x0 + k3, y0 + j3) j4 = h * g(t0 + h, x0 + k3, y0 + j3) x = x0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6 y = y0 + (j1 + 2 * j2 + 2 * j3 + j4) / 6 Cells(i + 9, 1) = i Cells(i + 9, 2) = t0 Cells(i + 9, 3) = x0 t0 = t0 + h x0 = x y0 = y Next i ’ 近似解曲線を散布図で描く With ActiveSheet.Shapes.AddChart.Chart .ChartType = xlXYScatterSmoothNoMarkers ’ 平滑線付き散布図データマーカーなし .SetSourceData Range(Cells(8, 2), Cells(8 + cnt, 3)) ’ 作図範囲

End With End Sub

’ 計算結果の消去 Sub Erase_Click()

Range(Cells(9, 1), Cells(9 + cnt, 3)).Clear ’ セルに格納したデータを消去 If FL = 1 Then ActiveSheet.ChartObjects("1").Activate ’ グラフの消去 Selection.Delete FL = 0 End If End Sub ロトカ・ヴォルテラ方程式

2

階常微分方程式を解くのに

1

階連立常微分方程式に変換し，ルンゲ・クッタ法で その解き方を見てきました。そこで次に最初から

1

階の連立常微分方程式で記述される系の振る舞いをルンゲ・ クッタ法で解いてみましょう。生物の捕食者と被食者（食う

/

食われる）の関係による個体数の変動を表す数理モ デルとして有名なロトカ・ヴォルテラ10の方程式を取り上げます。















dx

dt

= ax − cxy

dy

dt

= −by + dxy

(2.35)

第

1

式は被食者の個体数増殖速度，第

2

式は捕食者の個体数増殖速度を表しています。

a, b, c, d

は実数のパラメー ターで

₍

a

：

x

の繁殖率

,

b

：自然減少，死亡率

c : y

の捕食率

,

d

：

y

の繁殖率

(2.36)

ここでは

a = 1, b = 0.5, c = 0.1, d = 0.02

とし，初期条件は

x(0) = 100, y(0) = 7

と設定しました。 ●

Lotka-Volterra

方程式の

VBA

ログラム   ＜ソースコード＞ Option Explicit Private FL As Integer ’ グラフ消去のフラグ Function f(a, c, t, x, y)

10_{Alfred James Lotka}

f = a * x - c * x * y End Function Function g(b, d, t, x, y) g = -b * y + d * x * y End Function ’**ルンゲ・クッタ法（ロトカ・ヴォルテラ方程式）*********** Sub ルンゲクッタ法 ()

Dim t0 As Single, x0 As Single, y0 As Single Dim t As Single, x As Single, y As Single Dim i As Integer, cnt As Integer, h As Single

Dim k1 As Single, k2 As Single, k3 As Single, k4 As Single Dim j1 As Single, j2 As Single, j3 As Single, j4 As Single Dim a As Single, b As Single, c As Single, d As Single ’ パラメータの値 a = Range("H4").Value b = Range("H5").Value c = Range("H6").Value d = Range("H7").Value ’ 初期値 FL = 1 t0 = Range("C4").Value x0 = Range("C5").Value y0 = Range("C6").Value h = Range("E4").Value cnt = Range("E5").Value ’ ルンゲクッタ法逐次計算 For i = 0 To cnt k1 = h * f(a, c, t0, x0, y0) j1 = h * g(b, d, t0, x0, y0) k2 = h * f(a, c, t0 + h / 2, x0 + k1 / 2, y0 + j1 / 2) j2 = h * g(b, d, t0 + h / 2, x0 + k1 / 2, y0 + j1 / 2) k3 = h * f(a, c, t0 + h / 2, x0 + k2 / 2, y0 + j2 / 2) j3 = h * g(b, d, t0 + h / 2, x0 + k2 / 2, y0 + j2 / 2) k4 = h * f(a, c, t0 + h, x0 + k3, y0 + j3) j4 = h * g(b, d, t0 + h, x0 + k3, y0 + j3) x = x0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6 y = y0 + (j1 + 2 * j2 + 2 * j3 + j4) / 6 Cells(i + 9, 1) = i Cells(i + 9, 2) = t0 Cells(i + 9, 3) = x0 Cells(i + 9, 4) = y0 t0 = t0 + h x0 = x y0 = y Next i ’ 近似解曲線を散布図で描く

With ActiveSheet.Shapes.AddChart.Chart

.ChartType = xlXYScatterSmoothNoMarkers ’ 平滑線付き散布図データマーカーなし .SetSourceData Range(Cells(8, 2), Cells(8 + cnt, 4)) ’ 作図範囲

End With End Sub ’ 計算結果の消去 Sub Erase_Click() Dim n As Integer n = Range("E5").Value

Range(Cells(9, 1), Cells(9 + n, 4)).Clear ’ セルに格納したデータを消去 If FL = 1 Then ActiveSheet.ChartObjects("1").Activate ’ グラフの消去 Selection.Delete FL = 0 End If End Sub

2.3.3

3

階常微分方程式

3

階常微分方程式を取り上げます。これより高階常微分方程式への拡張も容易にできると思います。

d

3

_y

dx

3

= f

µ

x, y,

dy

dx

,

d

2

_y

dx

2

¶

,

（初期条件：

y(x

0

) = a

0

, y

0

(x

0

) = a

1

, y

00

(x

0

) = a

2

)

(2.37)

ルンゲ・クッタ法を適用するために

1

階の連立微分方程式に変換します。

dy

dx

= z,

d

2

_y

dx

2

=

dx

dx

= u

(2.38)

とおけば，

(2.37)

は次の

1

階連立常微分方程式に変換できます。

dy

dx

= z,

dz

dx

= u,

du

dx

= f (x, y, z, u)

(2.39)

計算のアルゴリズムを見通しよくする都合上，上式を次のように書き換えておきます。

dy

dx

= p(x, y, z, u),

dz

dx

= q(x, y, z, u),

du

dx

= r(x, y, z, u)

(2.40)

刻み幅を

h

とすると

(2.29)

でやったように各ステップの増分は























































































































































k

1

= p(x

i

, y

i

, z

i

, u

i

)h

j

1

= q(x

i

, y

i

, z

i

, u

i

)h

`

₁

= r(x

_i

, y

_i

, z

_i

, u

_i

)h

k

2

= p

µ

x

i

+

h

₂

, y

i

+

k

₂

1

, z

i

+

j

₂

1

, u

i

+

`

₂

1

¶

h

j

2

= q

µ

x

i

+

h

₂

, y

i

+

k

₂

1

, z

i

+

j

₂

1

, u

i

+

`

₂

1

¶

h

`

₂

= r

µ

x

_i

+

h

2

, y

i

+

k

1

2

, z

i

+

j

1

2

, u

i

+

`

1

2

¶

h

k

3

= p

µ

x

i

+

h

₂

, y

i

+

k

₂

2

, z

i

+

j

₂

2

, u

i

+

`

₂

2

¶

h

j

3

= q

µ

x

i

+

h

₂

, y

i

+

k

₂

2

, z

i

+

j

₂

2

, u

i

+

`

₂

2

¶

h

`

3

= r

µ

x

i

+

h

₂

, y

i

+

k

₂

2

, z

i

+

j

₂

2

, u

i

+

`

₂

2

¶

h

k

₄

= p(x

_i

+ h, y

_i

+ k

₃

, z

_i

+ j

₃

, u

_i

+ `

₃

)h

j

₄

= q(x

_i

+ h, y

_i

+ k

₃

, z

_i

+ j

₃

, u

_i

+ `

₃

)h

`

4

= r(x

i

+ h, y

i

+ k

3

, z

i

+ j

3

, u

i

+ `

3

)h

(i = 0, 1, 2, · · · , n − 1)

(2.41)

となり，各増分に重み付けして加重平均した値を最終増分とすれば

y

_i+1

= y

_i

+

1

6

(k

1

+ 2k

2

+ 2k

3

+ k

4

)

z

i+1

= z

i

+

1

₆

(j

1

+ 2j

2

+ 2j

2

+ j

4

)

u

_i+1

= u

₁

+

1

6

(`

1

+ 2`

2

+ 2`

2

+ `

4

)

(2.42)

が得られます。  それでは次の

3

階常微分方程式を解いてみましょう。

d

3

_y

dx

3

−

d

2

_y

dx

2

−

dy

dx

+ y = e

x _{（初期値：}

_{y(0) = −1, y}

0

_{(0) = 2, y}

00

_{(0) = 1.5}

_）

_(2.43)

(2.43)

の厳密解は

y = (0.25x

2

e

2x

+ xe

2x

− 1)e

−xです。

dy/dx = z, dz/dx = u

とおくと

(2.43)

は次の

1

階連立常 微分方程式に変換されます。

_





























dy

dx

= z = p(x, y, z, u)

dz

dx

= u = q(x, y, z, u)

du

dx

= e

x

_{+ u + z − y = r(x, y, z, u)}

(2.44)

刻み幅

h = 0.1

とした場合の計算結果を下図に示します。厳密解と高い精度で一致しています。 ●ルンゲ・クッタ法

3

階常微分方程式の

VBA

ログラム   ＜ソースコード＞ Option Explicit Private FL As Integer ’ グラフ消去のフラグ Private cnt As Integer Function p(x, y, z, u) p = z End Function Function q(x, y, z, u) q = u End Function Function r(x, y, z, u) r = Exp(x) + u + z - y End Function ’**ルンゲ・クッタ法（3 階常微分方程式）*********** Sub ルンゲクッタ法 ()

Dim x0 As Single, y0 As Single, z0 As Single, u0 As Single Dim x As Single, y As Single, z As Single, u As Single Dim i As Integer, h As Single

Dim j1 As Single, j2 As Single, j3 As Single, j4 As Single Dim l1 As Single, l2 As Single, l3 As Single, l4 As Single ’ 初期値 FL = 1 x0 = Range("C4").Value y0 = Range("C5").Value z0 = Range("C6").Value u0 = Range("C7").Value h = Range("E4").Value cnt = Range("E5").Value ’ ルンゲクッタ法逐次計算 For i = 0 To cnt k1 = h * p(x0, y0, z0, u0) j1 = h * q(x0, y0, z0, u0) l1 = h * r(x0, y0, z0, u0) k2 = h * p(x0 + h / 2, y0 + k1 / 2, z0 + j1 / 2, u0 + l1 / 2) j2 = h * q(x0 + h / 2, y0 + k1 / 2, z0 + j1 / 2, u0 + l1 / 2) l2 = h * r(x0 + h / 2, y0 + k1 / 2, z0 + j1 / 2, u0 + l1 / 2) k3 = h * p(x0 + h / 2, y0 + k2 / 2, z0 + j2 / 2, u0 + l2 / 2) j3 = h * q(x0 + h / 2, y0 + k2 / 2, z0 + j2 / 2, u0 + l2 / 2) l3 = h * r(x0 + h / 2, y0 + k2 / 2, z0 + j2 / 2, u0 + l2 / 2) k4 = h * p(x0 + h, y0 + k3, z0 + j3, u0 + l3) j4 = h * q(x0 + h, y0 + k3, z0 + j3, u0 + l3) l4 = h * r(x0 + h, y0 + k3, z0 + j3, u0 + l3) y = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6 z = z0 + (j1 + 2 * j2 + 2 * j3 + j4) / 6 u = u0 + (l1 + 2 * l2 + 2 * l3 + l4) / 6 Cells(i + 10, 1) = i Cells(i + 10, 2) = x0 Cells(i + 10, 3) = y0

Cells(i + 10, 4) = (0.25 * x0 ^ 2 * Exp(2 * x0) + x0 * Exp(2 * x0) - 1) * Exp(-x0) x0 = x0 + h y0 = y z0 = z u0 = u Next i ’ 近似解曲線を散布図で描く With ActiveSheet.Shapes.AddChart.Chart .ChartType = xlXYScatterSmoothNoMarkers ’ 平滑線付き散布図データマーカーなし .SetSourceData Range(Cells(9, 2), Cells(9 + cnt, 4)) ’ 作図範囲

End With End Sub

’ 計算結果の消去 Sub Erase_Click()

Range(Cells(10, 1), Cells(10 + cnt, 4)).Clear ’ セルに格納したデータを消去 If FL = 1 Then ActiveSheet.ChartObjects("1").Activate ’ グラフの消去 Selection.Delete FL = 0 End If End Sub