# Solving implicit function numerically and plotting the solution against a parameter

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:

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):

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.

• D[1 - e, e] seems a somewhat confusing way to write -1, and I wanted to check, is -1 correct? Commented Jul 5, 2023 at 15:40
• @MichaelE2, yes, -1 would be a simpler way to put it.
– ppp
Commented Jul 5, 2023 at 21:46
• @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.) Commented Jul 6, 2023 at 4:36
• @MichaelE2, this is what I get: shorturl.at/fBHOX
– ppp
Commented Jul 6, 2023 at 5:31

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


• 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?
– ppp
Commented Jul 7, 2023 at 7:13
• @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. Commented Jul 7, 2023 at 11:51
• 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!
– ppp
Commented Dec 20, 2023 at 1:15

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}


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


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

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

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]


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


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