Solving problem with ODE45

1 view (last 30 days)
Yi Li
Yi Li on 17 Oct 2021
Commented: Yi Li on 8 Nov 2021
clc;
clear all;
close all;
[t1,y1] = ode45(@vdp1,[0 20],[2; 0]);
t2(1,:)=t1(1,:);
y2(1,:)=y1(1,:);
figure
plot(t1,y1(:,1),'-o',t1,y1(:,2),'-o')
for i=2:size(t1)
[a,b]=ode45(@vdp1,[0 t1(i) 2*t1(i)],y2(i-1,:));
t2(i,:)=a(2,:);
y2(i,:)=b(2,:);
end
figure
plot(t2,y2(:,1),'-o',t2,y2(:,2),'-o')
Hey guys, when I am solving a function step by step with ODE45, the result is completely different from the direct solution. Can someone tell me why?
  2 Comments
Jan
Jan on 17 Oct 2021
Edited: Jan on 17 Oct 2021
Note: for i=2:size(t1) is fragile. size() replies a vector, but the colon operator uses the first element only. Better: for i=2:size(t1, 1) or for i=2:numel(t1).
Please provide vdp1 to allow use to run your code.
Yi Li
Yi Li on 17 Oct 2021
Edited: Yi Li on 17 Oct 2021
Thank you!
function dydt = vdp1(t,y)
%VDP1 Evaluate the van der Pol ODEs for mu = 1
%
% See also ODE113, ODE23, ODE45.
% Jacek Kierzenka and Lawrence F. Shampine
% Copyright 1984-2014 The MathWorks, Inc.
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

Sign in to comment.

Accepted Answer

Jan
Jan on 17 Oct 2021
In the loop you are calculating the integration from 0 to 2*t(i) with the initial value set to y2(i-1, :). But this is another integration, because you start from t=0 with another initial value. Instead of the time interval [0 t1(i) 2*t1(i)] the interval [t1(i-1), t1(i)] would create the same output:
for i = 2:numel(t1)
[a,b] = ode45(@vdp1, [t1(i-1), t1(i)], y2(i-1, :));
t2(i,:) = a(end, :);
y2(i,:) = b(end, :);
end
figure
plot(t2, y2(:, 1), '-o', t2, y2(:, 2), '-o');

More Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!