6

scipy の odeint を使用して 2 次 ODE を解こうとしています。私が抱えている問題は、単純化されたスニペットに見られるように、関数が暗黙的に二次項に結合されていることです (例のふりをした物理学を無視してください): 

import numpy as np
from scipy.integrate import odeint

def integral(y,t,F_l,mass):
    dydt = np.zeros_like(y)
    x, v = y
    F_r =  (((1-a)/3)**2 + (2*(1+a)/3)**2) * v # 'a' implicit 
    a  = (F_l - F_r)/mass

    dydt = [v, a]
return dydt


y0 = [0,5]
time = np.linspace(0.,10.,21)
F_lon = 100.
mass = 1000.

dydt = odeint(integral, y0, time, args=(F_lon,mass))

この場合、暗黙の変数を代数的に解くことが可能であることに気付きましたが、実際のシナリオでは、F_rと の評価aと代数操作の間に多くのロジックがあり、失敗します。

DAE は MATLAB のode15i関数を使用して解決できると思いますが、可能であればそのシナリオを回避しようとしています。

私の質問は - Python で暗黙の ODE 関数 (DAE) を解決する方法はありますか (できれば scipy)? そして、そうするために上記の問題を提起するより良い方法はありますか?

a最後の手段として、前の時間ステップから渡すことが許容される場合があります。dydt[1]各時間ステップの後に関数に戻すにはどうすればよいですか?

4

2 に答える 2

3

代数操作が失敗した場合は、たとえばfsolve各タイムステップで実行して、制約の数値解を求めることができます。

import sys
from numpy import linspace
from scipy.integrate import odeint
from scipy.optimize import fsolve

y0 = [0, 5]
time = linspace(0., 10., 1000)
F_lon = 10.
mass = 1000.

def F_r(a, v):
    return (((1 - a) / 3) ** 2 + (2 * (1 + a) / 3) ** 2) * v

def constraint(a, v):
    return (F_lon - F_r(a, v)) / mass - a

def integral(y, _):
    v = y[1]
    a, _, ier, mesg = fsolve(constraint, 0, args=[v, ], full_output=True)
    if ier != 1:
        print "I coudn't solve the algebraic constraint, error:\n\n", mesg
        sys.stdout.flush()
    return [v, a]

dydt = odeint(integral, y0, time)

明らかに、これは時間積分を遅くします。fsolve適切な解決策が見つかることを常に確認し、出力をフラッシュして、発生したときにそれを認識してシミュレーションを停止できるようにします。

前のタイムステップで変数の値を「キャッシュ」する方法については、デフォルト引数が関数定義でのみ計算されるという事実を利用できます。

from numpy import linspace
from scipy.integrate import odeint

#you can choose a better guess using fsolve instead of 0
def integral(y, _, F_l, M, cache=[0]):
    v, preva = y[1], cache[0]
    #use value for 'a' from the previous timestep
    F_r = (((1 - preva) / 3) ** 2 + (2 * (1 + preva) / 3) ** 2) * v 
    #calculate the new value
    a = (F_l - F_r) / M
    cache[0] = a
    return [v, a]

y0 = [0, 5]
time = linspace(0., 10., 1000)
F_lon = 100.
mass = 1000.

dydt = odeint(integral, y0, time, args=(F_lon, mass))

トリックが機能するためには、cacheパラメーターが変更可能でなければならないことに注意してください。それが、リストを使用する理由です。デフォルト引数の仕組みに慣れていない場合は、このリンクを参照してください。

2 つのコードは同じ結果を生成しないことに注意してください。数値の安定性と精度の両方について、前のタイムステップでの値の使用には十分注意する必要があります。ただし、2番目は明らかにはるかに高速です。

于 2014-05-10T11:08:10.760 に答える