0

私は 2 つのプログラムを持っています。1 つは N 結合 ODE を取り込むことができ、もう 1 つは 2 つの結合 ODE を使用します。2 つの同じ ODE を両方のコードに同じ時間範囲で入力すると、異なる答えが得られます。私は正しい答えを知っているので、N many プログラムが間違っていると推測できます。

2方程式専用のコードは次のとおりです。

# solve the coupled system dy/dt = f(y, t)
def f(y, t):
    """Returns the collections of first-order 
    coupled differential equations"""
    #v11i = y[0]
    #v22i = y[1]
    #v12i = y[2]
    print y[0]

    # the model equations
    f0 = dHs(tRel,vij)[0].subs(v12,y[2])
    f1 = dHs(tRel,vij)[3].subs(v12,y[2])
    f2 = dHs(tRel,vij)[1].expand().subs([(v11,y[0]),(v22,y[1]),(v12,y[2])])

    return [f0, f1, f2]


# Initial conditions for graphing
v110 = 6               
v220 = 6                 
v120 = 4                  
y0 = [v110, v220, v120]       # initial condition vector

sMesh  = np.linspace(0, 1, 10e3)   # time grid

# Solve the DE's
soln = odeint(f, y0, sMesh)

そして、ここにN方程式専用のものがあります：

def f(y, t):
    """Returns the derivative of H_s with initial
    values plugged in"""
    # the model equations 
    print y[0]

    for i in range (0,len(dh)):
        for j in range (0,len(y)):
            dh[i] = dh[i].subs(v[j],y[j])

    dhArray = []  
    for i in range(0,len(dh)):
          dhArray.append(dh[i])

    return dhArray

sMesh  = np.linspace(0, 1, 10e3)   # time grid
dh = dHsFunction(t, V_s).expand()
soln = odeint(f, v0, sMesh)

ここで、dHs(tRel,vij) = dHsFunction(t,V_s) つまり、まったく同じ ODE です。同様に y0 と v0 はまったく同じです。しかし、N many ケースで y[0] を出力すると、次の出力が得られます。

6.0
5.99999765602
5.99999531204
5.97655553477
5.95311575749
5.92967598021
5.69527820744
5.46088043467
5.2264826619
2.88250493418
0.53852720647
-1.80545052124
-25.2452277984
-48.6850050755
-72.1247823527
-306.522555124

次の 2 つの専用ケースとは対照的に:

6.0
5.99999765602
5.99999765602
5.99999531205
5.99999531205
5.98848712729
5.98848712125
5.97702879748
5.97702878476
5.96562028875
5.96562027486
5.91961750442
5.91961733611
5.93039037809
5.93039029335
5.89564277275
5.89564273736
5.86137647436
5.86137638807
5.82758984835

etc.

ここで、2 番目の結果は正しい結果であり、適切なグラフをグラフ化します。

より多くのコードが必要な場合やその他のことがあれば教えてください。ありがとう。

4

1 に答える 1

1

の 2 番目のバージョンfは、グローバル変数の値を変更しますdh。最初の呼び出しで値を代入すると、その後のすべての呼び出しで同じ値が使用されます。

関数内などで使用して回避してくださいdh_tmp = list(dh)。

于 2013-10-21T18:52:19.763 に答える