Comparisons of numerical solution methods
11 ビュー (過去 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!