2
$\begingroup$

Consider an objective function

v = ((i + lambda) Uem - (-1 + e) Uun)/(i (1 - e + i + lambda))

where

n = 1; beta = 0.6; W = 60; i=0.05;
L = (w/(beta (a e)^beta))^(1/(beta - 1));
lambda = ((1 - e) n L)/(1 - e n L);
Uem = w - 1/(w (1 - e));
b[e_] := b /. Solve[a W (1 + i)^(1/(1 - e) + 1/lambda - 1)- b \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(\*FractionBox[\(1\), \(lambda\)] - 1\)]\(\((1 + i)\)^j\)\) == b + a W, b][[1]];
Uun = b[e] - 1/b[e];

with the conditions of $0\leq e \leq 1$ and $U_{em} > U_{un}$.

I would like to do some plotting by taking the following steps.

Step 1: Find e that maximizes v which will be yielded as a function of w with the parameter a, i.e., $e = f(w; a)$.

I plotted the f function for different values of a = ${0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}$ using the following code:

brf[a_, w_] = ArgMax[{v, e >= 0, e <= 1, Uem > Uun}, e];
Plot[{brf[0.1, w] Boole[0.2 <= w <= 6], brf[0.2, w] Boole[0.5 <= w <= 6], brf[0.3, w] Boole[0.8 <= w <= 6], brf[0.4, w] Boole[1.1 <= w <= 6], brf[0.5, w] Boole[1.4 <= w <= 6], brf[0.6, w] Boole[1.7 <= w <= 6], brf[0.7, w] Boole[2 <= w <= 6], brf[0.8, w] Boole[2.3 <= w <= 6], brf[0.9, w] Boole[2.6 <= w <= 6]}, {w, 0, 6}, MaxRecursion -> 5, PlotLegends -> LineLegend[Table[a, {a, 0.1, 0.9, 0.1}], LegendLabel -> a]]

which produces enter image description here

Step 2: Find e and w as the solution of $\frac{\partial f}{\partial w}=\frac{e}{w}$ for varying values of a$\in [0,1]$. Graphically, this is equivalent to finding in the graph below the tangent point of a curve where a ray from the origin (straight dotted line) cross. (I produced this graph by manually adding the rays to the graph above.)

enter image description here

Step 3: Plot each of e and w of the tangent points against a$\in [0,1]$ on the x-axis.

The code for is:

Clear["Global`*"];
n = 1; beta = 0.6; W = 60;
L = (w/(beta (a e)^beta))^(1/(beta - 1));
lambda = ((1 - e) n L)/(1 - e n L);
Uem = w - 1/(w (1 - e));
b[e_] := b /. Solve[a^2 W (1 + i)^(1/(1 - e) + 1/lambda -1) - b \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(\*FractionBox[\(1\), \(lambda\)] - 1\)]\(\((1 + i)\)^j\)\) ==b + a^2 W, b][[1]];
Uun = b[e] - 1/b[e];
v = ((i + lambda) Uem - (-1 + e) Uun)/(i (1 - e + i + lambda));
f[w_, a_, i_] = ArgMax[{v, e >= 0, e <= 1, Uem > Uun}, e];
With[{i = 5/100}, sol[a_] := sol[a] = FindRoot[D[f[w, a, i], w] == e/w, {{e, 0.5}, {w, 0.5}}]; {estar[a_], wstar[a_]} := {e /. sol[a], w /. sol[a]}; estarplot = Plot[estar[a], {a, 0, 1}, AxesLabel -> {s, SuperStar[\[CurlyEpsilon]]}, PlotRange -> {{0, 1}, {0, 1}}]; wstarplot = Plot[wstar[a], {a, 0, 1}, AxesLabel -> {s, SuperStar[w]}, PlotRange -> {{0, 1}, {0, 10}}];]
e = estarplot;
w = wstarplot;
Style[Row[{e, w}], ImageSizeMultipliers -> {0.6, 0.6}]

And I get empty graphs along with a list of error messages as follows. Can anyone help please?

enter image description here

Edit: As Daniel Huber's answer below suggests, it seems the problem is related to FindRoot where w gets replaced by a number but I use w as a variable for differentiation. But consider the following code, which is almost the same method applied to a different objective function.

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

In this code, the issue raised by Dainel Huber does not emerge but yields plots nicely. Why is this so? What's the difference? In any case, does this mean that the issue is not what Daniel Huber mentioned?

$\endgroup$
7
  • 1
    $\begingroup$ Testing your functions with a few typical values may help you debug this code. For example, f[.5,.5,.5,5/100] generates error messages. I would check that the functions b[e] and f[...] are working properly for several values of the parameters before trying to use the functions in plots or elsewhere. $\endgroup$
    – LouisB
    Dec 20, 2023 at 1:56
  • $\begingroup$ ` FindRoot[D[f[e, w, 1, i], w] == e/w, {{e, 0.5}, {w, 0.5}}]` "e" in f[e...] should be a variable not a number. $\endgroup$ Dec 20, 2023 at 9:14
  • $\begingroup$ @LouisB, thanks for your comment. I already checked and confirmed that the functions b[e] and f[...] work. I edit the post to add this. Can you please verify? Thanks much! $\endgroup$
    – ppp
    Dec 20, 2023 at 14:02
  • $\begingroup$ @DanielHuber, thanks for your comment! But isn't e in f[e...] already a variable, not a number? $\endgroup$
    – ppp
    Dec 20, 2023 at 14:20
  • 1
    $\begingroup$ With the revised function and stated parameters I get a plot similar to the one in the question. When I evaluate f[5], I get 0.584, which says there is a peak in $v$ at $(e,w)=(0.584,5)$. But, when I evaluate Plot[v /. w -> 5, {e, .25, .75}] I don't see a peak. No peak in Plot3D[v, {e, .15, .95}, {w, 3, 12}] either. I suspect Maximize is returning a value on the boundary, but I haven't pursued it. By the way, you don't need to Solve to find $b$. Just push the symbols around a little to get $b=...$. $\endgroup$
    – LouisB
    Dec 20, 2023 at 14:52

