2011-06-28 32 views
5

基本上我試圖找到矩陣的特徵值,大約需要12個小時。當它結束時,它說它找不到所有的特徵向量(實際上幾乎沒有),我對它找到的那些特徵向量持懷疑態度。我所能做的只是發佈我的代碼,我希望有人能夠向我提出一些建議。我對mathematica並不是很有經驗,也許運行時間很慢,不好的結果與我有關,而不是mathematica的能力。感謝任何回覆的人,我真的很感激。使用mathematica計算特徵值的問題

cutoff = 500; (* set a cutoff for the infinite series *) 
numStates = cutoff + 1; (* set the number of excited states to be printed *) 
If[numStates > 10, numStates = 10]; 

    $RecursionLimit = cutoff + 256; (* Increase the recursion limit to allow for the specified cutoff *) 
(* set the mass of the constituent quarks *) 
m1 := mS; (* just supposed to be a constant *) 
m2 := 0; 

(* construct the hamiltonian *) 
h0[n_,m_] := 4 Min[n,m] * ((-1)^(n+m) * m1^2 + m2^2); 

v[0,m_] := 0; 
v[n_,0] := 0; 
v[n_,1] := (8/n) * ((1 + (-1)^(n + 1))/2); 
v[n_,m_] := v[n - 1, m - 1] * (m/(m - 1)) + (8 m/(n + m - 1))*((1 + (-1)^(n + m))/2); 

h[n_,m_] := h0[n,m] + v[n,m]; 

(* construct the matrix from the hamiltonian *) 
mat = Table[h[n,m], {n, 0, cutoff}, {m, 0, cutoff}] // FullSimplify; 

(* find the eigenvalues and eigenvectors, then reverse the order *) 
PrintTemporary["Finding the eigenvalues"]; 
{vals, vecs} = Eigensystem[N[mat]] // FullSimplify; 

$RecursionLimit = 256; (* Put the recursion limit back to the default *) 

還有一點我的代碼,但這是真正放緩的點。我應該提到的一件事是,如果我將m1和m2設置爲零,我沒有任何問題,但將m1設置爲常數會使所有事情都陷入困境。

+2

它可能是值得指出的是時間的顯著塊都花在構建矩陣了(即使有記憶化作爲蒂莫建議)。'RSolve'爲'v'的遞歸定義提供了一個明確的形式,雖然修改未確定的函數(通過初始條件)可能會因分支切割等原因而變得複雜。在任何情況下,如果您進一步擴展它,這可能是看着。 – acl

回答

9

你的問題是,常數mS仍然是象徵性的。這意味着Mathematica試圖解析求解特徵值而不是數值。如果你的問題允許你爲mS選擇一個數值,你應該這樣做。

的其他不相關的問題,你必須是你使用的是遞推公式,並要使用,例如,記憶化在下面一行

v[n_, m_] := v[n, m] = v[n - 1, m - 1]*(m/(m - 1)) 
        + (8 m/(n + m - 1))*((1 + (-1)^(n + m))/2); 

額外v[n, m] =存儲值對於給定nm因此,您不必每次在Table[]中調用h[n, m]時一路緩存到v[0,0]

有了這兩件事情照顧我的舊核心2二重奏不到一分鐘做特徵值。

+3

我真正想要的是能夠將mS保持爲常量,這樣當我得到一些解決方案時,我可以自己改變mS值(即,我真的想用mS來獲得解決方案)。但是,這肯定是有道理的,因爲它不再能夠找到數字解決方案。我想我可以從一開始就爲m1指定一個數值(我可以改變那裏而不是稍後)。無論如何,感謝回覆,遞歸技巧非常棒! – adhanlon

+0

如果你還沒有這樣做,看看你截斷= 5時得到了什麼。對於mS使用特定的數字可能是最好的選擇。 –

3

這是Timo回答的後續。我想顯示一個數字,所以我把它作爲答案而不是評論。

鑑於您想查找具有501 x 501個符號元素的矩陣的特徵值。 [順便說一句,你把它們稱爲常量,但這是一個誤用。常量是剛剛定義的,具有名稱的固定值。你在Timo的回答中描述的是一個象徵性的變量。]​​

很高興看到完全符號矩陣對特徵值計算的作用。這是針對2×2矩陣:

Array[f, {2, 2}] // Eigenvalues 

(* ==> 
{1/2 (f[1, 1]+f[2, 2]-Sqrt[f[1, 1]^2+4f[1, 2] f[2, 1]-2 f[1, 1] f[2, 2]+f[2, 2]^2]), 
1/2(f[1, 1]+f[2, 2]+Sqrt[f[1, 1]^2+4 f[1, 2] f[2, 1]-2 f[1, 1] f[2, 2]+f[2, 2]^2])} 
*) 

它佔用Array[f, {2, 2}] // Eigenvalues//ByteCount = 3384字節。這個爆炸非常快:7x7解決方案已經佔用了70 MB(需要幾分鐘的時間才能找到這個解決方案)。事實上,有矩陣尺寸和字節計數之間找到一個很好的關係:

enter image description here

擬合函數是:字節計數= E ^(2.2403067075863197 + 2.2617380321848457 X矩陣大小)。

正如您所看到的,在宇宙結束之前將不會發現501 x 501符號矩陣的特徵值。

[BTW什麼是矩陣的所有格形式?]

+1

好圖!查看問題的另一種方法是認識到求解n×n矩陣的特徵值等價於求解n階多項式的根,並且已知的事實是,對於所有的n = 5以上的多項式的根。因此,MMA可能開始爲變量「mS」的值組裝一個龐大的條件列表,在該列表中可以求解某些特徵值。 – Timo

+0

@Timo事實上,MMA返回'Root'對象。對於大小爲6 x 6的上述矩陣,此Root的第一個成員已具有32,331個元素。 –

+0

感謝您清除常量和符號變量之間的差異。我在評論中應該說的是,我想要一個符號變量的答案,但是由於你已經證明這是非常不可能的,所以我將不得不將它定義爲一個常量並按照我看到的那樣改變這個常數適合。謝謝! – adhanlon