● 講義資料
▼ 講義予定
• 常微分方程式の境界値問題の数値解法
• 連立一次方程式の数値解法
● 前回の講義のまとめ
★ Runge-Kutta法の誤差の影響と安定性
• p次の離散変数法でつくった数値解{xn}と,厳密解xとの誤差は,
|x(tn)−xn| ≤Chp
と評価できる. （これが p 次であることの定義である）しかし, 実際に計算しているのは, 数値計算スキームで定義された{xn} ではなく,浮動小数点演算の丸め誤差を含んだ値の列 {xn}であり,xn−1−xn−1 の仮定の下で,
|xn−xn| ∼Cδ
が成り立つ. (C はスキームと微分方程式に依存する定数である）したがって,実際の計算値 と厳密解の値には
|x(tn)−xn| ≤C(hp−nδ) が成り立つ. ここで,t= 1まで計算すると考え, 1 =nhとおくと,
|x(tn)−xn| ≤C(hp−δ h)
となる. したがって,hp+1∼δ 程度よりも小さなhをとったとき,計算値 xn の値には丸め 誤差が無視できない程度入り込んでいることがわかる. 特に,hはある程度以上小さくとって はいけないことがわかる.
• x′(t) =−x(t),x(0) = 1の解はx(t) =e−tであり,安定な解の代表例である. したがって,こ の方程式の数値解{xn} を求めたとき,|xn| →0 が成り立つことが,その数値計算スキーム が持つべき必要条件であると考えられる.
しかし,修正オイラー法xn+1=xn+ 2hf(xn)をこの方程式に適用すると,どのようなh >0 をとっても|xn| → ∞となっていた. （数値的不安定性）hを適切にとることにより,この ような現象がRunge-Kutta型公式については成り立たないことを示す必要がある.
• x′(t) =λx(t), x(0) = 1 の数値解{xhn} が|xhn| →0 を満たすようなz =λhの領域を,その 数値計算スキームに対する安定領域と呼ぶ.
– 前進オイラー法の場合,xn+1 =xn+λhxnであるので,xn+1= (1 +z)xn,xn= (1 +z)n となり,安定領域S は
S={z∈C:|1 +z|<1}
となる. したがって,h∈(0,2)を満たさない hをとると,x′ =−xに対する数値解は
|xn| →0を満たさないことがわかる.
– 一般のp段 p次陽的Runge-Kutta法の場合(p= 2, 3, 4), xn+1=
1 +z+· · ·+zp p!
xn
となる. （簡単な計算により確かめることができる.）したがって,安定領域S は S=
z∈C:
1 +z+· · ·+zp p!
<1
となる. この場合も, 安定領域の図を書けば,各 pに対して, あるhp >0が存在して, h∈(0, hp)を満たさない hをとるとx′ =−xに対する数値解は|xn| →0 を満たさな いことがわかる. （hp≥2であることは,計算からわかる.）
特に,hはある程度以上大きくとることはできないことがわかる.
• 以上により,陽的Runge-Kutta法については,刻み幅hはある適切な範囲でしか利用できな いことがわかる.
★ Runge-Kutta法と数値積分の関係
• Runge-Kutta法
ki =f(x0+h
s
X
j=1
aijkj), i= 1, . . . , s,
x1=x0+h
s
X
i=1
biki,
を x′(t) = f(t, x(t)) に適用することを考える. そのため X(t) = (t, x(t)), F(X(t)) = (1, f(t, x(t))すれば,微分方程式X′(t) =F(X(t))は,
d dt
t x(t)
!
= 1
f(t, x(t))
!
となる. よって, Runge-Kutta法は
ki=f(t0+hX
j=1
aij, x0+h
s
X
j=1
aijkj), i= 1, . . . , s,
x1=x0+h
s
X
i=1
biki,
と書くことができる.
• 特に, Runge-Kutta法を
x′(t) =f(t), x(0) = 0 に適用すると,
ki=f(t0+hX
j=1
aij), i= 1, . . . , s,
x1=x0+h
s
X
i=1
biki,
となるので,
xn+1=xn+h
s
X
i=1
bif(tn+ci), ci=
s
X
j=1
aij
となる. 一方,
x(tn+1)−x(tn) = Z tn+1
tn
f(s)ds であるので, Runge-Kutta法は,右辺の積分の近似値を
Z tn+1
tn
f(s)ds∼h
s
X
i=1
bif(tn+ci)
と計算していることがわかる.
前進オイラー法 この場合,
Z tn+1
tn
f(s)ds∼hf(tn)
となり,区間[tn, tn+1]上の積分を高さf(tn),幅hの長方形の面積で近似する計算式を 与えている.
後退オイラー法 この場合,
Z tn+1
tn
f(s)ds∼hf(tn+1)
となり,区間[tn, tn+1]上の積分を高さf(tn+1),幅 hの長方形の面積で近似する計算式 を与えている.
ホインの２次公式 この場合, Z tn+1
tn
f(s)ds∼ h
2(f(tn) +f(tn+1))
となり,区間[tn, tn+1]上の積分を, (tn,0), (tn, f(tn)), (tn+1, f(tn+1)), (tn+1,0)の４点 を頂点とする台形の面積で近似する計算式を与えている. これを台形公式とよぶ.
台形公式は,区間[tn, tn+1]上で,関数f を(tn, f(tn)), (tn+1, f(tn+1))を通る直線（１ 次関数）で近似して積分を求めたことに他ならない.
★ Symplectic法
• Runge-Kutta法は,広範囲の常微分方程式の初期値問題に適用可能な汎用的な方法であるが,
力学に由来するエネルギー保存則が成り立つ系では, Runge-Kutta法による数値解はエネル ギーを保存しない.
たとえば,x′′=−xの解に沿って,エネルギー E(t) =1
2((x(t))2+ (x′(t))2)
は一定値を持つ. この証明はE′(t) = 0を示す方法もあるが,方程式の両辺にx′ をかけると 0 =x′′(t)x′(t) +x(t)x′(t) =1
2 d
dt((x′(t))2+ (x(t))2)
となることからも証明できる.
しかし,前進オイラー法ではEn = (1/2)(x2n+vn2)は単調増加となり,古典的Runge-Kutta 法ではEn は単調減少となってしまい,そのような数値解は「よい」数値解と言うことは難 しい. そこで,エネルギーまたはそれに近いものが保存するような数値解を構成することが 望まれる.
• ここでは, Hamilton系と呼ばれる常微分方程式を考える.
未知関数を(q1, . . . , qN, p1, . . . , pn), qi: [0, T]−→R, pi: [0, T]−→R として, Hamiltonian H:R2N −→Rが与えられているととき,常微分方程式
dq dt =∂H
∂p, dp
dt =−∂H
∂q を考える.
x′′=−xに対応する問題としては,q=x,p=x′ と考え,H(q, p) = (1/2)(q2+p2)とすれば, dq
dt = ∂H
∂p =p, dp
dt =−∂H
∂q =−q
より d
dt dq dt = d
dtp=−q
となるので,x′′=−xを考えることと,H(q, p) = (1/2)(q2+p2)に対するHamilton系を考 えることは同じである.
また, HamlitonianH は解に沿って一定値をとる. すなわち, dH
dt = 0 が成り立つ. （実際に計算してみればよい.）
なお,系の運動エネルギーが(1/2)Pq2i で書けているとき, HamiltonianH は,系の全エネ ルギーを与える関数に一致する.
• Hamilton系の持つ重要な性質は, Hamiltonianが保存することのほかに,次のLiouville の 定理が成り立つことである.
φt:R2N −→R2N を初期条件(q(0), p(0))に対する解(q(t), p(t))によるt 秒後の位置とす る. （この写像を相流,または1-parameter 変換群と呼ぶ.）この時,φt は R2N の任意の 領域の面積を保存する. すなわち,φtの微分Dφtに対して,
|det(Dφt)|= 1 が成り立つ. （φtは面積保存写像となる.）
• このことから, Hamilton 系に対する数値解法として, (qn, pn)7→(qn+1, pn+1)が面積保存写 像となるような解法が望ましいことがわかる. そのような解法をSymplectic 解法と呼ぶ.
• HamiltonianH が
H(q, p) =T(p) +U(q) と分解できている場合に,以下の計算スキーム
qn+1 =qn+hdT dp(pn), pn+1=pn−hdU
dq(qn+1),
をSymplectic Euler 法とよび, 1次のSymplectic解法となっていることが証明できる.
実際,１段目の写像(qn, pn)7→(qn+1, pn)をS1h,２段目の写像(qn+1, pn)7→(qn+1, pn+1)を Sh2 と書くと,それぞれの微分は
DSh1= I ∂p∂q∂2T
0 I
! ,
DSh2= I 0
−∂p∂q∂2U I
! , となり,それぞれ
det(DS1h) = 1, det(DSh2) = 1,
となるので, Symplectic Euler法のスキーム Sh: (qn, pn)7→(qn+1, pn+1)は, det(DSh) = 1 を満たし,面積保存写像となる.
実際には,DSh1,DSh2,DSh は, Symplectic群
Sp(N) =
M ∈M2N(R) :MTJM =J , J = 0 I
−I 0
!
に入る. （なので, Symplectic解法と呼ぶ.）
Symplectic Euler法をH(q, p) = (1/2)(q2+p2)に適用すると, qn+1
pn+1
!
= 1 h
−h 1−h2
! qn
pn
!
となり,この行列の行列式はたしかに1となっている.
• p次のSymplectic 解法で計算すると, HamiltonianH は保存しないかわりに, ある関数H˜ が存在して,
H(q˜ n, pn) =一定値,
H(q˜ n, pn)−H(qn, pn) =O(hp) が成り立つ.
実際H(q, p) = (1/2)(q2+p2)の時には,
H˜(q, p) = (1/2)(q2+p2+hqp) となる.
• H(q, p) =T(p) +U(q)の時により高次のSymplectic解法を構成するには,上記のSh1,Sh2 を 合成して
S2bkhSc1kh· · ·Sb21hS1c1h
を考え, パラメータ{bj}, {cj} をp次となるようにきめればよい. 実際,以下のようにきめ る方法が知られている.
St¨olmer-Verlet (p= 2)
1 2
bi 0 1
ci 1/2 1/2 Ruth (p= 3)
1 2 3
bi 7/24 3/4 −1/24 ci 2/3 −2/3 −1
Sanz-Serna (p= 4)
1 2 3 4 5 6
bi 7/48 3/8 −1/48 −1/48 3/8 7/48
ci 1/3 −1/3 1 −1/3 1/3 0
• 陰的Runge-Kutta法も,条件式の取り方によってはSymplectic解法とすることができる.
実際,陰的Runge-Kutta法がSymplecticとなる条件は, biaij+bjaji−bibj = 0 が成り立つことである.
陰的Runge-Kutta法で, Symplectic解法となるものの例としては, 陰的中点法（１段２次Gauss-Legrandre法）
1/2 1/2 1
２段４次Gauss-Legrandre法 3−√
3 6
1 4
3−2√ 3 3 +√ 12
3 6
3 + 2√ 3 12
1 4
1/2 1/2
３段６次Gauss-Legrandre法 5−√
15 10
5 36
80−24√ 15 360
50−12√ 15 1 360
2
50 + 15√ 15 360
2 9
50−15√ 15 5 +√ 360
15 10
50 + 12√ 15 360
80 + 24√ 15 360
5 36
5/18 4/9 5/18
などが知られている.
★ プログラム例１
以下のプログラムは, 2×2 行列の和を計算するプログラムの例である.
/* MAX x MAX の正方行列２つの和を計算する */
#include <stdio.h>
#define MAX (2)
void print_matrix(double *a, int size) ;
void add_matrix(double *c, double *a, double *b, int size) ;
int main(int argc, char **argv) {
double a[MAX*MAX] = {1.0, 2.0, 2.0, 3.0} ; double b[MAX*MAX] = {-1.0, -2.0,
-2.0, -3.0} ; double c[MAX*MAX] ;
int i, j ;
print_matrix(a, MAX) ; print_matrix(b, MAX) ; add_matrix(c, a, b, MAX) ; print_matrix(c, MAX) ; return 0 ;
}
void print_matrix(double *a, int size) {
int i, j ;
for(i=0;i<size;i++) {
for(j=0;j<size;j++) printf("%+f ", a[i*size+j]) ; printf("\n") ;
}
printf("\n") ; return ; }
void add_matrix(double *c, double *a, double *b, int size) {
int i, j ;
for(i=0;i<size;i++) {
for(j=0;j<size;j++) c[i*size+j] = a[i*size+j] + b[i*size+j] ; }
return ; }
● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
1. （★★★）4×4 行列の積を計算するプログラムを書きなさい.
2. （★★★）4×4 行列の転置行列を計算するプログラムを書きなさい.
3. （★★★）連立一次方程式
3x+ 12y+ 9z= 3 2x+ 5 y+ 4z= 4
−x+ 3 y+ 2z= 5
をGauss-Jordanの消去法（掃出し法）で解くプログラムを書きなさい.
