常微分方程的数值实现与解析解有很大不同-解网

问：

我正在尝试解决自由落体的 ODE，包括空气阻力。

因此，我将我的 ODE 定义为：

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

在我看来，这应该正确地代表系统，因为

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

哪里。a=vdot

现在我使用显式欧拉方法求解这个常微分方程，如下所示：

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

它们似乎有很大不同，如下所示。

我哪里做错了？

python 数值方法微分方程

请不要或标记交叉发布。在这里 scicomp.stackexchange.com/questions/43370/....具有讽刺意味的是，我认为这个问题最不偏离主题的地方是 math.SE，因为它的数字内容相当基本，并且没有犯或观察到编码错误。因此，问题实际上在于参考或分析解决方案。

1赞 lastchance 10/18/2023

@Axel，您的解析解应为 v=（beta/k）tanh（beta.t）和 x=（1/k）ln（cosh（beta.t）），其中 beta=sqrt（m.g.k）。显然，有几种不同的写作方式。您给出的解析解是电阻力表达式中 k.v 而不是 k.v**2 的不同微分方程。

答：

1赞 JustLearning 10/21/2023 #1

您提供的解析解是力的解

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

请注意，没有平方。一旦你解决了这个问题，两种解决方案就会重合。v

上一个：追加以列出 JAX 的替代工作流

下一个：如何获得某些参数的最佳值，这些参数将为我提供最高速度

常微分方程的数值实现与解析解有很大不同

Numerical implementation of ODE differs largely from analytical solution

评论