数値シミュレーションの基礎オイラー法とルンゲ・クッタ法

(1)

数値シミュレーションの基礎オイラー法とルンゲ・クッタ法

黒田紘敏（理学研究院数学部門）

数理解析学特論A /フロンティア数理物質科学II 第7回

2020年6月25日

数理解析学特論A (第回) 数値シミュレーションの基礎年月日

(2)

Introduction

今日は常微分方程式について扱う．有名な微分方程式は，ニュートンの運動方程式

md²x(t)

d t² = F(x(t)) などがある．

常微分方程式の解の存在を示すことは一般的に難しくないが，その具体的な解表示を得ることは難しいことが多い．

(It is not so difficult to show the existence of a solution to an ordinary differential equation. However, it is generally more difficult to express the solutions specifically.)

(Example) 例えば，次の微分方程式

x^′(t)= x(t)+e^x(t)sinx(t), x(0)= 1. の解を x(t) =· · · の形で具体的に示すのは難しい（と思う）．

(3)

Contents

Aim

解の具体的な表示が得られないときに，PCによるシミュレーションでそのグラフの概形を描く方法を学ぶ．そのための数値計算の基本的な考え方と，オイラー法や自然科学分野で広く用いられているルンゲ・クッタ法と呼ばれる有名な計算法を理解する．

(When the solution cannot be shown speciﬁcally, learn how to draw the outline of the graph by numerical simulation on a PC. To understand the basic concepts of numerical computation and the famous Euler’s method and the Runge-Kutta method, which is widely used in the natural

sciences.)

1 Introduction

2 常微分方程式の解の一意性(Uniqueness of solutions to ODEs)

3 数値計算のアイデア(Idea of Numerical calculations)

4 オイラー法(Euler’s method)

5 ルンゲ・クッタ法(Runge–Kutta method)

(4)

ODE

今日は次のような形の常微分方程式(ODE)を扱う．

dx(t)

d t = f(t,x(t)) (t >0), x(0) = a

ここで，f(t,x)は既知の関数(given function)，実数 aは初期値(initial

data)である．変数を省略して

dx

d t = f(t,x) or x^′ = f(t,x) のようにも表す．

(Example)

f(t,x) = x =⇒ dx(t) d t = x(t) f(t,x) =−t x² =⇒ dx(t)

d t =−t x(t)²

(5)

Review

Example 1

Find the solution to the following ODE: dx

d t = x, x(0)= 1.

Answer

両辺を xで割ってから積分すれば

∫ 1

x dx d t d t =

∫

d t =⇒ log|x|= t+C

となる（Cは積分定数）．これを式変形すれば

|x|= e^t⁺^C =⇒ x(t)= ±e^Ce^t = Ae^t (A=±e^C)

となる．ここで，x(0) = A=1なので，求める解はx(t)= e^t となる．

(Remark) この解は時間大域的に t ≧0の範囲で存在する．

(This is the time global solution on[0,∞).)

(6)

Review

Example 2

d t =2t x², x(0)=1.

Answer

両辺を x²で割ってから積分すれば

∫ 1

x² dx d t d t =

∫

2t d t =⇒ −1

x = t²+C =⇒ x(t) =− 1 t²+C となる（Cは積分定数）．ここで，初期条件より定数を決定すれば

x(0) =−1

C =1 =⇒ C =−1 ∴ x(t) = 1 1− t² (Remark) この解は0≤ t< 1の範囲で存在する（解の爆発現象）．

(This solution blows up in ﬁnite time.)

(7)

Review

Problem 1

d t =−2t x, x(0)= 1.

(8)

Review (Appendix)

Example 3

Find the solution to the following ODE: dx d t = √

x, x(0)= 0.

Answer

x >0として，両辺を √

xで割ってから積分すれば

∫ 1

√x dx

d t d t =

∫

d t =⇒ 2√

x = t+C =⇒ x(t) = (t+C)² 4

となる．ここで x(0)=C²/4 =0より，x(t) = t²/4となる．

一方，x(t)= 0 (t≧ 0)は明らかに方程式をみたし，これも解である．

(Remark) 一般に無条件では解の一意性は成り立たない．

(In general, the uniqueness of the solutions does not hold.)

(9)

