1
$\begingroup$

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

  • 4
    Oh my. It would really have been useful to have pasted in stuff in `InputForm[]`; I know *Mathematica* enough to analyze snippets even without the program, but pasting in a whole notebook is sorta confusing in plaintext.2011-08-24

2 Answers 2

3

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}}
  • 0
    For this case is ok, but if I increase the order of matrices, I had again the complex solutions. See here http://pastebin.com/raw.php?i=2xxyD9vV2011-08-24
  • 0
    I suspect really bad conditioning; I'd try out `SchurDecomposition[]` (with the `Pivoting` option enabled if need be) on the original matrix and check for any 2-by-2 blocks...2011-08-24
  • 0
    @derdack Here is again is a loss of accuracy since `//N` uses floating points arithmetics. This will give real answers: `N[Solve[ Det[D0] == 0, x], 20]`.2011-08-24
  • 0
    Congratulations. Thank you very much Mr Andrew!2011-08-24
  • 0
    Can I decrease the time for computing because I used just one core for solving this determinant and when matrices are high order, I had a problem.2011-08-24
  • 0
    @derdack May be using `Rationalize` as suggested by Searke could help here.2011-08-24
  • 0
    I'll repeat, @derdack: if you want even more helpful advice than what has already been given, either paste in your problem in `InputForm[]` or upload the notebook you have to an appropriate filehost (e.g. [this](http://www.sendspace.com/)).2011-08-25
  • 0
    Ok J.M. so, can I use more cores for solving this, N[Solve[ Det[D0] == 0, x], 20] to obtain results faster. Parallelize command, or I must reformulate something for faster computing in matrix calculus. Reg.2011-08-25
  • 0
    Matrices are quickly computed here. It is finding roots of a polynomial with large coefficients that takes time.2011-08-25
  • 0
    Yes Andrew, you are right. Problem is polynomial with large coefficients. But can I accelerate N[Solve[ Det[D0...? It is not numerical calculation because of command Solve? Or improve something?2011-08-25
  • 0
    One can try `NSolve[Det[D0] == 0, x, WorkingPrecision -> 20]`, I don't know if it uses the same algorithm as `N[Solve[...]]`.2011-08-25
  • 0
    No, can not work, but thank you so much!2011-08-25
5

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.

  • 0
    You will probably want to use NSolve instead of Solve. NSolve gives out numeric answers.2011-08-23