# How to numerically solve an equation

I would like to solve for $$r_2$$ the following two (independent) equations:

$$x_{Min}(a,b,c,r_1,F)=x_{Max}(a,b,c,r_2,F)$$ (1)

and

$$x_{Min}(a,b,c,r_2,F)=x_{Max}(a,b,c,r_1,F)$$ (2)

where $$a$$, $$b$$, $$c$$, $$r_1$$ and $$F$$ are parameters such that $$a>b>0$$, $$c>0$$, $$r_1>0$$ and $$F\geq0$$ and where $$x_{Min}$$ and $$x_{Max}$$ are functions defined below.

Ideally, I would like to solve these equations algebraically, but I do not think this is possible. Hence, I am trying to solve them numerically for given values of parameters $$a$$, $$b$$, $$c$$ and $$r_1$$. I am interested in how the solution $$r_2$$ depends on $$F$$. Ultimately, I would like to plot $$r_2(F)$$.

Here is how the $$x_{Min}$$ and $$x_{Max}$$ functions are defined:

Consider the following function $$f$$ of $$x$$ depending on 4 other parameters

$$f(x)=(c-ax)(c+x(2rx-a))/(3bx^2)$$

For $$a>b>0$$, $$c>0$$, $$r>0$$, this function has a unique maximum for $$x>c/a$$. I solve algebraically for this maximum: $$x_{fmax}=argmax_{x>c/a} f$$ and $$maxf=f(x_{fmax})$$.

Then I solve for $$x_{Min}$$ and $$x_{Max}$$ such that for a given constant $$F\in[0,maxf]$$, $$f(x_{Min})=f(x_{Max})=maxf-F$$.

$$x_{Min}$$ and $$x_{Max}$$ are functions of the parameters $$a$$, $$b$$, $$c$$, $$r$$ and $$F$$. Then, for given $$r_1$$, I want to find $$r_2$$ such that

$$x_{Min}(a,b,c,r_1,F)=x_{Max}(a,b,c,r_2,F)$$ and $$x_{Min}(a,b,c,r_2,F)=x_{Max}(a,b,c,r_1,F)$$

When $$r_1=1$$, the first equation has a solution $$\in[1,\infty)$$ increasing in $$F$$ and the second one has a solution $$\in[0,1]$$ decreasing in $$F$$.

For the first equation, I have an (almost working) solution using FindRoot (see code below), but this solution is not perfect as the resulting plot is not continuous. This solution does not work at all for the second equation.

Here is my code which derives the two equations to solve from the initial function $$f$$ and then attempts to solve the two equations.

(*Initial function definition*)
f[a_, b_, c_, r_,
x_] = ((c - a x) (c + x (-a + 2 x r)))/(3 b x^2)

(*Solving for the maximum*)
derf[a_, b_, c_, r_, x_] = FullSimplify[D[f[a, b, c, r, x], x]];
Solve[derf[a, b, c, r, x] == 0, x]
(*I pick manually the relevant root from the three solutions*)
xfmax[a_, b_, c_, r_] = Simplify[(2 2^(1/3) 3^(1/6) (3 I + Sqrt[3]) a^2 c r +
2^(2/3) 3^(
1/
3) (1 -
I Sqrt[
3]) (a^2 c r (9 c r +
Sqrt[3] Sqrt[c r (-4 a^2 + 27 c r)]))^(2/3))/(12 a^(5/3)
r (c r (9 c r + Sqrt[3] Sqrt[c r (-4 a^2 + 27 c r)]))^(1/3)),Assumptions -> {c > 0, a > 0, b > 0, r > 0}];
(*Evaluate the maximum*)
maxf[a_, b_, c_, r_] = Simplify[f[a, b, c, r, xfmax[a, b, c, r]], Assumptions -> {c > 0, a > 0, b > 0, a > b, a > c}]

