3.9.7 微分方程式の数値解法

1.6.4で説明した NDSolveを使うと，微分方程式を数値解析的に解くことができる．また，単独の微分方程式でも連立の微分方程式でも NDSolveで解ける．さらに，常微分方程式でも偏微分方程式でも扱える．解こうとする系が微分方程式ならば，未知の関数 が方程式にいくつあっても構わないが．すべての は単一の独立変数 xのみに依存しているものとする．偏微分方程式の場合は，複数の独立変数を含んでいる．NDSolveは微分方程式と代数方程式が混ざった「微分代数方程式」も扱うことができる．

微分方程式の数値解法

NDSolveで求まる関数 の解は，InterpolatingFunctionオブジェクトを使って表される．このオブ ジェクトは，独立変数 xが xminから xmaxで指定される範囲にあるとき の近似を与える．

NDSolveでは解を求めるために反復手法が使われる．それは xのある値を始点とし，一連の刻みを取りながら xminから xmaxの範囲をくまなく探索していく．

探索を開始するには， NDSolveの式の中に と の導関数に関する適切な初期条件や境界条件が備わっている必要がある．これらの条件を与えることで， xの特定点における [x]，さらに導関数 '[x]の値を決定していく．少なくとも常微分方程式において，この条件は xのどの点でもよい．それは NDSolveが xminから xmaxの全範囲をくまなく自動的に探すからである．

yについて解く． y[0]の初期条件を与え， xは 0から 2の範囲で探索する．

