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

Question

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

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

附加约束是：

when θ_iθ_(i+1)<0, then θ′_(i+1)=0.9θ′_i,

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

向量上如果y（t）=（θ（t），ω（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）²−1. &书信电报； 0和0.9·（y₂−1） 否则。这仍然只是f对y的依赖，可以简单地实现编程}

例如，您可以定义：

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

虽然这可能是可行的，但由于f不是连续的（或具有连续的导数），您可能会面临问题，特别是如果您想使用更先进的带步长控制的积分器，而不是自制的积分器。 要避免这种情况，请将if–else构造替换为（尖锐的）乙状结肠。有关这方面的详细信息，请参见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

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