0
$\begingroup$

I want to get the solution to a 6th degree polynomial in a certain range.

I have the following code

DefiningExpression = Simplify[r^4 (-3 M + r)^2 - a^4 (-1 + x) (M + r)^2 + 2 a^2 r^2 (3 (-1 + x) M^2 - 2 M r - (-1 + x) r^2)]; sol = Solve[DefiningExpression == 0, r]; CloseBound = 2 M (1 + Cos[2/3 ArcCos[-a/M]]); FarBound = 2 M (1 + Cos[2/3 ArcCos[a/M]]); Plot[{Lookup[sol[[1]], r], Lookup[sol[[2]], r], Lookup[sol[[3]], r], Lookup[sol[[4]], r]}, {x, 0, 2}] 

I want to solve the DefiningExpression for $r$ in terms of $x$. Unfortunately this is a no-so-friendly 6th degree polynomial, so we won't find any analytic solutions. The additional requirement is that only the following values of $r$ are of interest: $$r \in [2 M (1 + \cos[2/3 \arccos[-a/M]]), 2 M (1 + \cos[2/3 \arccos[a/M]])].$$

When plotting this, we get that the solutions between these bounds are spread out between the different solutions of sol:

Here $a=0.7$, $M=1$.

Also the constraints on $a$ and $M$ are that $r>M>a>0$.

How can I get a solution of $r$ in terms of $x$ where $r$ lies within this range?

$\endgroup$
7
  • 5
    $\begingroup$ Please post the Mathematica code for your equation so that we don't have to type it in from scratch. Also, are there ranges of $a$ and $M$ that are of interest? $\endgroup$ Commented Nov 22, 2023 at 16:58
  • 1
    $\begingroup$ @JimB: You are completely right. This is not a good practice not to present a Mathematica code. However, a TeX code in the text can be seen by clicking the right button of a mouse and choosing Show Math As/ TeX Commands and can be copied. In the question under consideration the TeX code is very close to Mathematica code. $\endgroup$ Commented Nov 22, 2023 at 17:50
  • 2
    $\begingroup$ I can't believe I'm saying this: Maybe I'm more picky than you. I also like to be able to assume that the OP has a copy of Mathematica. $\endgroup$ Commented Nov 22, 2023 at 18:00
  • $\begingroup$ @JimB: The OP wrote " I'll test this out as soon as possible!" so one may draw the conclusion the OP has access to Mathematica. $\endgroup$ Commented Nov 22, 2023 at 18:22
  • $\begingroup$ @JimB I was in a rush and didn't know immediately how to write down mathematica code in stackexchange so I opted for LaTeX. I'm also new to this site of stack as you may see. I'm updating the post now. $\endgroup$ Commented Nov 23, 2023 at 8:33

2 Answers 2

1
$\begingroup$

Evaluating M and a, this can be done as follows.

M = 3; a = 1; Reduce[r^4 (-3 M + r)^2 - a^4 (-1 + x) (M + r)^2 + 2 a^2 r^2 (3 (-1 + x) M^2 - 2 M r - (-1 + x) r^2) == 0 && r >= 2 M (1 + Cos[2/3*ArcCos[-a/M]]) && r <= 2 M (1 + Cos[2/3*ArcCos[a/M]]), r, Reals] 

