使用Mathematica中的库函数求解Duffing方程（技术贴）

Duffing方程是非线性动力学中的经典模型，用于描述具有非线性恢复力的受迫振动系统。其一般形式为：

Duffing方程揭示了线性系统无法描述的真实世界复杂性，是理解非线性动力学的重要工具。

核心特点：

• 非线性恢复力

恢复力包含线性项（(\alpha x)）和立方项（(\beta x^3)），导致系统呈现非线性特性：

• 硬弹簧效应（(\beta > 0)）：刚度随位移增大而增强。
• 软弹簧效应（(\beta < 0)）：刚度随位移增大而减弱。
• 典型现象

• 双稳态：系统可能有两个稳定平衡态，出现跳跃现象（振幅随频率突跳）。
• 混沌行为：特定参数下，响应呈现对初值敏感的混沌轨迹。
• 次谐波与超谐波共振：非线性导致非整数倍频的响应。

• 物理意义

方程可模拟机械结构（如夹紧梁、弹簧振子）在周期驱动下的非线性振动，尤其在‘大变形’或‘材料非线性’显著时适用。

接下来，使用Mathematica做简单的数值求解。

数值模拟：

下面代码，定义了要研究的方程：

\[Epsilon] = 3/10; \[Omega] = 1; \[Gamma] = 15/100; deq = x''[t] + \[Gamma] x'[t] - x[t] + x[t]^3 == \[Epsilon] Cos[\[Omega] t];

对于这些参数，系统轨道其实是确定性混沌的。在Mathematica中，可以使用NDSolve来数值求解这个方程。这里可以绘制平面相图

pplot = ParametricPlot[ Evaluate[{x[t], x'[t]} /. NDSolve[{deq, x'[0] == x[0] == 0}, x, {t, 0, 25}]], {t, 0, 25}]

混沌轨迹通常对初始条件极为敏感。此示例展示了如何跟踪一系列初始条件下轨迹的变化。

一种方法是：初始时在直线上设置足够密集的间隔点，使曲线在特定时间内保持可分辨。这种方法相对简单但计算量大，且难以预判需要设置多少个初始点。

另一种方法是：随着轨迹演化自适应地添加点， 自动解析特征。具体流程是：先将所有轨迹推进一个时间步长，然后根据特定准则添加新点。本示例仅以点间距为判断准则，更复杂的标准（如曲率）也可使用，但单纯基于距离的准则更简洁。

最高效的求解方式是使用NDSolve，并通过向量初始条件描述一阶系统进行求解。

将Duffing方程改写成一阶方程组：

Clear[advance]; advance[{initx_, inity_}, t0_, t1_, opts___] := Block[{t, x, y}, {x[t1], y[t1]} /. First[NDSolve[{x'[t] == y[t], x[t0] == initx, y'[t] == \[Epsilon] Cos[\[Omega] t] + x[t] (1 - x[t]^2) - \[Gamma] y[t], y[t0] == inity}, {x, y}, {t, t0, t1}, opts]]];
v = {0., 0.}; \[CapitalDelta]t = .1; points = Table[v = advance[v, t - \[CapitalDelta]t, t], {t, \[CapitalDelta]t, 25, \[CapitalDelta]t}]; Show[pplot, Graphics[{Red, Point[points]}]]

NDSolve允许初始值为一个向量，因此函数advance可以同时处理多个轨道。注意：此时x和y坐标要分离，替换成 {{x1,x2,...},{y1,y2,...}}.

v = {Range[0., 1., .1], ConstantArray[0., {11}]}; \[CapitalDelta]t = .1; points = Table[Transpose[ v = advance[v, t - \[CapitalDelta]t, t]], {t, \[CapitalDelta]t, 1.5, \[CapitalDelta]t}]; ListLinePlot[Transpose[points]]

可见点间距已显著增大。细化方案的核心思路是：在间距过大前，于相邻点之间插入新点。

此函数定义了一个细化函数 refine，当给定点集 {p1,p2,…,pn} 时，将生成包含原有点的新点集，确保相邻点间最大间距不超过容差 tol，新增点通过线性插值定位。

