如何将Maple代码重写为Matlab，fsolve()和resolve()函数? [英] How to rewrite Maple code to Matlab, functions fsolve() and solve()?

问题描述

我有一个常微分方程(ODE)-Van Der Pol振荡器问题:

I have an ordinary differential equation (ODE) -- Van Der Pol oscillator Problem:

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

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

y(0)= 0，y'(0)= 0.5，[0,10]中的x​​，a = 0.025.

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

ODE在Maple Software上使用 fsolve()函数进行了求解.

The ODE was solved on Maple Software with the fsolve() function.

我需要在Matlab(R2013a版本)上重写代码.

I need to rewrite the code on Matlab (version 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中的问题?谁能给我一个思路，重写代码吗?

Question. Is it possible to solve the problem in Matlab R2013a? Can anyone please give an idea how to rewrite the code?

推荐答案

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

The way to find a numerical solution in matlab is to define

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

(或使用完整的函数声明进行定义)，然后使用数值求解器调用它

(or define it using the full function declaration) and call this with a numerical solver

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));

---

如果要构建自己的求解器，则应首先了解该方法.看来这是一种四阶方法，与线性多步方法相似，但并不完全等同于线性多步方法，就像用来计算分子动力学的Beeman-Schofield方法一样.

---

If you want to build your own solver, you should first understand the method. It seems that it is a 4th order method similar but not exactly equal to a linear multi-step method, like the Beeman-Schofield methods introduced to compute molecular dynamics.

每个步骤的输入是值 y(n)，y(n + 1)，输出减少为 y(n + 2)值将转移到下一步.在方法步骤内，在所有采样时间 t(n)，t(n +1)，t(n + 3/2)，t(n + 2).目的是获得浮点结果，因此没有必要用符号定义方程式，从而使系统适应 fsolve 的界面.

The input to each step are the values y(n), y(n+1), the output reduces to y(n+2) with no other values taken over to the next step. Inside the method step additional values are computed for y(n+3/2) and derivatives z at all sample times t(n), t(n+1), t(n+3/2), t(n+2). The aim is to get floating point results, thus it makes no sense to define the equations symbolically, thus adapt the system to the interface of fsolve.

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

The system that is solved in the Maple code can be simplified by observing that there are 6 equations, thus 6 unknown parameters. The method is based on a 5th degree polynomial with 6 coefficients as unknowns which interpolates the values, first and second derivatives of the solution. This gives 2 equations for the given values y0,y1 and 4 equations for the exact satisfaction of the differential equation at the interpolation points.

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

在该方法的步骤中执行类似的操作

and for the step of the method do something like

  %%%% ======== 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))，调用求解器端从解决方案数组中提取所需的值.

which set some predictor values (only O(h) exact), sets the tolerances for the solver so that the solver error is hopefully smaller than the discretization error which is O(h^6), calls the solver end extracts the wanted value from the solution array.

为使系统求解函数，首先将常量和变量提取为命名变量以提高可读性，计算给定点的函数值，然后返回离散方程式中的缺陷数组.

For the system to solve the function first extracts the constants and variables into named variables for better readability, computes the function values at the given point and then returns the array of the defects in the discretized equations.

  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 

这篇关于如何将Maple代码重写为Matlab，fsolve()和resolve()函数?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！