2 Answers 2

3
+150
$\begingroup$

You can use an interpolating function to simplify and speed-up the code:

B[a_] := B[a] = Interpolation[Table[{w, brf[a, w]}, {w, 3 a - .1, 6, .2}] // Quiet]

Then, the slope is determined by

wc[a_] := FindRoot[B[a]'[wc] wc == B[a][wc], {wc, 0.4 + 4.1 a}][[1, 2]]

such that

Show[
    Plot[Evaluate@Table[B[a][w] Boole[3 a - .1 < w], {a, .1, .9, .1}], {w, 0, 6}] // Quiet,
    Plot[Evaluate@Table[B[a]'[wc[a]] w, {a, .1, .9, .1}], {w, 0, 6}, PlotStyle -> Dashed]
]

which produces

enter image description here

--

The last step, if I understood correctly, is to plot $w_c(a)$ and $f(w_c(a))$ as a function of $a$. This is simple to do:

Table[{a, wc[a]}, {a, .1, .9, .05}] // ListPlot
Table[{a, B[a][wc[a]]}, {a, .1, .9, .05}] // ListPlot

which produces

enter image description here

and

enter image description here

You can clearly see the phase transition around $a~\sim0.21$.

$\endgroup$
18
  • 1
    $\begingroup$ @ppp 1) the function is not analytic at the bottom so the curving up is an artifact of the interpolation I am using. It is just numerical error. You can reduce it by adding more points, the code will run more slowly though. The rest of the curve is unaffected. 2) Just some parameters that I found work best. For example 3a-.1 is just your choice of lower bound on the Booles in your own code, .2,.5,.8,... for a=.1,.2,.3,.... 3) wc is the critical w, i.e., the point where $f'(w_c)=e/w_c$. The eq inside FindRoot is just this equation, $w_c$ is the point where the line meets the curve $\endgroup$ Jan 2 at 17:06
  • 1
    $\begingroup$ @ppp 1) There are artifacts in your graph too. It is possible to reduce the error by increasing the number of points before interpolating (namely, replace the .2 inside {w, 3 a - .1, 6, .2} by a smaller number), but it slows down the computation. Increase the number of points as much as you want, and as much as your patience allows. 2) yes, B is the same as $f$, the double bracket is syntactic sugar, formally equivalent to B[a, w]. $\endgroup$ Jan 3 at 21:13
  • 1
    $\begingroup$ @ppp Add Joined -> True to ListPlot[...]. Result. $\endgroup$ Jan 8 at 18:36
  • 1
    $\begingroup$ @ppp Works fine on my laptop, see i.imgur.com/eXn7NeO.png Did you try on a fresh kernel? What version are you using? Can you run the code line by line and let me know which specific part makes it hang? $\endgroup$ Jan 9 at 15:36
  • 1
    $\begingroup$ I tried fresh kernel and it worked! Thanks so much! $\endgroup$
    – ppp
    Jan 12 at 17:34
2
$\begingroup$

Try:

sol[0.1]

and you get

General::ivar: 0.5` is not a valid variable.

enter image description here

Now change the limit of w in FindRoot to 0.4:

 FindRoot[D[f[e, w, a, i], w] == e/w, {{e, 0.5}, {w, 0.4}}]

and you get:

General::ivar: 0.4` is not a valid variable.

Therefore, the error is that w get replaced by a number, but you use w as a variable for differentiation.

The solution is to use another variable e.g. ww for differentiation and afterwards replace ww by w.

$\endgroup$
10
  • $\begingroup$ Thanks, Daniel! I think I could use ww for differentiation, i.e. to replace w with ww in f[] and f[]=e/w, but I'm not sure how to replace it by w afterwards. Can you please help me on this? $\endgroup$
    – ppp
    Dec 21, 2023 at 21:14
  • $\begingroup$ /. ww->w or in long form: ReplaceAll $\endgroup$ Dec 22, 2023 at 8:01
  • $\begingroup$ Thanks for your reply! I'm sorry but I'm struggling with how to incorporate ww and replacing it to w afterwards as you suggested. May I ask for some more details how to include these into my code? I really appreciate it! $\endgroup$
    – ppp
    Dec 22, 2023 at 14:28
  • $\begingroup$ By rewriting your code, I just saw an additional error. f does not depend on w. Therefore, you can not take the derivative of f relative to w. Neither does f depend on i , but this should not give an error. $\endgroup$ Dec 22, 2023 at 17:49
  • $\begingroup$ If I run v as defined in the code, I get a long output of v as a very complicated function of e, w, a, and i. Then, if we find e that maximizes v, wouldn't I get that e as a function of w, a, and i? $\endgroup$
    – ppp
    Dec 22, 2023 at 23:23

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.