Idea of Numerical calculations（数値シミュレーションの考え方）

dx(t)

d t = f(t,x(t)) (t >0), x(0) = a

Motivation

微分計算を四則演算で置き換える．

(Replacing differential calculations with four arithmetic operations.)

Reason

PCを使えば四則演算や反復計算については容易に実行できる．

(Four arithmetic operations and iterations can be easily performed in numerical simulations.)

h >0が十分小さければ，導関数の定義より dx(t)

d t = lim

h→0

x(t+h)−x(t)

h ≒ x(t+h)−x(t)

h for 0 < h≪1

と期待できる．

(10)

Discretization（離散化）

dx(t)

d t = f(t,x(t)) (t >0), x(0) = a

最初に時間変数tに関する離散化を行う．0< h≪ 1を1つ固定し t0 =0, tn+1 = tn+h (n=0,1,2, . . .)

とおく．当然ながら

t_n= nh (n=0,1,2, . . .) である．

t t

₁

t

₂

t

₃

t

₄

t

₅

t

₀

=0

h h h h h

(11)

Euler’s method

dx(t)

d t = f(t,x(t)) (t >0), x(0) = a x(t)を上の微分方程式の解とする．(Letx(t)be a solution.) x(t)^の t= tnにおける微分係数の定義より

dx

d t(t_n) = lim

h→0

x(tn+h)−x(tn)

h ≒ x(tn+h)−x(tn) h

となる．さらに，tn+1 = tn+hなので，微分方程式と合わせて f(tn,x(tn))= dx

d t(tn)≒ x(t_n₊₁)−x(t_n) h が得られる．ゆえに

x(t_n₊₁) ≒ x(t_n)+hf(t_n,x(t_n)). という（導関数を含まない）近似式が得られた．

(12)

Euler’s method

近似式において，“≒^”を等式にして得られる漸化式

x0 = a, xn+1 = xn+hf(tn,xn) (n=0,1,2, . . .) により数列{x_n}^{を定める．}

(In the approximation, the sequence{xn}is determined by a reduction formula obtained by making≒an equal sign.)

このとき，第 n項x_nは未知関数x(t)の t= t_nにおける値x(t_n)に近いと考えられる．

(In this case, thenth termxnis considered to be an approximation of the valuex(t_n).)

Remark

f(t_n,x_n)が“t = t_nにおける接線の傾き”のような役割を果たしている．

(f(t_n,x_n)plays a role like the slope of the tangent at t = t_n.)

(13)

Euler’s method

x₀ = a, x_n₊₁ = x_n+hf(t_n,x_n) (n=0,1,2, . . .)

t t

1

t

2

t

3

t

4

t

5

t

0

=0

h h h h h

x

a x

₁

x(t

1

) x(t)

x=a+hf(0,a)

(Remark) 私たちには解 x(t)の本当のグラフはわからない．そこで，

x^′(0)の代わりに f(t0,x0) = f(0,a)の値を傾きとして採用している．

(14)

Euler’s method

x₀ = a, x_n₊₁ = x_n+hf(t_n,x_n) (n=0,1,2, . . .)

t t

1

t

2

t

3

t

4

t

5

t

0

=0

h h h h h

x

a x

₁

x(t

₂

) x(t) x

2

x=f(t

₁

,x

₁

)(t-t

₁

)+x

₁

(Remark) 私たちには真の値x(t₁)はわからない．そこで，x(t₁) ≒ x₁とみて f(t1,x1)の値を傾きとしている．

(15)

Euler’s method (Example)

dx(t)

d t = −2t x(t) (t >0), x(0)= a

上の微分方程式にオイラー法を適用すると，f(t,x) =−2t xより xn+1 = xn+hf(tn,xn)= xn−2tnxnh=(1−2nh²)xn

となる．例えば，初期値をa = x₀ =10，刻み幅をh=0.1とすれば x₁ =(1−2·0·0.1²)x₀ =10

x₂ =(1−2·1·0.1²)x₁ =9.8 x₃ =(1−2·2·0.1²)x₂ =9.408 x₄ =(1−2·3·0.1²)x₃ =8.84352 x₅ =(1−2·4·0.1²)x₄ =8.136038 と順次 xnが構成できる．

(16)

Euler’s method (Example)

(17)

Advantage

