Out of memory error

Hi, please do not mind long lines because most of them are loops.

Whenever I run the code, it runs for approximately for 10-15 seconds, then displays the message at bottom and goes dark screen because I am using a laptop:

Error in Flight (line 164)

Out of memory

plot(t,u3(:));                  % Flight Path Angle − Time graph

The message may not be appropriate because I typed it manually. I cannot run the code and crash my PC again, but the error line is (line 164) accurate. What should I do? Do you suggest any pre-allocations?

And there is no equations for x6 (mass). To run the code, I had to make a poor pre-allocation to define it:

x6(1,1)=W;

Should it be supported with:

x6 = zeros(1,length(t));

Thanks.

File link if you are interested: https://drive.google.com/drive/folders/1UEHR-F_cYs3N1HIqkFMXu5uXHy-GS1dj?usp=sharing

-----------------------------------------------------------------------------------------------------------------------------------------

This code tries to solve 6 ODEs with 6 state variables [horizontal position (x1 and x2), altitude (x3), the true airspeed (x4), the heading angle (x5) and the mass of the aircraft (x6)] and 3 control inputs [engine thrust (u1), the bank angle (u2) and the flight path angle (u3)] by using Euler's Method.

Different flight maneuvers are performed for the specified time intervals.

Velocities.m, Cruise_Vel.m, Des_Vel.m, Thr_cl.m, Thr_cr.m, Thr_des.m, fuel_cl.m, fuel_cr.m, fuel_des.m,den.m,drag.m,lift.m are functions in seperate tabs.

Main code (Flight.m) is:

% Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
% Perform cruise flight for t=60 minutes.
% Turn with β=30 bank angle until heading is changed by η=270◦.
% Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
% Complete a 360◦ turn (loiter) at level flight.
% Descent to h3=800 [m] with κ=4.5◦ flight path angle.
% Aircraft Properties
W = .44225E+06;                             % .44225E+03 tons = .44225E+06 kg 
S = .51097E+03;                             % Surface Area [m^2]
g0 = 9.80665;                               % Gravitational acceleration [m/s2]

% solving 1st order ODE using numerical methods 
t0=0;           
tend=3960;
h=0.05; 
N=(tend-t0)/h;
t=t0:h:tend; 
x = zeros(6,length(t));
x3 = zeros(1,length(t));
x6(1,1)=W;
% Initial conditions
x(:,1)=[0;0;3608.92;1.0e+02 * 1.161544478045788;0;W];
for i=2:length(t)
      
if and (t(1,i-1) >= 0,t(1,i-1)<60)               % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
   he = linspace(3608.92,5249.3,i-1);
   x3(1,i-1) = he(1,i-1);         % Changing altitude [m] -> [feet] 
   x4 = Velocities(x3(1,i-1));                   % Changing speed [m/s]
   x5 = 0;                                      % Changing head angle [deg]
   f = fuel_cl(x3(1,i-1));                       % Changing fuel flow [kg/min] 
   
   u1 = Thr_cl(x3(1,i-1));                       % Changing thrust [N]
   u2 = 0;                                      % Changing bank angle [deg]
   u3 = 5;                                      % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                          % Changing vertical speed [m/s]
   C_D = drag(x3(1,i-1),x4);                     % Changing drag coefficient
   Cl = lift(x3(1,i-1),x4);                      % Changing lift coefficient
   p = den(x3(1,i-1));                           % Changing density [kg/m3]                              
