ODE 的数值实现与解析解有很大不同

Question

我正在尝试求解包括空气阻力在内的自由落体的常微分方程。

因此我将 ODE 定义为：

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

我认为这应该正确地代表系统，因为

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

哪里

a=vdot

。

现在我使用显式欧拉方法求解此 ODE，如下所示：

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

乍一看，这似乎工作正常。然而，我将它与我在网上找到的 ODE 的解析解进行了绘制

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

而且它们似乎差异很大，如下所示。

我哪里出错了？

Answer 1

您提供的解析解就是力的解

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

注意没有平方

v

。一旦解决了这个问题，两个解决方案就会一致。