多重振り子　正準運動エネルギーをpythonで計算したい

import numpy as np  
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import matplotlib.animation as animation
import pandas as pd

df_param = pd.read_excel('python N重振り子 7.xlsx',header=1,sheet_name=None, engine='openpyxl')

df_param['Sheet1'].index = ['質量m(kg)', '長さL(m)','初速度0(m/s)']
df_param['Sheet1'] = df_param['Sheet1'].drop('個数',axis=1)
s=45
m1 = np.array([df_param['Sheet1'].iat[0, i] for i in range(len(df_param['Sheet1'].columns))]).astype(np.float64).round(2)
L1 = np.array([df_param['Sheet1'].iat[1, i] for i in range(len(df_param['Sheet1'].columns))]).astype(np.float64).round(2)
theta1 = np.array([s,s,s,s,s,s,s])
x1 = np.array([np.radians(i) for i in theta1]).astype(np.float64)
v1 = np.array([df_param['Sheet1'].iat[2, i] for i in range(len(df_param['Sheet1'].columns))]).astype(np.float64).round(2)

df_param['Sheet2'].index = ['質量m(kg)', '長さL(m)','初速度0(m/s)']
df_param['Sheet2'] = df_param['Sheet2'].drop('個数',axis=1)


m2 = np.array([df_param['Sheet2'].iat[0, i] for i in range(len(df_param['Sheet2'].columns))]).astype(np.float64).round(2)
L2 = np.array([df_param['Sheet2'].iat[1, i] for i in range(len(df_param['Sheet2'].columns))]).astype(np.float64).round(2)
theta2 = np.array([s,s,s,s,s,s,s+10**-10])
x2 = np.array([np.radians(i) for i in theta2]).astype(np.float64)
v2 = np.array([df_param['Sheet2'].iat[2, i] for i in range(len(df_param['Sheet2'].columns))]).astype(np.float64).round(2)

ms=[m1,m2]
Ls=[L1,L2]
thetas=[theta1,theta2]
xs=[x1,x2]
vs=[v1,v2]

#初期条件のコピー
xs0 = xs.copy()

print('m=',ms,type(ms),len(ms))
print('L=',Ls,type(Ls),len(Ls))
print('theta=',thetas,type(thetas),len(thetas))
print('x=',xs,type(xs),len(xs))
print('v=',vs,type(vs),len(vs))
print('列数=',len(df_param['Sheet1'].columns))

#初期位置の確認
n = len(df_param['Sheet1'].columns)
x_ini_cor1 = np.zeros(n,dtype=np.float64)
y_ini_cor1 = np.zeros(n,dtype=np.float64)

def ini_cor_func1(j):
    if j == 0:
        x_ini_cor1[j] = L1[j]*np.sin(x1[j])
        y_ini_cor1[j] = -L1[j]*np.cos(x1[j])
    else:
        x_ini_cor1[j] = L1[j]*np.sin(x1[j]) + x_ini_cor1[j-1]
        y_ini_cor1[j] = -L1[j]*np.cos(x1[j])+ y_ini_cor1[j-1]

    return x_ini_cor1[j], y_ini_cor1[j]

for j in range(n):
    x_ini_cor1[j] , y_ini_cor1[j] = ini_cor_func1(j)

x_ini_cor2 = np.zeros(n,dtype=np.float64)
y_ini_cor2 = np.zeros(n,dtype=np.float64)

def ini_cor_func2(j):
    if j == 0:
        x_ini_cor2[j] = L2[j]*np.sin(x2[j])
        y_ini_cor2[j] = -L2[j]*np.cos(x2[j])
    else:
        x_ini_cor2[j] = L2[j]*np.sin(x2[j]) + x_ini_cor2[j-1]
        y_ini_cor2[j] = -L2[j]*np.cos(x2[j])+ y_ini_cor2[j-1]

    return x_ini_cor2[j], y_ini_cor2[j]

for j in range(n):
    x_ini_cor2[j] , y_ini_cor2[j] = ini_cor_func2(j)


xplot1_ = np.insert(x_ini_cor1,0,0)
yplot1_ = np.insert(y_ini_cor1,0,0)
xplot2_ = np.insert(x_ini_cor2,0,0)
yplot2_ = np.insert(y_ini_cor2,0,0)


init = 1
end = 1000
dt = 0.001
h = dt
loop = int(end/h)

n = len(df_param['Sheet2'].columns)
g = 9.8  

