WOLFRAM

Wolfram 语言与系统参考资料中心

杜芬方程的敏感度

探索杜芬方程对初始条件的敏感度.

为了简化接下来的计算，使用只有初始条件作为参数而其他值固定的 ParametricNDSolveValue：
Wolfram Language code: xp = ParametricNDSolveValue[{x''[t] + δ x'[t] + α x[t] + βx[t] ^ 3 == γ Cos[ω t], x[0] == a, x'[0] == b} /. {δ -> 0.15, α -> -1, β -> 1, γ -> 0.3, ω -> 1}, x, {t, 0, 100}, {a, b}]
关于 a 和 b 的 ParametricFunction 导数分别是关于 x 的初始值和它导数的敏感度：
Wolfram Language code: {xpa = Head[D[xp[a, b], a]], xpb = Head[D[xp[a, b], b]]}
显示开始于 {0,0} 的解的图线和敏感度的对数图线：
Wolfram Language code: GraphicsRow[{Plot[xp[0, 0][t], {t, 0, 100}], LogPlot[Abs[xpa[0, 0][t]], {t, 0, 100}], LogPlot[Abs[xpb[0, 0][t]], {t, 0, 100}]}, ImageSize -> Large]
巨大的敏感性表示这些参数有很小的变化也会导致解的大的偏差. 下面显示的是关于从原点的任一方向的小扰动的解：
Wolfram Language code: ϵ = 10 ^ -6; Plot[{xp[0, 0][t], xp[ϵ, 0][t], xp[0, ϵ][t]}, {t, 0, 100}, ImageSize -> Medium]
另一种可视化敏感度的方法是局部把它分解成在相平面上平行和垂直于轨迹的方向的分量.
定义一个函数，给出垂直于轨迹的向量作为时间的函数：
Wolfram Language code: perp[a_, b_][t_ ? NumberQ] := Module[{f = xp[a, b], dt, do}, dt = Normalize[{f'[t], f''[t]}]; do = Reverse[dt] * {1, -1} ];
现在定义函数，给出在平行和垂直方向关于 x 的初始值的敏感度的分量幅度：
Wolfram Language code: spar[a_, b_][t_ ? NumberQ] := Module[{s = xpa[a, b], f = xp[a, b]}, Abs[Normalize[{f'[t], f''[t]}].{s[t], s'[t]}]]
Wolfram Language code: sperp[a_, b_][t_ ? NumberQ] := Module[{s = xpa[a, b]}, Abs[perp[a, b][t].{s[t], s'[t]}]]
对于该方程，两个分量以可比较的速率增长.
分解分量可以在相平面上有很好的可视化. 改变敏感度的尺度可以查看不同的时间间隔：
Wolfram Language code: scale = 10 ^ -8;Show[ParametricPlot[{xp[0, 0][t], xp[0, 0]'[t]} + s sperp[0, 0][t] perp[0, 0][t], {t, 0, 100}, {s, -scale, scale}, Mesh -> None, PlotPoints -> {1000, 2}], ParametricPlot[{xp[0, 0][t], xp[0, 0]'[t]}, {t, 0, 25}, PlotStyle -> Red], ImageSize -> Medium]