Mathematica 8 can not solve Eigenvalue problem

Question

I cannot explain this eigenvalue problem. Both matrices are positive definite and symmetric, but some solutions are complex.

But if it is matrices lower order, solutions are real. As I increase the number of order the matrices, I got more and more complex solutions (and in same time matrices are positive definite and symmetric.) What is the problem? I try to increase precision, but nothing.

From here you can copy and paste in Mathematica worksheet to see.

http://pastebin.com/raw.php?i=hM6DBVk5

Regards

Answer 1

The problem here is that elements of the matrix $B$ are floating point numbers. And may be somewhere along the way loss of accuracy happens. If to make $B$ a matrix with rational elements, then the answers will be real. For example,

B = Round[B 10^30]/10^30;
D0 := -x*A + B;
Solve[Det[D0] == 0, x]//N

gives

{{x -> 14847.993867199206}, {x -> 148848.79011442242}, {x -> 
   609570.78149666}, {x -> 1.8969502426287285*^6}, 
   {x -> 8.338192814798254*^6}, {x -> 2.2248467905151808*^8}, {x -> 
   2.3271962509336177*^8}, {x -> 2.4808572853336135*^8}, 
   {x -> 2.746203941496255*^8}, {x -> 3.9842849588286966*^8}, {x -> 
   1978.2943756360853}, {x -> 55725.53921427025}, 
   {x -> 321863.65320953313}, {x -> 1.1383725875194645*^6}, {x -> 
   5.278252037812979*^6}, {x -> 2.1976933210891974*^8}, 
   {x -> 2.2745356114392006*^8}, {x -> 2.4019361153335622*^8}, {x -> 
   2.6138647573316988*^8}, {x -> 3.547388273684647*^8}}

Answer 2

I am not sure if this the appropriate place to be asking this question. In the future, I would try maybe linking to a CDF document or giving a summary of the problem you are solving so that people can weigh in without being Mathematica gurus.

The problem here is that you are doing floating point arithmetic. The results are not accurate because of it. Maybe the matrices have a bad condition number. Mathematica is fortunately capable of working with numbers symbolically which is what you will want to do here.

In particular notice that the matrix b is filled with floating point numbers. These can be converted into symbolically processable numbers with the Rationalize command.

After defining the matrix b evaluate:

B = Rationalize[B, 0]

B will be filled with rational numbers instead of floating point numbers.

The results from Solve are all now real numbers.