3

SciPy の integrate.odeint 関数を使用して、硬い ODE のシステムを統合しています。統合は簡単ではなく、時間がかかるため、対応する jacobian も使用しています。方程式を並べ替えることで、ヤコビアンを帯行列として定義できます。API ドキュメントに従って、muおよびmlパラメータで形状を定義したいと思います。残念ながら、ドキュメントは少しあいまいであるため、jacobian 関数を実装する方法を理解できませんでした。

odeintを呼び出す方法を確認するために、次の (ややばかげた) コードを使用しています。

from scipy.integrate import odeint
lmax = 5

def f1(y, t):
    ydot = np.zeros(lmax)
    for i in range(lmax):
        ydot[i] = y[i]**2-y[i]**3
    return ydot

def fulljac(y, t,):
    J = np.zeros((lmax, lmax))
    J[0,0]=-3*y[0]**2 + 2*y[0]
    J[1,1]=-3*y[1]**2 + 2*y[1]
    J[2,2]=-3*y[2]**2 + 2*y[2]
    J[3,3]=-3*y[3]**2 + 2*y[3]
    J[4,4]=-3*y[4]**2 + 2*y[4]
    return J

## initial conditions and output times
delta = 0.0001;
yini  = np.array([delta]*lmax)
times = np.linspace(0, 2/delta, 100)

y, infodict = odeint(f1, yini, times, Dfun=fulljac, full_output=1)
print("f1: nst: {0}, nfe: {1}, nje: {2}".format(infodict["nst"][-1],
                                                infodict["nfe"][-1],
                                                infodict["nje"][-1]))

完全な NxN ヤコビ行列を使用すると、統合は成功します。対角線のみを使用し、mu=0およびml=0を使用すると、積分も成功します。

バンド行列のユース ケースをテストするために、 mu=1とml= 1を使用して人工的な 3xN バンド ヤコビアンを作成しています。ここで、対角線からの導関数はすべてゼロです。これにより、ソルバーの奇妙な動作が発生します (非対角がゼロでない元の問題で見られるものと同様です)。

def bandjac(y, t):
    J = np.zeros((lmax, 3))
    J[0,1]=-3*y[0]**2 + 2*y[0]
    J[1,1]=-3*y[1]**2 + 2*y[1]
    J[2,1]=-3*y[2]**2 + 2*y[2]
    J[3,1]=-3*y[3]**2 + 2*y[3]
    J[4,1]=-3*y[4]**2 + 2*y[4]
    return J

y, infodict = odeint(f1, yini, times, Dfun=bandjac, full_output=1, mu=1, ml=1)
print("f1: nst: {0}, nfe: {1}, nje: {2}".format(infodict["nst"][-1],
                                                infodict["nfe"][-1],
                                                infodict["nje"][-1]))

SciPyのintegrate.odeintでbanded jacobianオプションを使用する適切な方法は何ですか?

4

2 に答える 2

1

完全を期すために、私は自分の質問に答えています。

Warren Weckesser が指摘したように、Scipy <0.14.0 には、odeint がバンド付きヤコビアンを処理する方法に関するバグがあります。

odeint の現在のドキュメントには次のように記載されています。

Dfun は、行にゼロ以外のバンドが含まれる行列を返す必要があります (最小の対角線から開始)。

これは間違っていると思います。代わりに、最も高い対角線から開始する必要があります。

次のコード スニペットは、Dfun が ( test_integrate.py 単体テストから派生した) jacobianを返す方法を示しています。

def func(y, t, c):
    return c.dot(y)

def jac(y, t, c):
    return c

def bjac_rows(y, t, c):
    return np.array([[   0,    75,       1,  0.2], # mu - upper 
                     [ -50,  -0.1,  -1e-04,  -40], # diagonal
                     [ 0.2,   0.1,    -0.2,   0]]) # lower - ml

c = array([[-50,    75,     0,   0],
            [0.2, -0.1,     1,   0],
            [0,    0.1, -1e-4,   0.2],
            [0,      0,   -0.2, -40]])

y0 = arange(4)

t = np.linspace(0, 50, 6)

# using the full jacobian
sol0, info0 = odeint(func, y0, t, args=(c,), full_output=True,Dfun=jac)
print("f1: nst: {0}, nfe: {1}, nje: {2}".format(info0["nst"][-1],
                                                info0["nfe"][-1],
                                                info0["nje"][-1]))

# using the row based banded jacobian
sol2, info2 = odeint(func, y0, t, args=(c,), full_output=True,Dfun=bjac_rows, ml=1, mu=1)
print("f1: nst: {0}, nfe: {1}, nje: {2}".format(info2["nst"][-1],
                                                info2["nfe"][-1],
                                                info2["nje"][-1]))

注: 転置されたバンド行列は、 col_deriv=Trueでは機能しないようです

于 2014-09-23T15:57:20.483 に答える
0

odeint、muおよびmlパラメータによって期待されるように、ヤコビアン行列をバンド形式に変換する関数を作成しました。入力は であると予想されますnp.array。最高のパフォーマンスを得るには、おそらく結果の行列を に変換してnp.arrayから、ヤコビアン関数にその転置を返させ、 を使用することをお勧めしますcol_der=True。現在の API ドキュメントとは対照的に、最初に最大の対角線から開始する必要があるという BrainGrylls の意見は正しいと思います。

def make_banded_jacobian(matrix):
    '''returns a banded jacobian list (in odeint's format), along with mu and ml parameters'''

    # first find the values of mu and ml
    dims = matrix.shape[0]
    assert dims == matrix.shape[1]
    mu = 0
    ml = 0

    for row in xrange(dims):
        for col in xrange(dims):
            if matrix[row][col] != 0:
                if col > row:
                    dif = col - row
                    mu = max(mu, dif)
                else:
                    dif = row - col
                    ml = max(ml, dif)

    banded = []

    for yoffset in xrange(-mu, ml+1):
        row = []

        for diag in xrange(dims):
            x_index = diag
            y_index = diag + yoffset

            if y_index < 0 or y_index >= dims:
                row.append(0.0)
            else:
                row.append(matrix[y_index][x_index])

        banded.append(row)

    return (banded, mu, ml)
于 2016-09-30T20:37:55.273 に答える