Simulating a free-fall project with ode45

Hallo everybody, now that I've the basic solution for my ode45-method, I've another three questions.

1. How do I add differntials of higher or lower orders in one of the other differntial-equation, like in my code. 'rho' is depended on the height h, but that is the first order of the ode, so how do I implement, the first order in the second order?

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);
       
       v = hv(2); % velocity, positive downwards
       
       dv = [
-v;                          % dh/dt
             g-1/(2*m)*c_w*A*rho_0 * exp(-g*M*h(1)/R*T)*v^2];  % dv/dt
end

2. Since my project needs to show all forms of motion-equations (acceleration, speed, height), how do I plot the acceleration, with is the derivate of the ode45-method

3. And since it is required to name a timespan, and especelly the end, is there a possibility to let the method end, even when you do not exactly know the time needed, to reach a height of 0, in my specific case?

Answer 1

[Mischa Kim]

Hallo Lukas,

1. The height is the first component of the state vector. I'll just call the state vector X is in the code below. So you can just back out h = X(1).
2. Acceleration...you have the expression in the dX vector, the second component. To compute the acceleration as a function of time simply compute this expression in your main function, the function that is calling Beschleunigung().
3. The best and most elegant way to do this is to use events functions. Run in the command window the built-in demo ballode and check out the code.

function dX = Beschleunigung(~, X) % X is the state vector X = [h; v]
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 = X(1); % position or height
v = X(2); % velocity, positive downwards

dX = [-v;                                             % dh/dt
       g-1/(2*m)*c_w*A*rho_0 * exp(-g*M*h/R*T)*v^2];  % dv/dt; T not yet defined?
end