Trying to run ODE45 with system of two equations but its not working, any advice ? Thanks

g_i=9.81;
Re=6378 ;
h=100;
T=100;

Ae = 0.9; % sq m exit nozzle area
Ar = 1.7;

mSecondStage = 111500; % kg
m0_propellant = 433100; %kg
mEmpty = 22200; % kg
m0_full = m0_propellant + mSecondStage + mEmpty; % kg
Cd_up = 0.3;
Cd_down = 0.82;
N = 9; % Number of engines
A = 3.7.^2 / 4; % Diameter of 3.7m
Thrust_vacuum = 8227000; % N
Thrust_ground = 7607000; % N
Isp_vacuum = 311; % s
Isp_ground = 282; % s
vjet = 3000;
TargetAltitude = 2000000; % m


y = zeros(3,1);
tspan=[0 5];

fname=@(t,y) [((T - D)/m) - g_i.* sin(y(2)); -(g_i./(y(1)) - y(1)./Re +y(3)) .* cos(y(2) ./ y(1)); y(1).*sin(y(2)) ]

YO = [1000; 1];
[t,Y]=ode45(fname,tspan,YO); % the problem is here 
plot(t,Y(:,3),'r')

function [pe, Ae, mdot] = SolveThrust(vjet, T_sea, T_vac, Isp_sea, Isp_vac, N,g_h)
% Vacuum first
mdot_vac = T_vac / (Isp_vac * g_h);
pe_Ae = T_vac - mdot_vac * vjet;


% Sea level
mdot_sea = T_sea / (Isp_sea * 9.81);
Ae = (T_sea - mdot_sea * vjet - pe_Ae) / -101325;
pe = pe_Ae / Ae;

% Per Engine
Ae = Ae / N;
mdot = (mdot_vac + mdot_sea) / (2 * N);
end

global atmosphere % it takes atmospheric conditions, I am trying to simulate rocket tflight trajectory
 

function [D,g_h, q] = getdragandaccngrav(Re, Cd, Ar, velabs )