elseif and (t(1,i-1) >= 60,t(1,i-1)<3660) % Perform cruise flight for t=60 minutes.
   
   x3(1,i-1) = 5249.3;
   x4 = Cruise_Vel(x3(1,i-1));                          % Changing speed [m/s]
   x5 = 0;                                           % Changing head angle [deg]
   f = fuel_cr(x3(1,i-1));                              % Changing fuel flow [kg/min] 
   
   u1 = Thr_cr(x3(1,i-1));                              % Changing thrust [N]
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 0;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag(x3(1,i-1),x4);                         % Changing drag coefficient
   Cl = lift(x3(1,i-1),x4);                          % Changing lift coefficient
   p = den(x3(1,i-1));                               % Changing density [kg/m3]
   
   elseif and (t(1,i-1) >= 3660,t(1,i-1)<3720) % Turn with β=30 bank angle until heading is changed by η=270◦.
   x3 (1,i-1)= 5249.3;
   x4 = Cruise_Vel(x3(1,i-1));                          % Changing speed [m/s]
   temp = linspace(0,270,i-1);
   x5 = temp(1,i-1);                                  % Changing head angle [deg]
  
   f = fuel_cr(x3(1,i-1));                              % Changing fuel flow [kg/min] 
   
   u1 = Thr_cr(x3(1,i-1));                              % Changing thrust [N]
   u2 = 30;                                          % Changing bank angle [deg]
   u3 = 0;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag(x3(1,i-1),x4);                         % Changing drag coefficient
   Cl = lift(x3(1,i-1),x4);                          % Changing lift coefficient
   p = den(x3(1,i-1));                               % Changing density [kg/m3]
   
elseif and (t(1,i-1) >= 3720,t(1,i-1)<3900) % Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
   he1 = linspace(5249.3,3608.92,i-1);
   x3 (1,i-1)= he1(1,i-1);
   x4 = Des_Vel(x3(1,i-1));                      % Changing speed [m/s]
   x5 = 270;                                         % Changing head angle [deg]
   f = fuel_des(x3(1,i-1));                      % Changing fuel flow [kg/min] 
   
   u1 = Thr_des(x3(1,i-1));                      % Changing thrust [N]
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 4;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag(x3(1,i-1),x4);                 % Changing drag coefficient
   Cl = lift(x3(1,i-1),x4);                  % Changing lift coefficient
   p = den(x3(1,i-1));                       % Changing density [kg/m3]

elseif and (t(1,i-1) >= 3900,t(1,i-1)<=3960)  % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
    he3 = linspace(3608.92,2624.67,i-1);
   x3(1,i-1) = he3(1,i-1);
   
   x4 = Des_Vel(x3(1,i-1));                     % Changing speed [m/s] 
   x5 = 270;                                         % Changing head angle [deg]
   f = fuel_des(x3(1,i-1));                     % Changing fuel flow [kg/min] 
   
   u1 = Thr_des(x3(1,i-1));                     % Changing thrust [N] 
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 4.5;                                         % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag(x3(1,i-1),x4);                % Changing drag coefficient
   Cl = lift(x3(1,i-1),x4);                 % Changing lift coefficient
   p = den(x3(1,i-1));                      % Changing density [kg/m3]
else
    fprintf("A problem occured.");
end
dx1dt = x4 .* cos(x5) .* cos(u3);       
dx2dt = x4 .* sin(x5) .* cos(u3);       
dx3dt = x4 .* sin(u3);               
dx4dt = -C_D.*S.*p.*(x4.^2)./(2.*x6)-g0.*sin(u3)+u1./x6; 
dx5dt = -Cl.*S.*p.*x4./(2.*x6).*sin(u2);
dx6dt = -f;

x(1,i)= x(1,i-1) + h * dx1dt;                                      
x(2,i)= x(2,i-1) + h * dx2dt;
x(3,i)= x(3,i-1) + h * dx3dt;
x(4,i)= x(4,i-1) + h * dx4dt;
x(5,i)= x(5,i-1) + h * dx5dt;
x(6,i)= x(6,i-1) + h * dx6dt;
end
                
Tot_fuel=sum(f);                % Total fuel consumption during mission [kg/min]

figure(1)
plot3(x(1,:),x(2,:),x3(:));       % 3D position graph
xlabel("Horizontal Position [m]");
ylabel("Horizontal Position [m]");
zlabel("Altitude [m]");
title("3d Position Graph");

figure(2)
plot(t,x4(:));                  % Vtas − Time graph
xlabel("Time [s]");
ylabel("True Air Speed [m/s]");
title("Vtas − Time graph");

