deviation in solutions (differential equations) using odeint vs runge-kutta-4th

I am modelling a Coupled Spring-Mass-System: Two objects with masses m1 and m2, and are coupled through springs with spring constants k1 and k2, with damping d1 and d2.

Method 1: taking a cookbook-script using ODEINT to solve the differential equations.

# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

def vectorfield(w, t, p):
    """
    Defines the differential equations for the coupled spring-mass system.

    Arguments:
        w :  vector of the state variables:
                  w = [x1,y1,x2,y2]
        t :  time
        p :  vector of the parameters:
                  p = [m1,m2,k1,k2,b1,b2]
    """
    x1, y1, x2, y2 = w
    m1, m2, k1, k2, b1, b2 = p

    f = [y1,
         (-d1 * y1 - k1 * x1  + k2 * (x2 - x1 )+ d2 * (y2-y1)) / m1,
         y2,
         (-d2 * (y2-y1) - k2 * (x2 - x1 )) / m2]

    return f

# Parameter values
tau=1.02 
f_1=0.16 
f_2=tau*f_1

m1 = 2000000 # mass1
m2 = 20000 # mass2

# Spring constants
k1 = m1*pow(2*np.pi*f_1,2)
k2 = m2*pow((2*np.pi*f_2),2)


# damping

d1 = (0.04/2/np.pi)*2*pow(k1*m1,0.5)



d_d1=6000
l_p=9.81/pow(2*np.pi*f_2,2)
b=0.3
d2=d_d1*pow((l_p-b)/l_p,2)

# Initial conditions
# x1 and x2 are the initial displacements; y1 and y2 are the initial velocities
x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0

# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 100.0
numpoints = 5001 

# Create the time samples for the output of the ODE solver.
# I use a large number of points, only because I want to make
# a plot of the solution that looks nice.
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

# Pack up the parameters and initial conditions:
p = [m1, m2, k1, k2, d1, d2]
w0 = [x1, y1, x2, y2]

# Call the ODE solver.
wsol = odeint(vectorfield, w0, t, args=(p,),
              atol=abserr, rtol=relerr)

#    convert zip object to tuple object
temp=tuple(zip(t, wsol[:,0],wsol[:,1],wsol[:,2],wsol[:,3]))
#    convert tulpe object to list then to dataframe
df=pd.DataFrame(list(map(list,temp)))

# Create plots with pre-defined labels.
fig, ax = plt.subplots()
ax.plot(df.loc[:,0], df.loc[:,1],label='displacement object1') 


legend = ax.legend(loc='upper right', shadow=None, fontsize='small')
ax.set_xlabel('time [s]', fontdict=None, labelpad=None, loc='center')
ax.set_ylabel('pos [m] or V [m/s]', fontdict=None, labelpad=None, loc='center')

Method 2: Using Runge-Kutta-4th Order (programmed by myself) to solve the differential equations

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib.pyplot as plt




# Parameter values
tau=1.02 
f_1=0.16 
f_2=tau*f_1

m1 = 2000000 # mass1
m2 = 20000 # mass2

# Spring constants
k1 = m1*pow(2*np.pi*f_1,2)
k2 = m2*pow((2*np.pi*f_2),2)


# damping
d1 = (0.04/2/np.pi)*2*pow(k1*m1,0.5)



d_d1=6000
l_p=9.81/pow(2*np.pi*f_2,2)
b=0.3
d2=d_d1*pow((l_p-b)/l_p,2)

def system1(x1,y1,x2,y2):
    return (-d1 * y1 - k1 * x1 + k2 * (x2 - x1 )+ d2 * (y2-y1)) / m1 

def system2(x1, y1, x2, y2):
    return (-d2 * (y2-y1) - k2 * (x2 - x1)) / m2


def runge_kutta_4_sys1(f, x1, y1, x2, y2, h):
    k1 = f(x1,y1,x2,y2)
    k2 = f(x1,y1+k1*h/2,x2,y2)
    k3 = f(x1,y1+k2*h/2,x2,y2)
    k4 = f(x1,y1+k3*h,x2,y2)
 
    

    temp1= y1
    temp2= y1+k1*h/2
    temp3= y1+k2*h/2
    temp4= y1+k3*h
    
    
    displace_solution = x1 + (temp1 + 2 * temp2 + 2 * temp3 + temp4)*h / 6
    velocity_solution = y1 + (k1 + 2 * k2 + 2 * k3 + k4)*h / 6
    acc_solution = f(x1,y1,x2,y2)
 
    return displace_solution, velocity_solution, acc_solution

def runge_kutta_4_sys2(f, x1, y1, x2, y2, h):
    k1 = f(x1,y1,x2,y2)
    k2 = f(x1,y1,x2,y2+k1*h/2)
    k3 = f(x1,y1,x2,y2+k2*h/2)
    k4 = f(x1,y1,x2,y2+k3*h)
 
    

    temp1= y2
    temp2= y2+k1*h/2
    temp3= y2+k2*h/2
    temp4= y2+k3*h
    
    
    displace_solution = x2 + (temp1 + 2 * temp2 + 2 * temp3 + temp4)*h / 6
    velocity_solution = y2 + (k1 + 2 * k2 + 2 * k3 + k4)*h / 6
    acc_solution = f(x1,y1,x2,y2)
 
    return displace_solution, velocity_solution, acc_solution

