2011-06-28 33 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