2
$\begingroup$

This post is a clarified version of the question in this one.

Consider an implicit function, $F(e,w,a,b,i,n)=0$, whose code is:

D[w - 1/((1 - e) w), e] ==  D[1 - e, e] ((-((-a (i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) + (-1 + a^2) (-1 + e)^2 w - a (-1 + e) (i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) w^2)/(a (-1 + e) i (1 - e + i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) w))) - (-((-a (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1))) + (-1 + a^2) (-1 +   e) (-1 + e - i) w -  a (-1 + e) (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(
    1 - n e (1/(b (a e)^(b)))^(1/(b - 1))) w^2)/( a (-1 + e) i (1 - e + i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) w))))

From the above implicit function, $e$ can be implicitly defined as a function of $w$ with the rest of the variables taken as constant.

My task takes three steps as follows.

Step 1: I numerically examined how this function $e=e(w)$ looks like for any given value of $a = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1$ along with the other parameters fixed at $n=1$, $b=0.7$, $i=0.1$. My code for this is:

Clear["Global`*"]
n = 1;
eqns[i_, b_] = {D[w - 1/((1 - e) w), e] == 
D[1 - e, e] ((-((-a (i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(
         1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) + (-1 + 
         a^2) (-1 + e)^2 w - 
      a (-1 + e) (i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(
         1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) w^2)/(
     a (-1 + e) i (1 - e + i + (
        n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(
        1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) w))) - (-((-a (
       n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(
       1 - n e (1/(b (a e)^(b)))^(1/(b - 1))) + (-1 + a^2) (-1 + 
         e) (-1 + e - i) w - 
      a (-1 + e) (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(
       1 - n e (1/(b (a e)^(b)))^(1/(b - 1))) w^2)/(
     a (-1 + e) i (1 - e + i + (
        n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(
        1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) w)))), 0 <= i <= 1, 0 <= a <= 1, 0 <= e <= 1, w >= 0, 0 <= b <= 1};
ContourPlot[Evaluate@ Table[Simplify@    eqns[1/10, 7/10][[1]], {a, {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1}}], {w, 0, 5}, {e, 0, 1}, PlotPoints -> 50, MaxRecursion -> 4, FrameLabel -> (Style[#, 14] & /@ {w, e}), RotateLabel -> False, PlotLegends -> Placed[LineLegend[{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1}, 
LegendLabel -> "a"], {.7, .4}]]

And the result is: enter image description here

Step 2: I would like to use this $e=e(w)$ function numerically found as above to obtain the solution for $e$ and $w$ that solve $\frac{\partial e}{\partial w}=\frac{e}{w}$. That is, the pair of solutions $(e^* , w^* )$ for each $a$ will be obtained at the tangent point between the $e=e(w)$ function and the ray from the origin. Graphically, it will look like this (In the figure, I indicated the solution pairs in red dot):

enter image description here

As is clearly described in the above figure, I will be able to obtain a unique $(e^* , w^* )$ for each $a$.

Step 3: Finally, what I would like to produce two separate plots, one plotting $e^* $ against $a$ and the other plotting $w^* $ against $a$. These will look something like below. Step 1 is already done. Now, I need help on Steps 2 & 3. (But I don't need to produce the second figure, the one with the red dots added. All I need to produce are the last two plots.) Please help.

enter image description here enter image description here

Improve this question
$\endgroup$
4
  • $\begingroup$ D[1 - e, e] seems a somewhat confusing way to write -1, and I wanted to check, is -1 correct? $\endgroup$
    – Michael E2
    Commented Jul 5, 2023 at 15:40
  • $\begingroup$ @MichaelE2, yes, -1 would be a simpler way to put it. $\endgroup$
    – ppp
    Commented Jul 5, 2023 at 21:46
  • $\begingroup$ @pppI get this, no errors: i.sstatic.net/ZP7RV.jpg -- What was the first error or two you got? (I'm using V13.2.1, if that matters.) $\endgroup$
    – Michael E2
    Commented Jul 6, 2023 at 4:36
  • $\begingroup$ @MichaelE2, this is what I get: shorturl.at/fBHOX $\endgroup$
    – ppp
    Commented Jul 6, 2023 at 5:31

3 Answers 3

Reset to default
5
$\begingroup$

This is clunky and can probably be optimized but it seems to work:

Define:

Clear["Global`*"]
f[e_, w_, a_, b_, i_, n_] = D[w - 1/((1 - e) w), e] - D[1 - e, e] ((-((-a (i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) + (-1 + a^2) (-1 + e)^2 w - a (-1 + e) (i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) w^2)/(a (-1 + e) i (1 - e + i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) w))) - (-((-a (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1))) + (-1 + a^2) (-1 +   e) (-1 + e - i) w -  a (-1 + e) (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(
    1 - n e (1/(b (a e)^(b)))^(1/(b - 1))) w^2)/( a (-1 + e) i (1 - e + i + (n (1 - e) (1/(b (a e)^(b)))^(1/(b - 1)))/(1 - n e (1/(b (a e)^(b)))^(1/(b - 1)))) w))))

In other words, the contour curves you have plotted are those for which $$f(e(w),w) = 0. \tag{1}$$ Along such a curve, you will also have $$ \frac{d}{dw} \left[f(e(w),w)\right] = 0 \quad \Rightarrow \quad \frac{\partial f}{\partial e} \frac{de}{dw} + \frac{\partial f}{\partial w} = 0 $$ and the requirement that $de/dw = e/w$ then leads us to the condition $$ e \frac{\partial f}{\partial e} + w \frac{\partial f}{\partial w} = 0. \tag{2} $$ Equations (1) and (2) can be solved numerically for a fixed value of $a$, and the results can be plotted as a function of $a$ as follows:

With[{i = 1/10, b = 7/10, n = 1},
 sol[a_] := sol[a] =
   FindRoot[{
     f[e, w, a, b, i, n] == 0, 
     e D[f[e, w, a, b, i, n], e] + w D[f[e, w, a, b, i, n], w] == 0}, 
     {{e, 0.5}, {w, 0.5}}];
 {estar[a_], wstar[a_]} := {e /. sol[a], w /. sol[a]};
 estarplot = 
  Plot[estar[a], {a, 0.1, 1}, AxesLabel -> {a, SuperStar[e]}];
 wstarplot = 
  Plot[wstar[a], {a, 0.1, 1}, AxesLabel -> {a, SuperStar[w]}];
 ]
estarplot
wstarplot

enter image description here

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ thanks for your help. It worked great. May I ask one more related question? The constant parameter values used are i=0.1, b=0.7, n=1. I'm also trying to use Manipulate to see how the results might change depending on the variations of these three parameters. Can you please help? $\endgroup$
    – ppp
    Commented Jul 7, 2023 at 7:13
  • $\begingroup$ @ppp I think that adapting Michael E2's method to use ParametricNDSolve (with i, b, and n as parameters) would be your best bet in that case. My method is really not well-adapted to the sort of rapid updating that Manipulate requires. $\endgroup$
    – Michael Seifert
    Commented Jul 7, 2023 at 11:51
  • $\begingroup$ I was struggling to applying your method here to a similar problem but without success as you can see in this post: mathematica.stackexchange.com/questions/295283/… May I sincerely ask for your help once more? Thank you in advance! $\endgroup$
    – ppp
    Commented Dec 20, 2023 at 1:15
6
$\begingroup$

Aother way, using NDSolve to construct $w^*$ and $e^*$:

The basic ODE system is obtained by differentiating the OP's equation (myEq) with respect to a and constraining the line through the origin to be perpendicular to the gradient of the equation (at $(w(a),e(a))$:

D[{ (* (w,e) to satisfy the given equation *)
    myEQ /. Equal -> Subtract,
    (* vector (w,e) perpendicular to gradient *)
    {w, e} . D[myEQ /. Equal -> Subtract, {{w, e}}]
    } /. {w -> w[a], e -> e[a]}
  , a] == 0

Code:

myEQ = First@eqns[i, b];

(* solve for an initial point *)
icEQ = Block[{n = 1, i = 1/10, b = 7/10, a = 1/2},
      {#, Dt[#, w]} &@myEQ
      ] /. {e -> e[w]} /. {e'[w] -> e[w]/w} /. {e[w] -> e} // Simplify;
icEQ = icEQ /. Equal -> Subtract // Simplify;
icSOL = NSolve[icEQ == 0 && 0 < e < 1 && 0 < w < 2, {e, w}]
(*  {{e -> 0.577473, w -> 1.6782}}  *)

{e\[FivePointedStar], w\[FivePointedStar]} = NDSolveValue[{
   ode = Solve[
      Block[{n = 1, i = 1/10, b = 7/10},
       D[{
           myEQ /. Equal -> Subtract,
           {e, w} . D[myEQ /. Equal -> Subtract, {{e, w}}]
           } /. {w -> w[a], e -> e[a]}
         , a] == 0
       ],
      {e'[a], w'[a]}
      ] /. Rule -> Equal,
   e[1/2] == (e /. First@icSOL), w[1/2] == (w /. First@icSOL)},
  {e, w}, {a, $MachineEpsilon, 1}]

(* cp = OP's ContourPlot[..] *)
Show[
  cp,
  Graphics[{
    PointSize@Medium,
    Point[
     Table[{w\[FivePointedStar][a], e\[FivePointedStar][a]}, {a, 0.1, 
       0.9, 0.1}]],
    Opacity[0.3],
    Line[
     Table[{{0, 0}, {4, 
        4 e\[FivePointedStar][a]/w\[FivePointedStar][a]}}, {a, 0.1, 1,
        0.1}]]
    }],
  ParametricPlot[{w\[FivePointedStar][a], e\[FivePointedStar][a]}, 
   Evaluate@Flatten@{a, e\[FivePointedStar]["Domain"]}, 
   PlotStyle -> Magenta, PlotRange -> All], Axes -> True,
  AspectRatio -> Automatic, ImageSize -> Large, 
  PlotRange -> {{0, 4}, {-0.2, 1.1}}
  ] /. {0.7`, 0.4`} -> {0.88`, 0.4`}

enter image description here

ListLinePlot[{e\[FivePointedStar], w\[FivePointedStar]}, 
 PlotLegends -> 
  Block[{e\[FivePointedStar], w\[FivePointedStar]}, 
   HoldForm /@ {e\[FivePointedStar], w\[FivePointedStar]}], 
 PlotRange -> All]

enter image description here

In response to comments: Large image of a session in V13.3.0 for Mac OS X ARM (64-bit) (June 3, 2023); code cut and pasted from MSE (the code for cp from the OP as indicated and the code from "Code" in this answer):

https://i.sstatic.net/fCKPo.jpg

Improve this answer
$\endgroup$
14
  • $\begingroup$ Michael, thanks! I tried to reproduce your result, but your code yields errors. Can you please confirm? $\endgroup$
    – ppp
    Commented Jul 6, 2023 at 4:15
  • $\begingroup$ I still don't reproduce your result. Can you please confirm your code? Thank you very much! $\endgroup$
    – ppp
    Commented Jul 18, 2023 at 10:04
  • $\begingroup$ "I still don't reproduce your result." NONE of it worked? $\endgroup$
    – Michael E2
    Commented Jul 18, 2023 at 23:22
  • $\begingroup$ correct, this is what I get: shorturl.at/fBHOX $\endgroup$
    – ppp
    Commented Jul 19, 2023 at 5:01
  • $\begingroup$ I would really appreciate if you could help me on this whenever you have time. Thank you so much in advance! $\endgroup$
    – ppp
    Commented Jul 22, 2023 at 7:45
4
$\begingroup$

This needs some work, that needs be done for every tangent point. As an example I show it for a==1. To start, we define a function of e,w and a:

f[e_, w_, a_] = Simplify[eqns[1/10, 7/10]][[1, 1]];

The next step is to create some points on the curve around the tangent point. For a==1 this will be points in the region 2<w<.2.8 and 0.4<e<0.7:

pts = Table[{w, e /. Evaluate@
     NSolve[{f[e, w, 1] == 0, 0.4 < e < 0.7}, e, Reals][[1]]}, {w, 2, 2.8, 0.01}]

Next we fit some suitable function through the points. I choose a 4th order poly:

poly = LinearModelFit[pts, {1, x, x^2, x^3, x^4}, x]

enter image description here

Finally we can solve the tangent equation:

xt = x /. Solve[{fun'[x] == poly[x]/x, 2 < x < 3}, x, Reals][[1]]

2.53014

Lets check if this worked by plotting the curves:

ContourPlot[Evaluate[f[e, w, 1] == 0], {w, 0, 5}, {e, 0, 1}, 
 Epilog -> {Point[{xt, fun[xt]}], Line[{{0, 0}, {xt, fun[xt]}}]}]

enter image description here

Now this has to be repeated for the other tangent points.

Improve this answer
$\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.