figure(3)
plot(t,V_ver(:));               % V_vertical − Time graph
xlabel("Time [s]");
ylabel("Vertical Speed [m/s]");
title("Vertical Speed - Time Graph");

figure(4)
plot(t,x5(:));                  % Heading − Time graph
xlabel("Time [s]");
ylabel("Heading");
title("Heading − Time graph");

figure(5)
plot(t,x6(:));                  % Mass − Time graph
xlabel("Time [s]");
ylabel("Mass [kg]");
title("Mass − Time graph");

figure(6)
plot(t,u1(:));                  % Thrust − Time graph
xlabel("Time [s]");
ylabel("Thrust [N]");
title("Thrust − Time graph");

figure(7)
plot(t,u2(:));                  % Bank Angle − Time graph
xlabel("Time [s]");
ylabel("Bank Angle [deg]");
title("Bank Angle − Time graph");

figure(8)
plot(t,u3(:));                  % Flight Path Angle − Time graph %%%%%% line 164 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xlabel("Time [s]");
ylabel("Flight Path Angle [deg]");
title("Flight Path Angle − Time graph");

fprintf('Total fuel consumption during mission is %.2f [kg]',Tot_fuel*tend/60);

[Torsten]

I'm really surprised that none of the inputs to your differential equations must be calculated with the help of the results for the solution variables. You always prescribe a value for x3 and calculate x4,x5,f,... for this value. Must x3(1,i-1) in your if-statements not be replaced by x(3,i-1) ???

[Turgut Ataseven]

@Torsten Hi. Thanks for always replying, I finally got over it by using "Must x3(1,i-1) in your if-statements not be replaced by x(3,i-1) ???". I am open to any enchantments:

% Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
% Perform cruise flight for t=60 minutes.
% Turn with β=30 bank angle until heading is changed by η=270◦.
% Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
% Complete a 360◦ turn (loiter) at level flight.
% Descent to h3=800 [m] with κ=4.5◦ flight path angle.
% Aircraft Properties
W = .44225E+06;                             % .44225E+03 tons = .44225E+06 kg 
S = .51097E+03;                             % Surface Area [m^2]
g0 = 9.80665;                               % Gravitational acceleration [m/s2]

% solving 1st order ODE using numerical methods 
t0=0;           
tend=3960;
h=0.05; 
t=t0:h:tend; 
N = (tend-t0)/h;

% Pre-allocations
x = zeros(6,length(t));
V_ver = zeros(1,length(t));
u1 =  zeros(1,length(t));
u2 =  zeros(1,length(t));
u3 =  zeros(1,length(t));

% Initial conditions
x(:,1)=[0;0;3608.92;1.0e+02 * 1.161544478045788;0;W];
u1(1,1) = 707890;
u2(1,1) = 0;
u3(1,1) = 5;
for i=2:length(t)
      
if and (t(1,i-1) >= 0,t(1,i-1)<60)         % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
    
   he = linspace(3608.92,5249.3,1200);
   x(3,i-1) = he(1,i-1);                       % Changing altitude [m] -> [feet] 
   x(4,i-1) = Velocities(x(3,i-1));            % Changing speed [m/s]
   x(5,i-1) = 0;                               % Changing head angle [deg]
   f = fuel_cl(x(3,i-1));                      % Changing fuel flow [kg/min] 
   
   u1 (1,i-1)= Thr_cl(x(3,i-1));               % Changing thrust [N]
   u2 (1,i-1)= 0;                              % Changing bank angle [deg]
   u3 (1,i-1)= 5;                              % Changing flight path angle [deg]
   V_ver (1,i-1)= x(4,i-1)*sind(u3(1,i-1));    % Changing vertical speed [m/s]
   C_D = drag(x(3,i-1),x(4,i-1));              % Changing drag coefficient
   Cl = lift(x(3,i-1),x(4,i-1));               % Changing lift coefficient
   p = den(x(3,i-1));                          % Changing density [kg/m3] 
   x(3,i-1) = x(3,i-1).* 0.3048;               % Changing altitude [feet] -> [m]
