数学原論 8 (高山)82005 年 7 月 5 日

差分化の復習

微分方程式の近似解を求める漸化式の導出方法（差分法) については前々回のプリントで説明した. 今回はこの方法を別の角度から説明しよう.
今 f (t) を未知関数とする微分方程式

f 0 (t) = a(t)f (t) + b(t),f (0) = cを解く問題を例として考える. ここで a(t), b(t) は時刻 t の関数であり, c は定数である. たとえば a(t) = −1, b(t) = | cos(t)|, c = 2 とすれば上の微分方程式は f 0 (t) =−f (t) + | cos(t)|, f (0) = 2 となる.
f (t) は時刻 t おけるある量, たとえば温度とか細菌の数とか物体の位置や電流の大きさを表すとすれば, 初期条件 f (0) = c での微分方程式を満たす関数 f (t) を求めることは, 時
刻 0 でのこれらの値から時刻 t での値を予想することにほかならない.
さて h = ∆t を微小な時間, たとえば h = 0.00001 とする. このとき微分の定義より f 0 (t)は (f (t + h) − f (t))/h にほぼ等しい. よってf (t + h) ≈ f (t) + hf 0 (t) = f (t) + (∆t)f 0 (t)である. 微分方程式の右辺を用いて f 0 (t) を書き換えると,

f (t + h) ≈ f (t) + h(a(t)f (t) + b(t))(1)となる. つまり f (t) の h = ∆t 秒後の値 f (t + h) は f (t) + h(a(t)f (t) + b(t)) に大体等しいということになる. つまり (1) は 微小時間 h 秒未来 の f の値を現在の f の値で表す式とみなせる.
(1) より t = 0 とすれば,f (h) ≈ f (0) + h(a(0)f (0) + b(0)) = c + h(a(0)c + b(0))
である. これで時刻 t = h での f (t) の近似値が求まった. 次に (1) で t = h とすれば,f (h + h) ≈ f (h) + h(a(h)f (h) + b(h)),
t = 2h とすれば
f (2h + h) ≈ f (2h) + h(a(2h)f (2h) + b(2h))となり, 時刻 t = 2h, t = 3h の f (t) の近似値が次々とさだまっていくこととなる. まとめ
ると漸化式
fk+1 = fk + h(a(hk)fk + b(hk)), f0 = c(2)で数列 fk を決めていけばそれが時刻 t = hk での f (t) の近似値となるのである.

9

連立方程式の差分化

問 9.1 上の考え方を用いて y1 (t), y2 (t) を未知関数とする次の連立の微分方程式の近似解を求めるプログラムを書きなさい. ((fk , gk ) の値をプロットしていくプログラムを書きな
さい. )y10 = (2 − y2 )y1 , y20 = (2y1 − 3)y2 , y1 (0) = 4, y2 (0) = 1.
参考: 種族 2 は種族 1 を食べる. y1 (t) は時刻 t における種族 1 の数. y2 (t) は時刻 t における種族 2 の数. と解釈する場合もある.
答え. h を微小な数とする. y1 (t) = f (t), y2 (t) = g(t) とおこう. このとき微分方程式よりf (t + h) − f (t) ≈ h(2 − g(t))f (t)g(t + h) − g(t) ≈ h(2f (t) − 3)g(t)である. よって,f (t + h) ≈ f (t) + h(2 − g(t))f (t)g(t + h) ≈ g(t) + h(2f (t) − 3)g(t)
この式の右辺は時刻 t での f, g の値できまる量であり, 左辺は時刻 t + h での f, g の値である. したがって漸化式

fk+1 = fk + h(2 − gk )fk
gk+1 = gk + h(2fk − 3)gk
f0 = 4,g0 = 1
で数列 fk , gk をきめていけば fk , gk は 時刻 hk での f (t), g(t) の近似値となる.

¶³

