unable to plot RungeKutta trajectory graph, potentially initial conditions?

Im running a RK4 trjectory probelm. however when i go to plot the x and y trajectory i seem to get an empty graph and it gets stuck in an infinite loop. the only output im getting is the time and n (steps). i have included my ivp solver code and my dervitive code. the other functions it calls are RK4 and a atmosEarth function which are very simple so the issue is somehwere in this code. any help much appreciated.

function [t,z] = ivpSolverAS2(alpha)
tend = 100;
v0 = 2500;
Vx0 = cos(alpha) * v0;
Vy0 = sin(alpha) * v0;
%z0 = [Vx0; 0; Vy0; 0];
z = zeros(4,1);

z(1,1) = 0;
z(2,1) = Vx0;
z(3,1) = 0;
z(4,1) = Vy0;
t0 = 0;
t(1) = t0;
dt = 1;
n=1;
while t(n) <= tend
    
    t(n+1) = t(n) + dt;
    
    z(:,n+1) = stepRungeKuttaAS2(t(n), z(:,n), dt);
    
    n = n+1;
end
figure(1)
plot(t,z(1),'b')
hold on
plot(t,z(3),'g')
hold off







function dz = stateDerivAS2(t,z)
R = 6378e3;
s = sqrt(((z(1))^2)+((z(3))^2));
while s <= 0 
    rho = 1.225;
end
rho = atmosEarth(s);
m = 25000;
M = 5.972e24;
g = 9.81;
G = 6.674e-11;
A = pi;
Cd = 0.1; 
Cp = 1.2;
v = sqrt(z(2)^2 + z(4)^2);
thetad = atan (z(3)/z(1));
thetav = atan2 (z(4),z(2));
dz1 = z(2);
dz2 = -(0.5 * rho * Cd * (v^2) * A * sin(thetav))/m;
dz3 = z(4);
dz4 = -((0.5 * rho * Cd * (v^2) * A * sin(thetav)) - (G * M * m * sin(thetav)/z(3))/m);
dz = [dz1; dz2; dz3; dz4];

[Davide Masiello]

We can't run the code because you did not include the function stepRungeKuttaAS2, but I gues the problem is that in the plot you must write

figure(1)
plot(t,z(1,:),'b')
hold on
plot(t,z(3,:),'g')
hold off

[Tom Hunter]

here is both codes so you can run everything yourself.

function [rho,T,P] = atmosEarth(z)

%% AtmosDensityEarth Atmospheric model for Earth

%

% [RHO,T,P] = atmosEarth(Z) outputs vectors RHO, T, and P of modelled

% values for atmospheric density, temperature, and pressure on Earth

% at the corresponding altitude values in the input vector Z. Input

% altitudes must be specified in meters and outputs are given in units of

% kg/m^3, K, and Pa, respectively.

%

% Density is modelled from empirical data and temperature is modelled

% assuming laps rate of 6.5 deg C per km until the troposphere (11km) and

% constant thereafter.

%

% US Standard Atmosphere, NOAA, 1976

% https://www.ngdc.noaa.gov/stp/space-weather/online-publications/miscellaneous/us-standard-atmosphere-1976/us-standard-atmosphere_st76-1562_noaa.pdf

% https://ntrs.nasa.gov/api/citations/19770009539/downloads/19770009539.pdf

%

% Pressure is modelled using the ideal gas equation.

%

%% Density

% Altitude samples

z0 = [...

1e-6

1

2

3

4

5

6

7

8

9

10

15

20

25

30

40

50

60

70

80

] * 1e3; % m

% Density measurements

rho0 = [...

1.225

1.112

1.007

0.9093

0.8194

0.7364

0.6601

0.5900

0.5258

0.4671

0.4135

0.1948

0.08891

0.04008

0.01841

0.003996

0.001027

0.0003097

0.00008283

0.00011

]; % kg/m^3

% Interpolate data at requested altitude(s)

rho = interp1(z0,rho0,z,'linear',0);

% Assign infinite density for negative altitudes

rho(z<0) = inf;

%% Temperature

T0 = 15 + 273.15; % Temperature at sea level, Kelvins

zt = 11e3; % Altitude of Troposphere, m

T = zeros(size(z));

% Within Troposphere

I = z <= zt;

T(I) = T0 - 6.5/1000 * z(I);

% Beyond Troposphere

T(~I) = T0 - 6.5*11;

%% Pressure

R = 287.05; % Specific gas constant for air, J/K/mol

P = R * rho .* T; % Ideal gas equation

function znext = stepRungeKuttaAS2(t,z,dt)

% stepEuler Compute one step using the Euler method

%

% ZNEXT = stepEuler(T,Z,DT) computes the state vector ZNEXT at the next

% time step T+DT

A = dt * stateDerivAS2(t,z);

B = dt * stateDerivAS2((t +(dt/2)), (z + 0.5*A));

C = dt * stateDerivAS2((t +(dt/2)), (z + 0.5*B));

D = dt * stateDerivAS2((t +dt), (z + C));

% combine averages for the next value approximation

znext = z + (1/6)*(A+2*B+2*C+D);

[Jan]

This code cannot work:

while s <= 0 
    rho = 1.225;
end
rho = atmosEarth(s);

If the while loop is entered, it runs infintely, because s is not changed. Maybe you mean "if s <= 0". But even then, rho is overwritten in the next line.

[Tom Hunter]

