How to | NDSolveの結果をプロットする方法
NDSolveは微分方程式を数値的に解く.解は,さまざまな方法で即座に使用できるような形で返される.典型的な使用例の1つに解のプロットを作成するということがある.
ode1 = {y''[x] + Sin[y[x]] == 0, y[0] == 2, y'[0] == 1};NDSolveで,上の方程式を第1引数,解かなくてはならない関数である
を第2引数,そして独立変数の範囲を第3引数とする:
sol = NDSolve[ode1, y, {x, -10, 10}]Plot[y[x] /. sol, {x, -10, 10}]解をその導関数(あるいは複数の従属変数)と一緒にプロットするというのはよくあることである.それぞれを別の色にしたい場合には,プロットしたいものをEvaluateで囲むとよい:
Plot[Evaluate[{y[x], y'[x], y''[x]} /. sol], {x, -10, 10}]2度および3度の自由度を持つシステムについては,位相面にプロットした方がより多くのことが分かることが多い.これはParametricPlotを使って行うことができる:
ParametricPlot[Evaluate[{y[x], y'[x]} /. sol], {x, -10, 10}]NDSolveは,方程式系を解くことができる:
s = NDSolve[{x'[t] == -y[t] - x[t] ^ 2, y'[t] == 2x[t] - y[t] ^ 3, x[0] == y[0] == 1}, {x, y}, {t, 20}]Plotを使ってその方程式系をプロットすることができる:
Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, 20}]ParametricPlotを使ってプロットすることもできる:
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 20}]Manipulateと一緒に位相面プロットを使うことは簡単で,この方法で初期条件をいろいろ変化させることができる:
Manipulate[
Module[{sol = NDSolve[{y''[x] + Sin[y[x]] == 0, y[0] == p[[1]], y'[0] == p[[2]]}, y, {x, 0, T}]},
ParametricPlot[Evaluate[{y[x], y'[x]} /. sol], {x, 0, T}, PlotRange -> {{-4, 4}, {-3, 3}}]],
{{p, {2, 1}}, Locator}, {{T, 10}, 0, 100}]Locatorを使うと,点をドラッグして初期条件を変更することができるようになる.パラメータTで問題が解かれる区間を制御することができる.
偏微分方程式については,いくつかの選択肢があることが多い.例えば,Wolframの非線形波動方程式 [もっと詳しく]を考えてみる:
wws = NDSolve[{Subscript[∂, t, t]u[t, x] == Subscript[∂, x, x]u[t, x] + (1 - u[t, x]^2) (1 + 2 u[t, x]), u[0, x] == E^-x^2, u^(1, 0)[0, x] == 0, u[t, -10] == u[t, 10]}, u, {t, 0, 10}, {x, -10, 10}]解の全体像をつかみたい場合は,Plot3Dを使うとよい:
Plot3D[u[t, x] /. wws, {t, 0, 10}, {x, -10, 10}]解の詳細についてよい情報が得られることが多いもう1つの方法は,DensityPlotを使う方法である:
DensityPlot[u[t, x] /. wws, {t, 0, 10}, {x, -10, 10}]以下のような時間発展方程式については,アニメーションを使うとよい直観が得られることが多い.一般に一番よい結果はListAnimateを使って得られる.まずすべてが同じPlotRangeで同じ時間間隔離れているプロットのリストを作成する:
plots = Table[Plot[u[t, x] /. wws, {x, -10, 10}, PlotRange -> {-2, 2}], {t, 0, 10, .25}];ListAnimate[plots]この場合は,ゼロの初期条件によって背景が変化するので,波動を見ることは少々難しい.この問題を取り除く簡単な方法は,これに対応する常微分方程式を解いてそれを引くという方法である:
bg = NDSolve[{Subscript[∂, t, t]u[t] == (1 - u[t]^2) (1 + 2 u[t]), u[0] == 0, u'[0] == 0}, u, {t, 0, 10}]plots = Table[Plot[(u[t, x] /. wws) - (u[t] /. bg), {x, -10, 10}, PlotRange -> {-2, 2}], {t, 0, 10, .25}];ListAnimate[plots]