(*Solving for the two intersections*)
Solve[CC - F == f[a, b, c, r, x], x]
(*I pick manually the two relevant roots*)
x1[a_, b_, c_, r_, F_, CC_] =Simplify[-((-a^2 + 3 b CC - 3 b F - 2 c r)/(
6 a r)) - (12 a^2 c r - (-a^2 + 3 b CC - 3 b F - 2 c r)^2)/(3 2^(
2/3)
a r (2 a^6 - 18 a^4 b CC + 54 a^2 b^2 CC^2 - 54 b^3 CC^3 +
18 a^4 b F - 108 a^2 b^2 CC F + 162 b^3 CC^2 F +
54 a^2 b^2 F^2 - 162 b^3 CC F^2 + 54 b^3 F^3 - 24 a^4 c r +
36 a^2 b c CC r + 108 b^2 c CC^2 r - 36 a^2 b c F r -
216 b^2 c CC F r + 108 b^2 c F^2 r + 60 a^2 c^2 r^2 -
72 b c^2 CC r^2 + 72 b c^2 F r^2 +
16 c^3 r^3 + \[Sqrt]((2 a^6 - 18 a^4 b CC +
54 a^2 b^2 CC^2 - 54 b^3 CC^3 + 18 a^4 b F -
108 a^2 b^2 CC F + 162 b^3 CC^2 F + 54 a^2 b^2 F^2 -
162 b^3 CC F^2 + 54 b^3 F^3 - 24 a^4 c r +
36 a^2 b c CC r + 108 b^2 c CC^2 r - 36 a^2 b c F r -
216 b^2 c CC F r + 108 b^2 c F^2 r + 60 a^2 c^2 r^2 -
72 b c^2 CC r^2 + 72 b c^2 F r^2 + 16 c^3 r^3)^2 +
4 (12 a^2 c r - (-a^2 + 3 b CC - 3 b F - 2 c r)^2)^3))^(
1/3)) +
1/(6 2^(1/3)
a r) (2 a^6 - 18 a^4 b CC + 54 a^2 b^2 CC^2 - 54 b^3 CC^3 +
18 a^4 b F - 108 a^2 b^2 CC F + 162 b^3 CC^2 F +
54 a^2 b^2 F^2 - 162 b^3 CC F^2 + 54 b^3 F^3 - 24 a^4 c r +
36 a^2 b c CC r + 108 b^2 c CC^2 r - 36 a^2 b c F r -
216 b^2 c CC F r + 108 b^2 c F^2 r + 60 a^2 c^2 r^2 -
72 b c^2 CC r^2 + 72 b c^2 F r^2 +
16 c^3 r^3 + \[Sqrt]((2 a^6 - 18 a^4 b CC + 54 a^2 b^2 CC^2 -
54 b^3 CC^3 + 18 a^4 b F - 108 a^2 b^2 CC F +
162 b^3 CC^2 F + 54 a^2 b^2 F^2 - 162 b^3 CC F^2 +
54 b^3 F^3 - 24 a^4 c r + 36 a^2 b c CC r +
108 b^2 c CC^2 r - 36 a^2 b c F r - 216 b^2 c CC F r +
108 b^2 c F^2 r + 60 a^2 c^2 r^2 - 72 b c^2 CC r^2 +
72 b c^2 F r^2 + 16 c^3 r^3)^2 +
4 (12 a^2 c r - (-a^2 + 3 b CC - 3 b F - 2 c r)^2)^3))^(
1/3), Assumptions -> {c > 0, a > 0, b > 0, a > b, a > c}];