D = 0.5 * rho * Cd * Ar*  velabs.^2;
q = 0.5 * rho * velabs.^2;
go = 9.81;
g_h = go * ((Re)^2 / (Re + vz)^2 ;     
end 

figure(1);
plot(t, vz);
xlabel('Time');
ylabel("Altitude");
title("Altitude vs Time");

figure(2);
plot(t, vx);
xlabel('Time');
ylabel("Horizontal Position");
title("X vs Time");

figure(3);
plot(t, Fd);
xlabel('Time');
ylabel("Drag");
title("Fd vs Time");

[Jan]

The line "global atmosphere" seems to be orphaned.

Please explain "it is not working" with any details. We cannot guess, what's going wrong, but you do see this impportant detail on your screen already.

[SL]

Global Atmopshere is another mat file with atmopshere properties. But the real issue is that ODE is not executing the equations correctly, I am not able to get an outuput and I dont understand where I am going wrong.

[Torsten]

In your function definition, D and m are undefined.

Further you divide by y(1) which is 0 at time t=0.

[SL]

D and m are obtained from getdragandaccngrav function . y1 is acceleration therefore, initially zero or 1 m/s^2; and m can be any value for now, as I need to get this code running

[Torsten]

Matlab does not know how to get D and m if you just define the function in the way you did.

You will have to pass a function handle to ODE45 and calculate the constants in a respective function called "fname".

You divide by y(1) which is zero at the beginning, and MATLAB will punish that - whether the setting is senseful or not.

[James Tursa]

Please post the differential equations you are solving. E.g., an image of the equations would be acceptable. Then we can compare this with your code to determine correctness. In particular, if you are simulating the motion of a rocket in a 2D plane, then I would expect a 2nd order vector equation (F=ma where F and a are 2D vectors), which would yield the equivalent of four (4) 1st order scalar equations. Yet you only have three (3) 1st order scalar equations. Or if you are simply simulating a rocket going straight up on a line then the 1D problem would be a single 2nd order scalar equation which would yield two 1st order scalar equations. Again not three 1st order scalar equations. Something seems wrong here.

Answer 1

[Bjorn Gustavsson]

Well, your problem is that your function fname will return a 3x1 array - that correspods to three coupled first-order ODEs. Since you seem to want to integrate the rocket-equation you should have some other number of coupled ODEs. Also when you define a dynamic function, the variables that are not defined as input-variables will be taken as whatever the variables with those names in the workspace at creation-time. Since I cannot see what D is or how it is defined that might be a second problem.

g_i=9.81;
Re=6378 ;
h=100;
T=100;

Ae = 0.9; % sq m exit nozzle area
Ar = 1.7;

mSecondStage = 111500; % kg
m0_propellant = 433100; %kg
mEmpty = 22200; % kg
m0_full = m0_propellant + mSecondStage + mEmpty; % kg
Cd_up = 0.3;
Cd_down = 0.82;
N = 9; % Number of engines
A = 3.7.^2 / 4; % Diameter of 3.7m
Thrust_vacuum = 8227000; % N
Thrust_ground = 7607000; % N
Isp_vacuum = 311; % s
Isp_ground = 282; % s
vjet = 3000;
TargetAltitude = 2000000; % m


y = zeros(3,1);
tspan=[0 5];

fname=@(t,y) [((T - D)/m) - g_i.* sin(y(2));...
              -(g_i./(y(1)) - y(1)./Re +y(3)) .* cos(y(2) ./ y(1));
              y(1).*sin(y(2)) ]; % Here you define your ODE-function returning a 3-by-1 array

YO = [1000; 1];
[t,Y]=ode45(fname,tspan,YO); % the problem is here - no it is above

HTH

[SL]

D is defined below in the function;

[Bjorn Gustavsson]

Yes. But that is not how the dynamic function works. It will take the contents of the variable with that name in the current workspace, if it is not defined it doesn't help that there are some other function that returns a variable that is named something similar as an output-argument. You first have to call that function to create that variable, or if it is to be varying along the trajectory you'll have to include that as a call inside the dynamic function:

fname=@(t,y) [((T - getdragandaccngrav(Re, Cd, Ar, sqrt(sum(y(:).^2)))/m) - g_i.* sin(y(2));...
              -(g_i./(y(1)) - y(1)./Re +y(3)) .* cos(y(2) ./ y(1));
              y(1).*sin(y(2)) ]; % Here you define your ODE-function returning a 3-by-1

But you will still have the problem that you have a 3x1 ODE-definition and a 2x1 initial condition.

HTH

[Torsten]

... and m is unknown and has to be calculated from SolveThrust, I guess.

[Bjorn Gustavsson]

@Torsten: Ah missed that one.

@SL: It would be best if you defined the rocket-equation/equation of motion as a separate sub-function in this file - that way you will have a far easier way to call the other functions inside the ODE-function than with a dynamic function-definition. Far easier.

[SL]

@Bjorn Gustavsson I have three equations, one for acceleration, gamma Flight path angle and then height.

I will do that and ask for your advice.

Thanks very much,

[James Tursa]

Please post all of the equations and initial conditions. An image would be acceptable.

[SL]

[Bjorn Gustavsson]

So that is a system of four coupled ODEs. You should write a function returning the derivatives as a function of time. Something like this:

function dxyvgdt = rocket_eqs(t,xhvg,m)
% Description here
% Calling:
%   dxyvgdt = rocket_eqs(t,xhvg,m)
% Input variables:
%  t    - 
%  xhvg - 
%  m    - 
% Output variables:
%  dxyvgdt - 
% Example:
%  Just one single call to illustrate the return of the derivatives and how
%  to interpret them.
%  Then one use-case with ode45

% 0: Definition of physical constants:
R_E = 6370e3 + ; % Radius of Earth, deliberately broken
g = 9.82;        % (m/s^2) Gravitational acceleration at ground 
% Note the inclusion of units you use in the documentation!

% Extract the 4 different components
% for ease of reading - as per "programmes should be written for humans to
% read and computers to occassionally run" (paraphrased from one of the
% giants, Knuth(?) or someone at that level) 
x = xhvg(1);
h = xhvg(2);
v = xhvg(3);
gamma = xhvg(4);

% Caclulate the remaining terms of th ODEs
[T,D] = my_T_n_D_function(t,xhvg,m);

% Return the four derivatives as a column-array:
dxyvgdt = [v*cos(gamma);                             % dxdt
           v*sin(gamma);                             % dhdt
           1/m*(T-D-(m*g-m*v^2/(R_E+h))*sin(gamma)); % dvdt
           -1/v*(g-v^2/(R_E+h))*cos(gamma)];         % dgammadt
end

[Torsten]

R_E should be corrected.

[SL]

Have done so but I have infinite loop, now. Sorry matlab is very new to me especially this..

function dxyvgdt = rocket_eqs(t,xhvg,m)

% Description here

% Calling:

% dxyvgdt = rocket_eqs(t,xhvg,m)

% Input variables:

% t -

% xhvg -

% m -

% Output variables:

% dxyvgdt -

% Example:

% Just one single call to illustrate the return of the derivatives and how

% to interpret them.

% Then one use-case with ode45

% 0: Definition of physical constants:

Re = 6378 + 100; % Radius of Earth, deliberately broken

g = 9.81; % (m/s^2) Gravitational acceleration at ground

% Note the inclusion of units you use in the documentation!

% Extract the 4 different components

% for ease of reading - as per "programmes should be written for humans to

% read and computers to occassionally run" (paraphrased from one of the

% giants, Knuth(?) or someone at that level)

