# Numerical implementation of ODE differs largely from analytical solution

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.

Where did I go wrong here?

• 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.. Oct 18, 2023 at 7:11
• 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/…
– Axel
Oct 18, 2023 at 7:30
• 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. Oct 18, 2023 at 7:37
• 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/… Oct 18, 2023 at 7:57
• @WolfgangBangerth I'm Australian ;-)
– Axel
Oct 19, 2023 at 8:44

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.

• 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.
– Axel
Oct 18, 2023 at 8:16

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.

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