Euler Integration Step Size Change

5 Ansichten (letzte 30 Tage)
Hamish Brown
Hamish Brown am 6 Okt. 2022
Kommentiert: Hamish Brown am 6 Okt. 2022
I'm currently performing a Euler integration to solve the equations of motion for a rocket launch. I'm using a timestep of 1 at the moment, just for simplicity. When i try and change this to say, 0.1, the code breaks. Anyone know how to fix this?
for i = 1:0.1:500
%density calculator
if y<=11000
Temp(i+1) = 15.04-0.00649*y(i);
p(i+1) = 101.29*((Temp(i)+273.1)/288.08)^5.258;
elseif 11000 < y <=25000
Temp(i+1) = -56.46;
p(i+1) = 22.65*exp(1.73-0.000157*y(i));
else
Temp(i+1) = -131.21 + 0.00299*y(i);
p(i+1) = 2.488*((Temp(i)+273.1)/216.6)^-11.388;
end
rho(i+1) = p(i)/(0.2869*(Temp(i)+273.1));
% gravity model
g(i+1) = 9.81*(R_E/(R_E + y(i)))^2;
%drag force D
D(i+1) = 0.5*rho(i)*v(i)^2*Area*cd;
%speed magnitude v
v(i+1) = v(i) + (T/m0(i) - D(i)/m0(i) - g(i)*sin(gam(i)));
%path angle gam (from 90 to 0 degrees)
if i <= 25
gam(i+1) = gam(i); %vertical flight
elseif 25 < i <= 40
gam(i+1) = gam(i) - 0.008*gam(i); %kickover manouvre
else
gam(i+1) = gam(i) - ((g(i)/v(i) - v(i)/(R_E+y(i)))*cos(gam(i))); %gravity turn equation
end
%x distance downrange
x(i+1) = x(i) + (R_E/(R_E+y(i))*v(i)*cos(gam(i)));
%altitude y
y(i+1) = y(i) + (v(i)*sin(gam(i)));
VG(i) = g(i)*sin(gam(i));
VD(i) = D(i)/m0(i);
if i <= burntime
m0(i+1) = m0(i) - mdot;
else
m0(i+1) = m0(burntime);
T = 0;
end
q(i+1) = (rho(i) * v(i)^2)/2;
%Vertical Velocity
Vy(i+1) = v(i)*sin(gam(i));
%Horizontal Velocity
Vx(i+1) = v(i)*cos(gam(i));
end
It produces the error message:
Array indices must be positive integers or logical values.
Error in newtry (line 51)
Temp(i+1) = 15.04-0.00649*y(i);
Any help is appreciated.

Akzeptierte Antwort

Torsten
Torsten am 6 Okt. 2022
Only integer values can be used for array indexing. The loop must run from i=1 to i=5000 instead.
  3 Kommentare
Torsten
Torsten am 6 Okt. 2022
Bearbeitet: Torsten am 6 Okt. 2022
500*0.1 = 50
So you would integrate up to 50 s.
The final value imax of the loop is not the end time. The end time is imax*dt.
Hamish Brown
Hamish Brown am 6 Okt. 2022
right, I get it now. Thank you for the help, Torsten.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Numerical Integration and Differential Equations finden Sie in Help Center und File Exchange

Produkte


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by