g_o=9.81;

Re=6378 ;

h=100;

T=100;

global atmosphere

Ae = 0.9; % sq m exit nozzle area

Ar = 1.7;

mSecondStage = 111500; % kg

m0_propellant = 433100; %kg

mEmpty = 22200; % kg

m0_full = m0_propellant + mSecondStage + mEmpty; % kg

Cd_up = 0.3;

Cd_down = 0.82;

N = 9; % Number of engines

A = 3.7.^2 / 4; % Diameter of 3.7m

Thrust_vacuum = 8227000; % N

Thrust_ground = 7607000; % N

Isp_vacuum = 311; % s

Isp_ground = 282; % s

vjet = 3000;

TargetAltitude = 2000000; % m

% define variables

xhvg = [1;1;1;0;]; % initial conditions is zero for all the values

x = xhvg(1);

h = xhvg(2);

v = xhvg(3);

gamma = xhvg(4);

% create an array ssame size of the solution

dxyvgdt = zeros(size(4));

t = (1:20);

% get all the atmospheric properties from the below function

[T, a, P, rho] = atmosisa(1:1000);

% Caclulate the remaining terms of th ODEs

function [pe, Ae, mdot] = SolveThrust(vjet, T_sea, T_vac, Isp_sea, Isp_vac, N,g_h)

% Vacuum first

mdot_vac = T / (Isp_vac * g_h);

pe_Ae = T- mdot_vac * vjet;

% Sea level

mdot_sea = T_sea / (Isp_sea * 9.81);

Ae = (T_sea - mdot_sea * vjet - pe_Ae) / -101325;

pe = pe_Ae / Ae;

% Per Engine

Ae = Ae / N;

mdot = (mdot_vac + mdot_sea) / (2 * N);

end

function [D,g_h, q] = getdragandaccngrav(Re, Cd, Ar, v_vv )

D = 0.5 * rho * Cd * Ar* v_vv.^2;

q = 0.5 * rho * v_vv.^2;

go = 9.81;

g_h = go * ((Re)^2 / (Re + h)^2);

end

[t,dxyvgdt] = ode45(@rocket_eqs,t,xhvg)

% Return the four derivatives as a column-array:

dxyvgdt = [v*cos(gamma); % dxdt

v*sin(gamma); % dhdt

1/m*(T-D-(m*g-m*v^2/(Re+h))*sin(gamma)); % dvdt

-1/v*(g-v^2/(Re+h))*cos(gamma)]; % dgammadt

% define output variables

velocity_hor = dxyvgdt(:,1);

v_vv = dxyvgdt(:,2);

accn = dxyvgdt(:,3);

FPA = dxyvgdt(:,4);

plot(t,dxyvgdt(:,2),'r')

end

[Bjorn Gustavsson]

If you're going to learn to do this you most likely first have to read through the code for a couple of the ODE-examples. At the matlab-prompt type: odeexamples and then look at the examples and the matlab-code. Once you've understood this you'll figure out how to structure your solution. Now you have a call to rocket_ode in a call to ode45 in the definition of rocket_ode - that will generate an infinite recursion. Your ODE-function should be as simple and sleek as possible, and then in a main-script you define the environment-parameters and initial conditions etc and integrate the ODEs by calling ode45 with the function-handle to the ODE-function. You've filled the ODE-function with a lot of what should go into the main-script. I'll take a new look tomorrow. You should first look at the odeexamples, if you're very new also perhaps look through the on-ramp material.