In[1]:= NDSolve[{y'[x] == y[x], y[0] == 1}, y, {x, 0, 2}]

Out[1]=

xの探索範囲は 0から 2のままだが，初期値は y[3]について与えている．

In[2]:= NDSolve[{y'[x] == y[x], y[3] == 1}, y, {x, 0, 2}]

Out[2]=

簡単な境界値の問題を解く．

In[3]:= NDSolve[{y''[x] + x y[x] == 0, y[0] == 1, y[1] == -1},
y, {x, 0, 1}]

Out[3]=

NDSolveを使う際は， の解を完全に決定するのに十分な初期条件，または，境界条件を与えておく必要がある． DSolveを使って微分方程式の厳密解を求めるときには，初期条件が足りなくてもどうにかなる．これは，指定しなかった初期条件のため発生する自由度を明示化するため，積分定数 C[i]が自動的に挿入されるからである．ところが， NDSolveでは解は数値として与えなければならない．数値では余分に発生する自由度が表せない．このため，始めから解を決定するのに十分な初期条件，または，境界条件を与えなければならない．

階の微分方程式を解くとき，通常，最大で 次までの導関数について初期条件を与えるか，または， 点における境界条件を与える必要がある．

3階の微分方程式を解く．3階なので，最大で2次までの導関数を与えなければいけない．

In[4]:= NDSolve[
{ y'''[x] + 8 y''[x] + 17 y'[x] + 10 y[x] == 0,
y[0] == 6, y'[0] == -20, y''[0] == 84},
y, {x, 0, 1} ]

Out[4]=

得られた解をプロットする．

In[5]:= Plot[Evaluate[ y[x] /. % ], {x, 0, 1}]

Out[5]=

別の3階微分方程式を解く．境界条件を3点について与える．

In[6]:= NDSolve[
{ y'''[x] + Sin[x] == 0,
y[0] == 4, y[1] == 7, y[2] == 0 }, y, {x, 0, 2}]

Out[6]=

Mathematicaでは，線形結合した関数値や微分値なら何でも境界条件として使える．

In[7]:= NDSolve[{ y''[x] + y[x] == 12 x,
2 y[0] - y'[0] == -1, 2 y[1] + y'[1] == 9},
y, {x, 0, 1}]

Out[7]=

多くの場合，与える初期条件は必ず同一の x値，例えば を含む必要がある．このため， xminから xmaxの範囲は別に与えなくてもよくなる． xの範囲を x, として指定すると，自動的にからの範囲で有効な解が生成される．

0〜2の範囲で有効な解を求める．

In[8]:= NDSolve[{y'[x] == y[x], y[0] == 1}, y, {x, 2}]

Out[8]=

初期条件とする方程式はどんな種類のものでもよい．場合によっては，これらの方程式は複数の解を持つことがある．そのような場合は，それに対応した形で NDSolveが複数の解を生成する．

この初期条件は，複数の解を与える．

In[9]:= NDSolve[{y'[x]^2 - y[x]^2 == 0, y[0]^2 == 4},
y[x], {x, 1}]

Out[9]=

求まる解をすべてプロットする．

In[10]:= Plot[Evaluate[ y[x] /. % ], {x, 0, 1}]

Out[10]=

NDSolveを使って，結合された微分方程式の系を解くことも可能である．

結合した方程式の対を数値的に解く．

In[11]:= sol = NDSolve[
{x'[t] == -y[t] - x[t]^2, y'[t] == 2 x[t] - y[t],
x[0] == y[0] == 1}, {x, y}, {t, 10}]

Out[11]=

求まった yの解をプロットする．

In[12]:= Plot[Evaluate[y[t] /. sol], {t, 0, 10}]

Out[12]=

xと yの両方を使いパラメトリックプロッ トをする．

In[13]:= ParametricPlot[Evaluate[{x[t], y[t]} /. sol],
{t, 0, 10}, PlotRange -> All]

Out[13]=

微分方程式に未知の関数を使うとき，関数の名前として必ずしも単一のシンボルを与える必要はない．未知の関数が複数個あるとき等は，例えば y[i]といった関数名にするとよい．

5つの結合した微分方程式と初期条件を構築する．

In[14]:= eqns = Join[
Table[ y[i]'[x] == y[i-1][x] - y[i][x], {i, 2, 4} ],
{y[1]'[x] == -y[1][x], y[5]'[x] == y[4][x],
y[1][0] == 1},
Table[ y[i][0] == 0, {i, 2, 5}]
]

Out[14]=

この方程式を解く．

In[15]:= NDSolve[eqns, Table[y[i], {i, 5}], {x, 10}]

Out[15]=

解をプロットする．

In[16]:= Plot[ Evaluate[Table[y[i][x], {i, 5}] /. %],
{x, 0, 10} ]

Out[16]=

NDSolveは値がリストや配列である関数も扱える．y[0] == , , ... , ]のような初期条件を与えると，NDSolveはyを値が長さnのリストである関数だと仮定する．

これは4つの微分方程式が組み合された系の解を求める．

In[17]:= NDSolve[{y''[x] == -Table[Random[], {4}, {4}] . y[x],
y[0] == y'[0] == Table[1, {4}]}, y, {x, 0, 8}]

Out[17]=

これが解である．

In[18]:= With[{s = y[x] /. First[%]},
Plot[{s[[1]], s[[2]], s[[3]], s[[4]]}, {x, 0, 8},
PlotRange -> All]]

Out[18]=

NDSolveの特別オプション

NDSolveには方程式を解く多くのメソッドが備わっているが，基本的にはすべてのメソッドが独立変数 x内で一連のステップを取ることで解を求める．この手順は適応型であり，ステップ幅は常に適当であるように自動調整される．値が特定の領域で急激に変化するときは，解を見逃さないようにステップ幅が徐々に狭められるか，あるいはメソッドが変更される．

導関数に不連続点がある微分方程式を解いてみる．

In[19]:= NDSolve[
{y'[x] == If[x < 0, 1/(x-1), 1/(x+1)],
y[-5] == 5},
y, {x, -5, 5}]

Out[19]=

NDSolveの解の探索において， の近傍ではステップ幅が狭められ，折れ部でも正確に追跡ができている．

In[20]:= Plot[Evaluate[y[x] /. %], {x, -5, 5}]

Out[20]=

NDSolveの探索手順は適応型であるから，変数 xの変化率が大きく異なるいくつかの要素を持つ，いわゆる「扱いにくい」微分方程式を解くこともできる．

これらの方程式において， yは zより急激に変化する．

In[21]:= sol = NDSolve[
{y'[x] == -40 y[x], z'[x] == -z[x]/10,
y[0] == z[0] == 1},
{y, z}, {x, 0, 1}]

Out[21]=

それでも，両方とも NDSolveで正確に追跡できている．

In[22]:= Plot[Evaluate[{y[x], z[x]} /. sol], {x, 0, 1},
PlotRange -> All]

Out[22]=

NDSolveで使われる手順では，解を確実に追跡できるまでステップ幅を徐々に狭めている．ただし，この手順では，真の解に特異点が存在すると問題が起る．つまり，ステップ幅の細分化が限りなく行われてしまい，探索が終了しなくなってしまう．この問題を避けるため， NDSolveに許す最大ステップ数をオプションMaxStepsで指定しておく．常微分方程式の場合，このオプションはデフォルトで MaxSteps -> 10000になっている．

NDSolveは10000ステップに達した後に停止する．

In[23]:= NDSolve[{y'[x] == -1/x^2, y[-1] == -1}, y[x], {x, -1, 0}]

Out[23]=

確かに， において解は特異点を持つ．

In[24]:= Plot[Evaluate[y[x] /. %], {x, -1, 0}]

Out[24]=

解が滑らかなほとんどの関数では，MaxStepsのデフォルト値で十分間に合う．解が複雑な構造を取るようなときは， MaxStepsを大きく取っておくとよいだろう．また， MaxSteps -> Infinityと設定しておくとステップの上限なしで探索ができる．

ローレンツ(Lorenz)方程式の解をここまで求めるためには， MaxStepsのデフォルトによる制限を除いておく必要がある．

In[25]:= NDSolve[ {x'[t] == -3 (x[t] - y[t]),
y'[t] == -x[t] z[t] + 26.5 x[t] - y[t],
z'[t] == x[t] y[t] - z[t],
x[0] == z[0] == 0, y[0] == 1},
{x, y, z}, {t, 0, 200}, MaxSteps->Infinity ]

Out[25]=

求まった解を3Dパラメトリック形式でプロットする．

In[26]:= ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. %],
{t, 0, 200}, PlotPoints -> 10000]

Out[26]=

連立した微分方程式を解くとき， NDSolveではすべての方程式に適応するようにステップ幅が自動調整される．しかし，場合によっては，最初に取るステップ幅が大きすぎて解の重要な特徴を見逃しかねない．このような問題を避けるには，オプション StartingStepSizeを使い第1ステップの幅を特別に指定しておく．

微分代数方程式の数値解を求める

NDSolveに与える方程式すべてが導関数を含んでいる必要はない．単に代数的なだけでもよい．NDSolveを用いて多くの「微分代数方程式」を解くことができる．

これは微分代数方程式系を解く．

In[27]:= NDSolve[{x'[t] == y[t]^2 + x[t] y[t],
2 x[t]^2 + y[t]^2 == 1,
x[0] == 0, y[0] == 1}, {x, y}, {t, 0, 5}]

Out[27]=

これが解である．

In[28]:= Plot[Evaluate[{x[t], y[t]} /. %], {t, 0, 5}]

Out[28]=

偏微分方程式を数値的に解く

これは波動方程式の数値解を求める．結果は2次元の補間関数である．

In[29]:= NDSolve[{D[u[t, x], t, t] == D[u[t, x], x, x],
u[0, x] == Exp[-x^2], Derivative[1,0][u][0, x] == 0,
u[t, -6] == u[t, 6]}, u, {t, 0, 6}, {x, -6, 6}]

Out[29]=

結果のプロットを生成する．

In[30]:= Plot3D[Evaluate[u[t, x] /. First[%]],
{t, 0, 6}, {x, -6, 6}, PlotPoints->50]

Out[30]=

これで非線形波動方程式の数値解を求める．

In[31]:= NDSolve[
{D[u[t, x], t, t] ==
D[u[t, x], x, x] + (1 - u[t, x]^2)(1 + 2u[t, x]),
u[0, x] == Exp[-x^2], Derivative[1, 0][u][0, x] == 0,
u[t, -10] == u[t, 10]}, u, {t, 0, 10}, {x, -10, 10}]

Out[31]=

これは結果の3Dプロットである．

In[32]:= Plot3D[Evaluate[u[t, x] /. First[%]],
{t, 0, 10}, {x, -10, 10}, PlotPoints->80]

Out[32]=

これは結果の高解像度密度プロットである．

In[33]:= DensityPlot[Evaluate[u[10 - t, x] /. First[%%]],
{x, -10, 10}, {t, 0, 10},
PlotPoints -> 200, Mesh -> False]

Out[33]=

これはこの方程式の2+1次元におけるバージョンである．

In[34]:= eqn = D[u[t, x, y], t, t] == D[u[t, x, y], x, x] +
D[u[t, x, y], y, y]/2 + (1 - u[t, x, y]^2)(1 + 2u[t, x, y])

Out[34]=

これで方程式を解く．

In[35]:= NDSolve[{eqn, u[0, x, y] == Exp[-(x^2 + y^2)], u[t, -5, y] == u[t, 5, y],
u[t, x, -5] == u[t, x, 5], Derivative[1, 0, 0][u][0, x, y] == 0}, u,
{t, 0, 4}, {x, -5, 5}, {y, -5, 5}]

Out[35]=

これは解のプロットの配列を生成する．

In[36]:= Show[GraphicsArray[
Partition[
Table[ Plot3D[Evaluate[u[t, x, y] /. First[%]],
{x, -5, 5}, {y, -5, 5}, PlotRange -> All,
PlotPoints -> 100, Mesh -> False,
DisplayFunction -> Identity], {t, 1, 4}], 2]]]

Out[36]=