MATLAB: Simulating a free fall-project with ode45

MATLAB and Simulink Student Suiteode45

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

Best Answer

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
Related Question