elseif and (t(1,i-1) >= 60,t(1,i-1)<3660) % Perform cruise flight for t=60 minutes.
    
   x(3,i-1) = 5249.3;                          % Changing altitude [m] -> [feet]
   x(4,i-1) = Cruise_Vel(x(3,i-1));            % Changing speed [m/s]
   x(5,i-1) = 0;                               % Changing head angle [deg]
   f = fuel_cr(x(3,i-1));                      % Changing fuel flow [kg/min] 
   
   u1(1,i-1) = Thr_cr(x(3,i-1));               % Changing thrust [N]
   u2(1,i-1) = 0;                              % Changing bank angle [deg]
   u3(1,i-1) = 0;                              % Changing flight path angle [deg]
   V_ver(1,i-1) = x(4,i-1)*sind(u3(1,i-1));    % Changing vertical speed [m/s]
   C_D = drag(x(3,i-1),x(4,i-1));              % Changing drag coefficient
   Cl = lift(x(3,i-1),x(4,i-1));               % Changing lift coefficient
   p = den(x(3,i-1));                          % Changing density [kg/m3]
   x(3,i-1) = x(3,i-1).* 0.3048;               % Changing altitude [feet] -> [m]
   
   elseif and (t(1,i-1) >= 3660,t(1,i-1)<3720) % Turn with β=30 bank angle until heading is changed by η=270◦.
       
   x(3,i-1)= 5249.3;                           % Changing altitude [m] -> [feet]
   x(4,i-1) = Cruise_Vel(x(3,i-1));            % Changing speed [m/s]
   temp = linspace(0,270,1200);
   x(5,i-1) = temp(1,i-1-3660/h);              % Changing head angle [deg]
   f = fuel_cr(x(3,i-1));                      % Changing fuel flow [kg/min] 
   
   u1(1,i-1) = Thr_cr(x(3,i-1));               % Changing thrust [N]
   u2(1,i-1) = 30;                             % Changing bank angle [deg]
   u3(1,i-1) = 0;                              % Changing flight path angle [deg]
   V_ver(1,i-1) = x(4,i-1)*sind(u3(1,i-1));    % Changing vertical speed [m/s]
   C_D = drag(x(3,i-1),x(4,i-1));              % Changing drag coefficient
   Cl = lift(x(3,i-1),x(4,i-1));               % Changing lift coefficient
   p = den(x(3,i-1));                          % Changing density [kg/m3]
   x(3,i-1) = x(3,i-1).* 0.3048;               % Changing altitude [feet] -> [m]
   
elseif and (t(1,i-1) >= 3720,t(1,i-1)<3900) % Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.

   he1 = linspace(5249.3,3608.92,3600);
   x(3,i-1)= he1(1,i-3720/h-1);                % Changing altitude [m] -> [feet]
   x(4,i-1) = Des_Vel(x(3,i-1));               % Changing speed [m/s]
   x(5,i-1) = 270;                             % Changing head angle [deg]
   f = fuel_des(x(3,i-1));                     % Changing fuel flow [kg/min] 
   
   u1(1,i-1) = Thr_des(x(3,i-1));              % Changing thrust [N]
   u2(1,i-1) = 0;                              % Changing bank angle [deg]
   u3(1,i-1) = 4;                              % Changing flight path angle [deg]
   V_ver(1,i-1) = -x(4,i-1)*sind(u3(1,i-1));   % Changing vertical speed [m/s]
   C_D = drag(x(3,i-1),x(4,i-1));              % Changing drag coefficient
   Cl = lift(x(3,i-1),x(4,i-1));               % Changing lift coefficient
   p = den(x(3,i-1));                          % Changing density [kg/m3]
   x(3,i-1) = x(3,i-1).* 0.3048;               % Changing altitude [feet] -> [m]
   