# initial state
t = init

tpoints = np.arange(init, end , h)
xpoints1 = []
vpoints1 = []

# A = np.zeros((n,n),dtype=np.float64)
# B = np.zeros((n,n),dtype=np.float64)

E1 = -np.ones_like(x1)


def N_func1(t, x1, v1):
    A1 = np.zeros((n,n),dtype=np.float64)
    B1= np.zeros((n,n),dtype=np.float64)

    for i in range(n):
        for j in range(n):
            for k in range(max(i,j),n):
                A1[i][j] += m1[k]
                B1[i][j] += m1[k]
            if i == j:
                A1[i][j] *= L1[j]
                B1[i][j] *= g * np.sin(x1[i])
            else:
                A1[i][j] *= L1[j]*np.cos(x1[i]-x1[j])
                B1[i][j] *= L1[j]*v1[j]**2*np.sin(x1[i]-x1[j])


    #逆行列の計算
    inv_A1 = np.linalg.inv(A1)

    #inv_A*Bを計算
    inv_A1_B1 = np.dot(inv_A1, B1)

    F1 =  np.dot(inv_A1_B1, E1)

    return F1

xpoints1 = []
vpoints1 = []

#配列要素数の定義
j11 = np.zeros_like(v1)
k11 = np.zeros_like(x1)

j12 = np.zeros_like(v1)
k12 = np.zeros_like(x1) 

j13 = np.zeros_like(v1)
k13 = np.zeros_like(x1)

j14 = np.zeros_like(v1)
k14 = np.zeros_like(x1)


def RK1(t,x1,v1):

    vt1 = v1.copy()
    xt1 = x1.copy()
    xpoints1.append(xt1)
    vpoints1.append(vt1)


    j11 = N_func1(t, x1, v1) * h
    k11 = v1 * h

    j12 = N_func1(t + h / 2, x1 + k11 / 2, v1 + j11 / 2) * h
    k12 = (v1 + j11/ 2)* h

    j13 = N_func1(t + h / 2, x1 + k12 / 2, v1 + j12 / 2) * h
    k13 = (v1 + j12/ 2)* h

    j14 = N_func1(t + h, x1 + k13, v1 + j13)*h
    k14 = (v1 + j13)* h

    v1 += (j11 + 2*j12 + 2*j13 + j14)/6
    x1 += (k11 + 2*k12 + 2*k13 + k14)/6

    return x1,v1,xpoints1 ,vpoints1

E2 = -np.ones_like(x2)

def N_func2(t, x2, v2):
    A2 = np.zeros((n,n),dtype=np.float64)
    B2 = np.zeros((n,n),dtype=np.float64)

    for i in range(n):
        for j in range(n):
            for k in range(max(i,j),n):
                A2[i][j] += m2[k]
                B2[i][j] += m2[k]
            if i == j:
                A2[i][j] *= L2[j]
                B2[i][j] *= g * np.sin(x2[i])
            else:
                A2[i][j] *= L2[j]*np.cos(x2[i]-x2[j])
                B2[i][j] *= L2[j]*v2[j]**2*np.sin(x2[i]-x2[j])


    #逆行列の計算
    inv_A2 = np.linalg.inv(A2)

    #inv_A*Bを計算
    inv_A2_B2 = np.dot(inv_A2, B2)

    F2 =  np.dot(inv_A2_B2, E2)

    return F2

xpoints2 = []
vpoints2 = []

#配列要素数の定義
j21 = np.zeros_like(v2)
k21 = np.zeros_like(x2)

j22 = np.zeros_like(v2)
k22 = np.zeros_like(x2) 

j23 = np.zeros_like(v2)
k23 = np.zeros_like(x2)

j24 = np.zeros_like(v2)
k24 = np.zeros_like(x2)


def RK2(t,x2,v2):

    vt2 = v2.copy()
    xt2 = x2.copy()
    xpoints2.append(xt2)
    vpoints2.append(vt2)

    j21 = N_func2(t, x2, v2) * h
    k21 = v2 * h

    j22 = N_func2(t + h / 2, x2 + k21 / 2, v2 + j21 / 2) * h
    k22 = (v2 + j21/ 2)* h

    j23 = N_func2(t + h / 2, x2 + k22 / 2, v2 + j22 / 2) * h
    k23 = (v2 + j22/ 2)* h

    j24 = N_func2(t + h, x2 + k23, v2 + j23)*h
    k24 = (v2 + j23)* h

    v2 += (j21 + 2*j22 + 2*j23 + j24)/6
    x2 += (k21 + 2*k22 + 2*k23 + k24)/6

    return x2,v2,xpoints2 ,vpoints2


