RK4 Python Explanation

Question

我想使用我在网上找到的RK4的实现，但是我很难将自己的头缠绕在网上找到的实现。例如：

def rk4(f, x0, y0, x1, n):
    vx = [0] * (n + 1)
    vy = [0] * (n + 1)
    h = (x1 - x0) / float(n)
    vx[0] = x = x0
    vy[0] = y = y0
    for i in range(1, n + 1):
        k1 = h * f(x, y)
        k2 = h * f(x + 0.5 * h, y + 0.5 * k1)
        k3 = h * f(x + 0.5 * h, y + 0.5 * k2)
        k4 = h * f(x + h, y + k3)
        vx[i] = x = x0 + i * h
        vy[i] = y = y + (k1 + k2 + k2 + k3 + k3 + k4) / 6
return vx, vy

有人可以帮我了解参数到底是什么？如果可能的话，我希望有一个更一般的解释，但是，如果更具体的说明使我更容易解释，那么我将专门用于理想的弹簧系统。

Answer 1

您在这里要求参数：

def rk4(f, x0, y0, x1, n):
    ...
    return vx, vy

• f是ode函数，称为微分方程y'(x)=f(x,y(x))，
为def f(x,y)
• (x0,y0)是初始点和值，
• x1是集成间隔[x0,x1]
的末尾

• n是子间隔或集成步骤的数量

• vx,vx是计算的样本点，vy[k]近似y(vx[k])。

您无法将其用于弹簧系统，因为该代码仅适用于标量v。您需要将其更改为与 numpy合作进行矢量操作。

def rk4(func, x0, y0, x1, n):
    y0 = np.array(y0)
    f = lambda x,y: np.array(func(x,y))
    vx = [0] * (n + 1)
    vy = np.zeros( (n + 1,)+y0.shape)
    h = (x1 - x0) / float(n)
    vx[0] = x = x0
    vy[0] = y = y0[:]
    for i in range(1, n + 1):
        k1 = h * f(x, y)
        k2 = h * f(x + 0.5 * h, y + 0.5 * k1)
        k3 = h * f(x + 0.5 * h, y + 0.5 * k2)
        k4 = h * f(x + h, y + k3)
        vx[i] = x = x0 + i * h
        vy[i] = y = y + (k1 + 2*(k2 + k3) + k4) / 6
    return vx, vy