Simulating a free fall-project with ode45

Hallo everyone, I'm trying to simulate a free fall-project, from the stratosphere (about 50 km) to the ground. So, in order to conclude an appropriate solution, drag must be put in consideration. I've allready found a proper solution for the Euler-method, but not for the ode45-method. Furthermore I need to put the [vector-]solutions of the acceleration (the differnecial) for every single calculation-step, in the velocity terms, I mean you could realize this one easy, by putting a while-slope in the code, that the steps would repeat itself until completion. But now to my issue, I've tried many approches with the ode45-method, but non of those have worked out.

Here my code:

function [a] = Test2
        
        start=[0];
        
        tspan=[0 100];
        
        


        [t,A] = ode45(@Beschleunigung, tspan, start);
        
        %a(end+1) = - g + 1/(2*m)* c_w * A * rho(end)* (v(end))^2; 
        plot(t,A(:,1));
end        
        function dv = Beschleunigung(t, v)
        
        g = 9.81;
        m = 120;
        c_w = 0.28;
        A = 2.7;
        rho = 1.2041;
        
        dt = 0.5;
        
        hi = 40000;
        vi = 0;
        
        t = [0]; % Deklaration des Zeit-Arrays
        h = [hi]; % Deklaration des Höhe-Arrays
        v = [vi]; % Deklaration des Geschwindigkeits-Arrays
        a = [0]; % Deklaration des Beschleunigungs-Arrays
        
        dv = zeros(2,1);
        a=dv(1,end);
        dv(1,end+1) = g-(1/2*m)*c_w*A*rho*v^2;
        v(end+1)=v(end)+dv(1,end)*dt;
        h(end+1)=h(end)+v(end)*dt;
        t(end+1)=t(end)+dt;
        
        end

[Lukas Menzel]

And here are the errors of my programm

>> Test2
Error using odearguments (line 93)
BESCHLEUNIGUNG must return a column vector.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options,
  varargin);

Error in Test2 (line 10)
        [t,A] = ode45(@Beschleunigung, tspan, start);

Answer 1

[Alan Stevens]

More like this perhaps:

hi = 40000;     % Initial height
start=[hi 0];   % [Initial height, initial velocity]
tspan=[0 100];
[t,A] = ode45(@Beschleunigung, tspan, start);

%a(end+1) = - g + 1/(2*m)* c_w * A * rho(end)* (v(end))^2; 
subplot(2,1,1)
plot(t,A(:,1)),grid
xlabel('t'),ylabel('height')
subplot(2,1,2)
plot(t,A(:,2)),grid
xlabel('t'),ylabel('velocity')
       
function dv = Beschleunigung(~, hv)

        g = 9.81;
        m = 120;
        c_w = 0.28;
        A = 2.7;
        rho = 1.2041;
       
       v = hv(2); % velocity, positive downwards
       
       dv = [-v;                        % dh/dt
             g-1/(2*m)*c_w*A*rho*v^2];  % dv/dt

end

[Lukas Menzel]

That answers my question thank you, very much:D

[Lukas Menzel]

I have another question: How do I plot the actual differential (the acceleration)?

[Alan Stevens]

Probably the simplest way here is as follows:

g = 9.81;
hi = 40000;     % Initial height
start=[hi 0];   % [Initial height, initial velocity]
tspan=[0 100];
[t,A] = ode45(@Beschleunigung, tspan, start);
accn = [g; diff(A(:,2))./diff(t)];
%a(end+1) = - g + 1/(2*m)* c_w * A * rho(end)* (v(end))^2; 
subplot(3,1,1)
plot(t,A(:,1)),grid
xlabel('t'),ylabel('height')
subplot(3,1,2)
plot(t,A(:,2)),grid
xlabel('t'),ylabel('velocity')
subplot(3,1,3)
plot(t,accn),grid
xlabel('t'),ylabel('acceleration')

       
function dv = Beschleunigung(~, hv)

        g = 9.81;
        m = 120;
        c_w = 0.28;
        A = 2.7;
        rho = 1.2041;
       
       v = hv(2); % velocity, positive downwards
       
       dv = [-v;                        % dh/dt
             g-1/(2*m)*c_w*A*rho*v^2];  % dv/dt

end

[Lukas Menzel]

Thank you, very much for the additional support:D

[Lukas Menzel]

I've got one final question:

dv = [-v;                        % dh/dt
             g-1/(2*m)*c_w*A*rho_0 * exp(-g*M*h/R*T_EO)*v^2];  % dv/dt

Since the drag is dependent of the density and the density is dependent on the height, how can I implement the h, preferably with an additional function? I've tried it, but the programm hasn't computed it correctly.

g = 9.81;
hi = 40000;     % Initial height
start=[hi 0];   % [Initial height, initial velocity]
tspan=[0 100];
[t,A] = ode45(@Beschleunigung, tspan, start);
accn = [g; diff(A(:,2))./diff(t)];
%a(end+1) = - g + 1/(2*m)* c_w * A * rho(end)* (v(end))^2; 
subplot(3,1,1)
plot(t,A(:,1)),grid
xlabel('t'),ylabel('height')
subplot(3,1,2)
plot(t,A(:,2)),grid
xlabel('t'),ylabel('velocity')
subplot(3,1,3)
plot(t,accn),grid
xlabel('t'),ylabel('acceleration')
       
function dv = Beschleunigung(~, hv)
g    = 9.81
R    = 8.31446261;
M    = 0.0289586;
m    = 120;
c_w  = 0.28;
A    = 2.7;
p_0  = 101325;
T_EO = 293.15;

rho_0 = (p_0 * M) / (R * T_EO);
       
       h = hv(1);
       v = hv(2); % velocity, positive downwards
       
       dv = [-v;                        % dh/dt
             g-1/(2*m)*c_w*A*rho_0 * exp(-g*M*h/R*T_EO)*v^2];  % dv/dt
end

[Alan Stevens]

I suspect you should have

dv = [-v;                        % dh/dt
             g-1/(2*m)*c_w*A*rho_0 * exp(-g*M*h/(R*T_EO))*v^2];  % dv/dt

Notice the parentheses around (R*T_E0). The temperature should multiply R, with both in the denominator. The way you had written it, T_E0 was in the numerator.

[Lukas Menzel]

Yeah, thank you

[Lukas Menzel]

Oh and Merry Christmas