CODEWITHSUNDEEP
matlabode
CroCo
CroCo

Reputation: 5741

How to animate solution of ode with ode45 solver?

I'm trying to access the current solution the ode45 provides. All the examples I've seen, they solve the ode and then store the solution for further manipulation. While this is an alternative solution, but I need to access the solution at each iteration so I can plot them for animation purposes. This is the code for simple pendulum.

clear 
clc

theta0 = [pi/2; 0]; 
tspan = [0 10];
[t, theta] = ode45(@odeFun, tspan, theta0);
plot(t, theta(:,1));

function dydx = odeFun(t, theta)
    g = 9.8;
    l = 1;
    dydx = zeros(2, 1);
    dydx(1) = theta(2);
    dydx(2) = -g/l*theta(1);
end

My goal is achieve the following

clear 
clc

theta0 = [pi/2; 0]; 
t=0;
while t < 10
    [t, theta] = ode45(@odeFun, ..., ...);
    plot(t,theta)
    drawnow
    t = t + 0.1;
end

[t, theta] = ode45(@odeFun, tspan, theta0);
plot(t, theta(:,1));

function dydx = odeFun(t, theta)
    g = 9.8;
    l = 1;
    dydx = zeros(2, 1);
    dydx(1) = theta(2);
    dydx(2) = -g/l*theta(1);
end

Upvotes: 0

Views: 326

Answers (1)

AVK
AVK

Reputation: 2149

You can use the ode output function option. The ODE solver calls the output function after each successful time step. Here is a very simple example of using an output function:

clear 
clc
close all
theta0 = [pi/2; 0]; 
tspan = [0 10];
options = odeset('OutputFcn',@plotfun,'MaxStep',0.0001);
axes; hold on
[t, theta] = ode45(@odeFun, tspan, theta0, options);

function status = plotfun(t,y,flag)
    plot(t, y(1,:),'k-',t,y(2,:),'m-');
    drawnow
    status = 0;
end

function dydx = odeFun(t, theta)
    g = 9.8;
    l = 1;
    dydx = zeros(2, 1);
    dydx(1) = theta(2);
    dydx(2) = -g/l*theta(1);
end

The maximum integration step is set to 0.0001 so that the solution graphs are not drawn immediately.

Alternatively, you can calculate the solution by dividing the time interval into parts and using the end point of the previous part as the initial point for the next part:

clear 
clc
close all
axes; hold on
theta0 = [pi/2; 0]; 
tstart=0; dt= 0.1;
options= odeset('MaxStep',0.0001);
while tstart < 10
    [t, theta] = ode45(@odeFun,[tstart tstart+dt],theta0,options);
    plot(t,theta(:,1),'k-',t,theta(:,2),'m-');
    drawnow
    tstart = tstart + dt;
    theta0= theta(end,:);
end

function dydx = odeFun(t, theta)
    g = 9.8;
    l = 1;
    dydx = zeros(2, 1);
    dydx(1) = theta(2);
    dydx(2) = -g/l*theta(1);
end

And of course you can combine both approaches.

Upvotes: 2

Related Questions