x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0


h = 0.02
numpoints = 5000
time = 0

temp=0

df = pd.DataFrame(index=range(1+numpoints),columns=range(7))

df.iloc[0]=[time,x1,y1,temp,x2,y2,temp]

for i in range(numpoints):
    x1, y1, acc1 = runge_kutta_4_sys1(system1, x1, y1, x2, y2, h)
    x2, y2, acc2 = runge_kutta_4_sys2(system2, x1, y1, x2, y2, h)
    time=time+h
    df.iloc[i+1]=[time,x1,y1,acc1,x2,y2,acc2]
    df.iloc[i,3]=acc1
    df.iloc[i,6]=acc2

# Create plots with pre-defined labels.
fig, ax = plt.subplots()
ax.plot(df.loc[:,0], df.loc[:,1],label='displacement object1') 


legend = ax.legend(loc='upper right', shadow=None, fontsize='small')
ax.set_xlabel('time [s]', fontdict=None, labelpad=None, loc='center')
ax.set_ylabel('pos [m] or V [m/s]', fontdict=None, labelpad=None, loc='center')    

There is clearly deviation in the solutions between method 1 and method 2.

I assume the deviation is due to inaccurate Runge-Kutta-integration. Is my assumption correct? How to set or improve the accuracy for the Runge-Kutta-Methode?

I tried searching in internet to found the possible restricts and accuracy using runger-kutta method, without success.

As @lastchance first observed in the comments, your RK4 is not correct as it doesn't apply the method to all the four equations at the same time. I also don't understand the rationale behind the two velocity_solution equations.

Here's RK4 method applied to the system of four equations. For pedagogical reasons I add it explicitly, in all its "glorious" complexity. For slightly more complex systems, one would have to add abstractions that would allow the groups of k factors to be computed in loops.

def system1(x1, y1, x2, y2):
    return (-d1 * y1 - k1 * x1 + k2 * (x2 - x1) + d2 * (y2 - y1)) / m1

def system2(x1, y1, x2, y2):
    return (-d2 * (y2 - y1) - k2 * (x2 - x1)) / m2

def runge_kutta_4(f1, f2, x1, y1, x2, y2, h):
    k1x1 = y1
    k1x2 = y2
    k1y1 = f1(x1, y1, x2, y2)
    k1y2 = f2(x1, y1, x2, y2)

    k2x1 = y1 + h * k1y1 / 2
    k2x2 = y2 + h * k1y2 / 2
    k2y1 = f1(x1 + h * k1x1 / 2, y1 + h * k1y1 / 2, x2 + h * k1x2 / 2, y2 + h * k1y2 / 2)
    k2y2 = f2(x1 + h * k1x1 / 2, y1 + h * k1y1 / 2, x2 + h * k1x2 / 2, y2 + h * k1y2 / 2)

    k3x1 = y1 + h * k2y1 / 2
    k3x2 = y2 + h * k2y2 / 2
    k3y1 = f1(x1 + h * k2x1 / 2, y1 + h * k2y1 / 2, x2 + h * k2x2 / 2, y2 + h * k2y2 / 2)
    k3y2 = f2(x1 + h * k2x1 / 2, y1 + h * k2y1 / 2, x2 + h * k2x2 / 2, y2 + h * k2y2 / 2)

    k4x1 = y1 + h * k3y1
    k4x2 = y2 + h * k3y2
    k4y1 = f1(x1 + h * k3x1, y1 + h * k3y1, x2 + h * k3x2, y2 + h * k3y2)
    k4y2 = f2(x1 + h * k3x1, y1 + h * k3y1, x2 + h * k3x2, y2 + h * k3y2)

    x1_next = x1 + h * (k1x1 + 2 * k2x1 + 2 * k3x1 + k4x1) / 6
    x2_next = x2 + h * (k1x2 + 2 * k2x2 + 2 * k3x2 + k4x2) / 6
    y1_next = y1 + h * (k1y1 + 2 * k2y1 + 2 * k3y1 + k4y1) / 6
    y2_next = y2 + h * (k1y2 + 2 * k2y2 + 2 * k3y2 + k4y2) / 6

    acc1_next = f1(x1_next, y1_next, x2_next, y2_next)
    acc2_next = f2(x1_next, y1_next, x2_next, y2_next)

    return x1_next, x2_next, y1_next, y2_next, acc1_next, acc2_next

x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0

h = 0.02
numpoints = 5000
time = 0

temp1 = system1(x1, y1, x2, y2)
temp2 = system1(x1, y1, x2, y2)

df = pd.DataFrame(index=range(1 + numpoints), columns=range(7))

df.iloc[0] = [time, x1, y1, temp1, x2, y2, temp2]

for i in range(numpoints):
    x1, x2, y1, y2, acc1, acc2 = runge_kutta_4(system1, system2, x1, y1, x2, y2, h)
    time = time + h
    df.iloc[i + 1] = [time, x1, y1, acc1, x2, y2, acc2]
    df.iloc[i, 3] = acc1
    df.iloc[i, 6] = acc2