オイラー法によるコンピューターシミュレーションは難しくない．実際，

オイラー法は x_nから x_n₊₁を求めるのに四則演算や代入計算を繰り返すだけだからである．微分方程式の近似解を見つけるために，微分や積分の計算がまったく現れないことには注意してほしい．

(The sequencex_nis determined sequentially, and differential and integral calculations are not used to ﬁnd approximate solutions to the differential equations.)

Problem

微分方程式の解 x(t)がわからないときに，求めた近似値x_nと真の値 x(tn)の誤差を見積もることはできるか？もし可能ならばどのようにすればよいか？

(When the solutionx(t)to a differential equation is not known, can we estimate the error between the approximate value x_nand the true value x(tn)? If so, how can we do it?)

(18)

オイラー法の誤差評価(error estimation) x(t)を t= t_nのまわりでTaylor展開すれば

x(t_n₊₁) = x(t_n+h)

= x(t_n)+x^′(t_n)h+ x^′′(tn)

2 h²+ x^′′′(tn)

6 h³+· · ·

= x(t_n)+x^′(t_n)h+O(h²) (h→0) となる．

(Remark)

|h| ≪ 1とする(ex. h=10⁻²)．n≥2, a_n, 0のとき a₂h²+a₃h³+a₄h⁴+· · ·+anhⁿ

= h²(a₂+a₃h+a₄h²+· · ·+anhⁿ⁻²)

≒ h²(a₂+0)= a₂h²

このようにa₂h²+a₃h³+a₄h⁴+· · · ^はhについて少なくとも2次の微小量なのでO(h²) (h →0)と表す．

(19)

Consideration of Error

dx(t)

d t = f(t,x(t)) (t >0), x(0) = a

Euler’s method

x₀ = a, xn+1 = xn+hf(tn,xn) (n=0,1,2, . . .) 上の微分方程式の解x(t)^を t= tnのまわりでTaylor展開すれば

x(t_n₊₁) = x(t_n)+x^′(t_n)h+O(h²)

= x(t_n)+hf(t_n,x(t_n))+O(h²) (h→0)

となる．よって，漸化式で1ステップ進むのにO(h²)の項を無視しているので，その程度の誤差は生じうる．

もし x(t)のt =2での値を近似計算しようと思えば，刻み幅がhならば2/hステップの計算が必要となる．誤差は累積するから，この場合の誤差は大雑把に見ればO(h²)×2/h=O(h)程度となる．これよりオイラー法は精度よい計算法とはいえない．どう改良すればよいか？

(20)

Error estimate of Euler’s method

Assumption (Lipschitz condition)

ある定数L >0が存在して

|f(t,x)− f(t,x)˜ | ≤ L|x−x˜|

がすべての t,x,x˜ について成り立つ．これはリプシッツ条件と呼ばれる．

x₀ = a, x_n₊₁ = x_n+hf(t_n,x_n) (n=0,1,2, . . .) Taylor展開の式より

x(tn+1) = x(tn)+hf(

tn,x(tn))+O(h²) (h→0) なので，漸化式を引いて誤差を評価すれば

|x(tn+1)−xn+1| ≤ |x(tn)−xn|+h|f(

tn,x(tn))− f(tn,xn)|+O(h²)

≤ (1+Lh)|x(tn)−xn|+Ah² となる．ただし，Aはある正の定数．

(21)

Error estimate of Euler’s method

これより，x(t₀) = a= x₀と合わせて

|x(tn)−xn| ≤ (1+Lh)|x(tn−1)−xn−1|+ Ah²

≤ (1+Lh)²|x(tn−2)−xn−2|+ Ah²+(1+Lh)Ah²

≤ · · ·

≤ (1+Lh)ⁿ|x(t₀)−x₀|

+ Ah²{1+(1+Lh)+· · ·+(1+Lh)ⁿ⁻¹}

= 0+ Ah²(1+Lh)ⁿ−1

(1+Lh)−1 = Ah(1+Lh)ⁿ−1 L

が成り立つ．そこで，C= A/L^{とおけば，不等式}1+x≤e^x^により

|x(t_n)−x_n| ≤Ch{(1+Lh)ⁿ−1} ≤Ch(e^Ltⁿ−1) とわかる．従って，大域的な誤差はO(h)と評価できる．

