# Comparisons of numerical solution methods

11 views (last 30 days)
Enrico Incardona on 5 Dec 2022
Answered: Enrico Incardona on 8 Dec 2022
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.
Jiri Hajek on 6 Dec 2022
...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.

Sam Chak on 6 Dec 2022
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

Enrico Incardona on 8 Dec 2022
Hi everybody, thank you all for your answers; it really helped, I appreciate.
Have a nice day.
Enrico