2016-08-07 19 views
1

私は、SymphyでLagrangiansからODEを生成し、NumpyとScipyで数値積分するなど、いくつかのプロセスを自動化しようとしています。 最後に完全なコード。solve()として常微分方程式を生成した結果として、私は次のようなSympy式で辞書を得る:odeint()のsympy式が派生エラーを返します

{Derivative(lambda1(t), t): (y(t) + 1)/(x(t)*y(t)), 
Derivative(z(t), t): x(t), 
Derivative(x(t), t): y(t)*z(t), 
Derivative(y(t), t): -x(t)*z(t) 
} 

これから私はscipyのダウンロードからodeint()で微分方程式のシステムを統合したいと思います。このために、として紹介するために、def Field(Q,t):の中の辞書から式(例えば、lambdify)を抽出する必要があります。私は困難と実行のはここです:

私が最初に

def Equ2(nQ,t,Q,Field): 
x1,y1,z1,lamb1 = nQ 
dQ =[] 
for f in Q: 
    dQ.append(lambdify(Q, Field[f.diff(t)],'numpy')(x1,y1,z1,lamb1)) 
return dQ[0:len(nQ)] 

を試みたが、それは2つの引数を持つフィールドを取り、私はのオプションarga=()でそれを渡そうとしましたように、odeint()に行くことができませんodeint()、私(ロング)エラーを与える:

ValueError     Traceback (most recent call last) 
<ipython-input-20-63f086b8a252> in Equ2(nQ, t, Q, Field) 
     20  dQ =[] 
     21  for f in Q: 
---> 22   dQ.append(lambdify(Q, Field[f.diff(t)],'numpy')(x1,y1,z1,lamb1)) 
     23  return dQ[0:len(Q)-1] 

[...] 
ValueError: 
Can't calculate 1st derivative wrt 14.0430379424125. 

をだから私は基本的に同じことをしようとしたが、ループせずに、

def Equ1(nQ,t): 
x1,y1,z1,lamb1 = nQ 
dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy')(x1,y1,z1,lamb1) 
dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy')(x1,y1,z1,lamb1) 
dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy')(x1,y1,z1,lamb1) 
dlam = lambdify((x,y,z,lam[0]), field[lam[0].diff(t)],'numpy')(x1,y1,z1,lamb1) 
return [dx,dy,dz] 

と(と思う)同じ問題があります。私は単純にしようとした場合

ValueError        Traceback (most recent call last) 
<ipython-input-20-63f086b8a252> in Equ1(nQ, t) 
     9 def Equ1(nQ,t): 
    10  x1,y1,z1,lamb1 = nQ 
---> 11  dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy')(x1,y1,z1,lamb1) 
    12  dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy')(x1,y1,z1,lamb1) 
    13  dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy')(x1,y1,z1,lamb1) 

[...] 
ValueError: 
Can't calculate 1st derivative wrt 17.6326726993661. 

を:

def Equ0(nQ,t): 
x,y,z,lamb = nQ 
dx = y*z 
dy = -x*z 
dz = x 
dlam = (y+1.)/(x*y) 
return [dx,dy,dz] 

統合がうまく動作します。また、私がEquX()という関数を、同様の引数を使って呼び出すと、odeint()の中に入るとうまく動作します。

FULL CODE

from sympy import * 
from sympy.physics.mechanics import dynamicsymbols 
from numpy import linspace, sin, cos 
from scipy.integrate import odeint 

t = Symbol('t') 
x = Function('x')(t) 
y = Function('y')(t) 
z = Function('z')(t) 

lam = dynamicsymbols('lambda1:{0}'.format(5)) 
f = x.diff(t)- y*z 
eq = Matrix([x.diff(t) - lam[0].diff(t)*y*x*z+z, 
     y.diff(t) +x*z, 
     z.diff(t)-x 
     ]) 

field = solve(list(eq)+[f],[x.diff(t),y.diff(t),z.diff(t),lam[0].diff(t)]) 


def Equ0(nQ,t): 
    x,y,z,lamb = nQ 
    dx = y*z 
    dy = -x*z 
    dz = x 
    dlam = (y+1.)/(x*y) 
    return [dx,dy,dz] 

