Thanks to "Lastchance" for your hint. And Thanks to "Kikon" for the detailed explanation.
@Kikon: I have modified your script to see the difference between synchronous and sequential RK4 integration of system 1 &2.
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(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
# runge kutta 4th integration,
'''
x1: system 1 displacement
y1: system1 velocity
acc1: system 1 acceleration
x2: system 2 displacement
y2: system2 velocity
acc2: system 2 acceleration
h: time step
'''
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]
# runge kutta 4th integration - test with sequential integration of system 1 and system 2
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_seq = pd.DataFrame(index=range(1 + numpoints), columns=range(7))
df_seq.iloc[0] = [time, x1, y1, temp1, x2, y2, temp2]
for i in range(numpoints):
x1, temp_x2, y1, temp_y2, acc1, temp_acc2 = runge_kutta_4(system1, system2, x1, y1, x2, y2, h)
temp_x1, x2, temp_y1, y2, temp_acc1, acc2 = runge_kutta_4(system1, system2, x1, y1, x2, y2, h)
time = time + h
df_seq.iloc[i + 1] = [time, x1, y1, acc1, x2, y2, acc2]
# Create plots with pre-defined labels.
fig, ax = plt.subplots()
ax.plot(df.loc[:,0], df.loc[:,1],label='displacement TT_synch')
ax.plot(df_seq.loc[:,0], df.loc[:,1],'--',label='displacement TT_seq')
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]', fontdict=None, labelpad=None, loc='center')
the result shows that there is no difference between synchron and sequential calculation.
synchron vs sequential
Therefore I have got below conclusions:
- the deviation I asked at the beginning ("odeint vs runge-kutta-4th") is due to incorrect differential equations
- the synchronous and sequential calculation for system 1 &2 using RK4 have no impact on calculation results.
What is your opinion?