100 SET WINDOW -1,10,-1,10
110 DRAW axes
120 LET f=4130 LET g=1140 LET h=0.0001
150 FOR k=0 TO 200000
160
LET f2 = f+h*(2-g)*f
170
LET g2 = g+h*(2*f-3)*g
180
SET COLOR 0190
PLOT POINTS: f,g
200
SET COLOR 4210
PLOT POINTS: f2,g2
220
LET f = f2230
LET g = g2240 NEXT k250 ENDµ´

http://www.math.kobe-u.ac.jp/ ˜ taka/2005/prayp.bas よりダウンロードできる.
問 9.2 初期条件 f0 = 4, g0 = 1 を変更して種族 1, 2 が絶滅したり絶滅寸前になる初期条
件を実験的に決めよ.

10

カオス的力学系

解の関数が非常に複雑な形を描くような微分方程式があるということは 19 世紀から 20世紀前半の数学者により理論的にある程度解明されていたが, その理論は難解で科学者一般の認識となることはなかった.
グラフィックを容易に表示できる計算機と差分法による微分方程式の解法がこの状況を一変させることとなる. 1970 年代以降のことである. 解の関数が初期値に敏感に依存し, 非常に複雑な形を描くような微分方程式は カオス的な微分方程式 と呼ばれている. 現在はカオス的な方程式を含むいわゆる 複雑系 の考え方が提唱され多くの科学者の研究対象となっている.
ここではカオス的な微分方程式の代表の一つである ローレンツ方程式 の解のグラフを描
くプログラムを実行してその複雑な解を見てみよう.
課題 1(1) 次のプログラムを入力して実行し, 解のグラフを観察せよ.
(2) 180, 190, 200 行で初期条件を設定している. これらの値を変えて解のグラフがどのように変わるか調べよ. 解のグラフは 8 の字を描く. 8 の下の部分を何度か回ってから 8 の上の部分へ移り上の部分を何度かまわる. これを繰り返す. 回る回数を数えよ. 動きが早すぎるときは 230 行目と 240 行目の間にwait delay 0.1
命令をいれるとよい. 0.1 は 0.1 秒停止することを意味する.
¶³

100 ! lorentz.bas. Solving p1’ = -a p1 + a p2,
110 ! p2’ = -p1 p3 + b p1 - p2
120 ! p3’ = p1 p2 + c p3
130 SET WINDOW -25,25,-25,25
140 DRAW axes
150 LET a=10160 LET b=20170 LET c=2.66
180 LET p1=0190 LET p2 = 3
200 LET p3 = 0210 LET dt = 0.004
220 LET t = 0230 FOR t=0 TO 50 STEP dt
240
LET q1 = p1+dt*(-a*p1+a*p2)
250
LET q2 = p2+dt*(-p1*p3+b*p1-p2)
260
LET q3 = p3+dt*(p1*p2-c*p3)
262
SET COLOR 0264
PLOT POINTS: p1,p2266
SET COLOR 4270
PLOT POINTS: q1,q2280
LET p1 = q1
290
LET p2 = q2300
LET p3 = q3310 NEXT t
320 ENDµ´

ローレンツ方程式は 3 つの未知関数 p1 (t), p2 (t), p3 (t) についての次の連立微分方程式である. a, b, c には適当な数字をいれる.

p01 = −ap1 + ap2 ,p02 = −p1 p3 + bp1 − p2 ,p03 = p1 p2 + cp3
上のプログラムでは (p1 (t), p2 (t)) をプロットしている.
http://www.math.kobe-u.ac.jp/ ˜ taka/2004/lorentz2.txt よりダウンロードできる.
差分化の解説 240 行目を変形すると

q1 − p1= −a ∗ p1 + a ∗ p2dtである. dt を微小な数 (上のプログラムでは 0.004) としたとき, 変数 q1 に p1 (t + dt), 変数 p1 に p1 (t), 変数 q2 に p2 (t + dt), 変数 p2 に p2 (t), 変数 q3 に p3 (t + dt), 変数 p3 にp3 (t), が対応する. したがって 240 行目は−ap1 (t) + ap2 (t) =p1 (t + dt) − p1 (t)' p01 (t)dt(一つ目の微分方程式) にほかならない. 同様に 250, 260 行目はそれぞれ 2 つ目, 3 つ目の微分方程式にほかならない.

