1
$\begingroup$

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

Solution of 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:

Solution of Equation 2

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}]

enter image description here

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

Improve this question
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Apr 18 at 1:15
  • $\begingroup$ Thank you for your comment! I tried to reformulate my question. I hope it is clearer now. $\endgroup$
    – Lednacek
    Apr 19 at 14:00

1 Answer 1

Reset to default
2
$\begingroup$

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}]

enter image description here

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}]

enter image description here

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]

enter image description here

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

enter image description here

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

enter image description here

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}]

enter image description here

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]

enter image description here

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

Improve this answer
$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$
    – Lednacek
    Apr 20 at 13:50
  • $\begingroup$ @Lednacek You are right. I am correcting the typo now. Thanks. $\endgroup$
    – bbgodfrey
    Apr 20 at 13:59
  • $\begingroup$ 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. $\endgroup$
    – Lednacek
    Apr 20 at 15:41
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Apr 20 at 15:49
  • $\begingroup$ Thank you very much!!! This is awesome! $\endgroup$
    – Lednacek
    Apr 20 at 16:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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