Access current value ode45

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

[James Tursa]

You need to show more of your code. That small snippet you currently have posted doesn't help us much.

[Aarron Sheppard]

I've edited my post to include my function

[Steven Lord]

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]

Good spot, thank you!

Answer 1

[James Tursa]

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.

[Aarron Sheppard]

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]

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]

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]

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).