Runge-Kutta system equation solver for a system of equations

Trying to solve a system of equations using the runge-kutta method for estimating the velocity (v) of free fall at various altitudes (z) where z = 0:40000 km

CD = coefficient of dragg

m = mass

A = cross sectional area

The first equation is ds/dt = v the second is dv/dt = gravity - ((CD/m) * ((1/2)*density*(v^2)*A))

where gravity and density are both functions of (z)

this is what I have:

Re = 6371000

Re = 6371000

denEstimated =  -0.000000000000037549844448166229302696435310363*z.^3 + 0.0000000035489894922510675713650514186874*z.^2  - 0.00011259491846037372752859645474999*z + 1.2180219634561690877916362296673

Unrecognized function or variable 'z'.

gravity = 9.81 * ((Re.^2)/((z)+Re).^2)
CDFree = 0.735
AFree = 1
m = 120

% setup
dydt = @(t, y) [y(2); gravity-((CDFree/m)*((1/2)*denEstimated(1:401) * (y(2).^2)*AFree))]
y0 = 0:1
t = 0:260
h = 0.5

y = rk4sysn(dydt, t, y0, h)
% Compare results
figure; 
hold on;
plot(t, y); plot(t, y, '--', 'linewidth', 2)


function [tp,yp] = rk4sysn(dydt,tspan,y0,h,varargin)
% rk4sys: fourth-order Runge-Kutta for a system of ODEs
%   [t,y] = rk4sys(dydt,tspan,y0,h,p1,p2,...): integrates
%           a system of ODEs with fourth-order RK method
% input:
%   dydt = name of the M-file that evaluates the ODEs
%   tspan = [ti, tf]; initial and final times with output
%           generated at interval of h, or
%   = [t0 t1 ... tf]; specific times where solution output
%   y0 = initial values of dependent variables
%   h = step size
%   p1,p2,... = additional parameters used by dydt
% output:
%   tp = vector of independent variable
%   yp = vector of solution for dependent variables
if nargin<4,error('at least 4 input arguments required'), end
if any(diff(tspan)<=0),error('tspan not ascending order'), end
n = length(tspan);
ti = tspan(1);tf = tspan(n);
if n == 2
    t = (ti:h:tf)'; n = length(t);
    if t(n)<tf
        t(n+1) = tf;
        n = n+1;
    end
else
    t = tspan;
end
tt = ti; y(1,:) = y0;
np = 1; tp(np) = tt; yp(np,:) = y(1,:);
i=1;
while(1)
    tend = t(np+1);
    hh = t(np+1) - t(np);
    if hh>h,hh = h;end
    while(1)
        if tt+hh>tend,hh = tend-tt;end
        k1 = dydt(tt,y(i,:),varargin{:})'
        ymid = y(i,:) + k1.*hh./2;
        k2 = dydt(tt+hh/2,ymid,varargin{:})';
        ymid = y(i,:) + k2*hh/2;
        k3 = dydt(tt+hh/2,ymid,varargin{:})';
        yend = y(i,:) + k3*hh;
        k4 = dydt(tt+hh,yend,varargin{:})';
        phi = (k1+2*(k2+k3)+k4)/6;
        y(i+1,:) = y(i,:) + phi*hh;
        tt = tt+hh;
        i=i+1;
        if tt>=tend,break,end
    end
    np = np+1; tp(np) = tt; yp(np,:) = y(i,:);
    if tt>=tf,break,end
end
end

Answer 1

[Walter Roberson]

Are you trying to run the ode on a whole series of elevations simultaneously, or does one of your state variables represent the evalation?

gravity = 9.81 * ((Re.^2)/((z)+Re).^2)

If you are trying to run a whole series of elevations simultaneously, then here z would need to be the list of elevations. Your value for Re appears to be Earth radius in metres, so your z would have to be a list of elevations in metres that you wanted to process simultaneously, and you would need to change the / to ./

But if one of your state variables is elevation, then you would need gravity to be a function of an input evalation, not a constant based upon a list of pre-determined elevation. You would need to use an anonymous function, and in the ode you would need to pass the elevation state variable to the anonymous function.

denEstimated =  -0.000000000000037549844448166229302696435310363*z.^3 + 0.0000000035489894922510675713650514186874*z.^2  - 0.00011259491846037372752859645474999*z + 1.2180219634561690877916362296673

z is not defined. Again, are you wanting to process a list of pre-defined elevations, or are you wanting to process according to a state variable?

dydt = @(t, y) [y(2); gravity-((CDFree/m)*((1/2)*denEstimated(1:401) * (y(2).^2)*AFree))]

Are you indexing a vector denEstimated at locations 1:401 ? Or are you invoking a function with input [1, 2, 3, ... 401] ? Either way you are probably expecting that the result would be length 401, which would tend to imply that you are running a number of different elevations simultaneously. If denEstimated is made into a function, then are the values 1 to 401 in the proper units for the function? Is that a density estimate for a height in metres like the gravity needs? If so then does it need metres above the surface of the Earth or does it need metres relative to the centre of the Earth ?

[David Wonus]

I have a maxtrix of elevations of length 401 from 0:40000m both gravity and density of air are dependent on the elevations (z) these elevations are relative to the surface such that 0 = sea level

[David Wonus]

tspan = [0 260];
y0 = [0; 0]; % initial conditions for s and v
h = 0.5; % step size
[tp,yp] = rk4sys(@myODEFree,tspan,y0,h);

z = yp(:,1)/1000; % extract solution for altitude and convert to km
v = yp(:,2)*3.6; % extract solution for velocity and convert to km/h
plot(z,v)
ylabel('Velocity (km/h)')
xlabel("Altitude Descent (km)")

tspan = [260 560];
y0 = [2567; 0]; % initial conditions for s and v
h = 0.5; % step size
[tp,yp] = rk4sys(@myODEPara,tspan,y0,h);

z = yp(:,1)/1000; % extract solution for altitude and convert to km
v = yp(:,2)*3.6; % extract solution for velocity and convert to km/h
plot(z,v)
ylabel('Velocity (km/h)')
xlabel("Altitude Descent (km)")

function dydt = myODEFree(t,y)
      % Constants
    CD = 0.735; % Drag coefficient
    m = 120; % kg
    A = 1; % m^2
    z = y(1);
    v = y(2);
    g = 9.81 * (6371000^2)./((6371000+z).^2);
    density = -0.000000000000037549844448166229302696435310363*z.^3 + 0.0000000035489894922510675713650514186874*z.^2 - 0.00011259491846037372752859645474999*z + 1.2180219634561690877916362296673;
    dydt = zeros(2,1);
    dydt(1) = v;
    dydt(2) = g - ((CD/m)*((1/2)*density*v^2*A));
end
% Define a function that evaluates the ODEs
function dydt = myODEPara(t,y)
      % Constants
    CD = 1.75; % Drag coefficient
    m = 120; % kg
    A = 25; % m^2
    z = y(1);
    v = y(2);
    g = 9.81 * (6371000^2)./((6371000+z).^2);
    density = -0.000000000000037549844448166229302696435310363*z.^3 + 0.0000000035489894922510675713650514186874*z.^2 - 0.00011259491846037372752859645474999*z + 1.2180219634561690877916362296673;
    dydt = zeros(2,1);
    dydt(1) = v;
    dydt(2) = g - ((CD/m)*((1/2)*density*v^2*A));
end

%I am still using the rk4sys code, I no longer have errors but it seems
%like something is wrong as theoretically the velocity should increase
%until reaching terminal velocity until eventually decreasing through the
%atmosphere as air resistance increases with lower altitude