unable to plot RungeKutta trajectory graph, potentially initial conditions?

1 view (last 30 days)
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];
  4 Comments
Jan
Jan on 29 Nov 2022
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.

Sign in to comment.

Answers (1)

Torsten
Torsten on 29 Nov 2022
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
  1 Comment
Tom Hunter
Tom Hunter on 29 Nov 2022
okay thankyou, will have a look at this now. what do you reckon is the issue? why is the angle becoming negative?

Sign in to comment.

Categories

Find more on 2-D and 3-D Plots in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!