x2[a_, b_, c_, r_, F_, CC_] = Simplify[-((-a^2 + 3 b CC - 3 b F - 2 c r)/(
6 a r)) + ((1 +
I Sqrt[
3]) (12 a^2 c r - (-a^2 + 3 b CC - 3 b F - 2 c r)^2))/(6 2^(
2/3)a r (2 a^6 - 18 a^4 b CC + 54 a^2 b^2 CC^2 - 54 b^3 CC^3 +
18 a^4 b F - 108 a^2 b^2 CC F + 162 b^3 CC^2 F +
54 a^2 b^2 F^2 - 162 b^3 CC F^2 + 54 b^3 F^3 - 24 a^4 c r +
36 a^2 b c CC r + 108 b^2 c CC^2 r - 36 a^2 b c F r -
216 b^2 c CC F r + 108 b^2 c F^2 r + 60 a^2 c^2 r^2 - 72 b c^2 CC r^2 + 72 b c^2 F r^2 +
16 c^3 r^3 + \[Sqrt]((2 a^6 - 18 a^4 b CC +
54 a^2 b^2 CC^2 - 54 b^3 CC^3 + 18 a^4 b F -
108 a^2 b^2 CC F + 162 b^3 CC^2 F + 54 a^2 b^2 F^2 -
162 b^3 CC F^2 + 54 b^3 F^3 - 24 a^4 c r +
36 a^2 b c CC r + 108 b^2 c CC^2 r - 36 a^2 b c F r -
216 b^2 c CC F r + 108 b^2 c F^2 r + 60 a^2 c^2 r^2 -
72 b c^2 CC r^2 + 72 b c^2 F r^2 + 16 c^3 r^3)^2 +
4 (12 a^2 c r - (-a^2 + 3 b CC - 3 b F - 2 c r)^2)^3))^(
1/3)) -
1/(12 2^(1/3)
a r) (1 - I Sqrt[3]) (2 a^6 - 18 a^4 b CC + 54 a^2 b^2 CC^2 -
54 b^3 CC^3 + 18 a^4 b F - 108 a^2 b^2 CC F + 162 b^3 CC^2 F +
54 a^2 b^2 F^2 - 162 b^3 CC F^2 + 54 b^3 F^3 - 24 a^4 c r +
36 a^2 b c CC r + 108 b^2 c CC^2 r - 36 a^2 b c F r -
216 b^2 c CC F r + 108 b^2 c F^2 r + 60 a^2 c^2 r^2 -
72 b c^2 CC r^2 + 72 b c^2 F r^2 +
16 c^3 r^3 + \[Sqrt]((2 a^6 - 18 a^4 b CC + 54 a^2 b^2 CC^2 -
54 b^3 CC^3 + 18 a^4 b F - 108 a^2 b^2 CC F +
162 b^3 CC^2 F + 54 a^2 b^2 F^2 - 162 b^3 CC F^2 +
54 b^3 F^3 - 24 a^4 c r + 36 a^2 b c CC r +
108 b^2 c CC^2 r - 36 a^2 b c F r - 216 b^2 c CC F r +
108 b^2 c F^2 r + 60 a^2 c^2 r^2 - 72 b c^2 CC r^2 +
72 b c^2 F r^2 + 16 c^3 r^3)^2 +
4 (12 a^2 c r - (-a^2 + 3 b CC - 3 b F - 2 c r)^2)^3))^(
1/3)];

(*Definitions of xMin and xMax*)
xMin[a_, b_, c_, r_, F_] = x2[a, b, c, r, F, maxf[a, b, c, r]];
xMax[a_, b_, c_, r_, F_] = x1[a, b, c, r, F, maxf[a, b, c, r]];

(*Attempting to solve the two equations with a=10, b=c=r1=1*)
(*Equation 1*)
NH1[F_] := NH1[F] = r2 /. FindRoot[xMin[10, 1, 1, 1, F] == xMax[10, 1, 1, r2, F], {r2, 1.001}]
Plot[{NH1[F]}, {F, 0, 2}]

