Creating vector in ODE Function to plot graphs

2022 7 9
1 个回答
回答已采纳
更新时间：2022 7 10
10 次查看（30 天）

0 个投票

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

请先登录，再回答此问题。

 采纳的回答

Star Strider
Star Strider 2022-7-10
Ran in:

0 个投票

Ran in:
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
.

7 个评论
显示 5更早的评论 隐藏 5更早的评论

Leon Waschkowski
Leon Waschkowski 2022-7-10
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
Star Strider 2022-7-10
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
Leon Waschkowski 2022-7-10
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
Torsten 2022-7-10
Ran in:
Ran in:
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
Torsten 2022-7-10
Ran in:
Ran in:
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
Leon Waschkowski 2022-7-10
Thank you very much, you saved my day.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Assembly 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by