如何将 Maple 代码重写为 Matlab，函数 fsolve() 和 solve()？

Question

我有一个常微分方程（ODE）——范德波尔振荡器问题：

y''-2a (1-y^2)y'+y=0,

y(0)=0, y'(0)=0.5, x in [0,10], a =0.025.

ODE 是在 Maple Software 上使用 fsolve() 函数求解的。

我需要在 Matlab（版本 R2013a）上重写代码。

我的尝试如下。

n = 0
h = 0.1
syms z0; % z in point 0
syms y0; % z in point 0
f0 = 2 * 0.025 * z0 - 2 * 0.025 * ((y0)^2) * z0 - y0 
syms z1; % z in point 1
syms y1; % z in point 1
f1 = 2 * 0.025 * z1 - 2 * 0.025 * ((y1)^2) * z1 - y1
syms z2; % z in point 2
syms y2; % z in point 2
f2 = 2 * 0.025 * z2 - 2 * 0.025 * ((y2)^2) * z2 - y2
syms z32; % z in point 3/2
syms y32; % z in point 3/2
f32 = 2 * 0.025 * z32 - 2 * 0.025 * ((y32)^2) * z32 - y32
syms z3; % z in point 3
syms y3; % z in point 3
syms p;
f3 = 2 * p * (1-(y3)^2) * z3 - y3
syms z4; % z in point 4
syms y4; % z in point 4
f4 = 2 * p * (1-(y4)^2) * z4 - y4
syms z72; % z in point 7/2
syms y72; % z in point 7/2
f72 = 2 * p * (1-(y3)^2) * z72 - y72

[c1,c2,c3,c4]=solve(y0, z0, y2 - (-y0 + 2 * y1 + (1/12) * h^2 * f0 + 5/6 * h^2 * f1 + (1/12) * h^2 * f2), y32 - (-1/2 * y0 + 3/2 * y1 + 1/24 *h^2 *f0 + 13/32 * h^2 *f1 - 5/48 * h^2 * f32 + 1/32 * h^2 * f2), -360 * h * z0 - (89 * h^2 * f0 + 186 * h^2 * f1 + 33 * h^2 * f2 - 128 * h^2 * f32 + 360 * y0 -360 * y1))

我有警告：

Warning: 4 equations in 1 variables. 
> In C:\Program Files\MATLAB\R2013a\toolbox\symbolic\symbolic\symengine.p>symengine at 56
  In mupadengine.mupadengine>mupadengine.evalin at 97
  In mupadengine.mupadengine>mupadengine.feval at 150
  In solve at 170 
Warning: Explicit solution could not be found. 
> In solve at 179 

问题。 是否可以在 Matlab R2013a 中解决问题？谁能告诉我如何重写代码？

Answer 1

在matlab中求数值解的方法是定义

a=0.025;
VDP_ODE = @(t,u) [ u(2), -u(1)+2*a*u(2)*(1-u(1)^2) ]; % u = [ y, z ]

（或使用完整的function 声明定义它）并使用数值求解器调用它

h=0.1;
t = 0:h:10;
y0 = 0; z0 = 0.5;
u0 = [y0, z0 ];
[t,u] = ode45(VDP_ODE, t, u0);
figure(1); plot(t,u(:,1));
figure(2); plot(u(:,1),u(:,2));

---

如果您想构建自己的求解器，您应该首先了解方法。似乎它是一种类似于但不完全等于线性多步方法的 4 阶方法，例如用于计算分子动力学的 Beeman-Schofield 方法。

每个步骤的输入是值y(n), y(n+1)，输出减少到y(n+2)，没有其他值被带到下一步。在方法步骤中，在所有采样时间 t(n), t(n+1), t(n+3/2), t(n+2) 计算 y(n+3/2) 和导数 z 的附加值。目的是得到浮点结果，因此符号定义方程没有意义，从而使系统适应fsolve的接口。

在 Maple 代码中求解的系统可以通过观察有 6 个方程来简化，因此有 6 个未知参数。该方法基于以 6 个系数作为未知数的 5 次多项式，它对解的值、一阶和二阶导数进行插值。这给出了给定值y0,y1 的 2 个方程和 4 个方程，用于精确满足插值点处的微分方程。

Y = [ y0, step0(y0,z0) ]
for n=1:N-2
    Y(n+2,:)=stepN(Y(n,:), Y(n+1,:), h)
end

对于方法的步骤做一些类似的事情

  %%%% ======== General step ========
  %% fit a degree 5 polynomial to values, 1st and 2nd derivatives
  %% at t0, t1, t1h=t1+0.5*h, t2
  %% p(s) = sum c_k s^k/k!,  c0=y1,  y(t1+s) ~ p(s)
  %% eqn y0 = p(-h)
  %% def z0 = p'(-h)
  %% eqn f(y0,z0) = p''(-h)
  %% def z1 = p'(0)
  %% eqn f(y1,z1) = p''(0)
  %% eqn y2 = p(h), etc
  %%
  %% state vector for non-lin solver is u = [y2, c1, c2, ...c5 ]


  function y2 = stepN(y0,y1,h)
    zz = (y1-y0)/h;
    predict = [ y1+h*zz, zz, f(y1,zz), 0, 0, 0];
    options = optimset("TolX", 1e-1*h^6, "TolFun", 1e-2*h^6);
    u = fsolve(@(u) systemN(u,[y0,y1], h), predict, options);
    y2 = u(1);
  end

设置一些预测值（仅 O(h) 精确），设置求解器的容差，以便求解器误差希望小于离散化误差 O(h^6)，调用求解器端从解决方案数组。

对于求解函数的系统，首先将常数和变量提取到命名变量中以获得更好的可读性，计算给定点的函数值，然后返回离散方程中的缺陷数组。

  function eqn = systemN(u,init,h)
    [ y0, y1 ] = num2cell(init){:};
    [ y2, c1, c2, c3, c4, c5 ] = num2cell(u){:};
    p = @(h) y1 + h*(c1+h/2*(c2+h/3*(c3+h/4*(c4+h/5*c5))));
    dp = @(h) c1+h*(c2+h/2*(c3+h/3*(c4+h/4*c5)));
    d2p = @(h) c2+h*(c3+h/2*(c4+h/3*c5));
    z0 = dp(-h); z1 = c1;
    y1h = p(0.5*h); z1h = dp(0.5*h);
    z2 = dp(h);

    eqn = [ y2-p(h), 
            y0-p(-h), 
            f(y0,z0)-d2p(-h), 
            f(y1,z1)-c2, 
            f(y1h,z1h)-d2p(0.5*h),
            f(y2,z2)-d2p(h) ];
  end