movements of the descent vehicle in the atmosphere of Earth

I need help writting a Matlab script file, In the matlab I need to write the equations of motion of the space capsule that are as follows:

and I need the graph of time (x axis) vs speed (y axis).

[James Tursa]

What have you done so far? What specific problems are you having with your code? Are you required to write your own RK4 code and not use ode45?

[Ameer Hamza]

Is it necessary to write your own RK-4 code, or can you use ode45?

[danny crizalid]

is indifferent, but I want to do it for RK4 code

Answer 1

[James Tursa]

You have four differential equations involving four variables: V, theta, H, and L. You have to solve these simultaneously since they depend on each other. You can't just pick the one you want and try to solve it by itself like this:

f = @(t,v) v*sind(theta);

You need an f that takes in a 4-element state vector containing V, theta, H, and L, and produces a 4-element derivative vector containing dV/dt, dtheta/dt, dH/dt, and dL/dt. Then that is used in your RK4 looping code.

[danny crizalid]

I understand but there I have doubts I do not know how to correctly enter all the differential equations in the script.

[James Tursa]

Define a 4-element state vector y as

y(1) = V

y(2) = theta

y(3) = H

y(4) = L

Then your function would look something like this:

function dydt = myderiv(t,y,g,Px,R3,K)
V = y(1);
theta = y(2);
H = y(3);
L = y(4);
p = _____; % calculate atmosphere density here, a function of H
dVdt = _____; % you fill this in
dthetadt = _____; % you fill this in
dHdt = _____; % you fill this in
dLdt = _____; % you fill this in
dydt = [dVdt;dthetadt;dHdt;dLdt];
end

The parts for you to fill in are from your differential equations and your atmosphere model.

Then in you RK4 code you would define the function handle as

f = @(t,y) myderiv(t,y,g,Px,R3,K);

[Walter Roberson]

R=R3+h;

You have not defined h at that point.

[danny crizalid]

i have this error

Error: File: myderiv.m Line: 117 Column: 1

Function definitions are not permitted in this context.

[danny crizalid]

if i already saw my error i changed it was h = H

Error: File: myderiv.m Line: 117 Column: 1

Function definitions are not permitted in this context.

[Walter Roberson]

Which MATLAB release are you using? In particular are you using R2016a or earlier, or are you using R2016b or later?

[danny crizalid]

I'm using R2014a

[Walter Roberson]

in your release you need to start by renaming your script to something different than myderive.m such as myrk4.m

Once you have done that, cut the function myderive() out of the script and put it into myderive.m

That is, in your release scripts cannot have functions in them, and also you currently have a clash between naming your script myderive but also wanting to define a function named myderive: your script and function need to be named differently.

[danny crizalid]

I also thought that that was the problem and I tried it but the problem keeps insisting me, this message comes out on the matlab in the function part.

but I can't understand why, here I have another very separate code using function and I don't have any problem.

%%%Equation of Motion in respect to dv/dt AKA vdot
function vdot = f(v, t)
  %Set Constants 
  M= 90;              % Mass in Kg
  rho = 1.225;        % Air Density in kg/m^3
  gradient = 0;       % Incline kept as zero for simplicity atm
  Rolling = 0.003;    % Rolling Resistance
  g = 9.81;           % Gravitational acceleration
  DragArea = 0.05;    % Drag Area in m^2
  FEngine = 80;       % Engine Force, kept  constant for simplicity atm
  
    %%Equation of motion in respect to dv/dt
    vdot = (1/M)*(FEngine - 0.5*rho*DragArea*v^2) - (g*(Rolling + sin(gradient)));
end

[James Tursa]

You only have 4 states, not 6 states. K and t are not states. K is a constant and t is the independent variable. Also you need to pass in all constants that your derivative function needs, so you need to pass in rho0 and beta also. So your deriv function should look something like this:

function dydt = myderiv(t,y,g,Px,R3,K,rho0,beta) % added rho0 and beta
V = y(1);
theta = y(2);
H = y(3);
L = y(4);
% got rid of y(5) and y(6) which don't exist
p = rho0*exp(-H*beta);                  %atmosphere density a function of H
R = _____; % you need to fill this in
dVdt = -(g.*sind(theta))-(p.*(V.^2)/(2.*Px)); 
dthetadt = (p*V*K)/(2*Px)+(V^2-g*R)/(R*V)*cosd(theta); % Need V here, not v
dHdt = V*sind(theta);
dLdt = R3*V*cosd(theta)./(R);
dydt = [dVdt;dthetadt;dHdt;dLdt];
end

Then your function handle needs to include rho0 and beta also:

f = @(t,y) myderiv(t,y,g,Px,R3,K,rho0,beta);

