Python 中具有附加约束的 ODE 积分的龙格-库塔方法

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%。

向量 if y(t) = (θ( t), ω(t)), 则方程为ẏ = f(t,y), 其中f( t,y) = (y₂(t ), −by₂(t) − cos(y₁(t ))).

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

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) = −by2(t)− cos(y₁(t)) 如果 θ_i²−1=y₁(t)²−1<0 和 0.9·(y2−1) 否则。这仍然只是 f 对 y, 的依赖，可以简单地通过编程实现。

例如，您可以定义：

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

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

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

即先计算新值再应用条件。时间步长应该足够小，角度只改变几度。对于较大的时间步长，您必须拆分包含零交叉的步长，使用一些求根方法（例如正割法）来找到更准确的根时间。

Python 中具有附加约束的 ODE 积分的龙格-库塔方法

Runge–Kutta methods for ODE integration in Python with additional constraints

python

math

physics

differential-equations

runge-kutta