2重振り子の運動エネルギーの時間平均を積分して求めたい
解決したいこと
2重振り子の各質点の運動エネルギー $K_{i}(t)\equiv mv_{i}(t)^{2}/2$ の時間平均
$\overline{K_{i}(t)}\equiv(1/t)\int_{0}^{t}K_{i}(t’)dt’$ を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重振り子 2.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])
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)
ms=[m1]
Ls=[L1]
thetas=[theta1]
xs=[x1]
vs=[v1]
#初期条件のコピー
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)
xplot1_ = np.insert(x_ini_cor1,0,0)
yplot1_ = np.insert(y_ini_cor1,0,0)
init = 1
end = 100
dt = 0.01
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
for t in range(len(tpoints)):
# for t in range(2):
x1, v1, xpoints1, vpoints1= RK1(t,x1,v1)
xpoints1 = np.array(xpoints1)
vpoints1 = np.array(vpoints1)
vtips1=[]
for v in vpoints1:
vtips1.append(v[0])
vtips2=[]
for v in vpoints1:
vtips2.append(v[1])
xtips1=[]
for x1 in xpoints1:
xtips1.append(x1[0])
xtips2=[]
for x1 in xpoints1:
xtips2.append(x1[1])
vtips1=np.array(vtips1)
vtips2=np.array(vtips2)
xtips1=np.array(xtips1)
xtips2=np.array(xtips2)
#vtips3=np.array(vtips3)
K1=0.5*(vtips1)**2
U1=-g*np.cos(xtips1)
K2=0.5*(vtips1)**2+0.5*(vtips2)**2+vtips1*vtips2*np.cos(xtips1-xtips2)
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))
自分で試したこと
どんな方法を用いても、以下のエラーが発生してしまいました。
object of type 'int' has no len()
0
