1
$\begingroup$

I want to compute the solution to an equation with a parameter. It turns out that this solution is a sextic polynomial and mathematica outputs the solutions in a root form. Here's my code

Δ[r] = r^2 + a^2 - 2*M*r;
λ[r] = a + r/a*(r - 2*Δ[r]/(r - M));
η[r] = r^3*(4*a^2*M - r*(r - 3*M)^2)/(a^2*(r - M)^2);

tempsol = r /. NSolve[p == (λ[r]^2/η[r])^-1 /. {a -> 0.99, M -> 1}, r];

Plot[tempsol, {p, -2, 13}]

enter image description here

From this I want to extract the solutions for positive $p$ which are larger than $1$, as I need these for other expressions. These are basically the two continuous parts in orange and blue. However, close to $p=0$, mathematica jumps between the solutions because of the ordering NSolve uses. Is there a way to extract the solutions "continuously"?

I already had a look at this link. However, the solution here just helps for plotting.

Edit: Using the strategy in the first answer does not seem to help, as in my further computation, the region where the solution jumps cannot be plotted.

HTilde = Simplify[(2 a L M r[L, U] + Sqrt[r[L, U] (a^2 + r[L, U] (-2 M + r[L, U])) (U r[L, U]^3 - a^2 (L^2 - U) (2 M + r[L, U]))])/(r[L, U]^4 + a^2 r[L, U] (2 M + r[L, U])) /. r[L, U] -> r /. L -> \[Lambda] H /. U -> H^2 (\[Eta] + \[Lambda]^2) /. H -> Abs[H]]

\[Delta]\[Alpha]1 = -\[Lambda] (1 - HTilde / Abs[H]) /. r -> sol[[1]] /. p -> \[Eta]/\[Lambda]^2 /. {\[Eta] -> \[Beta]^2, \[Lambda] -> -\\[Alpha]} /. {a -> 0.99, M -> 1};

