How to rewrite Maple code to Matlab, functions fsolve() and solve()?

Question

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

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

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

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

I need to rewrite the code on Matlab (version R2013a).

My attempt is below.

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

I have warnings:

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 

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

Answer 1

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

---

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.

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.

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

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