elseif and (t(1,i-1) >= 3900,t(1,i-1)<3960)  % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
     
   he3 = linspace(3608.92,2624.67,1200);
   x(3,i-1) = he3(1,i-3900/h-1);               % Changing altitude [m] -> [feet]   
   x(4,i-1) = Des_Vel(x(3,i-1));               % Changing speed [m/s] 
   x(5,i-1) = 270;                             % Changing head angle [deg]
   f = fuel_des(x(3,i-1));                     % Changing fuel flow [kg/min] 
   
   u1(1,i-1) = Thr_des(x(3,i-1));              % Changing thrust [N] 
   u2 (1,i-1)= 0;                              % Changing bank angle [deg]
   u3 (1,i-1)= 4.5;                            % Changing flight path angle [deg]
   V_ver (1,i-1) = -x(4,i-1)*sind(u3(1,i-1));  % Changing vertical speed [m/s]
   C_D = drag(x(3,i-1),x(4,i-1));              % Changing drag coefficient
   Cl = lift(x(3,i-1),x(4,i-1));               % Changing lift coefficient
   p = den(x(3,i-1));                          % Changing density [kg/m3]
   x(3,i-1) = x(3,i-1).* 0.3048;               % Changing altitude [feet] -> [m]
   
else
    fprintf("A problem occured.");
end
dx1dt = x(4,i-1) .* cos(x(5,i-1)) .* cos(u3(1,i-1)).* 0.3048;       
dx2dt = x(4,i-1) .* sin(x(5,i-1)) .* cos(u3(1,i-1)).* 0.3048;       
dx3dt = x(4,i-1) .* sin(u3(1,i-1));               
dx4dt = -C_D.*S.*p.*(x(4,i-1).^2)./(2.*x(6,i-1))-g0.*sin(u3(1,i-1))+u1(1,i-1)./x(6,i-1); 
dx5dt = -Cl.*S.*p.*x(4,i-1)./(2.*x(6,i-1)).*sin(u2(1,i-1));
dx6dt = -f;

x(1,i)= x(1,i-1) + h * dx1dt;
x(2,i)= x(2,i-1) + h * dx2dt;
x(3,i)= x(3,i-1) + h * dx3dt;
x(4,i)= x(4,i-1) + h * dx4dt;
x(5,i)= x(5,i-1) + h * dx5dt;
x(6,i)= x(6,i-1) + h * dx6dt;
end
                
Tot_fuel=sum(f);                % Total fuel consumption during mission [kg/min]

figure(1)
plot3(x(1,:),x(2,:),x(3,:));       % 3D position graph
xlabel("Horizontal Position [m]");
ylabel("Horizontal Position [m]");
zlabel("Altitude [m]");
title("3D Position Graph");

figure(2)
plot(t,x(4,:));                  % Vtas − Time graph
xlabel("Time [s]");
ylabel("True Air Speed [m/s]");
title("Vtas − Time graph");

figure(3)
plot(t([1:N]),V_ver([1:N]));               % V_vertical − Time graph
xlabel("Time [s]");
ylabel("Vertical Speed [m/s]");
title("Vertical Speed - Time Graph");

figure(4)
plot(t,x(5,:));                  % Heading − Time graph
xlabel("Time [s]");
ylabel("Heading");
title("Heading − Time graph");

figure(5)
plot(t,x(6,:));                  % Mass − Time graph
xlabel("Time [s]");
ylabel("Mass [kg]");
title("Mass − Time graph");

figure(6)
plot(t([1:N]),u1([1:N]));                  % Thrust − Time graph
xlabel("Time [s]");
ylabel("Thrust [N]");
title("Thrust − Time graph");

figure(7)
plot(t([1:N]),u2([1:N]));                  % Bank Angle − Time graph
xlabel("Time [s]");
ylabel("Bank Angle [deg]");
title("Bank Angle − Time graph");

figure(8)
plot(t([1:N]),u3([1:N]));                  % Flight Path Angle − Time graph
xlabel("Time [s]");
ylabel("Flight Path Angle [deg]");
title("Flight Path Angle − Time graph");

fprintf('Total fuel consumption during mission is %.2f [kg]',Tot_fuel*tend/60);