DensityPlot[\[Delta]\[Alpha]1, {\[Alpha], -8, 0}, {\[Beta], 0, 8},PlotLegends -> Automatic, ColorFunction -> ColorData[{"SolarColors", "Reverse"}], PlotPoints -> 10, MeshFunctions -> {#3 &, #3 &}, Mesh -> {Range[-10, 10, 0.5], Range[-10, 10, 0.5]}, MeshStyle -> Opacity[0.3,Black], ClippingStyle -> Automatic]

The plot with one region clearly not being evaluated is given below. enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ Are you OK with piecewise definitions? $\endgroup$ Commented Apr 24 at 20:58
  • $\begingroup$ as long as it gives correct results for small p when plugging the solution into subsequent expressions $\endgroup$ Commented Apr 25 at 8:58

2 Answers 2

2
$\begingroup$
$Version

(* "14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023)" *)

Clear["Global`*"]

Δ[r] = r^2 + a^2 - 2*M*r;
λ[r] = a + r/a*(r - 2*Δ[r]/(r - M));
η[r] = r^3*(4*a^2*M - r*(r - 3*M)^2)/(a^2*(r - M)^2);

tempsol = 
  SolveValues[{p == (λ[r]^2/η[r])^-1 /. {a -> 99/100, 
      M -> 1}, p > 0, r > 1}, r];

sol = {Piecewise[List @@@ tempsol[[{1, 3}]]], 
   Piecewise[List @@@ tempsol[[{2, 4}]]]};

Plot[sol, {p, 0, 13}]

enter image description here

EDIT:

HTilde = Simplify[(2  a  L  M  r[L, U] + 
          Sqrt[r[L, 
             U]  (a^2 + r[L, U]  (-2  M + r[L, U]))  (U  r[L, U]^3 - 
              a^2  (L^2 - U)  (2  M + r[L, U]))])/(r[L, U]^4 + 
          a^2  r[L, U]  (2  M + r[L, U])) /. r[L, U] -> r /. 
      L -> λ  H /. U -> H^2  (η + λ^2) /. H -> Abs[H]];

δα1 = -λ  (1 - HTilde/Abs[H]) /. r -> sol[[1]] /. 
     p -> η/λ^2 /. {η -> β^2, λ -> -α} /. {a -> 99/100, M -> 1};

DensityPlot[δα1, {α, -8, 0}, {β, 0, 8},
 PlotLegends -> Automatic,
 ColorFunction -> ColorData[{"SolarColors", "Reverse"}],
 PlotPoints -> 25,
 MaxRecursion -> 3,
 MeshFunctions -> {#3 &, #3 &},
 Mesh -> {Range[-10, 10, 0.5], Range[-10, 10, 0.5]},
 MeshStyle -> Opacity[0.3, Black],
 ClippingStyle -> Automatic,
 Exclusions -> None]

enter image description here

$\endgroup$
7
  • $\begingroup$ I unfortunately only have mathematica 12, which doesn't have SolveValues as a function. Is there another way to do this? $\endgroup$ Commented Apr 25 at 8:53
  • $\begingroup$ r /. Solve[…] will give the same result. $\endgroup$
    – Bob Hanlon
    Commented Apr 25 at 14:57
  • $\begingroup$ This again works nicely to plot these parts, but using them in further calculations doesn't seem to work. See my edit above. $\endgroup$ Commented Apr 25 at 15:13
  • $\begingroup$ Works for me. See edit. $\endgroup$
    – Bob Hanlon
    Commented Apr 25 at 15:59
  • $\begingroup$ I'm sorry but I copy-pasted your code and I still get the same unevaluated region. Could it be an issue with the version of mathematica I am using? $\endgroup$ Commented Apr 25 at 16:41
0
$\begingroup$

One way to go at this is to determine where these root functions cross. That can be done by finding zeros of the discriminant polynomial. For that we first need to get a polynomial. I use exact values so that Together will behave.

delta[r] = r^2 + a^2 - 2*M*r;
lambda[r] = a + r/a*(r - 2*delta[r]/(r - M));
nu[r] = r^3*(4*a^2*M - r*(r - 3*M)^2)/(a^2*(r - M)^2);
poly = Numerator[
  Together[(lambda[r]^2/nu[r])^-1 - p /. {a -> 99/100, M -> 1}]]

(* Out[201]= -96059601 p - 192119202 p r + 492000399 p r^2 + 
 392040000 r^3 + 392040000 p r^3 - 900000000 r^4 - 1096020000 p r^4 + 
 600000000 r^5 + 600000000 p r^5 - 100000000 r^6 - 100000000 p r^6 *)

Next we compute the solutions and also the zeros of the discriminant of poly with respect to r.

rts = Solve[poly == 0, r];
rtsV = r /. rts;
dpoly = Discriminant[poly, r]

(* Out[197]= \
1393789027415502113641834611781759672320000000000000000000000000000 \
(-3557763000000000000 p^2 + 9626512203100000000 p^3 + 
   15625948562800000000 p^4 - 3472793096363850000 p^5 + 
   4110172748122553347 p^6 + 10024639204186403347 p^7) *)

Determine values of p where the solution is 1.

oneSol = Solve[poly == 0 /. r -> 1, p]
rtsV /. oneSol // N

(* Out[204]= {{p -> -(10000/199)}}

Out[205]= {{-0.441971, -0.385892, 0.962187, 1., 2.43284 - 0.231795 I, 
  2.43284 + 0.231795 I}} *)

So there is only one point where a root function goes from less to greater than 1, and it is the fourth of the root functions.

Find the zeros of the discriminant and then choose evaluation points between them to see which root functions might be greater than 1 in the respective regions (between discriminant zeros, that is).

pCrossings = 
  Sort[Union@
    Select[SolveValues[dpoly == 0, p], Element[#, Reals] &]];
testvals = 
  Join[{pCrossings[[1]] - 1}, 
   Most[pCrossings] + 
    Differences[pCrossings]/2, {pCrossings[[-1]] + 1}];
pCrossings // N
testvals // N

(* Out[210]= {-1., -0.964972, 0., 0.261643}

Out[211]= {-2., -0.982486, -0.482486, 0.130821, 1.26164} *)

Now we evaluate each root function at these test points, and keep the ones in each region that exceed 1 at the test points. We show the root numbers for each region.

rootvals = 
  Table[Select[rtsV, 
    With[{tval = (# /. p -> testval)}, 
      Element[tval, Reals] && tval >= 1] &], {testval, testvals}];
Map[#[[All, -1]] &, rootvals]

(* Out[229]= {{4}, {5, 6}, {3, 4}, {3, 4}, {1, 2}} *)

Recall that the fourth root is only valid for p>=-10000/199. With this knowledge, we determine that these are the desired roots for the regions indicated below. I remark that one could merge the third and fourth intervals since the relevant root numbers are the same.

Partition[Join[p /. oneSol, pCrossings, {Infinity}], 2, 1] // N

(* Out[232]= {{-50.2513, -1.}, {-1., -0.964972}, {-0.964972, 0.}, {0., 
  0.261643}, {0.261643, \[Infinity]}}

Also note that the second interval on the left is quite small. If you want to visualize the root jumping taking place around there, run this.

tempsol = 
  r /. NSolve[p == (lambda[r]^2/nu[r])^-1 /. {a -> 0.99, M -> 1}, r];
Plot[tempsol, {p, -1.05, -.92}, PlotRange -> {-4, 20}]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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