(*Equation 2*)
NH2[F_] := NH2[F] = r2 /. FindRoot[xMin[10, 1, 1, r2, F] == xMax[10, 1, 1, 1, F{r2, 0.9}]
Plot[{NH2[F]}, {F, 0, 2}]



Here are the plots that I get from my solution: From Equation 1

The solution should be a continuous increasing function of $$F$$ starting from 1 at $$F=0$$. So the plot is not perfect. Equation 2 is even worse:

For some reason, there seems to be a problem in Equation 2: even just plotting the relation as a function of $$r_2$$ for $$F=1$$ does not seem to work

Plot[xMin[10, 1, 1, r2, 1] - xMax[10, 1, 1, 1, 1], {r2, 0, 1}]


Could someone please suggest a better method to solve for $$r_{2}(F)$$ in the two equations?

• Please state simply what you are trying to accomplish. Describing your goal in terms of code that, as you say, does not give the desired result is not helpful. Commented Apr 18, 2022 at 1:15
• Thank you for your comment! I tried to reformulate my question. I hope it is clearer now. Commented Apr 19, 2022 at 14:00

To begin, find the maximum value of f subject to the constraints in the question.

f = ((c - a x) (c + x (-a + 2 x r)))/(3 b x^2);
MaxValue[{f, a > b > 0, c > 0, r > 0, x > c/a}, x]
(* Piecewise[{{(a^2 + 2*c*r - b*Root[16*a^4*c*r + 27*a^2*c^2*r^2 -
18*a^2*b*c*r*#1 - a^2*b^2*#1^2 + b^3*#1^3 & , 2])/(3*b),
r > 0 && a > 0 && 0 < c <= a^2/(8 r) && 0 < b < a},
{0, r > 0 && a > 0 && c > a^2/(8*r) && 0 < b < a}}, -Infinity] *)


Note that this calculation yields an additional constraint, c <= a^2/(8 r) or, equivalently, r <= a^2/(8 c). To illustrate, for the parameters in the question,

Plot[f /. {a -> 10, b -> 1, c -> 1, r -> 1}, {x, .1, 5}]


As r is increased, fmax decreases until it drops below the axis at r = 12.5, violating the requirement from the question that fmax be positive.

Plot[f /. {a -> 10, b -> 1, c -> 1, r -> 12.5}, {x, .1, .3}]


Next, compute and plot the values of xmin and xmax, as described in the question.

MaxValue[{f, x > c/a} /. {a -> 10, b -> 1, c -> 1, r -> 1}, x] // N;
Transpose[Chop[Table[SolveValues[% - F == f /. {a -> 10, b -> 1, c -> 1, r -> 1}, x],
{F, 0, %, .2}], 10^-7]] // Rest;
ListLinePlot[%, DataRange -> {0, %%}, AxesLabel -> {F, "xmin,xmax"},
PlotStyle -> Black]


Performing this same calculation for other values of r, say {1, 2, 4, 6, 8.5}, and superimposing the plots with Show yields

and with, say, {10^-8, 10^-5, 10^-3, 10^-2, .05, .2, 1}

Intersections of the red curves with the black curves provides qualitative solutions to the first two equations in the question. To obtain quantitive solutions, define

s[F_, rr_] := Module[{t, ff = f /. {a -> 10, b -> 1, c -> 1, r -> rr},
coa = 1/10}, t = MaxValue[{ff, x > coa}, x] // N;
SolveValues[t - F == ff, x] // Rest]


which provides {xmin, xmzx}. Then,

s1[F_] := Module[{tst}, tst = First[s[F, 1]];
r /. FindRoot[Last[s[F, r]] == tst, {r, 5, 1, 12.5},
Evaluated -> False]]


solves the first equation in the question.

Plot[s1[F], {F, 0.001, 3.2}, MaxRecursion -> 5, AxesLabel -> {F, r2}]


and

s2[F_] := Module[{tst}, tst = Last[s[F, 1]];
r /. FindRoot[First[s[F, r]] == tst, {r, .5, 0, 1},
Evaluated -> False]]


solves the second.

Plot[s2[F], {F, 0.001, 3.5},  MaxRecursion -> 5,
AxesLabel -> {F, r2}, PlotRange -> All]


Solutions do not exist for larger values of F than those in these plots. Each plot requires several minutes to compute

• Thank you SOOO much, this is amazing! I have two quick questions on your solution: (1) In the definition of s1[f] and s2[f] you are using the s[F,rr] function, right? So it should be FindRoot[Last[s[F, r]] and FindRoot[First[s[F, r]], correct? Commented Apr 20, 2022 at 13:50
• @Lednacek You are right. I am correcting the typo now. Thanks. Commented Apr 20, 2022 at 13:59
• Thank you! My 2nd question is: the definition of s[F,rr] assigns particular values to the parameters a, b, and c. Is there an easy way to postpone this assignment for the plotting stage? Ultimately, what I would like to do is to use something like Manipulate and have sliding buttons with which I can show how the plot of s1[F] and s2[F] varies with these parameters. I know this will be tricky because varying these parameters will change the optimisation constraints, but I would like to be able to visualise this at least on a common subrange for which the constraints are satisfied. Commented Apr 20, 2022 at 15:41
• @Lednacek Use s[F_, aa_, bb_, cc_, rr_] := Module[{t, ff = f /. {a -> aa, b -> bb, c -> cc, r -> rr}}, t = MaxValue[{ff, x > cc/aa}, x] // N; SolveValues[t - F == ff, x] // Rest], and pass aa, bb, cc through s1 and s2 as arguments. Note that different values for these constants may require changing {r, 5, 1, 12.5} in Findroot. Also, the range of values for which F is meaningful may change. too. Commented Apr 20, 2022 at 15:49
• Thank you very much!!! This is awesome! Commented Apr 20, 2022 at 16:49