Comparisons of numerical solution methods

Hello,

I solved the Van Der Pol equation (for m=1) with different numerical methods (ode 45, forward Euler, RK2 and RK4) and I get this graph. However, I can't understand why the curves are getting more and more offset from each other in the course of time. It is noticeable that, taking the ode45 curve as a reference, the forward euler curve shifts positively while the RK2 and RK4 curves shift negatively. Also, when measuring the period of each curve it seems to be constant! I am really confused...

Thank you in advance for your clarifications and have a nice day.

[Torsten]

There are good and bad codes ...

[Sam Chak]

The foundation of the these solvers is based on the Taylor series expansion.

Higher-order terms or derivatives lead to better accuracy in the solution but also become increasingly more complicated at the same time.

The accuracy of the solution in lower-order solvers like Euler's method can be compensated by having smaller step size. However, smaller step size demands large computational effort to obtain acceptable error levels.

[Davide Masiello]

Hi @Enrico Incardona, you should post your full code in order to enable us to help.

[Jiri Hajek]

...why the curves are getting more and more offset from each other in the course of time...

The short answer is: due to errors accumulating during the solution process in each of the solutions. Each ODE solver does make an error due to approximations used by its numerical method. But each solver makes a different error, that's why the solutions gradually differ from each other. An exception to this would be a process that converges towards a limit value. There, all solvers should converge to the same value, with largers differences somewhere between the initial time and infinity.

Answer 1

[Sam Chak]

Hi @Enrico Incardona,

Here is the comparison between the solution by ode45 and the solution by Forward_Euler.

As mentioned previously, choosing a smaller step size improves the accuracy of the solution for the Euler's method.

h      = 0.01;      % I think you used step size 0.2
tStart = 0;
tFinal = 10;
tspan  = tStart:h:tFinal;
y0     = [1; 0];
mu     = 1;
odefcn = @(t, y) vdp(t, y, mu);

% Using ode45 solver
[t, yODE45] = ode45(odefcn, tspan, y0);

% Using Forward Euler solver
yEuler = ForEuler(odefcn, tspan, y0);

% Plot to compare the solutions
plot(t, [yODE45(:,1)'; yEuler(1,:)]', 'linewidth', 1.5), grid on
title('Solutions of Van der Pol Equation (\mu = 1)');
xlabel('Time t');
ylabel('Solution y');
legend('y_{ode45}','y_{Euler}')

% Forward Euler
function u = ForEuler(f, x, u0)
    u(:,1) = u0;
    h      = x(2) - x(1); n = length(x);
    for i  = 1:n-1,
        u(:,i+1) = u(:,i) + h*f(x(i), u(:,i));
    end
end

% Van der Pol oscillator
function dydt = vdp(t, y, mu)
    dydt = [y(2); mu*(1 - y(1)^2)*y(2) - y(1)];
end 

Answer 2

[Enrico Incardona]

Hi everybody, thank you all for your answers; it really helped, I appreciate.

Have a nice day.

Enrico