2017-12-17 11 views
0

強制高調波発振器の微分方程式は、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() 

答えて

0

あなたが合格Fの値が生成された配列であるので、あなたがエラーを取得します。

numpyまたはscipyの補間関数を使用して、実際の関数を配列から外します。次に、時刻tでその関数を評価します。あるいは、強制的に強制的にスカラーの区分的に定義された関数として実装するt

はまた、あなたがodeinttに与えるサンプリング時間のリストがodeintはODE機能diffを呼び出すする時間を行うには(ほとんど)何も持っていないことに注意してください。独自の固定ステップ法を実装する必要があることを制御したい場合は、おそらくRunge-Kuttaではなく、いくつかのマルチステップ法を使用します。

関連する問題