Your initialization code needs to be for 4-element state vectors, e.g.,

%Set Initial Conditions
y0(1) = V;
y0(2) = theta;
y0(3) = H;
y0(4) = L;
%--------------------------------------------------------------------------
% Allocate state vector matrix and initialize first state
n = numel(t);
y = zeros(4,n);
y(:,1) = y0;

And your RK4 code needs to handle 4-element state vectors as well, so the y(j) scalars you currently have should be y(:,j) 4-element column vectors instead. E.g.,

z = t(2) - t(1); % stepsize
for j=1:n-1
    k1 = f(t(j)      , y(:,j)         );
    k2 = f(t(j) + z/2, y(:,j) + k1*z/2);
    k3 = f(t(j) + z/2, y(:,j) + k2*z/2);
    k4 = f(t(j) + z  , y(:,j) + k3*z  );
    y(:,j+1) = y(:,j) + (z/6)*(k1 + 2*k2 + 2*k3 + k4);
end

[Walter Roberson]

file myrk4.m

%==========================================================================
% format long
clear all
clc
%==========================================================================
%Input data
rho0 = 1.225;               %density at sea level (kg/m3)
V = 7800;                   %speed at re-entry into the atmosphere
H = 100e+3:-5000:0;         %altitude range (meters)
t = 0:45:900;               %time range
Gca = 8060;                 %mass
Cxa = 1.314;                %resistance coefficient
Sp  = 16;                   %capsule midsection area
theta = -2;                 %angle of entry atmosphere
L=1.68;                     %rango de superficie
K=0.5;                      %aerodynamic coefficient 
miu = 3.986e+14;            %gravitational parameter of the earth
beta=1.395e-04;             %atmospheric density gradient
R3= 6.371*1000000;          %earth radius
R=R3+H;
g=miu./R.^2;
Px=Gca/(Cxa*Sp);
%--------------------------------------------------------------------------
%Set Initial Conditions
y(1) = V;
y(2) = theta;
y(3) = H;
y(4) = L;
y(5) = K;
y(6) = t;
%--------------------------------------------------------------------------
%RUNGE KUTTA 4 ORDER
f = @(t,y) myderiv(t,y,g,Px,R3,K);
 for j=1:n-1
     k1 = f(t(j),y(j));
     k2 = f(t(j)+z/2,y(j)+k1*z/2);
     k3 = f(t(j)+z/2,y(j)+k2*z/2);
     k4 = f(t(j)+z,y(j)+k3*z);
     y(j+1)=v(j)+(z/6)*(k1+2*k2+2*k3+k4);
 end
 
plot(t,y ,'b-', 'LineWidth',2);
grid on

file myderiv.m

%--------------------------------------------------------------------------
%Differential equations of motion
% function dydt = myderiv(t,y,g,Px,R3,L,K)
function dydt = myderiv(t,y,g,Px,R3,K)
V = y(1);
theta = y(2);
H = y(3);
L = y(4);
K = y(5);
t = y(6);
p = rho0*exp(-H*beta);                  %atmosphere density a function of H
dVdt = -(g.*sind(theta))-(p.*(V.^2)/(2.*Px)); 
dthetadt = (p*v*K)/(2*Px)+(v^2-g*R)/(R*V)*cosd(theta);
dHdt = V*sind(theta);
dLdt = R3*V*cosd(theta)./(R);
dydt = [dVdt;dthetadt;dHdt;dLdt];
end

This will, however, fail because your H is a vector of values, but you do y(3) = H which tries to store the entire vector into a single location.

[Walter Roberson]

In your release, you cannot have a function line inside a script file.

R2016b is the first release that permits function to be defined inside scripts.

This is why I specifically showed which code belongs in which file.

[danny crizalid]

I already tried saving myderiv.m and myrkt4.m separately, but I have this problem:

I do not understand why?

H = 100e+3:-5000:0;         %altitude range (meters)
z=45;
t = 0:45:900;               %time range

H value 1x21 double

t value 1x21 double

[James Tursa]

Please show us your current code. You should have two separate files. One script file that sets constants and executes the RK4 code, and another derivative function file.

Also, I did not notice the large stepsize z that you have. That simply won't work. You need the t vector to have a much smaller stepsize. E.g.,

z = 0.5;
t = 0:z:900;

And H should not be a vector that you initialize. H should just be a single scalar number that is the initial value. Your RK4 code will be calculating the other H values dynamically as part of the simulation. E.g.,

H = 100000; % initial height (m)

[James Tursa]

You are going to get garbage out if you continue to use a huge stepsize. You need to drop the stepsize as I suggested. E.g.,

z = 0.5;
t = 0:z:900;