11

2 階の微分方程式の差分化

たとえば, 単振動の方程式

d2y+y =0dt2を数値的に解くことを考えてみよう.

¶³

物理既修者向け解説:
解く前に, この方程式の物理的意味を復習しておこう. 高校物理の最初の基本公式は質量 × 加速度 = 力なる関係式である. この関係式を Newton の運動方程式という. 物体の時刻 t における位置を q(t) とおくと, 速度は q 0 (t), 加速度は q 00 (t) である.
1 次元的に単振動をするバネについた 質量 1 の物体 W を考えよう. 時刻 t における,W の位置を y(t) とすることにする. ただし, バネが自然な長さにあるとき y = 0 とする. フックの法則によると, 自然な長さから y だけのびた (負のときはちじんだとみなす) とき物体 W にかかる力は −ky である. ここで k はバネできまる定数. ここでk = 1 と仮定して Newton の運動方程式を適用すると, 単振動の方程式

y 00 + y = 0を得る.
µ´

y(0) = 1, y 0 (0) = 0 を初期条件として y(t) = cos(t) がこの方程式の解であるが, この方程式を数値解法で解こう. 数値解法の利点は, cos や sin で解を書けないときでも, 微分方程式の近似解がわかることである. 式 (??) より, h が十分小さいとき,y(t + h) − 2y(t) + y(t − h)+ y(t)h2は大体 0 に等しい. したがって,y(t + h) = 2y(t) − y(t − h) − h2 y(t)が近似的になりたつとしてよいであろう. この式は時刻 t − h と t の y の値で, すこし先の時刻 t + h の y の値を表す式である. また y 0 (0) は y(h)−y(0)にほぼ等しい. y 0 (0) = 0 なhので, y(h) = y(0) が近似的に成り立つとしてよいであろう. したがって, 漸化式

yk+2 = 2yk+1 − yk − h2 hk+1 , y0 = 1, y1 = 1(3)を満たす数列を決めることにより, 解の近似をもとめることが可能であると予想できる.
つまり yk は hk 秒での y の値を近似していると予想される. このように解の近似を求める方法を 差分法 とよぶ. 漸化式 (3) を元の微分方程式の 差分化, または差分スキームとよぶ.
¶³! 振動の方程式 y’’+y = 0
LET x1=1
LET x2=1
LET dt=0.01 ! これは十分小さい数なら何でもよい.
! y(k+2), y(k+1), y(k) が x3, x2, x1 に対応.
! h が dt に対応.
FOR t=0 TO 3 STEP dt
LET x3 = 2*x2-x1-dt*dt*x2
PRINT t,x1
PLOT LINES: t,x1;
LET x1=x2
LET x2=x3
NEXT tENDµ´

上のプログラムで ! で始まる行は注釈行 (コメント行) であり, プログラムの実行には関係ない. ! から始まる注釈は行の途中から書きはじめてもよい. この場合は ! から行末までが注釈となる. 適宜注釈行を書くことにより, 人間がプログラムを読む助けとなる.
さて, これで微分方程式を近似的に解く問題が, 漸化式をみたす数列を求める問題になったのであるが, このような近似解が本当の解に収束するかとか, 全くことなる解しかとらえられない不安定現象が起こる場合があるとかの議論をやらないといけない. これはより上級の話題である.
注意: 漸化式は常に微小な数 h を含む. この h は十分小さければどんな数でも構わない.
たとえば 0.01 とか 0.001 など. 近似解が本当の解に収束するという議論では, h を小さくしていけば差分法で求めた近似解が本当の解に収束して行く (どんどん近くなる) ということを証明する.

問 11.1
h が十分小さいとき次の近似公式が成り立つ.

