用于 Python 中 ODE 集成的 Runge–Kutta 方法，具有附加约束

Question

我有一个关于使用 RK4 求解二阶微分方程的问题，考虑到一阶导数的附加约束。我正在做here 所示的示例，并进行了一些修改

θ′′(t) + b θ′(t) + c sin(θ(t)) = 0

附加约束是：

当θ_i θ_(i+1)θ′(i+1)=0.9θ′_i,

其中 i 是 t 的步数，i+1 是 i 之后的一步。在现实世界中，它说当位移方向反转时，它的速度会降低到 90%。

如果y(t) = (θ(t), ω em>(t))，则方程为 ẏ = f(t,y)，其中 f(t,y) = (y₂(t), -by₂(t) - cos(y₁(t)))。

在这个问题中，如果我想在 ω 或 θ′(t) （它们是同一件事）上添加约束，我应该如何修改代码？ 这是我的代码不起作用。附加条件使 θ' 不连续。以下“自制”解决方案无法正确更新 θ′。我是 Python 新手，这是我的第一篇 StackOverflow 帖子。非常感谢任何指导。

def rungekutta4(f, y0, t, args=()):
    n = len(t)
    y = np.zeros((n, len(y0)))
    y[0] = y0
    for i in range(n - 1):
        h = t[i+1] - t[i]
        if y[i][0]*y[i+1][0]<0:
            k1 = f(y[i], t[i], *args)
            k2 = f(y[i] + k1 * h / 2., t[i] + h / 2., *args)
            k3 = f(y[i] + k2 * h / 2., t[i] + h / 2., *args)
            k4 = f(y[i] + k3 * h, t[i] + h, *args)
            y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)*0.9
        else:
            k1 = f(y[i], t[i], *args)
            k2 = f(y[i] + k1 * h / 2., t[i] + h / 2., *args)
            k3 = f(y[i] + k2 * h / 2., t[i] + h / 2., *args)
            k4 = f(y[i] + k3 * h, t[i] + h, *args)
            y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)        
    return y

Answer 1

除非我完全误解了你，否则你想要的是 f: 中的简单大小写区分，在数学上，你有 f₂(t,y) = -by₂(t)-cos(y₁(t )) 如果 θ_i²−1=y₁(t)² -1y₂-1) 否则。这仍然只是 f 对 y 的依赖， 可以简单地以编程方式实现。

例如，您可以定义：

def f(y):
    θ = y[0]
    ω = y[1]
    return [
        θ,
        -b*ω-cos(θ) if θ**2<1 else 0.9*(ω-1)
    ]

虽然这可能会正常工作，但您可能会因为 f 不连续（或具有连续导数）而遇到问题，特别是如果您想使用更高级的步长积分器控制而不是您自制的。 为避免这种情况，请将if–else 构造替换为（尖锐的）sigmoid。有关这方面的更多详细信息，请参阅this answer of mine。

Answer 2

在当前的公式中，考虑到摆锤每次经过垂直位置时，它的速度会降低 10%，这可以近似地安排为

    for i in range(n - 1):
        h = t[i+1] - t[i]
        y[i+1] = RK4step(f,t[i],y[i],h, args)
        if y[i+1,0]*y[i,0] < 0: y[i+1,1] *= 0.9
    return y

也就是说，首先计算新值，然后应用条件。时间步长应该足够小，以使角度仅改变几度。对于较大的时间步长，您必须拆分包含过零的步长，使用一些求根方法（如割线法）来找到更准确的根时间。