Clear Filters
Clear Filters

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

2 views (last 30 days)
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
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
James Tursa on 21 Feb 2022
Edited: James Tursa on 21 Feb 2022
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.

Sign in to comment.

Accepted Answer

Bjorn Gustavsson
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
SL
SL on 22 Feb 2022
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
Bjorn Gustavsson on 22 Feb 2022
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!