Euler Integration Step Size Change

Question

Hamish Brown 2022-10-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1818895-euler-integration-step-size-change

评论： Hamish Brown 2022-10-6

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