(22)

Taylor expansion of a solutionx(t)

x(t)の t= t_nのまわりのTaylor展開で2次の項まで残せば x(t_n₊₁) = x(t_n+h)= x(t_n)+x^′(t_n)h+ x^′′(tn)

2 h²+ x^′′′(tn)

6 h³+· · ·

= x(t_n)+ f(t_n,x(t_n))h+ x^′′(t_n)

2 h²+O(h³) となる．これを基本アイデアとした計算法がホイン法である．

実際には，x^′′(t_n)を刻みの両端の微分係数を用いて近似し，オイラー法を援用して次のように漸化式を定める．

x(t_n₊₁) = x(t_n)+ h

2{x^′(t_n)+x^′(t_n+ h)}+O(h³)

= x(t_n)+ h

2{f(t_n,x(t_n))+ f(t_n₊₁,x(t_n₊₁))}+O(h³)

≒ x(t_n)+ h 2{f(

t_n,x(t_n))+ f(

t_n₊₁,x(t_n)+hf(

t_n,x(t_n)))}

(23)

Heun’s method（ホイン法）

ホイン法は通常次のように書かれる．

Heun’s method







k₁ = f(t_n,x_n)

k₂ = f(t_n+h,x_n+hk₁) xn+1 = xn+h k₁+k₂

2

t xn

tn

xn+hk1

k2

xn+1

(k1+k2)/2 k1

x k2

tn+1

(24)

The Runge-Kutta method（ルンゲ・クッタ法）

The Runge-Kutta method

x₀ = a, t_n= nhand forn=0,1,2, . . .











k₁ = f(t_n,x_n) k2 = f(tn+ h

2,xn+ h 2k1) k3 = f(tn+ h

2,xn+ h 2k2) k₄ = f(t_n+h,x_n+hk₃)

x_n₊₁ = x_n+h k₁+2k₂+2k₃+k₄ 6

ルンゲ・クッタ法は4次の項まで残した計算法で，1ステップの誤差は O(h⁵)程度である．そのため，大域的な誤差はO(h⁴)程度となる．

(25)

The Runge-Kutta method（ルンゲ・クッタ法）

t

_n

t

_n+1

k

₁

k

₂

h/2 h/2

x

_n

+hk

₃

x

_n+1

k

₄

(k

₁

+2k

₂

+2k

₃

+k

₄

)/6

k

₃

x

_n

+hk

₂

/2

x

_n

+hk

₁

/2 x

_n

(26)

The Runge-Kutta method（example）

(27)

Conclusion

Euler’s method Heun’s method

The Runge-Kutta method

数値計算法は今日紹介した以外にも多数あるので，精度とコストのバランスを考えなければならない．個人的には自然科学分野においてはルンゲ・クッタ法が広く使われている印象はある．

微分方程式の解の存在や一意性，および近似による誤差についてはシミュレーション結果だけではわからない．そのため，数値計算の妥当性を保証するための数学の重要性が増しつつあるように思う．

(Since there are many other numerical methods besides the ones introduced today, a balance between accuracy and computational complexity must be considered. I have the impression that the

Runge-Kutta method is widely used in the natural sciences.The numerical simulation results alone do not tell us whether the solution to the

differential equation exists or not, the uniqueness of the solution, or the error of the approximation. Therefore, the importance of mathematics in ensuring the validity of numerical calculations seems to be increasing.)

(28)

References

数値シミュレーションの文献は多数あるが，例えば

E. Kreyszig著,田村義保訳,技術者のための高等数学5 数値解析

（原書第8版）,培風館, 2003.

伊理正夫・藤野和建,数値計算の常識,共立出版, 1985.

山崎郭滋,偏微分方程式の数値解法入門,森北出版, 1993.

高見穎郎・河村哲也,東京大学基礎工学双書偏微分方程式の差分解法,東京大学出版会, 1994.

Erwin Kreyszig, Advanced Engineering Mathematics (10th Edition), Wiley, 2011.

Kendall E. Atkinson, Weimin Han, and David E. Stewart, Numerical Solution of Ordinary Differential Equations (1st Edition), Wiley, 2009.

数値シミュレーションの基礎 オイラー法とルンゲ・クッタ法