(x == 0 && r == Root[3 + # - 9 #^2 + #^3& , 3, 0]) || (0 < x < Root[-4096576 + 3799632 # + 272724 #^2 + 26407 #^3& , 1, 0] && (r == Root[9 - 9 x + (6 - 6 x) #1 + (-53 + 53 x) #1^2 - 12 #1^3 + (83 - 2 x) #1^4 - 18 #1^5 + #1^6 &, 1] || r == Root[ 9 - 9 x + (6 - 6 x) #1 + (-53 + 53 x) #1^2 - 12 #1^3 + (83 - 2 x) #1^4 - 18 #1^5 + #1^6 &, 2])) || (Root[-4096576 + 3799632 # + 272724 #^2 + 26407 #^3& , 1, 0] <= x < 1 && (r == Root[9 - 9 x + (6 - 6 x) #1 + (-53 + 53 x) #1^2 - 12 #1^3 + (83 - 2 x) #1^4 - 18 #1^5 + #1^6 &, 3] || r == Root[ 9 - 9 x + (6 - 6 x) #1 + (-53 + 53 x) #1^2 - 12 #1^3 + (83 - 2 x) #1^4 - 18 #1^5 + #1^6 &, 4])) || (x == 1 && (r == Root[-12 + 81 # - 18 #^2 + #^3& , 2, 0] || r == Root[-12 + 81 # - 18 #^2 + #^3& , 3, 0]))

Addition.For M = Pi; a = 1.; and with Chop the result is

(x == 0 && r == 9.28055) || (0 < x < 0.999592 && (r == Root[2.75795*10^23 - 2.75795*10^23 x + (1.75577*10^23 - 1.75577*10^23 x) #1 + (-1.62683*10^24 + 1.62683*10^24 x) #1^2 - 3.51153*10^23 #1^3 + (2.53804*10^24 - 5.58878*10^22 x) #1^4 - 5.2673*10^23 #1^5 + 2.79439*10^22 #1^6 &, 1] || r == Root[ 2.75795*10^23 - 2.75795*10^23 x + (1.75577*10^23 - 1.75577*10^23 x) #1 + (-1.62683*10^24 + 1.62683*10^24 x) #1^2 - 3.51153*10^23 #1^3 + (2.53804*10^24 - 5.58878*10^22 x) #1^4 - 5.2673*10^23 #1^5 + 2.79439*10^22 #1^6 &, 2])) || (0.999592 <= x <= 1. && (r == Root[2.75795*10^23 - 2.75795*10^23 x + (1.75577*10^23 - 1.75577*10^23 x) #1 + (-1.62683*10^24 + 1.62683*10^24 x) #1^2 - 3.51153*10^23 #1^3 + (2.53804*10^24 - 5.58878*10^22 x) #1^4 - 5.2673*10^23 #1^5 + 2.79439*10^22 #1^6 &, 3] || r == Root[ 2.75795*10^23 - 2.75795*10^23 x + (1.75577*10^23 - 1.75577*10^23 x) #1 + (-1.62683*10^24 + 1.62683*10^24 x) #1^2 - 3.51153*10^23 #1^3 + (2.53804*10^24 - 5.58878*10^22 x) #1^4 - 5.2673*10^23 #1^5 + 2.79439*10^22 #1^6 &, 4])) || (x == 1. && r == 10.5178)

10^24 does not make a good impression at me.

$\endgroup$
3
  • $\begingroup$ Thanks! I'll test this out as soon as possible! $\endgroup$ Commented Nov 22, 2023 at 17:00
  • $\begingroup$ It seems to give the correct answer, but plugging this into RegionPlot doesn't seem to work though... $\endgroup$ Commented Nov 23, 2023 at 9:21
  • $\begingroup$ @Geigercounter: Plot by parts. Plot[Root[2.75795*10^23-2.75795*10^23 x+(1.75577*10^23-1.75577*10^23 x) #1+(-1.62683*10^24+1.62683*10^24 x) #1^2-3.51153*10^23 #1^3+(2.53804*10^24-5.58878*10^22 x) #1^4-5.2673*10^23 #1^5+2.79439*10^22 #1^6&,1],{x,0,0.999592}] works well to me. $\endgroup$ Commented Nov 23, 2023 at 11:34
1
$\begingroup$

It is advantageous to solve the expression for x first and introduce new parameters a -> \[Alpha] M, r -> \[Rho] M:

solx = x /. Solve[DefiningExpression == 0, x][[1]] /. {a -> \[Alpha] M, r -> \[Rho] M} // Simplify[#, M > 0] & 

$x=\frac{\left(\alpha ^2 (\rho +1)+(\rho -3) \rho ^2\right)^2}{2 \alpha ^2 \rho ^2 \left(\rho ^2-3\right)+\alpha ^4 (\rho +1)^2}$

The constraints transform too

cond1=\[Rho] > 1 > \[Alpha] > 0(* from r>M>a>0 *) cond2=2 (1 + Cos[2/3*ArcCos[-\[Alpha]]]) <= \[Rho] <=2 (1 + Cos[2/3*ArcCos[\[Alpha]]]) 

There is only one parameter 1>\[Alpha]>0 left!

 plot[\[Alpha]_?NumericQ] := ParametricPlot[ Evaluate[{\[Rho], ((-3 + \[Rho]) \[Rho]^2 + \[Alpha]^2 (1 + \ \[Rho]))^2/(\[Alpha]^4 (1 + \[Rho])^2 + 2 \[Alpha]^2 \[Rho]^2 (-3 + \[Rho]^2))}], {\[Rho], 1, 5}, {x, 0, 1}, RegionFunction -> Function[{\[Rho], x}, \[Rho] > 1 > \[Alpha] > 0 && 2 (1 + Cos[(2 ArcCos[-\[Alpha]])/3]) <= \[Rho] <= 2 (1 + Cos[(2 ArcCos[\[Alpha]])/3])], FrameLabel -> {r/M, x}, PlotLabel -> "a/M=" <> ToString[\[Alpha]]] 

examplary simulation \[Alpha]=a/M=0.7

plot[0.7] 

enter image description here

$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.