0 个投票

Hi guys, I'm currently working on an assignment where I need to use ODE45. I have 4 equations to solve, which are linked to each other. The relevant section of the function i'm giving to ode45 looks somewhat like this:
[t,x]=ode45(@(t,x) rocketEqns(A,T,g0,m0,cd,rho0,me,htower,hscale,t),[t0 tbo],Initvals);
function[Variables]=rocketEqns(A,T,g0,m0,cd,rho0,me,htower,hscale,t)
%it is important that the rocket remains at a pitch angle of 89.85 degrees
%until past the tower height, ie alt>htower;
Variables=[0 0 0 0]'; %h,v,gamma,x
g=g0/((1+(Variables(1)/6378000))^2);
rho=rho0*exp((-1*Variables(1))/hscale);
q=0.5*rho*(Variables(2)^2);
D=q*A*cd;
if Variables(1)<htower
Variables(3)=0;
Variables(2)=(T-D-((m0-(me*t)))*g*sin(Variables(3))/(m0-(me*t)));
Variables(1)=Variables(2);
Variables(4)=0;
else
Variables(3)=(-g*cosd(Variables(3))/Variables(2));
Variables(2)=(T-D-((m0-(me*t)))*g*sin(Variables(3))/(m0-(me*t)));
Variables(1)=Variables(2)*sind(Variables(3));
Variables(4)=Variables(2)*cosd(Variables(3));
end
end
ie if the 1st variable is less than 90 then the second variable doesn't change, if the first variable is greater than 90 then the second variable changes by whatever is calculated inside the function. What I think is actually happening is that this is evaluating whether the rate of change of variable(1) is less than 90.
How do i access the current value? I tried getting it from the output but it gave me an error.
Thanks in advance

4 个评论
显示 2更早的评论 隐藏 2更早的评论

James Tursa
James Tursa 2018-3-23
You need to show more of your code. That small snippet you currently have posted doesn't help us much.
Aarron Sheppard
Aarron Sheppard 2018-3-23
I've edited my post to include my function
Steven Lord
Steven Lord 2018-3-23
In addition to the other comments or answers, this section of your code is suspicious.
Variables(2)=(T-D-((m0-(me*t)))*g*sin(Variables(3))/(m0-(me*t)));
Variables(1)=Variables(2)*sind(Variables(3));
Is Variables(3) an angle in degrees (in which case both lines should use sind) or an angle in radians (in which case both should use sin)?
Aarron Sheppard
Aarron Sheppard 2018-3-24
Good spot, thank you!

请先登录，再进行评论。

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

 采纳的回答

James Tursa
James Tursa 2018-3-23
编辑：James Tursa 2018-3-23

0 个投票

Some issues with your code:
1) You don't pass the state vector into your derivative routine! E.g. in this line:
[t,x]=ode45(@(t,x) rocketEqns(A,T,g0,m0,cd,rho0,me,htower,hscale,t),[t0 tbo],Initvals);
The state vector being passed from ode45 into your derivative routine is x (from the @(t,x) part), but x does not appear in the input argument list in your rocketEqns function! You need to pass that in. E.g.,
[t,x]=ode45(@(t,x) rocketEqns(A,T,g0,m0,cd,rho0,me,htower,hscale,t,x),[t0 tbo],Initvals);
2) You appear to be using Variables as both the state vector input and the derivative output in your function. E.g.,
function[Variables]=rocketEqns(A,T,g0,m0,cd,rho0,me,htower,hscale,t)
%it is important that the rocket remains at a pitch angle of 89.85 degrees
%until past the tower height, ie alt>htower;
Variables=[0 0 0 0]'; %h,v,gamma,x
g=g0/((1+(Variables(1)/6378000))^2);
rho=rho0*exp((-1*Variables(1))/hscale);
q=0.5*rho*(Variables(2)^2);
So, you list Variables as the output (meaning it will be the derivative of your state vector), and you set Variables=[0 0 0 0]' to begin with. All well and good. But then you turn around and start using it on the right-hand-side as if it were the incoming state vector, which it is not! To fix this, you need to use the incoming state vector (which you will fix from part 1 of my post) on the right hand side. E.g., something like this:
function[Variables]=rocketEqns(A,T,g0,m0,cd,rho0,me,htower,hscale,t,x) % <-- added the x
%it is important that the rocket remains at a pitch angle of 89.85 degrees
%until past the tower height, ie alt>htower;
Variables=[0 0 0 0]'; %h,v,gamma,x
g=g0/((1+(x(1)/6378000))^2); % <-- using incoming state vector x on rhs
rho=rho0*exp((-1*x(1))/hscale); % <-- using incoming state vector x on rhs
q=0.5*rho*(x(2)^2); % <-- using incoming state vector x on rhs
:
etc
You don't have to use the variable name x for the state vector, of course, if that name means something else to you in your problem context. You can use any variable name you like for the state vector.

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Aarron Sheppard
Aarron Sheppard 2018-3-23
Fantastic, thank you! as a follow up question, my equation
rho=rho0*exp((-1*x(1))/hscale);
should be a value calculated, not a derivative (the same applies for a few others) do i need to calculate these outside of the function and pass it in, or will it be ok where it is?
James Tursa
James Tursa 2018-3-23
I don't understand the question. That line calculates a value of rho based on the incoming state vector, which appears to me to be an atmospheric density calculation, which then gets used with the velocity to produce a drag force. It is fine to do this calculation inside your derivative routine. In fact, that is where it must be since it is a function of the current state vector. Is that what you are asking?
Btw, to make your drag force calculation more generic so that it covers the case of a rocket body falling down as well as going up, change this line:
q=0.5*rho*(Variables(2)^2);
to be this instead:
q=0.5*rho*abs(Variables(2))*Variables(2);
This will allow the drag force to always be in the opposite direction of the current velocity, regardless of whether it is going up or down.
Aarron Sheppard
Aarron Sheppard 2018-3-23
yes, that's what i was asking, thank you; im still very new to ode functions. Wouldn't i be using x(2) instead of Variables(2)? Thank you again
James Tursa
James Tursa 2018-3-23
Yes. I haven't looked over your code in much detail, but in general all of the Variables(etc) stuff on the right hand side should probably be x(etc).

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 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