Creating vector in ODE Function to plot graphs

Hi, i want to plot the acceleration over time but i dont know how i can create a vector for the acceleration dvdt and the time t. There are just single outputs for them when i run the code.

tRange = [0 120];
v0 = 6000;
[tSol,vSol] = ode45(@fallingbody,tRange,v0);
plot(tSol,vSol)

function dvdt = fallingbody(t,v)
    D = 3.3e-5;
    g = 3.72;
    dvdt = g - D*v^2
    
%plot(t,dvdt)
end

Answer 1

[Star Strider]

Probably the easiest way to return the derivative is to use a loop to calculate it from the original function using ‘tSol’ and the solved values for the velocity, ‘vSol’ —

tRange = [0 120];
v0 = 6000;
[tSol,vSol] = ode45(@fallingbody,tRange,v0);

for k = 1:numel(tSol)
    dvdt(k,:) = fallingbody(tSol(k),vSol(k));
end

figure
plot(tSol,vSol,  tSol,dvdt)
grid
legend('Velocity','Acceleration', 'Location','best')

function dvdt = fallingbody(t,v)
    D = 3.3e-5;
    g = 3.72;
    dvdt = g - D*v^2;
end

.

[Leon Waschkowski]

Thank you for your answer. Is there a way to create a vector in the function not only for the derivative but also for example for 'D'? Im trying to transfer this solution to a rocket trajectory simulation where i want to plot the Thrust over time.

[Star Strider]

My pleasure!

This is the only way that I am aware of that it is possible to calculate the derivative using the original function. An alternative would be to use the gradient function on the result, for example:

dvdt = gradient(vSol) ./ gradient(tSol);

Since ‘D’ is a constant, just define it again in the calling workspace.

[Leon Waschkowski]

For a constant like 'D' it works. In my example the thrust is time dependent (if loop) and when i try to create a vector it just creates a vector with the last value of the thrust. In this example i get a array with all zeros because at the end at t =120 i have a thrust of 0 and the phase t<=10 is not taken into account. I hope i described my problem understandably.

global thrust

tRange = [0 120];
v0 = 6000;
[tSol,vSol] = ode45(@fallingbody,tRange,v0);

for k = 1:numel(tSol)
    T(k,:) = thrust;
end

figure
plot(tSol,T)


function dvdt = fallingbody(t,v)

global thrust

if t <= 10
    thrust = 5;
else
    thrust = 0;
end

g = 3.72;
dvdt = g - thrust*v^2;
   
end

[Torsten]

tRange = [0 120];
v0 = 6000;
[tSol,vSol] = ode45(@fallingbody,tRange,v0);
for k = 1:numel(tSol)
    T(k,:) = Thrust(tSol(k));
end
figure
plot(tSol,T)

function dvdt = fallingbody(t,v)
  thrust = Thrust(t);
  g = 3.72;
  dvdt = g - thrust*v^2;
end

function thrust = Thrust(t)
  if t <= 10
    thrust = 5;
  else
    thrust = 0;
  end
end

[Torsten]

Or:

tRange = [0 120];
v0 = 6000;
[tSol,vSol] = ode45(@fallingbody,tRange,v0);
for k = 1:numel(tSol)
    [~,T(k,:)] = fallingbody(tSol(k),vSol(k,:));
end
figure
plot(tSol,T)

function [dvdt,thrust] = fallingbody(t,v)
  if t <= 10
    thrust = 5;
  else
    thrust = 0;
  end
  g = 3.72;
  dvdt = g - thrust*v^2;
end

[Leon Waschkowski]

Thank you very much, you saved my day.

[Star Strider]

@Torsten — Thank you!