RK4 Python Explanation

Question

I'd like to use an implementation of RK4 I found online for something, but I'm having a bit of difficulty wrapping my head around the implementations I have found online. For example:

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

Could someone please help me understand what exactly the parameters are? If possible, I'd like a more general explanation, but, if being more specific makes it easier to explain, I'm going to be using it specifically for an ideal spring system.

Why are you using some undocumented code you found online instead of a documented module like [this](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp) or [this](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#hnw93)? — Wrzlprmft, Dec 08 '17 at 10:09

Lutz Lehmann · Answer 1 · 2017-12-08T18:02:59.113

You are asking for the parameters here:

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

f is the ODE function, declared as def f(x,y) for the differential equation y'(x)=f(x,y(x)),
(x0,y0) is the initial point and value,
x1 is the end of the integration interval [x0,x1]
n is the number of sub-intervals or integration steps
vx,vx are the computed sample points, vy[k] approximates y(vx[k]).

You can not use this for the spring system, as that code only works for a scalar v. You would need to change it to work with numpy for vector operations.

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

RK4 Python Explanation

1 Answers1

Linked