如何用 Python 求解一个非常系数的微分方程？

Question

假设有一个关于弹力绳长度的方程式（用 x 表示），它取决于物体的质量，例如一个玩家（用 m 表示）。

假设弹力绳的自然长度为30米，即起始位置为x(0)=-30。

弹力绳的长度公式为：

    x''(m) = g + b/m*x(m) -a1/m*x'(m) - a2*|x'(m)|*x'(m)

其中g, a1, a2是常量； b 是阶跃函数：b = -k（另一个常数）当 x<0 和 b = 0 当 x>=0.

import numpy as np
from scipy.integrate import odeint

g = 9.8
a1, a2 = 0.6, 0.8
k = 20
b = [-k, 0]

def dx_dt(x, m):
    if x[0]>=0:
        return [x[1], g+b[0]/m*x[0]-a1/m*x[1]-a2/m*np.abs(x[1])*x[1]]
    else:
        return [x[1], g+b[1]/m*x[0]-a1/m*x[1]-a2/m*np.abs(x[1])*x[1]]

init = [[-30, 0], [-40, 0.0001]]

m = np.linspace(1, 100, 10000)

fig, ax = plt.subplots(1, 2, sharey=True, figsize=(6, 4))

for i in range(len(init)):
    xs = odeint(dx_dt, init[i], m)
    ax[i].plot(m, xs[:, 0], 'r-')
    ax[i].set_xlabel('mass (m)')
    ax[i].set_ylabel('length (x)')
    ax[i].set_xlim(xmin, xmax)
    ax[i].set_ylim(ymin, ymax)
    ax[i].set_title('init={}'.format(init[i]))

正确答案应该是正弦曲线

但是上面代码的结果是

代码有问题吗？

Answer 1

将长度坐标 x 更改为指向上方，没有跳线的绳索在位置 0 处静止，以便 x<0 绳索的行为类似于 spring。然后重力也指向下方。为此修改后的 ODE 函数是

def dx_dt(x, t, m):
    acc = -g-a1/m*x[1]-a2/m*np.abs(x[1])*x[1]
    if x[0] < 0: acc -= k/m*x[0]
    return [x[1], acc]

绘制 3 个不同的质量

for i,ini in enumerate(init):
    for m in masses: 
        xs = odeint(dx_dt, ini, t, args=(m,))
        ax[i].plot(t, xs[:, 0], label="m=%.2f"%m)

然后给图片

例如，如果第一种情况下的地面高度为-80m，则总高度为110m，则跳线的质量必须小于90kg。对于更精确的陈述，使用内插法和外推法或数值求解器来找到 x'(t)=0 的第一次和当时 x(t)=ground 的临界质量。

显然跳线在任何情况下都不会重新进入 "free-fall" 阶段 x>0。

如何用 Python 求解一个非常系数的微分方程？

How to solve a differential equation with a non-constant coefficient using Python ？

python

scipy

ode

differential-equations