ode45 error "Undefined function or variable"

I am modeling the position of a rocket launch using ode45. I was given all of the constant values, free body diagrams, etc, and came up with an equation for acceleration, a = (T - Fd - W)/m (where T is thrust, Fd is drag force as a function of velocity, W is the weight, and m is the mass). I am tasked with using ode45 to solve for the velocity and position, however I keep getting the error "Undefined function or variable v." I defined v as the first derivative of position, drdt(1), and I am unsure of why I am getting the error. My code is below:

In my first file I have:

%define constants
Cd = .08; %coeff of drag
m0 = 74100; %initial mass
mp = 52000; %mass of the propellant
tburn = 132; %burn time
dia = 2.4; %diameter of rocket
tspan = 1:132; %burn time vector
Thrust = 1058.697; %thrust force
mdot = 393.939; %rate of fuel consumption

G = 6.67408E-11; %gravitation constant
Me = 5.972E24; %mass of earth
Re = 6.371E6; %radius of earth

Area = (pi/4)*(dia^2); %cross sectional area of rocket

[t,v] = ode45(@(t,v) boost(m0,Cd,G,Me,Re,Area,tspan,Thrust),[0 132],[Re;0])

Then I have a function file for my boost function

function drdt = boost(m0,Cd,G,Me,Re,Area,tspan,Thrust)

pos = @(t) 5.0177.*t.^2;

r = pos(tspan); %position above earths surface
R = r + Re; %distance to center of earth
m = m0 - 393.939.*tspan; %mass as a function of time
Fg = (G.*Me.*m)/(R.^2); %gravity as a function of position
W = m*Fg; %weight
[T, a, P, rho] = atmosisa(r) %command to get atmosphere data

if r > 80000
    Fd = 0;
else
    Fd = (1/2).*Cd.*Area.*rho.*v.^2; %drag force calculation below 80,000
end

drdt(1) = v;
drdt(2) = (Thrust - Fd - W)./m;

Can anyone help me troubleshoot my code?

Answer 1

[James Tursa]

Here is your derivative function signature:

function drdt = boost(m0,Cd,G,Me,Re,Area,tspan,Thrust)

As you can see, there is no "v" variable in the input argument list. Hence the error when you try to use "v" in the function body.

That being said, your entire setup looks suspect. Your current derivative function signature is this:

@(t,v) boost(m0,Cd,G,Me,Re,Area,tspan,Thrust)

And then inside your derivative function you have this:

pos = @(t) 5.0177.*t.^2;
r = pos(tspan); %position above earths surface

This looks backwards to me. It appears that you are trying to calculate the position of the rocket inside your derivative function based on time. That is not how you should be doing this. The ode45 routine will pass in the current state of the system which will include the position and velocity. It is those values that you should be using in your derivative body code.

So, let's just make up a name for the current state vector and call it "y" (since that matches what the doc uses for examples). We could use any variable name for this, but "y" will hopefully make it easier for you to compare your code to the doc. y(1) is position and y(2) is velocity.

Then let's designate the initial state vector as y0 (again to match the doc although we could use any name).

So we have:

y0 = [Re;0]; % Start on Earth surface with 0 velocity.

Then your derivative function signature needs to be something like this:

@(t,y) boost(t,y,m0,Cd,G,Me,Re,Area,tspan,Thrust)  % <-- added t and y

Then inside your boost function you would have something like this:

function drdt = boost(t,y,m0,Cd,G,Me,Re,Area,tspan,Thrust)
v = y(2); % velocity
r = y(1) - Re; % position above earths surface
R = r + Re; % distance to center of earth
m = m0 - 393.939*t; % mass as a function of time

I can't run your code because I don't have all the necessary pieces, but make these changes and give it a try. E.g., call with

y0 = [Re;0]; % Start on Earth surface with 0 velocity.
[T,Y] = ode45(@(t,y) boost(t,y,m0,Cd,G,Me,Re,Area,tspan,Thrust),[0 132],y0)

A few more comments:

At some point this will all become invalid when the mass of your rocket drops below a reasonable value. So your simulation will only be valid up until that point. E.g., maybe cut the Thrust to 0 at some point in time, and freeze the mass at that point also. If you do this and you are not at escape velocity, then the rocket will eventually fall back to earth. In that case the drag force will be in the opposite direction. To make your drag force generic to cover this case, change the v.^2 calculation to v*abs(v). That way it will be valid regardless of the direction the rocket is currently traveling.

[Jacob D Dunahue]

I made the suggested changes, and it produced a new error:

In an assignment  A(:) = B, the number of elements in A and B must be the same.

Error in boost (line 20)
drdt(2) = (Thrust - Fd - W)./m

My new code, with the changes looks like this:

%define constants
Cd = .08; %coeff of drag
m0 = 74100; %initial mass
mp = 52000; %mass of the propellant
tburn = 132; %burn time
dia = 2.4; %diameter of rocket
tspan = 1:132; %burn time vector
Thrust = 1058.697; %thrust force
mdot = 393.939; %rate of fuel consumption

G = 6.67408E-11; %gravitation constant
Me = 5.972E24; %mass of earth
Re = 6.371E6; %radius of earth