def Equ1(nQ,t): 
    x1,y1,z1,lamb1 = nQ 
    dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy')(x1,y1,z1,lamb1) 
    dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy')(x1,y1,z1,lamb1) 
    dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy')(x1,y1,z1,lamb1) 
    dlam = lambdify((x,y,z,lam[0]), field[lam[0].diff(t)],'numpy')(x1,y1,z1,lamb1) 
    return [dx,dy,dz] 


def Equ2(nQ,t,Q,Field): 
    x1,y1,z1,lamb1 = nQ 
    dQ =[] 
    for f in Q: 
     dQ.append(lambdify(Q, Field[f.diff(t)],'numpy')(x1,y1,z1,lamb1)) 
    return dQ[0:len(Q)-1] 

q = [x,y,z,lam[0]] 
nq = [1,2,3,4] 
time=linspace(0,10,10) 

### This line works just fine: 
print Equ0(nq,t), Equ1(nq,t), Equ2(nq,t,q,field) #They give the same output 

sol0 = odeint(Equ0,nq,time) 
sol1 = odeint(Equ1,nq,time) #Errors here 
sol2 = odeint(Equ2,nq,time,args=(q,field)) #And here 

最後に完全なエラー:

--------------------------------------------------------------------------- 
ValueError        Traceback (most recent call last) 
<ipython-input-20-63f086b8a252> in Equ1(nQ, t) 
    9 def Equ1(nQ,t): 
10  x1,y1,z1,lamb1 = nQ 
---> 11  dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy')(x1,y1,z1,lamb1) 
12  dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy')(x1,y1,z1,lamb1) 
13  dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy')(x1,y1,z1,lamb1) 

/usr/local/lib/python2.7/dist-packages/sympy/core/expr.pyc in diff(self, *symbols, **assumptions) 
2864   new_symbols = list(map(sympify, symbols)) # e.g. x, 2, y, z--------------------------------------------------------------------------- 
ValueError        Traceback (most recent call last) 
<ipython-input-20-63f086b8a252> in Equ1(nQ, t) 
    9 def Equ1(nQ,t): 
10  x1,y1,z1,lamb1 = nQ 
---> 11  dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy')(x1,y1,z1,lamb1) 
12  dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy')(x1,y1,z1,lamb1) 
13  dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy')(x1,y1,z1,lamb1) 

    /usr/local/lib/python2.7/dist-packages/sympy/core/expr.pyc in diff(self, *symbols, **assumptions) 
    2864   new_symbols = list(map(sympify, symbols)) # e.g. x, 2, y, z 
2865   assumptions.setdefault("evaluate", True) 
---> 2866   return Derivative(self, *new_symbols, **assumptions) 
2867 
2868  ########################################################################### 

/usr/local/lib/python2.7/dist-packages/sympy/core/function.pyc in __new__(cls, expr, *variables, **assumptions) 
1068     ordinal = 'st' if last_digit == 1 else 'nd' if last_digit == 2 else 'rd' if last_digit == 3 else 'th' 
1069     raise ValueError(filldedent(''' 
---> 1070     Can\'t calculate %s%s derivative wrt %s.''' % (count, ordinal, v))) 
1071 
1072    if all_zero and not count == 0: 

ValueError: 
Can't calculate 1st derivative wrt 0.0. 

2865   assumptions.setdefault("evaluate", True) 
---> 2866   return Derivative(self, *new_symbols, **assumptions) 
2867 
2868  ########################################################################### 

/usr/local/lib/python2.7/dist-packages/sympy/core/function.pyc in __new__(cls, expr, *variables, **assumptions) 
1068     ordinal = 'st' if last_digit == 1 else 'nd' if last_digit == 2 else 'rd' if last_digit == 3 else 'th' 
1069     raise ValueError(filldedent(''' 
---> 1070     Can\'t calculate %s%s derivative wrt %s.''' % (count, ordinal, v))) 
1071 
1072    if all_zero and not count == 0: 

ValueError: 
Can't calculate 1st derivative wrt 0.0. 

TL; DRは、いくつかの派生エラーは、私が外odeintを再現カントodeint()内の表示されます()カスタムと作成された関数。

答えて

3

tをシンボルとして使用する場合は、関数宣言で浮動小数点数としてtを宣言しないでください。浮動小数点tを別の名前のsまたはttに置き換えてみてください。

関連する問題