2
$\begingroup$

I am trying to solve the ODE of a free fall including air resistance.

I therefore defined my ODE as:

def f(v, g, k, m):
    return g - k/m * v**2

which in my opinion should represent the system correctly because

m*a = m*g -k*v**2

where a=vdot.

Now I solve this ODE using the explicit Euler method like this:

h = 0.1
t = np.arange(0, 1000 + h, h)
v0 = 0

g = 9.81
k = 0.1
m = 1.

# Explicit Euler Method
v_num = np.zeros(len(t))
v_num[0] = 0

x_num = np.zeros(len(t))
x_num[0] = 100

for i in range(0, len(t) - 1):
    v_num[i + 1] = v_num[i] + h*f(v_num[i], g, k, m)
    x_num[i + 1] = x_num[i] - v_num[i + 1] * h

On first glance this seems to work fine. However I plotted it against the analytical solution of the ODE that I found online

v_ana = m*g*(1.-np.exp(-k/m*t))/k

And they seem to differ largely, as shown below.

Chart of numerical and analytical solution of ODE. Analytical solution reaches much higher terminal velocity

Where did I go wrong here?

$\endgroup$
7
  • $\begingroup$ Welcome to scicomp! If I where you, I'd skip the air resistance for a moment to see if the time-stepping is correct (compare to analytical solution). Then add the air-resistance term. I'm suspicius of the fact that the air-resistance term still contains the mass m.. $\endgroup$
    – MPIchael
    Oct 18, 2023 at 7:11
  • $\begingroup$ Hi @MPIchael thanks for your comment! Do you mean in the differential equation m*a = m*g -k*v**2? But how would it get cancelled out of that equation? The analytical solution I just took from wethestudy.com/mathematics/… $\endgroup$
    – Axel
    Oct 18, 2023 at 7:30
  • 2
    $\begingroup$ In the link with your analytical solution, the differential equation is m*a = m*g -k*v but you use the term -k*v**2 to account for air resistance. $\endgroup$ Oct 18, 2023 at 7:37
  • $\begingroup$ OP's analytical solution doesn't match the ODE. For the quadratic velocity drag force, the analytical solution is derived here: philosophicalmath.wordpress.com/2017/10/21/… $\endgroup$
    – cos_theta
    Oct 18, 2023 at 7:57
  • 2
    $\begingroup$ @WolfgangBangerth I'm Australian ;-) $\endgroup$
    – Axel
    Oct 19, 2023 at 8:44

2 Answers 2

4
$\begingroup$

From the link that you have posted in the comment, I guess your function definition is not correct.

def f(v, g, k, m):
    return g - k/m * v

Hope the above should fix the issue, but also when using explicit time stepping, be wary to choose the time step such that it is stable.

$\endgroup$
1
  • $\begingroup$ Thank you! I didn't pay attention that the motion equation looks different for low Reynolds numbers. So basically I implemented the low Reynolds analytical solution and the high Reynolds numerical solution. $\endgroup$
    – Axel
    Oct 18, 2023 at 8:16
7
$\begingroup$

The equation for the square resistance can be easily solved by remembering the pattern $$(e^y)''=e^y\,(y''+y'^2)$$ and considering $u=\exp(-cx)$. This gives $\dot u=cvu$ (with $v=-\dot x$) and $$\ddot u=cu(\dot v+cv^2).$$ With $c=\frac{k}{m}$ (and $\dot u(0)=cu(0)v(0)=0$) one finally gets $$ \ddot u = cgu, ~~~ u(t)=A\cosh(\sqrt{cg}t),\\ v(t)=\frac{\dot u}{cu}=\sqrt{\frac gc}\tanh(\sqrt{cg} t), \\ x(t)=-\frac1c\ln(u(t))=\frac1c[-\ln(A)-\ln(\cosh(\sqrt{cg}t))] $$ For large $t$ the velocity stabilizes at $\sqrt{\frac gc}=\sqrt{98.1}$ and the height drops linearly accordingly with speed about $10 m/s$.

More precisely, for large $t$ $$ x(t)\approx x(0)+\frac{\ln2}{c}-\sqrt{\frac gc}t, $$ which corrects the $t=0$ intercept of the asymptotic line to account for the initially slower fall.

$\endgroup$
2
  • 2
    $\begingroup$ +1 for such a delightful substitution which is much faster than the partial fraction expansion that is often seen. $\endgroup$ Oct 18, 2023 at 10:04
  • 1
    $\begingroup$ It comes from considering the equation as Riccati-equation for $v$. $\endgroup$ Oct 18, 2023 at 10:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.