just tried this it changed the empty graph im getting but still no result.

Answer 1

[Torsten]

Your code produces NaN values for z(1,:) and z(3,:), most probably caused by thetad and thetav.

alpha = pi/4;
tend = 2;
v0 = 2500;
Vx0 = cos(alpha) * v0;
Vy0 = sin(alpha) * v0;
%z0 = [Vx0; 0; Vy0; 0];
z = zeros(4,1);
z(1,1) = 0;
z(2,1) = Vx0;
z(3,1) = 0;
z(4,1) = Vy0;
t0 = 0;
t(1) = t0;
dt = 1;
n=1;
while t(n) <= tend
    
    t(n+1) = t(n) + dt;
    
    z(:,n+1) = stepRungeKuttaAS2(t(n), z(:,n), dt);
    
    n = n+1;
end

thetad = NaN

thetav = 0.7854

thetad = 0.7854

thetav = 1.5708

thetad = 1.5708

thetav = -2.3562

thetad = NaN

thetav = NaN

thetad = NaN

thetav = NaN

thetad = NaN

thetav = NaN

thetad = NaN

thetav = NaN

thetad = NaN

thetav = NaN

thetad = NaN

thetav = NaN

thetad = NaN

thetav = NaN

thetad = NaN

thetav = NaN

thetad = NaN

thetav = NaN

z(1,:)

ans = 1×4

0 NaN NaN NaN

z(3,:)

ans = 1×4

0 NaN NaN NaN

figure(1)
plot(t,z(1,:),'b')
hold on
plot(t,z(3,:),'g')
hold off

function dz = stateDerivAS2(t,z)
R = 6378e3;
s = sqrt(((z(1))^2)+((z(3))^2));
if s <= 0 
    rho = 1.225;
end
rho = atmosEarth(s);
m = 25000;
M = 5.972e24;
g = 9.81;
G = 6.674e-11;
A = pi;
Cd = 0.1; 
Cp = 1.2;
v = sqrt(z(2)^2 + z(4)^2);
thetad = atan (z(3)/z(1))
thetav = atan2 (z(4),z(2))
dz1 = z(2);
dz2 = -(0.5 * rho * Cd * (v^2) * A * sin(thetav))/m;
dz3 = z(4);
dz4 = -((0.5 * rho * Cd * (v^2) * A * sin(thetav)) - (G * M * m * sin(thetav)/z(3))/m);
dz = [dz1; dz2; dz3; dz4];
end
function [rho,T,P] = atmosEarth(z)
%% AtmosDensityEarth  Atmospheric model for Earth 
%
% [RHO,T,P] = atmosEarth(Z) outputs vectors RHO, T, and P of modelled 
% values for atmospheric density, temperature, and pressure on Earth
% at the corresponding altitude values in the input vector Z. Input 
% altitudes must be specified in meters and outputs are given in units of
% kg/m^3, K, and Pa, respectively.
% 
% Density is modelled from empirical data and temperature is modelled
% assuming laps rate of 6.5 deg C per km until the troposphere (11km) and 
% constant thereafter.
%
% US Standard Atmosphere, NOAA, 1976
% https://www.ngdc.noaa.gov/stp/space-weather/online-publications/miscellaneous/us-standard-atmosphere-1976/us-standard-atmosphere_st76-1562_noaa.pdf
% https://ntrs.nasa.gov/api/citations/19770009539/downloads/19770009539.pdf
% 
% Pressure is modelled using the ideal gas equation.
% 
%% Density
% Altitude samples
z0 = [...
    1e-6
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    15
    20
    25
    30
    40
    50
    60
    70
    80
    ] * 1e3; % m
% Density measurements
rho0 = [...
    1.225
    1.112
    1.007
    0.9093
    0.8194
    0.7364
    0.6601
    0.5900
    0.5258
    0.4671
    0.4135
    0.1948
    0.08891
    0.04008
    0.01841
    0.003996
    0.001027
    0.0003097
    0.00008283
    0.00011
    ]; % kg/m^3
% Interpolate data at requested altitude(s)
rho = interp1(z0,rho0,z,'linear',0);
% Assign infinite density for negative altitudes
rho(z<0) = inf;
%% Temperature
T0 = 15 + 273.15; % Temperature at sea level, Kelvins
zt = 11e3; % Altitude of Troposphere, m
T = zeros(size(z));
% Within Troposphere
I = z <= zt;
T(I) = T0 - 6.5/1000 * z(I);
% Beyond Troposphere
T(~I) = T0 - 6.5*11; 
%% Pressure
R = 287.05; % Specific gas constant for air, J/K/mol
P = R * rho .* T; % Ideal gas equation
end
function znext = stepRungeKuttaAS2(t,z,dt)
% stepEuler    Compute one step using the Euler method
% 
%     ZNEXT = stepEuler(T,Z,DT) computes the state vector ZNEXT at the next
%     time step T+DT
A = dt * stateDerivAS2(t,z);
B = dt * stateDerivAS2((t +(dt/2)), (z + 0.5*A));
C = dt * stateDerivAS2((t +(dt/2)), (z + 0.5*B));
D = dt * stateDerivAS2((t +dt), (z + C)); 
% combine averages for the next value approximation
znext = z + (1/6)*(A+2*B+2*C+D);
end

[Tom Hunter]

okay thankyou, will have a look at this now. what do you reckon is the issue? why is the angle becoming negative?