5.21 微分方程式を解く DSolve[], NDSolve[]

常微分方程式を解くコマンド DSolve[] があります。

一つの使い方は、次のように y[x] について解くというものです。

DSolve[y'[x]==a y[x],y[x],x]

結果は {{y[x]->Exp[a x]C[1]}} となります。 ここから何かするには、演算子 /. を使うことになるでしょう。 次の例では、初期値問題を解いて、結果をグラフに表し、 (少々強引に) 関数 x を上書き定義しています。

s=DSolve[{x''[t]==-x[t], x[0]==1, x'[0]==2},x[t],t]
Plot[x[t] /. s, {t,-10,10}]
x[t_]:=Evaluate[x[t] /. s[[1]]]	← 関数 x を定義
Plot[x[t], {t,-10,10}]
Remove[x,s,t]

連立の微分方程式を解くことも出来ます。

s=DSolve[{x'[t] == -y[t], y'[t]==2x[t], x[0]==1, y[0]==0},{x[t],y[t]},t]
ParametricPlot[{x[t],y[t]}/.s,{t,0,2Pi}]
Remove[s,t,x,y]

もう一つのやり方は、 y[x] でなく、y について解くと言うものです。 結果は (説明は省略しますが) 「純関数」 {{y->Function[{x},2Exp[a x]]}}になる。 つまり y に /. したときに、 y[] という関数が出来る、 つまり y'[x] や y[1] などの式が使える。

s=DSolve[{y'[x]==a y[x],y[0]==2},y,x]
y[x] /. s	← この段階では差がない
y'[x]-a y[x] /. s	← こういうことが出来る(解であることの確かめ)
Remove[a,s,x,x$,y]

初期値問題も解けます。

s=DSolve[{x''[t] == - omega^2 x[t], x[0] == x0, x'[0] == v0},x,t]
Remove[omega,s,t,t0,t$,v0,x,x0]

微分方程式は、 いわゆる式変形では解けない場合が多いことは良く知られています。 数値的に近似解を求める方法が良く知られていますが、 NDSolve[] は数値的に微分方程式を解くコマンドで、 結果は離散的な点での関数値を近似的に求めたものになります。 計算していない点での値は補間により簡略計算することになります。 こういうものを、 Mathematica では InterpolatingFunction と呼んでいるようです。) NDSolve[] の結果をグラフとして描く方法に注目して下さい。

振り子の方程式を解く

s=NDSolve[{y''[x]==-Sin[y[x]],y[0]==1,y'[0]==0}, y, {x,0,10}]
y[1.5] /. s
Plot[Evaluate[y[x] /. s], {x,0,10}]
Remove[s,x,y]