refine[dist_, {px_, py_}] := Module[{dx = Differences[px], dy = Differences[py], k, s, rx, ry}, s = MapThread[Function[k = Max[1, Ceiling[Sqrt[(#1^2 + #2^2)]/dist]]; N[Range[0, k - 1]]/N[k]], {dx, dy}]; rx = Apply[Join, Append[MapThread[#1 + #2*#3 &, {Drop[px, -1], dx, s}], px[[{-1}]]]]; ry = Apply[Join, Append[MapThread[#1 + #2*#3 &, {Drop[py, -1], dy, s}], py[[{-1}]]]]; {rx, ry}]

由于推进函数（advance）允许处理任意数量的点，因此推进给定初始条件集的流程为：1)使用 refine 函数生成间距足够密集的点集; 2)推进所有点的轨迹演化;3)重复上述过程.对于初始条件直线上的点，此算法通过持续细化并推进解，将相邻点间距严格控制在 0.025 以内，从而实现对非线性系统轨迹的高精度追踪。

v = {{0., 1.}, {0., 0.}}; dist = 0.025; \[CapitalDelta]t = 0.1; points = Table[v = refine[dist, v]; Transpose[ v = advance[v, t - \[CapitalDelta]t, t]], {t, \[CapitalDelta]t, 1.5, \[CapitalDelta]t}]; ListPlot[points]

当需要计算更大时间范围的解时，耗时可能急剧增加。由于相邻初始点的轨迹会呈指数级发散，新增点数将随时间呈指数增长，这会导致计算复杂度快速攀升。

优化计算效率的首要步骤是定位时间消耗瓶颈。 该计算过程仅涉及两个核心环节：1) 细化过程（Refinement）：动态插入新点以维持采样密度;2) 推进解的过程（Advancing）：数值积分求解轨迹演化

通过绘制轨迹演化过程中各步骤的耗时分布图，可精确识别计算资源的集中消耗点。例如：若细化耗时占比超过70%，需优化插值算法或降低容差阈值; 若推进过程占主导，可考虑采用更高效的ODE求解器（如并行化处理）。

可以看到：有一半的时间是用来每一步的细分上面. 因此改善此部分计算速度，将会很大地提升计算效率。

在Wolfram语言中，细化函数的加速优化面临挑战，因其处理的数组长度动态变化，既无法编译优化也无法转换为紧凑数组（packed arrays）。因此，针对该环节创建外部库函数（LibraryFunction）是合理的性能提升策略。

加载LibraryFunction中的refine函数：

refinelf = LibraryFunctionLoad["demo_numerical", "refine", {Real, {Real, 2}}, {Real, 2}];

在多数情况下，使用NDSolve内置的先进自适应步长方法进行微分方程数值求解更为高效。然而，本例中存在两个关键因素，使得编写简单的固定步长求解器更具优势：(1) 局部误差的指数增长特性无论采用何种方法，局部误差最终会呈指数级增长，导致数值解本质上仅能 跟踪真实解的轨迹。只要确保局部误差始终不超出阈值，长期来看，过度精细的误差控制对结果精度提升有限。(2) 步长与细化步骤的协同性对于适当高阶的方法（如经典四阶龙格-库塔法），ODE求解器的步长与细化步骤的间隔相匹配。通过 交替执行单步求解与细化操作，可在避免额外计算开销的前提下快速获得结果。尽管采用低阶方法（如二阶）可进一步节省时间，但需权衡精度损失风险。

advancelf = LibraryFunctionLoad["demo_numerical", "duffing_crk4", {{Real, _}, Real, Real}, {Real, _}]

经过优化，单步推进的库函数（LibraryFunction）显著提升了计算速度。由于各子步骤耗时过短（几乎可忽略），未单独绘制其耗时曲线。但当模拟时间足够长时，绘制子步骤耗时分布可发现， 微分方程求解仍是主要时间消耗环节。

在迭代过程高度优化后，现已能执行更大规模计算。通过从较短初始线开始演化，观察其逐步拉伸、折叠的序列图，可直观展示相邻初始条件在混沌系统中的指数发散与相空间形变过程。

AbsoluteTiming[v = {{0., 0.01}, {0., 0.}}; dist = 0.025; \[CapitalDelta]t = 0.1; \[CapitalDelta]T = 1.0; Tend = 60; plots = Table[Do[ v = advancelf[refinelf[dist, v], t - \[CapitalDelta]t, t], {t, T - 1 + \[CapitalDelta]t, T, \[CapitalDelta]t}]; points = Transpose[v]; Graphics[ Line[points, VertexColors -> Map[Hue, N[Range[Length[points]]]/Length[points]]], Axes -> True, PlotRange -> {{-1.5, 1.5}, {-1, 1}}], {T, \[CapitalDelta]T, Tend, \[CapitalDelta]T}];] Take[plots, {10, 60, 10}]