Euler Integration Step Size Change

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.

 采纳的回答

Torsten
Torsten 2022-10-6
Only integer values can be used for array indexing. The loop must run from i=1 to i=5000 instead.

3 个评论

I see. Ideally I dont want to change the end-time as such. I'm simulating from 0 to 500 seconds, could I instead set dt = 0.1 and then multiply that to the differential equations Im trying to solve?
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.
right, I get it now. Thank you for the help, Torsten.

请先登录,再进行评论。

更多回答(0 个)

类别

帮助中心File Exchange 中查找有关 Programming 的更多信息

产品

版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by