Comparisons of numerical solution methods
6 次查看(过去 30 天)
显示 更早的评论
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 个评论
Davide Masiello
2022-12-6
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
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
0 个评论
更多回答(1 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!