微分方程式と物理
森 真マルサスの微分方程式
人口は現在の人口に比例して増加するということでたてられたモ デルです． xは人口を表すとすると，その変化(微分)は人口xに比例すると 考え dx dt = ax が成り立つと考えました．この方程式は 1 x dx = a dt と変形してから，両辺を積分すると log x = at + c 書き直して x (t) = eatC (C = ec) と解けます．こういう方程式を変数分離形と言います．解はa > 0 ならば，t → ∞で無限に増えていき，a < 0ならば，t → ∞で0 に収束します．マルサスは食料は線形(at + bの形)でしか増産で きないので，人類はいずれ滅びると予想したのです．これがマル サスの人口論と言われ，真偽はともかく，多くの議論を生みまし た．数学が社会科学に用いられた初めての例ではないでしょうか．
コンピュータで求める方法
それには，微分の定義に戻る必要があります． dx dt = limh→0 x (t + h)− x(t) h ですから，hが十分小さければ x (t + h)− x(t) h はx (t)の時刻tでの微分に近いと考えられます．マルサスの方程 式は x (t + h)− x(t) h + ax(t) つまり x (t + h)+ x(t) + ahx(t) が成り立ちます．時刻を0, h, 2h, 3h, . . .と考えて，次々に代入し ていけば解が近似できます．Mathematica
で解く
a=1;h=0.1;
For[i = 1; x = 1; kinji ={{0, x}};, i <= 10, i++, x += a*x*h; AppendTo[kinji,{i*h, x}]];
ListPlot[kinji,PlotJoined−>True]
解の図
微分方程式 dx dt = x Figure: 初期値 1 と−1 のときのオイラー法による近似 0.2 0.4 0.6 0.8 1 -2 -1 1 2問題
人類が絶滅するのは嫌だというので，人口が増えると，増加が鈍 る項を加えて，マルサスの方程式を改良した人口モデルのロジス ティック方程式 dx dt = ax− bx 2 の解をオイラー法で近似してください．問題点
オイラー法では，より正確にするためにはhを小さくすればいい のですが，そうすると，計算の回数が増えます．コンピュータは 四捨五入をしますから，回数が増えると誤差が増える(丸め誤差) 危険があります．これを回避するには，もっと近似のよい式を使 う必要があります． Figure: h = 0.1のときの，真の解との差 (赤が真の解) 0.2 0.4 0.6 0.8 1 1.25 1.5 1.75 2 2.25 2.5 2.752
次元の線形微分方程式
2次元の微分方程式で dx dt = ax + by dy dt = cx + dy の形のものを，線形微分方程式といいます．行列による表現
dx dt = ax + by dy dt = cx + dy は A = ( a b c d ) とおくと d dt ( x y ) = A ( x y ) と変形できることにまず注目しましょう．対角化
Aが対角化P−1AP = ( λ 0 0 µ ) できれば， ( X Y ) = P−1 ( x y ) とお くと d dt ( X Y ) = d dtP −1(x y ) = P−1A ( x y ) = P−1AP ( X Y ) により， dX dt = λX dY dt = µY と2つのマルサスの微分方程式にわかれますので，解くことがで きます．対角化できる微分方程式
dx dt = 7 3x− 2 3y dy dt = − 1 3x + 8 3y を解こう．ここでは A = ( 7 3 − 2 3 −1 3 8 3 ) とおいて対角化すればよい．固有値，固有ベクトル
A ={{7/3, −2/3}, {−1/3, 8/3}} とおいて，固有値と固有ベクトルを求める． Eigenvalues[A] Q=Eigenvectors[A] Qは固有ベクトルが横に並んでいるので，縦に並べ直して P=Transpose[Q]対角化
Inverse[P].A.P とすれば，対角行列 ( 3 0 0 2 ) を得る．すなわち ( X Y ) = P−1 ( x y ) とおけば dX dt = 3X dY dt = 2Y になる．問題点
もちろん対角化できない場合も困りますが，それ以上に固有値が 複素数になってしまう固有値が複素数になる微分方程式 dx dt = x− y dy dt = x + y では上の方法が使えません． こんな場合でも，オイラー法で微分方程式を近似することができ ます．Mathematica
による解法
h = 0.1;
For[i = 1; x = 0.1; y = 0; kinji ={{x, y}}, i <= 67, i++, x += h(x - y); y += h(x + y); AppendTo[kinji,{x, y}]];
ListPlot[kinji, PlotJoined−> True, Ticks −> None]
この解法には誤りがある
方法に誤りはなさそうですが，新しいxをyを求めるのに使って しまっています． x += h(x -y); y += h(x + y); この場合には，xやyは大きくは変わりませんので，大した問題 にはなりませんが，正しくは，一旦，xをしまっておいて，全部 求めてからもとへ戻せばいいのです．Mathematica
による正しい解法
xを新しくしたxはuとしてしまっておいて，新しいyを求める
ときにはとっておいたxを使えばいいわけです．その後，uをx
に戻します．h = 0.1;
For[i = 1; x = 0.1; y = 0; kinji ={{x, y}}, i <= 67, i++,
u=x+ h(x - y); y += h(x + y);x=u; AppendTo[kinji, {x, y}]]; ListPlot[kinji, PlotJoined−> True, Ticks −> None]
問題
A = ( 1 2 1 1 ) の固有値と固有ベクトルを求め，さらにAを対角化してくださ い．さらにこの行列に対応する微分方程式 dx dt = x + 2y dy dt = x + y の解をオイラー法で求めてください．物理と微分方程式
ニュートンの力学は I 等速直線運動 I F = ma，力は加速度に比例する I 作用反作用の法則 というたった3つの原理に基づいています． これと万有引力の法則 F = GmM r2 からりんごの落下から太陽と惑星の方程式まで導きだしたのです．落下の方程式
物体は下に向かって，mg の力で引っ張られます． md 2y dt2 =−mg これは dy dt = v dv dt = −g と速度vも入れると簡単に解けて初期位置y0，初期速度v0y とお けば v (t) = v0y − gt y (t) = y0+ v0t− 1 2gt 2 となる．放物線
x方向も考えると，水平には力が働いていないので d2x dt2 = 0 になります．x方向の初期位置x0，初期速度v0xとすると x (t) = x0+ v0xt となります，これからtを消去すると y (t) = y0+ v0y v0x(x (t)− x0)− 1 2g ( x (t)− x0 v0x )2 と放物線が現れます．Mathematica
で描く
x , y方向それぞれに速度があるので4次元になります．そのうち，
x , yだけをとっておいて，図にします．g=1;h = 0.1; For[i = 1; x = 0; y = 0;vx=0.5;vy=20; kinji ={{x, y}}, i <= 100, i++,ux=x+ h*vx; vx += 0; uy=y+h*vy;vy+=-g*y; x=ux;y=uy; AppendTo[kinji,{x, y}]];
バネの方程式
バネの力は自然な状態からの長さの変化に比例する力が働きます (フックの法則)． mdx 2 dt2 =−kx これも dx dt = v dv dt = −kx と速度vも入れると，見通しが良くなりますが落下の方程式のよ うに解くことはできません．しかし，勘がよければ， x (t) = A sin(Bt + C )のような形をしていることに気がつくはず で，そうなると解は x (t) = A sin(√kt + C ) v (t) = A cos(√kt + C ) であることがわかり，A, Cは初期位置と初期速度から求めること ができます．ふりこの方程式
これまでの方程式は線形でしたから，線形代数で解くことができ ました．しかし，ふりこでは md 2x dt2 =−mg sin x の形をしています．xが0に近ければsin x∼ xとみることができ てバネの方程式と同じになります．これがふりこが振れ幅によら ず一定のリズムをする理由です．でも正しい解はこれとは違うは ずです．このような方程式を非線形の方程式といいます．解けな いものはコンピュータで近似解を求めるしか方法はありません．問題
ふりこの方程式を解くプログラムを作ってさまざまな初期値につ いて解の挙動を調べてください． 予想通りの行動はしてくれません．これはオイラー法があまり良 い近似でないために時間がたつと正確な解から遠ざかってしまう からです．問題
dx dt = −10x + 10y dy dt = 29x− y − xz dz dt = − 8 3z + xy の解をxy 平面に射影した図を作ってください．ローレンツの方程式
この方程式は気象の長期予報について研究していたローレンツが 考えたものです．この解をみて，彼は長期予報は不可能であると 結論をしました．解の形が蝶に似ていることから，北京で蝶が羽 ばたきをすると，ニューヨークに嵐が起きるなどと象徴的に言わ れています．ローレンツ方程式を
Mathematica
で解く
h = 0.01;
For[i = 1; x = 0.1; y = 0.1; z = 0.1; kinji ={{x, y}}, i < 5000, i++, u=x + h(-10x + 10y); v=y + h(28x - y - x*z); z += h(-8z/3 + x*y); x=u; y=v; AppendTo[kinji,{x, y}] ]
ローレンツ方程式の解
-15 -10 -5 5 10 15 20 -20 -10 10 20ファンデルポルの微分方程式
dx dt = y dy dt = εy (1− x 2)− x は周期解にまつわりついて行く図になります．ファンデルポルの微分方程式
h = 0.1; m = 2;
For[i = 1; x = 0.1; y = 0; kinji ={{x, y}}, i <=1000, i++,
u=x+ h*y;
y += h(m*(1 -x2)*y - x); x=u
AppendTo[kinji,{x, y}]];
ファンデルポルの微分方程式
コンピュータによる解法
オイラー法はx (t + h)をx (t) + x′(t)hで近似する方法です．これ はx (t + h)のテイラー展開 x (t + h) = x (t) + x′(t)h +x ′′(t) 2 h 2+· · · の初めの2項だけを用いたものです．もっとテイラー展開を長く 一致させればより正確な近似が得られるはずです．方程式を y′ = f (x , y ) とします．これを時間刻みh毎に初期値y (0)から導く．xをT ま でなら, 0, h, 2h, . . .とT /h回繰り返すことにしましょう．このm 回目をymで表しましょう．xm = mhです． l段階法ymを得るのに，ymからym−l+1を用います．この方法に は初期値として，y0 = y (0)の他にy1, . . . , yl−1までを別の方法で 求めることが必要になります． ym+1を求めるのにym+1′ = f (xm+1, ym+1)が必要な場合を陰的公 式，必要としないものを陽的公式といいます．陰的公式ではym+1 を解く必要があります．
オイラー法(陽的公式) ym+1= ym+ hym′ (ym′ = f (xm, ym)) 台形公式(陰的公式) ym+1 = ym+ h 2(y ′ m+ ym+1′ ) h ym′ オイラー法 h ym′ ym+1′ 台形法
オイラー法ではxm< x < xm+1の点では y (x ) = ym+ (x− xm)ym′ 台形公式では ym+1 = ym+ h 2y ′ m+ h 2y ′ m+1 = ym+ h 2y ′ m+ h 2(y ′ m+ hym′′ +· · · ) = ym+ hym′ + h2 2 y ′′ m+· · · なので y (x ) = ym+ (x− xm)ym′ + (x− xm) 2 y ′′ m+· · · となります．台形公式はテイラー展開が2次まで一致しているこ とがわかります．しかし，陰的公式なのでym+1′ = f (xm+1, ym+1) なので，ym+1を別の方法で近似するか，式を解かなければなりま せん．
