Numerical implementation of ODE differs largely from analytical solution

Question

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?

Answer 1

The analytic solution you provide is the solution for the force

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

Note the absence of a squared v. Once you fix that, both solutions coincide.