1
$\begingroup$

Defining two functions via NArgMax and plotting, there is an intersection and the Plot function works. But trying to solve for the intersection with various methods from Mathematica Stackexchange generates various errors. How to solve an equation system in which the functions are defined by the outputs of NArgMax?

MWE

    Clear[pl, ph, pil, pih, plbr, phbr]
pil[pl_, ph_] := pl*(ph - 2*pl);
pih[pl_, ph_] := ph*(1 - ph + pl);
plbr[ph_] = plbr[ph_?NumericQ] := NArgMax[pil[pl, ph], pl]
phbr[pl_] = phbr[pl_?NumericQ] := NArgMax[pih[pl, ph], ph]
NSolve[{plbr[ph] == pl, phbr[pl] == ph}, {pl, ph}]
FindRoot[{plbr[ph] == pl, 
  phbr[pl] == ph}, {{pl, 0.5, 0, 1}, {ph, 0.9, 0, 2}}]

This question and this one show a workaround:

    plt3 = Show[plt1, plt2]
Point@
 Graphics`Mesh`FindIntersections@plt3

which finds some intersections, fails to find some and adds some wrong ones.

My full code with various solution attempts:

    Clear[a, b, c, e, ql, qh, pl, ph, cl, ch, cdf, pil, pih, plbr]
$Assumptions = 
  0 <= a < b < c && 0 < e < 1 && 0 < ql < qh && 0 < pl < ph < c*qh && 
   pl < c*ql && 0 < cl < ch < ph && cl < pl;
cdf[x_] := 
  Piecewise[{{0, x < a}, {e*(x - a)/(b - a), 
     a <= x < b}, {e + (1 - e)*(x - b)/(c - b), b <= x < c}, {1, 
     x >= c}}];
pil[pl_, ph_] := (pl - cl)*(cdf[(ph - pl)/(qh - ql)] - cdf[pl/ql]);
pih[pl_, ph_] := (ph - ch)*(1 - cdf[(ph - pl)/(qh - ql)]);
a = 0.001; b = 0.5; c = 1.0; e = 0.9; cl = 0.; ch = 0.5;
ql = 1; qh = 1.6; qh2 = 1.8;
plbr[ph_] = plbr[ph_?NumericQ] := NArgMax[pil[pl, ph], pl]
phbr[pl_] = phbr[pl_?NumericQ] := NArgMax[pih[pl, ph], ph]
(*NSolve[{plbr[ph]\[Equal]pl,phbr[pl]\[Equal]ph},{pl,ph}]
FindRoot[{plbr[ph]\[Equal]pl,phbr[pl]\[Equal]ph},{{pl,(cl+c*ql)/2,cl,\
c*ql},{ph,(ch+c*qh)/2,ch,c*qh}}](*"The function value {-0.5`+Null} is \
not a list of numbers with dimensions"...*)
FindRoot[{plbr[phbr[pl]]\
\[Equal]pl},{{pl,(cl+c*ql)/2,cl,c*ql},{ph,(ch+c*qh)/2,ch,c*qh}}]*)

plbr[ph_?NumericQ] := NArgMax[pil[pl, ph], pl]
phbr[pl_?NumericQ] := NArgMax[pih[pl, ph], ph]
Reap[odsoln = 
   NDSolve[Evaluate[{yl'[ph] == plbr'[ph], yh'[pl] == phbr'[pl],(*yl[
      ch]\[Equal]plbr[ch]-ch,yh[cl]\[Equal]phbr[cl]-cl,*)
      yl[c*qh] == plbr[c*qh] - c*qh, yh[c*ql] == phbr[c*ql] - c*ql, 
      WhenEvent[{yl[ph] == 0, yh[pl] == 0}, Sow[{pl, ph}]]}], {yl[ph],
      yh[pl]}, {{pl, cl, c*ql}, {ph, ch, 
      c*qh}}];](*[[2,1]]*)
