2012-02-25 97 views
3

我很難解決光學bloch方程,這是一個具有複數值的一階ODE系統。我發現scipy可能會解決這樣的系統,但是他們的網頁提供的信息太少,我很難理解它。scipy中的複雜ODE系統

我有8再加一階微分方程,我應該產生一個功能,如:

def derv(y): 
    compute the time dervative of elements in y 
    return answers as an array 

然後做complex_ode(derv)

我的問題是:

  1. 我y是不是列表但矩陣,我怎麼能給一個電流輸出 適合complex_ode()?
  2. complex_ode()需要一個雅可比,我不知道如何開始構建一個 和它應該是什麼類型?
  3. 我應該如何將初始條件放在正常的ode和 時間linspace中?

這是SciPy的的complex_ode鏈接: http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.complex_ode.html

誰能給我提供更多資料,這樣我可以瞭解更多一點。

回答

5

我認爲我們至少可以將您指向正確的方向。光學 bloch方程是一個在科學界很好理解的問題,雖然不是由我:-),所以已經有關於這個特定問題的因特網 的解決方案。

http://massey.dur.ac.uk/jdp/code.html

然而,爲了滿足您的需求,你說使用complex_ode,我想 是好的,但我認爲只是普通的scipy.integrate.ode會工作得很好,以及 根據自己的文檔:

from scipy import eye 
from scipy.integrate import ode 

y0, t0 = [1.0j, 2.0], 0 

def f(t, y, arg1): 
    return [1j*arg1*y[0] + y[1], -arg1*y[1]**2] 
def jac(t, y, arg1): 
    return [[1j*arg1, 1], [0, -arg1*2*y[1]]] 
r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True) 
r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0) 
t1 = 10 
dt = 1 
while r.successful() and r.t < t1: 
    r.integrate(r.t+dt) 
    print r.t, r.y 

您也有一箇舊的更成熟,更好的 記錄功能的好處。我很驚訝你有8個,而不是9再加ODE的,但我敢肯定 你明白這一點比我更好的是,你是正確的,你的函數 的形式應該ydot = f(t,y),你叫def derv()的,但你 需要確保你的功能至少需要兩個參數 ,比如derv(t,y)。如果你的y在矩陣中,沒問題!只是「重塑」它 的derv(t,y)功能,像這樣:

Y = numpy.reshape(y,(num_rows,num_cols)); 

只要num_rows*num_cols = 8,你ODE的數量,你應該罰款。然後 在你的計算中使用矩陣。當您完成所有操作,只是一定要返回 像一個載體,而不是一個矩陣:

out = numpy.reshape(Y,(8,1)); 

不是必需的雅可比,但它可能會允許計算來更快速地進行 。如果你不知道如何計算這一點,你可能想諮詢 維基百科或微積分課本。這非常簡單,但可能很耗時。

至於初始條件,你應該已經知道這些應該 是,無論是複雜或實值。只要你選擇的值是在合理範圍內 ,應該沒有多大關係。

+0

非常感謝烏拉圭回合極大的答案,必須有成本ü相當多的時間。你給的鏈接絕對是完美的,我正在尋找。是的,對,應該是9 :)好吧,有時間到今天,如果我在某個地方遇到困難,我可以回到你身邊嗎?再次感謝。 – user1233157 2012-02-26 11:58:13