How to| 验证 NDSolve 的结果
对于大多数微分方程,由 NDSolve 给出的解是相当准确的. 然而,由于它的结果基于数值抽样和误差估计,偶尔可能会出现明显的误差. 当您需要确保解的质量时,对解做些基本检查是是一个好主意.
对 WorkingPrecision 高于默认 MachinePrecision 时计算得到的解进行比较往往是一种检查结果的有效方法.
解不同,存储解的数据的点的位置可能也会不同,导致在这些点处出现差异. 使用 InterpolationOrder->All 可以使这种差异保持在解的数值误差范围之内.
使用 NDSolve 和默认的 MachinePrecision 计算解:
lowsol = NDSolve[{y'[x] == 1 / (2 y[x]), y[0] == 1 / 10}, y, {x, 0, 1}, InterpolationOrder -> All];使用 NDSolve 和 WorkingPrecision->22 计算解:
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}]