清空筛选条件
清空筛选条件

Using ode45 for high altitude balloon trajectory: need some variables to update every iteration and need to plot altitude vs time.

4 次查看（过去 30 天）
Kellie Matthews
Kellie Matthews 2017-7-28
编辑： Torsten 2017-7-31
I'm trying to predict the trajectory of a high altitude balloon. This is the second order ODE I am using: (mass+cb*rho*vol)*z"=g(rho*vol-mass)-.5*rho*Cd*z'*|z'|*Ca
cb, Cd, g, and mass are constants. density[rho],volume, (approximate) surface area of a circle [Ca] all change with altitude.
I used this link to help me set up my second order ODE as a function: http://www.math.purdue.edu/~shen/cs614/projects/ode45.pdf
function xp=F(t,x) %xp = x prime or x'
xp=zeros(2,1); %output must be column vector
xp(1)=x(2);
xp(2)=(g*(RhoA*vol-mass)-.5*RhoA*realCD*x(2)*abs(x(2))*Ca)/(mass+cb*RhoA*vol);
end
This is what I need for my atmospheric properties that rely on altitude:
if (z <= 11000) %Meters (Troposhpere)
temp = 15.04 - 0.00649*z;
tempK = temp + 273.15;
p = 101.29*((temp+273.1)/(288.08)).^5.256; %kPa
elseif (z > 11000 && z < 25000) %Meters (Lower Stratosphere)
temp = -56.46;
tempK = temp + 273.15;
p = 22.65*exp(1.73-0.000157*z); %kPa
else %Upper Stratosphere
temp = -131.21 + 0.00299*z;
tempK = temp + 273.15;
p = 2.488 * ((temp+273.1)/216.6).^-11.388; %kPa
end
dTempK = abs(tempK - oldTempK);
RhoA = (p/(.2869*tempK));
Wg = Mb.*(1000*p).*vol/(r.*tempK);
radius = ((3/(4*pi))*vol).^(1/3);
Ca = pi*radius.^2;
old_z = z;
[t,x]=ode45('F',[0,tf],[0,0]);
hold on
plot(t,x(:,1))
z=x(i,1);
dz = z - old_z; %this is the change in altitude from the last second
dVol = (r/(p*Mb))*(Wg*dTempK/dt)*dt + (RhoA/p)*(vol)*dz;
vol = vol + dVol;

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

采纳的回答

Ari
Ari 2017-7-28
编辑：Ari 2017-7-28
or variables that change with time or state you should put their calculations inside the function xp. In your case, your states x seem to be [z;z']. Set z = x(1) in the beginning of the function and calculate the variable parameters before you calculate xp(2). It seems you will run into a problem trying to access the z value of the previous timestep (old_z) unless you use a persistent variable. You can try the following.
function xp = F(t,x)
persistent old_z;
z = x(1);
dz = z - old_z;
old_z = z; % set the value of old_z for next timestep
% calculate variable parameters
...
xp = zeros(2,1);
xp(1) = x(2);
xp(2) = (g*(RhoA*vol-mass)-.5*RhoA*realCD*x(2)*abs(x(2))*Ca)/(mass+cb*RhoA*vol);
end
The persisent variable will remain in memory between calls to the function.
  4 个评论
显示 2更早的评论
Kellie Matthews
Kellie Matthews 2017-7-31
I'm sorry, but this didn't work for me. After xp(2)=%ascent eq.% I have
dz = z - old_z; %this is the change in altitude from the last second
dVol = (r/(p*Mb))*(Wg*dTempK/dt)*dt + (RhoA/p)*(vol)*dz;
vol = vol + dVol;
However, dVol and vol are underlined in red and matlab says "The value assigned to variable '..' might be unused" Does this mean that the volume is not updated throughout the code? I tried the if statement you suggested and my graph came out exactly the same as before, with no apparent decrease in altitude caused by a decrease in volume.
This is what I tried with an if statement:
dz = z - old_z; %this is the change in altitude from the last second
if t <= 10800
dVol = (r/(p*Mb))*(Wg*dTempK/dt)*dt + (RhoA/p)*(vol)*dz;
vol = vol + dVol;
else
vol = vol - 10;
end
The variables dVol and vol are still underlined in red when I try this.
Torsten
Torsten 2017-7-31
编辑：Torsten 2017-7-31
I wonder why you don't solve an additional (third) ODE for "vol" together with the two ODEs for height and velocity:
dVol/dt = (r/(p*Mb))*(Wg*dTempK/dt) + (RhoA/p)*(Vol)*dz/dt ;
This way, you overcome all the problems from above.
Best wishes
Torsten.

请先登录，再进行评论。

更多回答（0 个）

类别

在 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!

Translated by