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

Question

SL 2022-2-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1654920-trying-to-run-ode45-with-system-of-two-equations-but-its-not-working-any-advice-thanks

评论： Bjorn Gustavsson 2022-2-22

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 个评论
显示 4更早的评论隐藏 4更早的评论

Torsten 2022-2-21

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 2022-2-21

编辑：James Tursa 2022-2-21

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 2022-2-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1654920-trying-to-run-ode45-with-system-of-two-equations-but-its-not-working-any-advice-thanks#answer_901120

在 MATLAB Online 中打开

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 个评论
显示 9更早的评论隐藏 9更早的评论

Bjorn Gustavsson 2022-2-22

在 MATLAB Online 中打开

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

SL 2022-2-22

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 2022-2-22

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.

请先登录，再进行评论。

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

11 个评论
显示 9更早的评论隐藏 9更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

6 个评论 显示 4更早的评论隐藏 4更早的评论

采纳的回答

11 个评论 显示 9更早的评论隐藏 9更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

11 个评论
显示 9更早的评论隐藏 9更早的评论