Issue with computing eigenvalues using mathematica

Question

Basically I'm trying to find the eigenvalues for matrix, and it takes about 12 hours. When it finishes, it says it couldn't find all the eigenvectors (actually barely any), and I'm skeptical about the ones it did find. All I can really do is post my code, and I'm hoping that someone might be able to make some suggestions to me. I'm not very experienced with mathematica and maybe the slow run time and the bad results has something to do with me and not mathematica's abilities. Thanks to anyone that replies, I really appreciate it.

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 *)

There is a bit more of my code, but this is the point where it is really slowing down. Something I should definitely mention, is that if I set both m1 and m2 to be zero, I don't really have any issues, but setting m1 to a constant makes everything go to hell.

Insanitary

it is probably worth pointing out that a significant chunk of time is spent building the matrix up (even with memoization as Timo suggested). RSolve gives an explicit form for your recursive definition of v, although fixing the undetermined function (via your initial conditions) may be complicated by branch cuts etc. In any case, if you scale this further, this may be something to look at. – Gaal

Answer 1

Your problem is that the constant mS remains symbolic. This means that Mathematica is trying to analytically solve for the eigenvalues instead of numerically. If your problem allows you to choose a numerical value for mS you should do so.

The other, unrelated, problem you have is that you are using a recursive formula and you want to use, e.g., memoization in the following line

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

The extra v[n, m] = stores the value for a given n and m so you don't have to recurse all the way to v[0,0] every time h[n, m] is called in Table[].

With those two things taken care of my old core 2 duo takes less than a minute to do the eigenvalues.

Esoterica

What I really want is to be able to keep mS as a constant, so that when I get some solutions I can vary the mS value myself (i.e. I'd really like to get the solution in terms of mS). But, that definitely makes sense that it is no longer able to find a numerical solution because of that. I think I can live with specifying a numerical value for m1 from the beginning (I could just vary that there instead of later). Anyway, thanks for the reply, that recursive trick is very great! – Insanitary

If you haven't already done so, take a look at what you get with cutoff=5. Using specific numbers for mS is probably the way to go. – Endostosis

Answer 2

This is a follow-up to Timo's answer. I'd like to show a figure so I place it as an answer instead of a comment.

Given that you want to find the Eigenvalues of a matrix that has 501 x 501 symbolic elements. [BTW you call them constants, but that's a misnomer. Constants are just defined, fixed values with a name. What you describe in your comment at Timo's answer is a symbolic variable.]

It's good to see what a fully symbolic matrix does for Eigenvalue calculations. This is for a 2 x 2 matrix:

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])}   
*)

It takes up Array[f, {2, 2}] // Eigenvalues//ByteCount = 3384 bytes. This explodes pretty fast: a 7x7 solution already takes up 70 MB (it takes several minutes to find this one). In fact, there is a nice relation to be found between matrix size and byte count:

The fitted function is: byte count =E^(2.2403067075863197 + 2.2617380321848457 x matrix size).

As you can see, a 501 x 501 symbolic matrix's Eigenvalues won't be found before the end of the universe.

[BTW what is the possessive form of matrix?]

Anorthic

Nice graph! Another way of looking at the problem is to realize that solving for the eigenvalues of an nxn matrix is equivalent to solving for the root of a n'th order polynomial and it is a known fact that there are no general solutions for all the roots of a polynomial for n=5 and up. Thus MMA probably starts assembling a humongous list of conditions for the values of the variable mS under which some of the eigenvalues can be solved for. – Esoterica

@Esoterica In fact MMA returns Root objects. For the above matrix with size 6 x 6, the first member of this Root already has 32,331 elements. – Anorthic

Thanks for clearing up the difference between a constant and symbolic variable. What I should have said in my comment is that I want an answer in terms of the symbolic variable, but that since you've shown that is quite impossible I am going to have to define it as a constant and vary that constant as I see fit. Thanks! – Insanitary