Area = (pi/4)*(dia^2); %cross sectional area of rocket

y0 = [Re;0];

[T, Y] = ode45(@(t,y) boost(t,y,m0,Cd,G,Me,Re,Area,tspan,Thrust),[0 132],y0)

and boost:

function drdt = boost(t,y,m0,Cd,G,Me,Re,Area,tspan,Thrust)

v = y(2);
r = y(1) - Re; %position above earths surface
R = r + Re; %distance to center of earth
m = m0 - 393.939.*tspan; %mass as a function of time
Fg = (G.*Me.*m)/(R.^2); %gravity as a function of position
W = m.*Fg; %weight
[T, a, P, rho] = atmosisa(r) %command to get atmosphere data

if r > 80000
    Fd = 0;
else
    Fd = (1/2).*Cd.*Area.*rho.*v.^2; %drag force calculation below 80,000
end

drdt(1) = v;
drdt(2) = (Thrust - Fd - W)./m

[James Tursa]

You did not replace the mass calculation with what I had given you. You still have this:

m = m0 - 393.939.*tspan; %mass as a function of time

When I gave you this:

m = m0 - 393.939*t; % mass as a function of time

Note that you are still using tspan, a time vector, whereas I use the single scalar t which is passed in from ode45.

Also read the extra info I have posted above about m not being reasonable at some point in time, and replacing v.^2 with v*abs(v).

Answer 2

[Walter Roberson]

You have

[t,v] = ode45(@(t,v) boost(m0,Cd,G,Me,Re,Area,tspan,Thrust),[0 132],[Re;0])

this does not pass v into the boost() function. Your declaration of function boost does not expect v as a parameter either.

When you have a function that uses a variable that is defined in a calling routine, MATLAB does not search upwards through the call tree to try to find a variable with that name. All variables must be either passed in as parameters or else lexically scoped (nested functions with shared variables.)

Also, you have

tspan = 1:132; %burn time vector

and you pass tspan into boost() . Inside boost() you have

pos = @(t) 5.0177.*t.^2;
r = pos(tspan); %position above earths surface

tspan is the vector of times, not the current time, so r is going to be a vector.

You then have

if r > 80000
    Fd = 0;
else
    Fd = (1/2).*Cd.*Area.*rho.*v.^2; %drag force calculation below 80,000
end

which is going to compare that vector r to the constant 80000. In MATLAB, when you test a vector, the result is considered true only if all elements are considered true, equivalent of

 if all(r > 80000)

which is almost certain to be false because the position vector will have a mix of above and below the value.

Perhaps you need to change your code to calculate r based upon the current time, rather than upon the tspan: if you can manage to make r a scalar instead of a vector then the test would make more sense.

If you do manage to make r a scalar rather than a vector, then you need to make a correction to the constant, as it is in imperial units whereas your mass of earth and radius of earth are in metric units.

The next problem after that is that you have a discontinuity in your calculations. Fd is increasing with the square of velocity (which increases as long as the thrust is going) but suddenly at 80000 chains you cut Fd to 0. The ode* routines cannot handle discontinuities in the functions. At every point of discontinuity you must stop the ode45() integration; you can then restart with another ode45() call that uses the output of the first call as the initial conditions for the second call.

Because your code is set up so that it does not change Fd at a fixed time, you cannot just break the calls up by tspan (the easy way): you need to program an event function to detect the boundary and signal for termination of the call. I would, however, put it to you that it would not take much work at all to transform the calculation of the termination of effective Fd from something calculated inside the routine over into something that could be calculated in advance -- at least if you follow the logic that you are currently using to calculate r. (A more accurate simulation would not calculate it based upon time, as the pos() function makes assumptions about converting time to height that should instead be the result of modeling -- for example with higher mass the climb rate would be slower -- so a more accurate simulation would use an event function anyhow.)

[James Tursa]

"... Perhaps you need to change your code to calculate r based upon the current time ..."

No. You need to get r from the current state vector.

"... Fd is increasing with the square of velocity ..."

Fd is proportional to the square of the velocity, but there is also a rho factor that comes from the atmosphere. At some point in time rho will be decreasing as a function of altitude faster than v^2 increases from that point on, so Fd will eventually decrease over time.

"... The ode* routines cannot handle discontinuities in the functions ..."

For the purposes of the simulation, that discontinuity in the derivative function might be able to be overcome by simply using tight tolerances on the relative error to manage integrating through that point gracefully. E.g.,

opts = odeset('RelTol',1e-10);
[T,Y] = ode45(@(t,y) boost(t,y,m0,Cd,G,Me,Re,Area,tspan,Thrust),[0 132],y0,opts);

E.g., compare these results where the derivative has a discontinuity at t=5:

[T,Y] = ode45(@(t,y)(t<=5)*1,[0 10],0);
figure;plot(T,Y); grid on
opts = odeset('RelTol',1e-10);
[T,Y] = ode45(@(t,y)(t<=5)*1,[0 10],0,opts);
figure;plot(T,Y); grid on

By specifying a tight relative tolerance, that corner has been managed gracefully. I suspect this rocket simulation will benefit similarly.