Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I want to solve the equation Det[matrix] == 0 which came from the consistency of solutions of a set of 8 equations. This matrix has two constants λ and μ and two variables $P$ and $q$. I want to get the plot of $P$ as a function of $q$, but the result of secular equation is giving a large output, and NSolve is not able to give a solution showing the error.

Solve::tdep: The equations appear to involve the variables to be solved for in an essentially non-algebraic way.

I am pasting the input form for the matrix here where I have already put the values for λ and μ. Please tell how to proceed. Ask for any clarifications. Please edit the title if you can make it more clear.

    {1, 1, 0, 0, -1, -1, -1, -1}, {0, 0, 1, 1, -1, 1, 1, -1}, {P, -P, -2, -2, 
    -Sqrt[1   +  P^2 + Sqrt[1 + 2*P^2]], 
    Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]], -Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]], 
    Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]]}, 
    {-2, -2, P, -P, -Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]], -Sqrt[
       1 + P^2 + Sqrt[1 + 2*P^2]], Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]], Sqrt[1 + P^2 -Sqrt[1 + 2*P^2]]}, 
   {E^((0.5*I)*P), E^((-0.5*I)*P), 0, 
    0, -E^(I*(-0.5*Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]] + 
     q)), -E^(I*(0.5*Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]] + 
     q)), -E^(I*(-0.5*Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]] + 
     q)), -E^(I*(0.5*Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]] + q))}, 
  {0, 0, E^((0.5*I)*P), 
   E^((-0.5*I)*P), -E^(I*(-0.5*Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]] + q)), 
   E^(I*(0.5*Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]] + q)), 
  E^(I*(-0.5*Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]] + 
    q)), -E^(I*(0.5*Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]] + q))}, 
  {E^((0.5*I)*P)*P, -(P/E^((0.5*I)*P)), -2*E^((0.5*I)*P), -2/
  E^((0.5*I)*P), -(E^(I*(-0.5*Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]] + q))*
  Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]]), 
 E^(I*(0.5*Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]] + q))*
 Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]], 
 -(E^(I*(-0.5*Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]] + q))*
  Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]]), 
 E^(I*(0.5*Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]] + q))*
 Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]]}, 
 {-2*E^((0.5*I)*P), -2/E^((0.5*I)*P), 
 P, -(P/E^((0.5*I)*
     P)), -(E^(I*(-0.5*Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]] + q))*
  Sqrt[1 + P^2 + 
    Sqrt[1 + 
     2*P^2]]), -(E^(I*(0.5*Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]] + q))*
 Sqrt[1 + P^2 + Sqrt[1 + 2*P^2]]), 
 E^(I*(-0.5*Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]] + q))*
  Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]], 
 E^(I*(0.5*Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]] + q))*
 Sqrt[1 + P^2 - Sqrt[1 + 2*P^2]]}}      
share|improve this question

closed as off-topic by belisarius, m_goldberg, Yves Klett, Michael E2, Jens Jul 8 '13 at 2:41

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question cannot be answered without additional information. Questions on problems in code must describe the specific problem and include valid code to reproduce it. Any data used for programming examples should be embedded in the question or code to generate the (fake) data must be included." – belisarius, m_goldberg
If this question can be reworded to fit the rules in the help center, please edit the question.

Try CharacteristicPolynomial[matrix, variable] reference.wolfram.com/mathematica/ref/… –  Sumit Jul 7 '13 at 15:06
What is a secular equation? You might try and find a smaller example where you have the same problem. Can you get the same error in a 2 by 2 matrix? I'll bet you can. Do this and provide the code for us, and maybe someone can help. –  bill s Jul 7 '13 at 15:54
@bills Re: secular equation. It's an equation not solvable inside a church –  belisarius Jul 7 '13 at 16:02
This question appears to be off-topic because it is too localized –  Yves Klett Jul 7 '13 at 18:56
add comment

1 Answer

This is more of a long comment. It's certainly not a solution.

The error message is because P appears inside exponential functions and outside them in algebraic functions such as Sqrt. Thus the equation is transcendental (non-algebraic). So you should try FindRoot instead of Solve or NSolve.

Further, since the complex exponential is involved, my guess is that for a given q there are infinitely many complex roots P, so defining P as a function of q needs clarification. Also, having so many roots, it may be difficult to apply FindRoot.

Here's a rough idea of what's going on. I stored the determinant in mydet. This shows the complex roots P (real and imaginary parts from -100 to 100), for q between 0 and 10. It's probably not very illuminating, except to show there are many roots for a given q. Since mydet is such a complicated function it took some time to evaluate. To speed things up, I compiled it (mydcf) and memoized it (myfn).