dyy(T + h) − y(T )(T ) 'dthd2 yy(T + h) − 2y(T ) + y(T − h)(T ) 'dt2h2
この近似公式を用いて次の微分方程式を差分化しなさい. (微分方程式の解を近似する数列の漸化式を求めよ.) 数列の漸化式を用いて数列の値を決めて行くプログラムを書きな
さい.
y 00 + 3y = 0, y(0) = 1, y 0 (0) = 0.
さてこのあたりで読者には徹底的に復習をすることをすすめる. (1) 微分方程式 y 00 +y = 0,y(0) = 1, y 0 (0) = 0 の差分化を独力で導けるだろうか? (2) 上のプログラムを自力で書けるだろうか? (1), (2) ができたら次へ進むとよいであろう. 孔子が論語の中でいっているように, “学びて時にこれを習う. またよろこばしからずや” である.

3. 第 3 回の講義では, 物体の落下を表す微分方程式 y 00 = −9.8 について学んだ. (a) 次のプログラムを入力, 実行せよ. どのようなグラフとなるか? (b) 物体が速度に比例する抵抗を受ける場合, 方程式は y 00 = −9.8 − ay 0 となる (a は適当な定数). これを解くプログ
ラムに変更せよ.

¶³100 ! 方程式 y’’=-9.8, y(0)=100, y’(0)=0
110 SET WINDOW -1,5,-10,110
120 DRAW axes
130 LET y1 = 100
140 LET y2 = 100
150 LET dt=0.001
160 ! y(k+1), y(k), y(k-1) が y3, y2, y1 に対応. h が dt に対応.
170 FOR t=0 TO 5 STEP dt
180
LET y3 = 2*y2 - y1 -9.8*dt*dt
190
PLOT LINES: t,y1
200
LET y1 = y2210
LET y2 = y3220 NEXT t
230 ENDµ´

この問題も, 上と同じ方針で解ける. つまり,y 00 (t) = −9.8 − ay 0 (t)を とりあえず差分化してみればよい. 差分化の計算を記す.

y 00 (t) 'よりy(t + h) − 2y(t) + y(t − h) 0y(t + h) − y(t) 00, y (t) ', y (t) = −9.8 − ay 0 (t)2hhy(t + h) − 2y(t) + y(t − h)y(t + h) − y(t)' −9.8 − a2hh

両辺に h2 を掛けて,y(t + h) − 2y(t) + y(t − h) ' −9.8h2 − ah [y(t + h) − y(t)]移項すると(1 + ah)y(t + h) ' (2 + ah)y(t) − y(t − h) − 9.8h2y(t + h) を y3, y(t) をy2, y(t − h) を y1 と置き換え, ' を = とおけば,(1 + a ∗ h) ∗ y3 = (2 + a ∗ h) ∗ y2 − y1 − 9.8 ∗ h2
したがって 180 行を次のように書き換えればよい.
¶³180 LET y3 = ((2+a*dt)*y2- y1 - 9.8*dt*dt)/(1+a*dt)µ´

12Debug の仕方

12.1

print 文をはさむ

12.2

自分が計算機となったつもりになって下のメモリに格納されたデータを書き変えていく

漸化式

fk+2 = fk+1 + fk ,f 1 = f2 = 1

できまる数列 (フィボナッチ数列) の 100 項まで計算して表示するプログラム.
10 a=120 b=130 for k=1 to 100
40
print a50
c = a+b60
a = b70b = c80 next k

アドレス

変数の名前

内容

内容

kabc1
1122123

内容

内容

内容

1. “内容” は 50 行目の実行が終了した時の各変数の値.
問題: 上の図の空白をうめよ.
問題: 他のプログラムについても同様な図を作成してプログラムの実行を追いなさい.
参考:1. アドレス (メモリの番地) は Basic のプログラムを読むときは必要ないが, C のプログラムを読むときは必須. 特にポインターを利用するとき.
2. C のプログラムを読むときは, 変数のサイズ, つまり 2 進数何桁の数を格納できるかにも注意.

dim を用いた実技試験問題の解答例. jitugi1b.bas
V = LdI/dt (コイル), I = CdV /dt (condenser), V = RI (抵抗). LC 共振回路.