– 初期値問題 – 数値計算講義第８回微分方程式の解法

数値計算講義 第８ 回

微分方程式の解法

初期値問題

カ ーネン コ

金 子

晃

[email protected] [email protected]

http://www.kanenko.com/

常微分方程式の初期値問題

１ 階常微分方程式の初期値問題

( dx

dt =f(t, x), t₀ ≤t≤t_max, x(t₀) =c

は， 区間

を

等分し て

N

^{と する と き ，} 前進差分近似

h =f(t, x(t)),

すなわち

から 得ら れる 反復公式

x₁ =x₀+hf(t₀, x₀),

· · · ·

x_n+1=x_n+hf(t₀+nh, x_n)

によ り ， 容易に解く こ と ができ る ．

の折れ線法

f c

c

( , )

t x

h 2 h 0

(

図は

の場合

連立常微分方程式の場合

dt =f(t, x, y), dy

dt =g(t, x, y), t₀≤t≤t_max, x(t0) =a, y(t0) =b

は， 次のよ う に離散化する ：

x₁ =x₀+hf(t₀, x₀, y₀), y₁=y₀+hg(t₀, x₀, y₀)

· · · ·

x_n+1=x_n+hf(t₀+nh, x_n, y_n), y_n+1=y_n+hg(t₀+nh, x_n, y_n),

２ 階単独常微分方程式の場合





 d²x

dt² =f t, x,dx dt

, t₀≤t≤t_max, x(t₀) =a, x^′(t₀) =b

は， そのま ま 離散化せず， １ 階連立常微分方程式

dt =y, dy

dt =f(t, x, y), t₀≤t≤t_max, x(t₀) =a, y(t₀) =b

に直し て から 上の離散化を 施すのが普通．

実装例

単独方程式の場合， 解法の核心部分はたっ たこ れだけ ：

時間の分割数を 設定

T=TMIN ! T

を 初期時刻に設定

X=C ! X

を 初期値に設定

DO 100 I=1, N

X=X+F(T,X)*H !

次の時刻の

を 現在の傾き から 計算

T=T+H !

時刻を 進める

100 CONTINUE

連立方程式の場合も ち ょっ と 変わる だけ ：

H=(TMAX-TMIN)/N !

時間の分割数を 設定

T=TMIN ! T

を 初期時刻に設定

X=A ! X

を 初期値に設定

Y=B ! Y

を 初期値に設定

DO 100 I=1, N

XS=X+H*F(T,X,Y) !

次の時刻の

を 計算し 保存

次の時刻の

を 計算

X=XS !

次の時刻の

を 保存値から 戻す

100 CONTINUE

モデルプロ グラ ムは， 結果を 直接描画する ため，

を リ ン ク し

ループ中の計算値を

指令で直接描画し て いる ．

gnuplot

で後から 描画する と き は

ループ中で計算し た

の値を

フ ァ イ ルに書き 出し

後で一括描画する

注意 ダミ ー変数

を 使う のを 忘れないよ う に ！

Y=Y+H*G(T,X,Y)

と 書く と ，

の折れ線法の正し い実装ではなく なる ．

間違え る と ど う なる か単振動の連立化方程式で観察し て みよ う

．

の場合，

ステッ プ分を 描画する と

տ

正し い

法

誤っ た方法

ど う も 話がう ま すぎる と ， 後者の軌道を

に拡大し て ，

でそれぞれ

ステッ プやっ て みる と

一周当たり の誤差のサイ ズは正し い実装と 同じ 程度だが，

軌道を 発散さ せないよ う な不思議な保存性を 持つ ！

Euler-Cauchy

法の誤差評価

|f(t, x)| ≤M₀

及び，

条件

|f(t, x)−f(s, y)| ≤m₁|t−s|+M₁|x−y|

を 仮定し ， １ ステッ プ当たり の近似解と 真の解と の差の拡大を 見る ． 真の解

は， 積分方程式

x(t_n+1) =x(t_n) + Z tn+1

tn

f(t, x(t))dt

を 満たし て いる ．

他方

近似解は，

x_n+1 =x_n+hf(t_n, x_n) =x_n+ Z _tn+1

tn

f(t_n, x_n)dt

を 満たし て いる ． よ っ て 差は，

|x(t_n+1)−x_n+1|

≤ |x(tn)−xn|+

Z tn+1

tⁿ

f(t, x(t))dt− Z tn+1

tⁿ

f(tn, xn)dt

≤ |x(t_n)−x_n|+ Z tn+1

tn

|f(t, x(t))−f(t_n, x_n)|dt

こ こ で

条件等を 用いて ， 積分記号下を

|f(t, x(t))−f(t_n, x_n)| ≤m₁(t−t_n) +M₁|x(t)−x_n|

≤m₁(t−t_n) +M₁|x(t_n)−x_n|+M₁|x(t)−x(t_n)|

と 分解する と ， 最後の項に

|x(t)−x(t_n)|=

Z _t

tn

f(s, x(s))ds

≤M₀(t−t_n)

(

こ こ に

は考え て いる 領域での

の上限

を 用いて ， 結局

|f(t, x(t))−f(t_n, x_n)| ≤(m₁+M₀M₁)(t−t_n) +M₁|x(t_n)−x_n|

よ っ て ， 最終的に

|x(tn+1)−xn+1|

≤ |x(t_n)−x_n|+m₁+M₀M₁

2 (t_n+1−t_n)²+M₁(t_n+1−t_n)|x(t_n)−x_n| i.e.

|x(t_n+1)−x_n+1| ≤(1 +M₁h)|x(t_n)−x_n|+ m₁+M₀M₁

2 h²

と いう １ 段分に対する 誤差の伝播公式が得ら れた．

こ れを

について 反復適用する ． ある いは，

|x(t_N)−x_N| ≤(1 +M₁h)|x(t_N−1)−x_N−1|+m₁+M₀M₁

2 h²

(1+M₁h)|x(t_N−1)−x_N−1| ≤(1+M₁h)²|x(t_N−2)−x_N−2|+(1+M₁h)m1+M0M1

2 h²

. . . .

(1+M₁h)^N−1|x(t₁)−x₁| ≤(1+M₁h)^N|x(t₀)−x₀|+(1+M₁h)^N−1m₁+M₀M₁

2 h²

を 総和し て ， 最後の項に等比級数の和の公式を 使う と

|x(tN)−xN| ≤(1 +M1h)^N|x(t0)−x0|+(1 +M₁h)^N −1 M₁h

m₁+M₀M₁

2 h²

こ こ で

は

によ ら ぬ定数，

ま た

よ り ，

|x(t_N)−x_N| ≤e^M1T|x(t₀)−x₀|+m₁+M₀M₁ 2M₁ e^M1Th

が得ら れ， １ 次の近似である こ と が確定し た．

普通は， 初期値には誤差は無いも のと し ， 右辺の第１ 項を 省く ．

初期値に誤差が有る と き は， 各段の計算にも 丸め誤差を 入れる のが妥当で，

こ の場合は上の評価は

|x(t_N)−x_N| ≤e^M1T|x(t₀)−x₀|+e^M1T M₁

m₁+M₀M₁

2 h+ ε

h

と なる ．

常微分方程式に対する 種々 の問題

初期値問題





 d²u

dx² +q(x)u=f(x), u(0) =u0, du

dx(0) =u1

離散化はも っ と も 簡単で

逐次代入で解ける

境界値問題





 d²u

dx² +q(x)u=f(x), u(0) = 0, u(1) = 0

離散化は連立一次方程式を 解く こ と

固有値問題





 d²u

dx² +q(x)u=λu, u(0) = 0, u(1) = 0

離散化は行列の固有値問題と なる

が有限行列のと き

Au=λu+f

が一意に解ける

が行列

の固有値でない

境界条件付き の線型常微分作用素は無限次元の行列だが

同じ こ と が成り 立つ

２ 次の

法

Euler

法は， 右辺が未知函数

を 含ま ないと き ， すなわち

単なる 定積分のと き は，

近似和に対応し ， はなはだ近似度が悪い．

こ れを 台形公式に相当する も のに改良し たのが， ２ 次の

法：

u_n=hf(t_n, x_n),

vn=hf(tn+1, xn+un), x_n+1=x_n+1

2(u_n+v_n)

時刻

における 傾き で予測し た到達点での傾き で進み直し たも のと 平均する と 覚え て おく と よ い．

実装の核心部分は

U=H*F(T,X) V=H*F(T+H,X+U) X=X+(U+V)/2

f( ,x ) t x

tn t_n+h=t_n+1 t

n u_n

n n

v_n

f (t_n+h,x_n+u_n)

４ 次の

法

こ れは定積分の

公式に相当する も ので，

普通

法と 呼べばこ れを 指す．

x_n

から

の更新アルゴリ ズムは次の通り ：

v_n=hf(t_n+h/2, x_n+u_n/2), w_n=hf(t_n+h/2, x_n+v_n/2), z_n=hf(t_n+h, x_n+w_n),

xn+1=xn+un/6 +vn/3 +wn/3 +zn/6

実装する と き は，

は共通の臨時変数を 使え ばよ いので

XS=X+U/6

U=H*F(T+H2,X+U/2) XS=XS+U/3

U=H*F(T+H2,X+U/2) XS=XS+U/3

U=H*F(T+H,X+U) X=XS+U/6

f( ,x ) t x

tn t_{n+ /2}h t_n+h=t_n+1 t

n u_n

n n

f (tn+ /2h ,xn+u_n/2) v_n

f (tn+h_/2,xn+v_n/2) w_n f (t_n+h,x_n+w_n)

z_n

※

が予め代入さ れて いる

少し でも 計算を 速く する けち な手法

Runge-Kutta

法の近似のオ ーダー

Runge-Kutta

法の打ち 切り 誤差の評価は非常に面倒で，

書かれて いる 書物はほと んど 無い．

よ く 書かれて いる のはオーダーだけの評価だが， それも ４ 次の公式は大変．

こ こ では２ 次の公式に対し て ， オーダー評価を 示す．

u_n=hf(t_n, x_n),

v_n=hf(t_n+h, x_n+u_n)

=hf(t_n, x_n) +h²f_t(t_n, x_n) +hu_nf_x(t_n, x_n) +O(h³)

=hf(t_n, x_n) +h²f_t(t_n, x_n) +h²f(t_n, x_n)f_x(t_n, x_n) +O(h³)

∴

2(un+vn)

=x_n+hf(t_n, x_n) +h²

2 f_t(t_n.x_n) +h²

2 f(t_n.x_n)f_x(t_n.x_n) +O(h³).

他方，

を 真の解と し て

f(s, x(s)) =f(t_n, x(t_n)) + (s−t_n)f_t(t_n, x(t_n))

+ (x(s)−x(t_n))f_x(t_n, x(t_n)) +O((s−t_n)²)

=f(t_n, x(t_n)) + (s−t_n)f_t(t_n, x(t_n))

+ (s−t_n)f(t_n, x(t_n))f_x(t_n, x(t_n)) +O((s−t_n)²)

を 積分し て ，

x(t_n+1) =x(t_n) +hf(t_n, x(t_n)) +h²

2 f_t(t_n, x(t_n)) + h²

2 f(t_n, x(t_n))f_x(t_n, x(t_n)) +O(h³)

よ っ て 両者の差を と り

の中の

を

で置き 換え る と

|x_n+1−x(t_n+1)| ≤ |x_n−x(t_n)|+M h|x_n−x(t_n)|+O(h³) 12

≤(1 +M h)|x_n−x(t_n)|+O(h³)

こ れを

について 書き 並べ

|x_N −x(t_N)| ≤(1 +M h)|x_N−1−x(t_N−1)|+O(h³)

(1 +M h)|x_N−1−x(t_N−1)| ≤(1 +M h)²|x_N−2−x(t_N−2)|+ (1 +M h)O(h³)

· · · ·

(1 +M h)^N−1|x₁−x(t₁)| ≤(1 +M h)^N|x₀−x(t₀)|+ (1 +M h)^N−1O(h³)

辺々 加え る と

に注意し て

|x_N −x(t_N)| ≤(1 +M h)^N|x₀−x(t₀)|+ (1 +M h)^N −1 M h O(h³)

≤e^M(b−a)|x₀−x(t₀)|+ e^M(b−a) M O(h²)

を 得る ．

※ 以上の評価では

の

に関する

連続性を 使う ので，

例え ば

が有界なら

上で

と 書いた部分は

こ の定数を 使え ば

見積り 可能なある

定数によ り

で抑え ら れ

最後の

も

にでき る

４ 次の

公式に対し て も 同様にし て １ 段の誤差が

従っ て 全体と し て の誤差が

である こ と が理論的に示せる ．

教科書の第８ 章参照

平面の力学系

１ 階連立常微分方程式は， 力学系

と も 呼ばれる ． 粒子が方程式で記述さ れる 法則に従い， 時間と と も に位置を 変え る 様子

流れ

を 研究する 学問である ．

こ の言葉の起源は

の運動方程式：

直線上を 保存力の場

に従っ て 運動する 質点は

dt² =f(x)

と いう 方程式を 満たす． 運動量座標

dt

を 導入する と ， １ 階連立系：





 dx

dt = p m, dp

dt =f(x)

こ れは相空間

と 呼ばれる

空間の流れを 記述する

実用的な力学系の例： 惑星の運動

２ 体問題

．

Md²R

dt² =γM m r−R

|R−r|³, md²r

dt² =γM m R−r

|R−r|³

前者が太陽， 後者が惑星と すれば

Md²R 14

dt² =γM m r−R

|R−r|³, md²r

dt² =γM m R−r

|R−r|³

二つの方程式を 加え る と

d²

dt²(MR+mr) = 0

∴

dt(MR+mr) =

一定

重心

M +m

の等速度運動を 表現する 式

重心座標を 導入し

重心から の相対運動を 考察する ：

r_G=r−MR+mr

M+m = M

M+m(r−R)

と 置けば

md²r_G

dt² =md²r

dt² =γM m R−r

|R−r|³ =−γ M³m (M+m)²

r_G

|r_G|³

こ れは， 相空間が４ 次元内の流れだが， 空間座標だけを リ アルタ イ ムで 描画する と ２ Ｄ で表現でき る

(

更に， 解析的手段で本当の平面の力学系に帰着でき る ．

こ の連立常微分方程式を 手計算で解いて 実際に楕円軌道を 得る と

の感激を 追体験でき る ので

ぜひやっ て みよ う

Van der Pol

方程式 平面の力学系の実用的な代表例と し て

緩和振動を 記述する 有名な方程式

x^′′−ε(1−x²)x^′+x= 0,

一階連立系と し て は

x^′ =y,

y^′=−x+ε(1−x²)y

|x|>1

のと き はブレ ーキを かけ

のと き はアク セルを かける こ と によ り

ち ょ う ど 周期

の安定な軌道を 維持さ せる 仕組み

高次元の力学系

物理法則は一般に時間に依存し な いので，

力学系の右辺は

を 含ま ないのが普通．

こ のよ う な方程式系は自励系と 呼ばれる ． 自励系では， 解軌道は時刻の平行移動で不変．

従っ て 初期値問題の解の一意性によ り ， 二つの解軌道は交わら ない．

特に， 平面の力学系は軌道が互いに他を 分離する ので， 構造が簡単と なる ． こ れに対し ， ３ 次元以上になる と ， 一方の軌道が他方に巻き 付いたり でき ， 非常に複雑な軌道が現れる ．

高次元力学系で古来有名なのが３ 体問題で， 未だに分から ないこ と が多い．

つい数年前にも 新たな周期解

八の字解

が発見さ れたり し て いる ． ま た， 偏微分方程式の離散化で乱流のモデル等と し て 得ら れた方程式系 の中には， カオス的な動き を する も のがあり ， 決定論的カオス と 呼ばれる ．

(Newton

力学では， 初期位置と 初期速度を 指定すれば， 質点のその後の運動は

一意に決ま っ て し ま う

決定論的

．

それなのに， その振舞が複雑で， 予測不可能に見え る

カオス的

ので，

こ う 呼ばれる ．

プロ グラ ム例

を コ ン パイ ルし 実行し て

決定論的カオスの世界を 味わっ て みよ ．

Lorenz

の方程式

気象学者が考え 出し た， 気体の運動方程式の簡略化モデルで，

x^′=σ(y−x), y^′=x(R−z)−y, z^′ =xy−bz R¨ossler

の方程式は

それを 更に単純化し た

x^′ =−y−z, y^′ =x+ay, z^′ =bx+z(x−c)

x y z

x y z

Lorenz

アト ラ ク タ

アト ラ ク タ ど ち ら も 単なる

法で書かれて いる ．

こ れら はカオス的な現象が見え て 来る ま で軌道を 追跡する と ，

かなり の誤差が累積する はずだが， も っ と も ら し い軌道に見え る のはな ぜか ？ ど の軌道も アナロ グ的には似たよ う な 振舞を する ので， 軌道は変わっ て も 人間には感知でき ない．

ある いみではイ ン チキを 見せて いる

本日の講義内容の自習課題

初期値問題

( dx dt =x,

x(0) = 1

^{を 解析的に解け}

ま たこ の近似解を 求める

法のプロ グラ ム

を コ ン パイ ル・ 実行し ， 分割数を いろ いろ 変え て 近似のオーダーを 数値的に観察せよ ． でき れば丸め誤差の影響も 考慮し ， 最良の刻み幅

を 推測せよ ．

２ 次元力学系の

およ び

を コ ン パイ ル・ 実行し ，

種々 の

に対し て 挙動の違いを 観察せよ ．

3 Runge-Kutta

法のプロ グラ ム

を

コ ン パイ ル・ 実行し ， 分割数を いろ いろ 変え て 近似のオーダーを 調べよ ． でき れば丸め誤差の影響も 考慮し ， 最良の刻み幅

を 推測せよ ．

4 planet.f,vanderpol.f,lorenz.f,rossler.f,eight.c

を コ ン パイ ル・ 実行し ， 出力結果を 味わえ ．

sphereplanet.c, sphereeight.c

はそれぞれ

の

OpenGL

版です

は伊藤先生の『 コ ン ピュ ータ グラ フ ィッ ク ス』

で後期に学べま すが

一足早く

等でコ ン パイ ルし て 味わっ て みま し ょ う

(opengl-gcc

はシェ ルスク リ プト です

伊藤先生は

にし ていま す

本日の範囲の試験予想問題

問題

単振動の方程式を １ 階連立化し た

dt =f(x, y) =y, dy

dt =g(x, y) =−x

に対し

エネ ルギー保存則

を 示せ

解の具体的な函数形を 用いず

微分方程式から 導け

課題

単振動の方程式を １ 階連立化し た上記の方程式に， ダミ ー変数を 使い忘れた誤っ た実装を し た場合には， かなり 大き い

について も ， 軌道が無限に発散し ないと いう 現象を 数学的に正当化せよ ．

ま た正し い

法では

こ の量は時間ステッ プと と も に

ど う 変化する か述べよ

［ ヒ ン ト ： 第

ステッ プの近似解

に対し

量

が

によ ら ず一定と なる こ と を 計算で示せ

］

課題

振幅の大き な振り 子の運動方程式は，

dt² =−ksinx

と なり ， 初等関数では解く こ と ができ ない． こ の近似解を 求める 方法を 述べよ ．

初期条件と し て は， 振り 子を 振れ角

の位置で静かに離し た場合 を 仮定せよ ． 実際に動く プロ グラ ムま で書く には及ばない．

問題

常微分方程式の境界値問題

−d²u

dx² +q(x)u=f(x) (0< x <1), u(0) =u(1) = 0

を 離散化し て 解く と き ， ど のよ う な計算が必要と なる か？ ただし ，

は既知の関数と する ．

問題

法について 述べよ ． ま た， こ れを 最も 簡単な微分 方程式

に適用する と ，

に対する ど んな積分公式と なる か？

問題

微分方程式

dx =y

の初期条件

を 満たす解の

における 近似値を

メ ッ シュ 幅

について

法で

計算し た値

すなわち ，

法の１ 段だけの計算値

を 示せ．