# from ipykernel import kernelapp as app

for t in range(len(tpoints)):
# for t in range(2):
    x1, v1, xpoints1, vpoints1= RK1(t,x1,v1)


for t in range(len(tpoints)):
# for t in range(2):
    x2, v2, xpoints2, vpoints2= RK2(t,x2,v2)

xpoints1 = np.array(xpoints1)
vpoints1 = np.array(vpoints1)
xpoints2 = np.array(xpoints2)
vpoints2 = np.array(vpoints2)


vtips1=[]
for v in vpoints1:
    vtips1.append(v[0])
vtips2=[]
for v in vpoints1:
    vtips2.append(v[1])
vtips3=[]
for v in vpoints1:
    vtips3.append(v[2])
vtips4=[]
for v in vpoints1:
    vtips4.append(v[3])
vtips5=[]
for v in vpoints1:
    vtips5.append(v[4])
vtips6=[]
for v in vpoints1:
    vtips6.append(v[5])
vtips7=[]
for v in vpoints1:
    vtips7.append(v[6])
xtips1=[]
for x1 in xpoints1:
    xtips1.append(x1[0])
xtips2=[]
for x1 in xpoints1:
    xtips2.append(x1[1])
xtips3=[]
for x1 in xpoints1:
    xtips3.append(x1[2])
xtips4=[]
for x1 in xpoints1:
    xtips4.append(x1[3])
xtips5=[]
for x1 in xpoints1:
    xtips5.append(x1[4])
xtips6=[]
for x1 in xpoints1:
    xtips6.append(x1[5])
xtips7=[]
for x1 in xpoints1:
    xtips7.append(x1[6])



vtips1=np.array(vtips1)
vtips2=np.array(vtips2)
vtips3=np.array(vtips3)
vtips4=np.array(vtips4)
vtips5=np.array(vtips5)
vtips6=np.array(vtips6)
vtips7=np.array(vtips7)
xtips1=np.array(xtips1)
xtips2=np.array(xtips2)
xtips3=np.array(xtips3)
xtips4=np.array(xtips4)
xtips5=np.array(xtips5)
vtips6=np.array(vtips6)
vtips7=np.array(vtips7)



K1=0.5*(vtips1)**2
U1=-g*np.cos(xtips1)
#plt.plot(tpoints,K1)
#plt.plot(tpoints,U1)


K2=0.5*(((vtips1*np.cos(xtips1)+vtips2*np.cos(xtips2))**2)
         +(vtips1*np.sin(xtips1)+vtips2*np.sin(xtips2))**2)
U2=-g*(np.cos(xtips1)+np.cos(xtips2))
##plt.plot(tpoints,K2)
#plt.plot(tpoints,U2)


K3=0.5*(((vtips1*np.cos(xtips1)+vtips2*np.cos(xtips2)+vtips3*np.cos(xtips3))**2)
        +(vtips1*np.sin(xtips1)+vtips2*np.sin(xtips2)+vtips3*np.sin(xtips3))**2)
U3=-g*(np.cos(xtips1)+np.cos(xtips2)+np.cos(xtips3))
#plt.plot(tpoints,K3)
#plt.plot(tpoints,U3)

K4=0.5*(((vtips1*np.cos(xtips1)+vtips2*np.cos(xtips2)+vtips3*np.cos(xtips3)+vtips4*np.cos(xtips4))**2)
        +(vtips1*np.sin(xtips1)+vtips2*np.sin(xtips2)+vtips3*np.sin(xtips3)+vtips4*np.sin(xtips4))**2)
U4=-g*(np.cos(xtips1)+np.cos(xtips2)+np.cos(xtips3)+np.cos(xtips4))
#plt.plot(tpoints,K4)
#plt.plot(tpoints,U4)

K5=0.5*(((vtips1*np.cos(xtips1)+vtips2*np.cos(xtips2)+vtips3*np.cos(xtips3)+vtips4*np.cos(xtips4)+vtips5*np.cos(xtips5))**2)
        +(vtips1*np.sin(xtips1)+vtips2*np.sin(xtips2)+vtips3*np.sin(xtips3)+vtips4*np.sin(xtips4)+vtips5*np.sin(xtips5))**2)
