Why is my implementation of the fourth-order Runge-Kutta method not working?

Question

After searching in this site and on my reference book, I found out I have no idea why my code is not working.

I made a fourth-order Runge-Kutta implementation for the mass-spring system (with amortization), as my professor showed us in class. However, the resulting grafic is quite weird, as you can see.

The code I ended up writing is this:

#! /usr/bin/env python3
#-*- coding: utf-8 -*-

def f(data, t, x1, v1):
    from math import cos

    F = data["F"]
    c = data["c"]
    k = data["k"]
    m = data["m"]
    omega = data["omega"]

    return( [v1, (F*cos(omega*t) - c*v1 - k*x1)/m] )

def run(data = {}):
    xi, vi, ti = [data["u1"]], [data["v1"]], [data["t_ini"]]
    dt = data["dt"]

    while ti[-1] <= data["t_end"]:
        xn = xi[-1]
        vn = vi[-1]
        tn = ti[-1]

        K1 = f(data, t = tn, x1 = xn, v1 = vn)
        K1 = [dt*K1[i] for i in range(len(K1))]

        K2 = f(data, t = tn + 0.5*dt, x1 = xn + 0.5*K1[0], v1 = vn + 0.5*K1[1])
        K2 = [dt*K2[i] for i in range(len(K2))]

        K3 = f(data, t = tn + 0.5*dt, x1 = xn + 0.5*K2[0], v1 = vn + 0.5*K2[1])
        K3 = [dt*K3[i] for i in range(len(K3))]

        K4 = f(data, t = tn + dt, x1 = xn + K3[0], v1 = vn + K3[1])
        K4 = [dt*K4[i] for i in range(len(K4))]

        xn = xn + (K1[0] + 2*K2[0] + 2*K3[0] + K4[0])/6
        vn = xn + (K1[1] + 2*K2[1] + 2*K3[1] + K4[1])/6

        ti.append(tn+dt)
        xi.append(xn)
        vi.append(vn)

    return(ti, xi, vi)

This is imported by the main.py file, which only contains the GUI and plot parts of the program, and the function was deduced in class, so I belive the error is in the Runge-Kutta itself. (Probably it's something stupid that I messed up.)

I tried switching the K's in "xn" and "vn", brute-forcing "F" and "c" values in f(), re-writing everything and writing each element of each K by hand (as K11, K12, K21, and so on), but it only gives an exponential result. Also, switching the return of f() for a numpy array didn't help anything.

I found some questions for the RK4 method, but I couldn't solve this problem, nor understand what's wrong. I have some knoledge of the method, but this is actually my first implementation of it, so please, any help is welcome.

I'm using the Anaconda distribution for python3, if it matters.

Answer 1

Could it be this line?

vn = xn + (K1[1] + 2*K2[1] + 2*K3[1] + K4[1])/6

Should it be:

vn = vn + (K1[1] + 2*K2[1] + 2*K3[1] + K4[1])/6

Answer 2

def rk4(f, xvinit, Tmax, N):
T = np.linspace(0, Tmax, N+1)
xv = np.zeros( (len(T), len(xvinit)) )
xv[0] = xvinit
h = Tmax / N
for i in range(N):
    k1 = f(xv[i])
    k2 = f(xv[i] + h/2.0*k1)
    k3 = f(xv[i] + h/2.0*k2)
    k4 = f(xv[i] + h*k3)
    xv[i+1] = xv[i] + h/6.0 *( k1 + 2*k2 + 2*k3 + k4)
return T, xv

This is what you can use for making a function with the rk4 algorithm