(* mydet = Det[matrix] *)
mydcf = Compile[{{x, _Real}, {y, _Real}, {q, _Real}},
   Evaluate[mydet /. P -> x + I y]
myfn[x_?NumericQ, y_?NumericQ, z_?NumericQ] := myfn[x, y, z] = mydcf[x, y, z]

foo = 0;
ContourPlot3D[Re[myfn[x, y, q]] == 0,
 {x, -20, 20}, {y, -20, 20}, {q, 0, 10},
 MaxRecursion -> 1, PlotPoints -> 9, ContourStyle -> None, BoundaryStyle -> None, 
 MeshFunctions -> {Function[{x, y, q}, Im[myfn[x, y, q]]]}, 
 Mesh -> {{0}}, MeshStyle -> Directive[Thickness[Medium], Darker@Blue], 
 EvaluationMonitor :> foo++]

Mathematica graphics

I noticed, and it's probably trivial, is that both P == 0 and P == -1 is a solution for all q. (The first two columns of the matrix are identical if I substitute P -> 0 or P -> -1.)

The only other thing I noticed is that there are quite a few common subexpressions in your matrix. Perhaps you can take advantage of that.

Response to comments

The imaginary part of P is given by the horizontal axis of the image on the left. One can by the vertical line through the middle that there appears to be a locus of real roots.

Here is how I use FindRoot, for an initial guess P0 for P and a fixed q == q0:

findP[P0_?NumericQ, q0_?NumericQ] := FindRoot[mydet /. q -> q0, {P, P0}]

There are two problems. One is that the gradient is complex, so it is unlikely that the real solutions will be found, unless it's a "trivial" root 0 or -1.

D[mydet, {{q, P}}] /. {q -> 1., P -> -10.}
(* {-5868.51 + 7246.41 I, -1157.66 - 1397.88 I} *)

findP[-9., -5.]
(* {P -> -2.46743 + 1.31051 I} *)

Another is that it's not clear that real P form a continuous function of a real q (except for the trivial roots 0 and -1). We can adjust the previous plot to focus on the real roots P. The first (on left) shows the "roots" {P, q} "found" by ContourPlot; the shows where the real and imaginary parts of mydet are zero (Red and Blue respectively). I used quotes because they are numerical found by approximate means and may be errors. One can see that the first plot shows the roots as a discontinuous set of points, which is entirely consistent with the dimensions of the problem (two real equations Re[mydet] == 0, Im[mydet] == 0 in two real unknowns P, q). The continuous loci P == 0 and P == -1 arise from singularities; there appear to be no other such loci in the region plotted. In the second plot one can see that some of the roots found correspond to where ContourPlot thinks the red and blue lines cross; these probably are real roots. Others correspond to where the red and blue lines appear to be tangent. Whether they are needs checking, and it seems unlikely they would be tangent at more than one point near each other as shown in the first plot.

rootsPlot = ContourPlot[Re[myfn[x, 0, q]] == 0,
 {x, -20, 20}, {q, -10, 10},
 PlotPoints -> 41, ContourStyle -> None, 
 MeshFunctions -> {Function[{x, q}, Im[myfn[x, 0, q]]]}, 
 Mesh -> {{0}}, MeshStyle -> Directive[Thickness[Medium], Darker@Blue], 
 FrameLabel -> {Re[P], q}]

ContourPlot[Re[myfn[x, 0, q]],
 {x, -20, 20}, {q, -10, 10},
 PlotPoints -> 41, Contours -> {0}, ContourStyle -> Red, 
 MeshFunctions -> {Function[{x, q}, Im[myfn[x, 0, q]]]}, 
 Mesh -> {{0}}, MeshStyle -> Directive[Thickness[Medium], Darker@Blue], 
 ContourShading -> None, FrameLabel -> {Re[P], q}]

Mathematica graphics

(Edit: I realized after posting that the labeling above is misleading. The value of P is real; there is no imaginary part.)

One can control-click on a Mac (right-click on Windows??) and use the "Get Coordinates" feature, to get an initial guess for P, q to feed to FindRoot.

findPQ[P0_?NumericQ, q0_?NumericQ] := FindRoot[{Re@mydet, Im[mydet]}, {P, P0}, {q, q0}]

The lower left root can be found this way:

findPQ[-10.72, -7.986]
(* {P -> -10.7133, q -> -7.96709} *)

You can get the approximate roots found in the first plot with the following. (One needs to examine the structure of the output rootsPlot to see what is there. In this case the mesh points consists of a single Point with a list of indices to the points in the GraphicsComplex as an argument.)

 First @ Cases[rootsPlot, 
   GraphicsComplex[pts_, g_, ___] :> 
    pts[[First @ Cases[g, Point[indices_] :> indices, Infinity]]], 
 {P_, q_} /; Chop[P] == 0 || Chop[P + 1] == 0]
share|improve this answer
Michael, thanks for the helpful answer/comment. Thanks for pointing it out that there are complex solutions to the equation. I now think that I am looking for only real solutions as they are the only relevant solutions to my problem. By saying 'P as a "function" of q', I should have been more careful and said that the aim is to get the possible values of P for each q. Thanks again. –  cleanplay Jul 8 '13 at 13:55
does FindRoot give the error "....is not a list of dimensions {1}" when there are multiple roots. How to solve the problem when there are multiple roots, like in this case. –  cleanplay Jul 8 '13 at 14:29
@user25957 See update. It appears very difficult to get FindRoot to find only real roots for this particular problem. Separating real and imaginary parts turns the problem into a real-number problem, but then you cannot search for P as a function of q (that is, for a specified value of q). –  Michael E2 Jul 8 '13 at 15:55
@user25957 One condition, the real part of myfn is zero (Red), is the contour line; the other condition, the imaginary part is zero (Darker@Blue), is the mesh line. In the first plot only the intersections are drawn, which is the only time both parts are zero, that is, that mydet is the complex number zero (for real P and q). Doesn't that seem right? –  Michael E2 Jul 9 '13 at 13:11
@user25957 It's a combination of passing an equation to ContourPlot, setting ContourStyle -> None to draw no contour line, and setting MeshFunction: it draws the "mesh" only on the equation; hence points instead of lines. In the second plot, a function, not an equation, is passed, and ContourPlot maps the whole domain, drawing mesh lines, even though the contour shading is turned off. –  Michael E2 Jul 9 '13 at 13:54
show 5 more comments

Not the answer you're looking for? Browse other questions tagged or ask your own question.