How to | NDSolveの結果をチェックする方法
ほとんどの微分方程式について,NDSolveで求めた結果は正確である.しかし,その結果は数値標本と誤差推定に基づいているため,重大な誤差が起ることが時にはある.解の質を確認したい場合は,解に対して基本的なチェックをいくつか行うとよい.
デフォルトのMachinePrecisionよりも高いWorkingPrecisionで計算した解と比べることは,結果をチェックする上で役に立つことが多い.
異なる解が異なる点における解のデータを保存するため,これらの点での差に繋がるという場合がある.この差が解の数値誤差より大きいものにならないようにするために,InterpolationOrder->Allを使う.
NDSolveをデフォルトのMachinePrecisionで使って,解を計算する:
lowsol = NDSolve[{y'[x] == 1 / (2 y[x]), y[0] == 1 / 10}, y, {x, 0, 1}, InterpolationOrder -> All];WorkingPrecision->22でNDSolveを使って,解を計算する:
highsol = NDSolve[{y'[x] == 1 / (2 y[x]), y[0] == 1 / 10}, y, {x, 0, 1}, WorkingPrecision -> 22, InterpolationOrder -> All];誤差はかなり小さいことが多いので,これらの誤差を対数尺度上で見ると便利である.RealExponent[x]は,実質的にはLog10[Abs[x]]と等しいが,ゼロにおいて特異点を持たないので,ある点においてゼロであるかもしれない差を見る場合にはよい選択である:
Plot[Evaluate[RealExponent[(y[x] /. lowsol) - (y[x] /. highsol)]], {x, 0, 1}]residual[x_] = y'[x] - 1 / (2 y[x]);NDSolveで使われる数値メソッドはどの点においても残差が小さくなるように設計されている.残差の対数をプロットすることができる:
Plot[Evaluate[RealExponent[{residual[x] /. lowsol, residual[x] /. highsol}]], {x, 0, 1}, PlotStyle -> {GrayLevel[0], RGBColor[1, 0, 0]}, AxesOrigin -> {0, 0}]これで分かるように,数値誤差はより高いWorkingPrecisionを使っている場合に著しく小さくなる.もちろんそのトレードオフは,より長い計算時間である.
ほとんどの場合,残差を小さくすると正確な数値解が得られるが,常にこれが真であるわけではない.その簡単な例としてDuffingの方程式がある:
sol = NDSolve[{x''[t] + 1 / 10 x'[t] + x[t] (x[t]^2 - 1) == 3 / 10 Cos[t], x[0] == x'[0] == 0}, x, {t, 0, 100}, InterpolationOrder -> All];hpsol = NDSolve[{x''[t] + 1 / 10 x'[t] + x[t] (x[t]^2 - 1) == 3 / 10 Cos[t], x[0] == x'[0] == 0}, x, {t, 0, 100}, WorkingPrecision -> 24, InterpolationOrder -> All];Plot[Evaluate[RealExponent[(x[t] /. sol) - (x[t] /. hpsol)]], {t, 0, 100}]