0

我有两个程序，一个可以接受 N 个耦合 ODE，一个可以使用 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 很多情况下打印 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.

其中第二个结果是正确的，并绘制了正确的图表。

如果需要更多代码或其他任何内容，请告诉我。谢谢。

4

1 回答 1

1

Your second version for f modifies the value of the global variable dh. On the first call, you substitute in values in it, and these same values are then used in all subsequent calls.

Avoid that by using e.g. dh_tmp = list(dh) inside the function.

于 2013-10-21T18:52:19.763 回答