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

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 Coordinate Reference Systems 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!

Translated by