With[{eql = plbr[ph] - pl, 
  eqh = phbr[pl] - ph}, 
 roots = Reap[
    NDSolve[{yl'[ph] == D[eql, ph], yh'[pl] == D[eqh, pl], 
       yl[0] == (eql /. ph -> 0), yh[0] == (eqh /. pl -> 0), 
       WhenEvent[eql == 0 && eqh == 0, Sow[{pl, ph}]]}, {}, {{pl, cl, 
        c*ql}, {ph, ch, c*qh}}];][[2, 1]]]
plbr[ph_] := NArgMax[pil[pl, ph], pl]
phbr[pl_] := NArgMax[pih[pl, ph], ph]
(*plbr[ph_?NumericQ]:=NArgMax[pil[pl,ph],pl]
plbr2[ph_?NumericQ]:=NArgMax[pil2[pl,ph],pl]*)
yl = 
 ImplicitRegion[{plbr[ph] == pl}, {pl, ph}]
yh = ImplicitRegion[{plbr[ph] == pl}, {pl, ph}];
Solve[{{pl, ph} \[Element] f1, {pl, ph} \[Element] f2}, {pl, ph}]

Also tried

 plbr[ph_] := 
 FunctionInterpolation[
  ArgMax[{pil[pl, ph], cl <= pl <= c*ql}, pl], {ph, ch, c*qh}]
phbr[pl_] := 
 FunctionInterpolation[
  ArgMax[{pih[pl, ph], ch <= ph <= c*qh}, ph], {pl, cl, c*ql}]
NSolve[{plbr[ph] == pl, phbr[pl] == ph}, {pl, 
  ph}](*With plbr[ph_?NumericQ] gives pl\[Rule]1.` Null,ph\[Rule]1.` \
Null*)

The plotting code that works:

    Clear[a, b, c, e, ql, qh, qh2, pl, ph, cl, ch, cdf, pil, pil2, pih, \
plbr, plbr2, plt1, plt2]
$Assumptions = 
  0 <= a < b < c && 0 < e < 1 && 0 < ql < qh < qh2 && 
   0 < pl < ph < c*qh && pl < c*ql && 0 < cl < ch < ph && cl < pl;
pil[pl_, ph_] := (pl - cl)*(cdf[(ph - pl)/(qh - ql)] - cdf[pl/ql]);
pil2[pl_, ph_] := (pl - cl)*(cdf[(ph - pl)/(qh2 - ql)] - cdf[pl/ql]);
pih[pl_, ph_] := (ph - ch)*(1 - cdf[(ph - pl)/(qh - ql)]);
pih2[pl_, ph_] := (ph - ch)*(1 - cdf[(ph - pl)/(qh2 - ql)]);
cdf[x_] := 
  Piecewise[{{0, x < a}, {e*(x - a)/(b - a), 
     a <= x < b}, {e + (1 - e)*(x - b)/(c - b), b <= x < c}, {1, 
     x >= c}}];
a = 0.001; b = 0.5; c = 1.0; e = 0.9; cl = 0; ch = 0.5;
ql = 1; qh = 1.6; qh2 = 1.8;
plbr[ph_] = plbr[ph_?NumericQ] := NArgMax[pil[pl, ph], pl]
plbr2[ph_] = plbr2[ph_?NumericQ] := NArgMax[pil2[pl, ph], pl]
phbr[pl_] = phbr[pl_?NumericQ] := NArgMax[pih[pl, ph], ph]
phbr2[pl_] = phbr2[pl_?NumericQ] := NArgMax[pih2[pl, ph], ph]

plt1 = ParametricPlot[{{plbr[ph], ph}, {plbr2[ph], ph}}, {ph, ch, 
    c*qh2}, PlotStyle -> {Orange, Blue}, 
   AxesLabel -> {Subscript[P, L], Subscript[P, H]}];
plt2 = ParametricPlot[{{pl, phbr[pl]}, {pl, phbr2[pl]}}, {pl, cl, 
    c*ql}, PlotStyle -> {Red, Green}];
Show[plt1, plt2]
$\endgroup$
0

Your Answer

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

Browse other questions tagged or ask your own question.