強制高調波発振器の微分方程式は、Mx '' + Lx '+(w^2)x = F(t)として与えられます。ここで、F(t)はソース項である。この問題を解決するために、関数diffで微分方程式を定義するコードを書きました。私はF(t)を与える別の関数 'generate_pulse'を書きました。odeintを使用してODEを解く際にソース項を渡す方法
次に、 'diff'関数を他のパラメータとともに呼び出すことによって微分方程式を解く 'odeint'を使用します。今度は、「diff」関数の中にF = 0(つまり、F(t)項を無視する)を指定すると、エラーメッセージは表示されません。私はF(t)を維持すると
F=0 #No error detected if I put F=0 here. Comment out this line to see the error
が、私はエラーメッセージを得る:「デフ」関数の内部を見ていてください「とValueError:順序で配列要素を設定する」を
この問題を解決するにはどうすればよいですか?
コード:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
import itertools
def diff(y, t,param):
M=param[0]
L=param[1]
w2=param[2]
F=param[3]
F=0 #No error detected if I put F=0 here. Comment out this line to see the error
x,v = y
dydt = [v, (F-L*v - w2*x)/M]
return dydt
def generate_pulse(t,Amp,init_delay,rtime,PW):
L=len(t)
t0=t[0:math.ceil(init_delay*L/100)]
t1=t[len(t0):len(t0)+math.ceil(rtime*L/100)]
t2=t[len(t0)+len(t1):len(t0)+len(t1)+math.ceil(PW*L/100)]
t3=t[len(t0)+len(t1)+len(t2):len(t0)+len(t1)+len(t2)+2*math.ceil(rtime*L/100)]
t4=t[len(t0)+len(t1)+len(t2)+len(t3):len(t0)+len(t1)+len(t2)+len(t3)+math.ceil(PW*L/100)]
t5=t[len(t0)+len(t1)+len(t2)+len(t3)+len(t4):len(t0)+len(t1)+len(t2)+len(t3)+len(t4)+math.ceil(rtime*L/100)]
t6=t[len(t0)+len(t1)+len(t2)+len(t3)+len(t4)+len(t5):]
s0=0*t0
s1=(Amp/(t1[-1]-t1[0]))*(t1-t1[0])
s2=np.full((1,len(t2)),Amp)
s2=list(itertools.chain.from_iterable(s2)) #The 'tuple' is converted into array
s3=-Amp/(t1[-1]-t1[0])*(t3-t3[0])+Amp
s4=np.full((1,len(t4)),-Amp)
s4=list(itertools.chain.from_iterable(s4)) #The 'tuple' is converted into array
s5=(Amp/(t5[-1]-t5[0]))*(t5-t5[0])-Amp
s6=np.full((1,len(t6)),0)
s6=list(itertools.chain.from_iterable(s6)) #The 'tuple' is converted into array
s=[s0,s1,s2,s3,s4,s5,s6]
s=list(itertools.chain.from_iterable(s))
return s
###############################################################################
# Main code from here
t = np.linspace(0, 30, 200)
y0 = [- 10, 0.0]
M=5
L = 1
w2 = 15.0
Amp=5
init_delay=10
rtime=10
PW=10
F=generate_pulse(t,Amp,init_delay,rtime,PW)
Param=([M,L,w2,F],) #Making the 'Param' a tuple. Because args of odeint takes tuple as argument.
sol = odeint(diff, y0, t, args=Param)
plt.plot(t, sol[:, 0], 'b', label='x(t)')
plt.plot(t,F,'g',label='Force(t)')
plt.legend(loc='best')
plt.show()