U5=-g*(np.cos(xtips1)+np.cos(xtips2)+np.cos(xtips3)+np.cos(xtips4)+np.cos(xtips5))
#plt.plot(tpoints,K5)
#plt.plot(tpoints,U5)

K6=0.5*(((vtips1*np.cos(xtips1)+vtips2*np.cos(xtips2)+vtips3*np.cos(xtips3)+vtips4*np.cos(xtips4)+vtips5*np.cos(xtips5)+vtips6*np.cos(xtips6))**2)
        +(vtips1*np.sin(xtips1)+vtips2*np.sin(xtips2)+vtips3*np.sin(xtips3)+vtips4*np.sin(xtips4)+vtips5*np.sin(xtips5)+vtips6*np.sin(xtips6))**2)
U6=-g*(np.cos(xtips1)+np.cos(xtips2)+np.cos(xtips3)+np.cos(xtips4)+np.cos(xtips5)+np.cos(xtips6))
#plt.plot(tpoints,K6)
#plt.plot(tpoints,U6)

K7=0.5*(((vtips1*np.cos(xtips1)+vtips2*np.cos(xtips2)+vtips3*np.cos(xtips3)+vtips4*np.cos(xtips4)+vtips5*np.cos(xtips5)+vtips6*np.cos(xtips6)+vtips7*np.cos(xtips7))**2)
        +(vtips1*np.sin(xtips1)+vtips2*np.sin(xtips2)+vtips3*np.sin(xtips3)+vtips4*np.sin(xtips4)+vtips5*np.sin(xtips5)+vtips6*np.sin(xtips6)+vtips7*np.sin(xtips7))**2)
U7=-g*(np.cos(xtips1)+np.cos(xtips2)+np.cos(xtips3)+np.cos(xtips4)+np.cos(xtips5)+np.cos(xtips6)+np.cos(xtips7))
#plt.plot(tpoints,K7)
#plt.plot(tpoints,U7)

theta_diff1=[]
for x in xpoints2-xpoints1:
    theta_diff1.append(x[-1])

E_all=K1+K2+K3+K4+K5+K6+K7+U1+U2+U3+U4+U5+U6+U7
E =(K1+U1)/E_all
E2=(K2+U2)/E_all
E3=(K3+U3)/E_all
E4=(K4+U4)/E_all
E5=(K5+U5)/E_all
E6=(K6+U6)/E_all
E7=(K7+U7)/E_all
fig, ax = plt.subplots()
ax.set_xlabel("time(s)")
ax.set_ylabel("Kinetic Energy(time averaged)")
ax.set_title("7 pendulum 45 degrees")

print(len(tpoints))
t=np.arange(1,999001,1)
print(len(t))
import itertools
import operator
l=[1]
for i in itertools.accumulate(l):
    plt.plot(tpoints,list(itertools.accumulate(K1))/t,label="1")
    plt.plot(tpoints,list(itertools.accumulate(K2))/t,label="2")
    plt.plot(tpoints,list(itertools.accumulate(K3))/t,label="3")
    plt.plot(tpoints,list(itertools.accumulate(K4))/t,label="4")
    plt.plot(tpoints,list(itertools.accumulate(K5))/t,label="5")
    plt.plot(tpoints,list(itertools.accumulate(K6))/t,label="6")
    plt.plot(tpoints,list(itertools.accumulate(K7))/t,label="7")
plt.legend()

fig, ax = plt.subplots()
ax.set_xlabel("time(s)")
ax.set_ylabel("Energy ratio (time averaged)")
ax.set_title("7 pendulum 45 degrees")

print(len(tpoints))
t=np.arange(1,999001,1)
print(len(t))
import itertools
import operator
l=[1]
for i in itertools.accumulate(l):
    plt.plot(tpoints,list(itertools.accumulate(E))/t,label="1")
    plt.plot(tpoints,list(itertools.accumulate(E2))/t,label="2")
    plt.plot(tpoints,list(itertools.accumulate(E3))/t,label="3")
    plt.plot(tpoints,list(itertools.accumulate(E4))/t,label="4")
    plt.plot(tpoints,list(itertools.accumulate(E5))/t,label="5")
    plt.plot(tpoints,list(itertools.accumulate(E6))/t,label="6")
    plt.plot(tpoints,list(itertools.accumulate(E7))/t,label="7")
plt.legend()