9. 微分方程式の数値解法（ルンゲ・クッタ法） 9.1.

(1)

9. 微分方程式の数値解法（ルンゲ・クッタ法）

9.1.

背景

全微分についての方程式が与えられたときに関数の形を求める問題を微分方程式という．微分の階数を頭につけて，「

n

階常微分方程式」と呼ばれることが多い．１階常微分方程式の形は以下の通りである．

) , ( t x f dt dx



(1)

これは，

x(t)

^の変化率

dx / dt

^が

f (t, x(t))

で与えられているときに

x(t)

^{を求める問題であ}

る．この解の一般解には積分定数がつく．特殊解を求めるためには，

t = 0

のときの

x(t )

^の値

である

x(0)

を与えなければならない．

同様にして，

n

階常微分方程式の形は以下の通りである．











  _

 1 1

2 2

, , , ,

,

n

n n

dt x d dt

x d dt dx x t f dt

x

d



これは，

x(t )

^の

^n階微分 d

ⁿ

x / dt

ⁿ ^が

f

で与えられているときに

x(t)

^{を求める問題である．}

ｎ階常特殊解を求めるためには，

t = 0

^のときの

x(t)

^{の値である}

x(0)

^{の他に}

x ¢ (0) ^,

¢¢

x (0) ^,…, x

^(n-1)

(0)

を与えなければならない．微分方程式は n 個の

1

階常微分方程式に分割で

きるので，1 階連立常微分方程式の解を求める方法がわかれば，

n

階常微分方程式の解を求めることができると考えてよいだろう．

実際問題としてみると，常微分方程式は解析解が求められる場合もあるし求められない場合もある．解析的に求めることができない微分方程式については，コンピュータを用いて数値解を求める．数値解とは

x(t)

^{定義域（すなわち}

t

の範囲）を適当に区切って（離散化して），それぞれの時刻における

x(t )

の値を求めることである．ここでは，この数値解の求め方を確認する．

この講義で紹介する数値解法は，差分法と Runge-Kutta 法である．差分法は単純なアルゴリズムであるが精度は悪い． Runge-Kutta 法は精度は良いがアルゴリズムは多少複雑になる．以下では定義域

0 £ t £ t

_max ^を

m

個の微小領域に等分割することを考える．このときの微小領域の幅を

Dt(= t

_max

/ m)

とおく．各種の数値解法は

t

^が

_D t

^{だけ進んだときに}

x(t + Dt)

^が

x(t )

よりもどれくらい大きくなるのかを見積もる方法を与えている．

9.2.

差分法のアルゴリズム

差分法は１階常微分方程式

f ( t , x ) dt

dx



(1)を Taylor

展開を用いて以下のように近似するものである

(2)

コンピュータプログラミング

9.3. Runge-Kutta

法のアルゴリズム

差分法は１階常微分方程式(1)を精度よくかつ効率よく計算する方法である． Runge-Kutta 法には幾つかのバリエーションがある．以下の計算式によって解く方法を４次の Runge-Kutta 法という．この方法は広く数値計算に用いられている方法である．

 

) ) ( , (

), 2 / ) ( , 2 / (

), 2 / ) ( , 2 / ( )),

( , (

2 2 6 1 ) ( ) (

3 4

2 3

1 2

1

4 3 2 1

k t x t t f t k

k t x t t f t k t x t f t k

k k k k t x t t x



















































9.4.

例（ロジスティック方程式）

ここでは，数値計算を行う例としてロジスティック方程式を考える．これは，成長速度が現在の値と現在の値の最大の値からの差の積に比例するというモデルを表している．

)

( N x

x r dt

dx

 

ここで，と

はそれぞれ定数である．

は

の成長の速さを決める定数であり，

は

が最終的に到達する値を表す定数である．

ここではロジスティック方程式の背景については詳しい解説を行わないが，その解は成長曲線と呼ばれ，さまざまなところに現れていることは覚えておくと良い．

9.5.

解析解のグラフを描く

はじめに，前回同様の方法でグラフを描き，解析解の形を確認しておく（これは例題だからできることであって，解析解がわからない問題では当然不可能）．この解と数値解とを較べて，数値解を計算する上で気になることを確認していくことにする．

 N x  e

^rNt

N

t

x   



1 1

) (

0

(3)

9.6.

差分法による実装

何通りか計算できるようにアクティブセルを基準として結果を示すようにする．今回は

f (t, x) = r x (N - x)

であることに注目する．

(4)

マクロが正しく動いているようであれば deltaT の値を変更して幾つかの場合を計算して比較を行うこと．すると，次のようなグラフが完成するはずである．

Sub fdLogistic()

Const N As Integer = 100 Const r As Double = 0.02 Const initialX As Double = 2

Const deltaT As Double = 0.05, maxT As Double = 6 Dim x As Double, newX As Double, t As Double Dim counter As Integer

ActiveCell.Value = "年"

ActiveCell.Offset(0, 1).Value = "人口(" & "DT=" & deltaT & ")"

ActiveCell.Offset(1, 0).Value = 0

ActiveCell.Offset(1, 1).Value = initialX x = initialX

counter = 2

For t = deltaT To maxT Step deltaT newX = x + deltaT * r * x * (N - x) ActiveCell.Offset(counter).Value = t

ActiveCell.Offset(counter, 1).Value = newX x = newX

counter = counter + 1 Next t

End Sub

(5)

DT

の値を小さくしたほうが解析解に近くなること，すなわち精度が良くなることがわかる．

9.7. Runge-Kutta

法による実装

差分法と同様に Runge-Kutta 法についてもチャレンジする．敢えて差分でエラーが出るような分割を採用してみる．

(6)

9. 微分方程式の数値解法（ルンゲ・クッタ法） 9.1.