Comparisons of numerical solution methods

11 ビュー (過去 30 日間)
Enrico Incardona
Enrico Incardona 2022 年 12 月 5 日
回答済み: Enrico Incardona 2022 年 12 月 8 日
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.
  4 件のコメント
Jiri Hajek
Jiri Hajek 2022 年 12 月 6 日
...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
Sam Chak 2022 年 12 月 6 日
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

その他の回答 (1 件)

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

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by