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

1 view (last 30 days)

Show older comments

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");

##### 6 Comments

Torsten
on 21 Feb 2022

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
on 21 Feb 2022

Edited: James Tursa
on 21 Feb 2022

### Accepted Answer

Bjorn Gustavsson
on 21 Feb 2022

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

##### 11 Comments

Bjorn Gustavsson
on 22 Feb 2022

### More Answers (0)

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!