微分方程式の数値解
（１）１階常微分方程式
前講で扱った数値積分
t
x = ∫ f(s) ds
a
は，微分方程式
dx/dt = f(t) , x(a) = 0
を満たす関数 x(t) を求めることであったと言ってもよい。微分方程式の言葉で言えば，変数分離型の最も簡単な場合である。一般の微分方程式では，右辺にも未知数 x が含まれており，
dx/dt = f(x, t)
の形になる。これを１階の常微分方程式という。(t, x) 平面上の各点で勾配の値 f(x, t) が与えられているから，出発点（初期値）を決めて次々にこの勾配の方向を繋いでいけば曲線が得られる。この曲線を表す関数 x(t) が，この微分方程式の解である。
（オイラー法） この解を数値的に求めるには，定積分のときと同様に，独立変数 t の変域を微小な幅 △t = h の区間に刻み
tn = t0 + nh
として
xn+1 = xn + h f(xn, tn)
で近似する。これを差分化という。右辺は常にそれまでにわかった量で表されているから，漸化式の形をしたこの差分方程式を次々に適用していけば，微分方程式の近似解を得ることができる。プログラムとしては単に
t = t0 + n*h
x = x + h*f(x, t)
と書くだけである。
差分による方法では，次に出てくる２つの方法を含めていずれも，xn+1 は xn と tn の関数を用いて
xn+1 = xn + h F(xn, tn)
の形に書ける。オイラー法ではもちろん，F(x, t) = f(x, t) である。１回の前進における誤差の目安として，tn における微分方程式の真の解 x(tn) を用いた差分解
x(tn) + h F(x(tn), tn)
と，tn+1 における真の解 x(tn+1) との差
| x(tn+1) - [x(tn) + h F(x(tn), tn)] |
あるいはこれを h で割った
δ = | [x(tn+1) - x(tn)]/h - F(x(tn), tn) |
を局所打ち切り誤差という。オイラー法では，テーラー展開を用いればδ〜h の程度であることがわかり，「１次近似」である。誤差は差分を繰り返すにしたがって次々に累積していき，前進を N 回行ったときの最大の誤差は Nhδ 程度であるから，δとしては h で割らない量を用いる方が便利かもしれない。
（ホイン法） オイラー法による行き先の点での勾配を求め，元の点に一度もどってこの勾配での差分を求め，平均値を採用する方法をホイン法といい，２次の近似となる。すなわち
f1 = f(xn, tn )
f2 = f(xn + h f1, tn + h)
として
xn+1 = xn + h [(f1 + f2)/2]
（ルンゲ-クッタ法） さらに
f1 = f(xn, tn )
f2 = f(xn + (h/2)f1, tn + h/2)
f3 = f(xn + (h/2)f2, tn + h/2)
f4 = f(xn + h f3, tn + h )
として
xn+1 = xn + h [(f1 + 2f2 + 2f3 + f4)/6]
とすれば，手間は素朴なオイラー法に比べておよそ４倍となるが，近似は４次で δ〜h4 となる。
微分方程式の右辺に未知数 x を含まない場合，すなわち f が t だけの関数の場合には，以上の方法はそれぞれ階段関数近似，台形公式，シンプソン公式に対応していることがわかるだろう。
（２）２階常微分方程式と連立１階微分方程式
おなじみの力学で現れる運動方程式は２階常微分方程式
d2x/dt2 = f(x, dx/dt, t)
である。（ただし１次元運動） これは，速度 v = dx/dt （または運動量 p = mv） を独立な未知数とすることにより
dx/dt = v
dv/dt = f(x, v, t)
と，未知数が２つの連立１階常微分方程式の形に書くことができる。解を決めるためには，２つの初期値 x(0), v(0) を与えなければならないことは言うまでもない。
連立方程式
dx/dt = f(x, y, t)
dy/dt = g(x, y, t)
にも前節と同じ方法を適用することができる。例えばルンゲ-クッタ法では
tn = t0 + n h
として
f1 = f(xn, yn, tn )
g1 = g(xn, yn, tn )
f2 = f(xn + (h/2)f1, yn + (h/2)g1, tn + h/2)
g2 = g(xn + (h/2)f1, yn + (h/2)g1, tn + h/2)
f3 = f(xn + (h/2)f2, yn + (h/2)g2, tn + h/2)
g3 = g(xn + (h/2)f2, yn + (h/2)g2, tn + h/2)
f4 = f(xn + h f3, yn + h g3, tn + h )
g4 = g(xn + h f3, yn + h g3, tn + h )
xn+1 = xn + h [(f1 + 2f2 + 2f3 + f4)/6]
yn+1 = yn + h [(g1 + 2g2 + 2g3 + g4)/6]
となる。プログラムで書くときにも，この通りの順番である。
未知数がもっとたくさんある多元連立常微分方程式では，配列（ベクトル）を用いる方が便利であろう。
dx[i]/dt = f(i, x[:], t) (i = 1, 2, ..., n)
以下では x[:] には配列要素 x[1], x[2], ... , x[n] を並べる。
または，Fortranの場合は配列引数 x(:) が書かれる。
に対して，プログラム風に書けば
t = t0 + n*h
f1[i] = f(i, x[:], t ) (i = 1, 2, ..., n)
z[i] = x[i] + (h/2)*f1[i] (i = 1, 2, ..., n)
f2[i] = f(i, z[:], t + h/2) (i = 1, 2, ..., n)
z[i] = x[i] + (h/2)*f2[i] (i = 1, 2, ..., n)
f3[i] = f(i, z[:], t + h/2) (i = 1, 2, ..., n)
z[i] = x[i] + h*f3[i] (i = 1, 2, ..., n)
f4[i] = f(i, z[:], t + h ) (i = 1, 2, ..., n)
x[i] = x[i] + h*(f1[i]+2*f2[i]+2*f3[i]+f4[i])/6 (i = 1, 2, ..., n)
となる。全行をまとめて i = 1, 2, ..., n と繰り返すのではなく，各行ごとに i = 1, 2, ..., n と繰り返して次の行へ行く。 一時的な変数 z を x と書いてしまったら間違いであることに注意せよ。
（３）差分力学系
微分方程式を差分化した場合，時間刻み幅 h をどの程度にとったらよいかを決める原理はなく、大きくとりすぎると振動して発散することさえ起きる。ある現象を探るとき，必ずしも法則を微分方程式で表現する必要はない。数値計算することを考えると、最初から差分方程式で考察しておく方が有利な場合もある。
例えば、１年で世代交代する昆虫の個体数の推移を想像してみよう。来年の個体数 x(n+1) は今年の個体数（親の数） x(n) だけで決まり
x(n+1) = F(x(n))
で与えられるとしよう。この写像関数 F(x) は、現在の個体数に比例して増えるプラス要素と、個体数が増えすぎると生息環境の悪化により生殖力が減じ、限界に達すると絶滅しまうマイナスの要素の重ねあわせであると考えられ、最も簡単なモデルとして
F(x) = A x ( 1 - x ) ( A≦４ ）
のような形で表現することができる。x=1 に達すれば絶滅するように規格化してある。この場合の差分方程式は、A＞1 のとき x≠0 の不動点（時間的に変動しない定常解） x=1 - 1/A をもち、この点での勾配 F'(1-1/A)=2-A の値から、A＜3 のとき安定、すなわち時間がたてばこの値に収束していくことがわかる。（ノートNo.3参照）
A＞３ では不安定で、この点の周りを振動しながら離れていくが、もう一つの不動点 x=0 も不安定であり、A≦４ としておけばx は 0＜x＜1 の間に制限されるから、不規則な振舞いをすることが予想される。特別な初期値をとれば、何世代かで元の値に戻る周期運動になるが、一般には不規則な無周期運動（カオス）となる。
（付１）差分解の安定性
ここまでに紹介した近似法の意味をもう一度，手で計算できる例で見てみよう。最も簡単な線形微分方程式
dx/dt = - x , x(0) = 1
の解は x(t) = e-t である。オイラー法では
xn+1 = xn - h xn = (1-h) xn
したがって
xn = (1-h)n = exp[n log(1-h)] = exp[-nh(1 + h/2 + h2/3 +...)]
厳密解との比は
xn/exp(-nh)=exp[n(log(1-h)+h)]=exp[-nh(h/2+h2/3+...)]〜1-th/2
となり，与えられた h に対して，どの程度の時間 t (= nh) まで有効かがわかる。
同様にしてホイン法では
xn = (1-h+h2/2)n = exp[n log(1-h+h2/2)]
xn/exp(-nh)=exp[-nh(h2/6+3h3/8+...)]〜1 - th2/6
ルンゲ-クッタ法では
xn = (1-h+h2/2-h3/6+h4/24)n = exp[n log(1-h+h2/2-h3/6+h4/24)]
xn/exp(-nh)=exp[-nh(h4/120+h5/144+...)]〜1 - th4/120
となる。
ここで，オイラー法の（手計算の）結果
xn = (1-h)n
において，公比が
|1 - h|＜1
を満たしていなければならないことは明かであろう。この例では h を
1 - h ＜ -1 すなわち h＞2
ほど大きく取ることは，まず考えられない。普通は素朴に考えても h＜１にとるからである。
しかしながら次のような例（過減衰振動）では，うっかりするとこの制約を見逃してしまうこともあり得る。
dx/dt = y , dy/dt = -16x - 10y , x(0) = 1 , y(0) = 0
厳密解は
x(t) = (4/3) e-2t - (1/3) e-8t
近似法も手で解くことができて，オイラー法では
xn = (4/3)(1-2h)n - (1/3)(1-8h)n
となり，収束条件は
|1-2h|＜1 , |1-8h|＜1 より 0＜h＜1/4
であって， h＞1/4 では差分操作ごとに激しく振動して発散してしまい，h＜1 だけでは不十分であることがわかる。この制約はホイン法，ルンゲ-クッタ法でもほぼ同じであり，近似の次数が上がることとは関係がないので注意する必要がある。
このように手計算で調べることのできるものでは予め近似の精度や収束条件を知ることができるが，実際に数値計算を必要とする例では一般にこうはいかない場合がほとんどである。残念ながら経験的にチェックするしか方法がないのが実状である。
（１）グラフに描いて変な振舞いがないか目で確かめる。
（２）微分解の漸近的な振舞いを調べておき比較する。
（３）h の値を何種類かとって結果が安定しているかどうか調べてみる。
など。
それでもなおかつ振動や発散が残る場合は、微分解そのものに不安定点が存在している可能性がある。そのような場合には，その近傍では近似の精度に十分な注意をはらうべきである。
最初から差分方程式で記述されるような現象では，（３）で見た例のように，不規則に振動する結果そのものが現象の本質を表していることもある。（カオス現象）
（付２）保存力学系
摩擦のない運動方程式のような保存量のある「保存力学系」では，差分化によって保存性が損なわれることは致命的である。例えば次のような１次元調和振動子系を考えてみる。
dx/dt = p
dp/dt = -x
この方程式は、あきらかにエネルギーに対応する保存量 x2+p2 をもつ。これをオイラー法で近似した場合
xn+1 = xn + h pn
pn+1 = -h xn + pn
したがって
xn+12 + pn+12 = (1 + h2) (xn2 + pn2)
となり，エネルギーは単調に増加し，保存則は損なわれる。この事情はルンゲ-クッタ法でも変わらず
xn+12 + pn+12 = (1 - h6/72 + h8/576) (xn2 + pn2)
で，オイラー法に比べれば遅いが今度は単調に減少する。
このような場合のために，保存性を尊重した差分法が工夫されている。（シンプレクティック法−−−早川先生の講義録を参照）
＜課題４−１＞
減衰振動の方程式
d2x/dt2 + 10 dx/dt +169 x = 0 初期条件 x(0)=1.2 , v(0)=0.0
の解を，t=0 から t=2.5 までオイラー法とルンゲ-クッタ法を用いて数値積分し，分かっている解
x(t) = exp(-5t) [ 0.5 sin(12t) + 1.2 cos(12t) ]
と比較せよ。数値的な比較は、差分化するときの時間幅 h より大きめの時間幅Δｔで定期的に行えばよい。グラフィックツール（パソコンならエクセルで可能）を使える人はグラフで表して比較してみよ。
＜課題４−２＞
ある化学反応系の２成分の濃度の変化を表現する非線形連立常微分方程式
dx/dt = A - Bx +x2y - x
dy/dt = Bx - x2y
ただし
A = 1.0 , B = 3.0
は，時間が十分たつと，ある（いびつな）周期的軌道に収束することが知られている。このような周期軌道をリミットサイクル（極限閉軌道）という。(0.0, 0.0) から出発した場合は閉軌道の外側から，また，不安定な定常点 (1.0, 3.0) のすぐ近くから出発した場合は内側から収束する。この二つの場合を比べながらリミットサイクルの（時間）周期を求めよ。座標データを眺めてほぼ周期的になっていることが確かめられたら，周期はｘ座標、ｙ座標の極大・極小などを用いてエイヤッと手計算で決めればよい。差分幅 h の程度の精度で求めることが可能であろう。ただし， 座標の出力数値データは膨大になるから印刷したり提出したりしないこと。できればリミットサイクルに近づいていく様子をグラフで描け。
＜課題４−３＞
(3)差分力学系で与えた例について，A＝4.0 のとき、適当な初期値（なるべく、いかにも有理数でないような値がよい）を与えて振舞いを調べてみよ。(0,1) 間を0.05刻みの区間に20等分しておき、x(n) の値の出現頻度のヒストグラムを描いてみよ。注 必ず「倍精度」で計算すること。
ほとんどの初期値に対して、この分布は N を差分回数（繰り返し回数）として
N/π[x(1-x)]1/2
で与えられることが知られている。このことを予め知った上でヒストグラムを設計してよい。