如何使用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是一个阶跃函数：当x <0时b = -k（另一个常数），当x> = 0时b = 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，软线的行为就像弹簧。然后重力也指向下方。为此，修改后的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的临界质量。

显然，在任何情况下，跳线都不会重新进入“自由落体”阶段x>0。