0
$\begingroup$

Using FindRoot on a smooth function of a parameter and a variable, the resulting solution should be a smooth function of the parameter, by the Implicit Function Theorem. However, Mathematica returns a very spiky unstable function. One way to fix the problem seems to avoid Newton's method. Are there others, in particular that work with more than one variable?

MWE:

Clear[cdf, s0, s, sol]
cdf[x_] := -2 (x - Sinh[1/2] + 
     Sinh[1/2 - x])/(-2 (1 - Sinh[1/2] + Sinh[1/2 - 1]))
s0[a_] := 
 s0[a] = x /. FindRoot[D[x*(cdf[a - x] - cdf[x]), x] == 0, {x, 0.5}]
Plot[s0[a], {a, 0, 1}]
(*A different root finding method does not cause instability.*)
s[a_] := s[a] = 
  x /. FindRoot[D[x*(cdf[a - x] - cdf[x]), x] == 0, {x, 0.01, 0.99}]
Plot[s[a], {a, 0, 1}]
(*The problem is worse with two variables*)
sol[c_] := 
 sol[c] = {p1, p2} /. 
   FindRoot[{D[p1*(cdf[p2 - p1] - cdf[p1]), p1] == 0, 
 D[(p2 - c)*(1 - cdf[p2 - p1]), p2] == 0}, {{p1, 0.2}, {p2, 0.5}}]
Plot[{sol[c][[1]], sol[c][[2]]}, {c, 0, 1}]
(*Using non-Newton methods of root finding does not help with two
variables.*)
sol[c_] := 
 sol[c] = {p1, p2} /. 
   FindRoot[{D[p1*(cdf[p2 - p1] - cdf[p1]), p1] == 0, 
     D[(p2 - c)*(1 - cdf[p2 - p1]), p2] == 0}, {{p1, 0.1, 0.8}, {p2, 
      0.2, 0.9}}]
Plot[{sol[c][[1]], sol[c][[2]]}, {c, 0, 1}]
$\endgroup$
1
$\begingroup$

FindRoot is sensitive to having a good starting value.

Clear[cdf, s0]

cdf[x_] := -2 (x - Sinh[1/2] + 
     Sinh[1/2 - x])/(-2 (1 - Sinh[1/2] + Sinh[1/2 - 1]))
s0[a_] := s0[a] = x /. FindRoot[D[x*(cdf[a - x] - cdf[x]), x] == 0, {x, 0.5}]

Plot[s0[a], {a, 0, 1}]

enter image description here

The fist part of the plot suggests that a lower starting value is needed, e.g., {x, 0.1}

Clear[cdf, s0]

cdf[x_] := -2 (x - Sinh[1/2] + 
     Sinh[1/2 - x])/(-2 (1 - Sinh[1/2] + Sinh[1/2 - 1]))
s0[a_] := s0[a] = x /. FindRoot[D[x*(cdf[a - x] - cdf[x]), x] == 0, {x, 0.1}]

Plot[s0[a], {a, 0, 1}]

enter image description here

Likewise for the two variable case:

Clear[sol]

sol[c_] := 
 sol[c] = {p1, p2} /. 
   FindRoot[{D[p1*(cdf[p2 - p1] - cdf[p1]), p1] == 0, 
     D[(p2 - c)*(1 - cdf[p2 - p1]), p2] == 0}, {{p1, 0.2}, {p2, 0.5}}]

Plot[{sol[c][[1]], sol[c][[2]]}, {c, 0, 1},
 PlotLegends -> Placed["Expressions", {0.75, 0.6}]]

enter image description here

The plot suggests better starting values might be {{p1, 0.1}, {p2, 0.9}}

Clear[sol]

sol[c_] := 
 sol[c] = {p1, p2} /. 
   FindRoot[{D[p1*(cdf[p2 - p1] - cdf[p1]), p1] == 0, 
     D[(p2 - c)*(1 - cdf[p2 - p1]), p2] == 0}, {{p1, 0.1}, {p2, 0.9}}]

Plot[{sol[c][[1]], sol[c][[2]]}, {c, 0, 1},
 PlotLegends -> Placed["Expressions", {0.75, 0.6}]]

enter image description here

|improve this answer|||||
$\endgroup$
3
$\begingroup$

I like using NDSolveValue for problems like this. Let's take your first example:

cdf[x_] := -2(x - Sinh[1/2] + Sinh[1/2 - x])/(-2 (1 - Sinh[1/2] + Sinh[1/2 - 1]))

eqn = D[x*(cdf[a - x] - cdf[x]), x] == 0;

To use NDSolveValue, we need an ODE and an initial condition. The initial condition can be obtained by solving the equation for a particular value of a:

x0 = x /. FindRoot[Evaluate[eqn /. a->0], {x, 1}]

2.98238*10^-16

The ODE is obtained by differentiating the equation with respect to the parameter a:

ode = D[eqn /. x->x[a], a];

Putting these together we have:

sol = NDSolveValue[{ode, x[0] == x0}, x, {a, 0, 1}];

Visualization:

Plot[sol[t], {t, 0, 1}]

enter image description here

For your second example we have:

eqn = {
    D[p1*(cdf[p2 - p1] - cdf[p1]), p1] == 0, 
    D[(p2 - c)*(1 - cdf[p2 - p1]), p2] == 0
};

Initial condition:

{p10, p20} = {p1, p2} /. FindRoot[Evaluate[eqn /. c->0], {p1, .1}, {p2, 1}]

{0.114562, 0.83039}

ODE:

deqn = D[eqn /. {p1 -> p1[c], p2 -> p2[c]}, c];

Using NDSolveValue:

sol = NDSolveValue[{deqn, p1[0] == p10, p2[0] == p20}, {p1, p2}, {c, 0, 1}];

Visualization:

Plot[Through @ sol[t], {t, 0, 1}]

enter image description here

|improve this answer|||||
$\endgroup$
2
$\begingroup$

ContourPlot shows you the solutions

gl[x_, a_]  = Evaluate[D[x*(cdf[a - x] - cdf[x]), x]] 
cp=ContourPlot[gl[x, a] == 0, {a, 0, 2}, {x, 0, 1/2},FrameLabel -> {a, x}]

enter image description here

The solutions can be evaluated using NDSolve

sol[a_] := NSolve[{ gl[x, a] == 0, 0 < x < 1}, x , Reals][[1]]

Show[{cp, ListPlot[Table[{a, x /. sol[a]}, {a, 0, 2, .1}]]}]

enter image description here

solF[a_] := FindRoot[{ gl[x, a] == 0 }, {x, .5, 0, 1} ] give the same results!

|improve this answer|||||
$\endgroup$
  • $\begingroup$ The solutions work for 1 variable, upvoted. solF[a_] := FindRoot puts some dots off the curve for me, but the curve is smooth. The NSolve...Show method is very slow for 2 variables, I have not obtained the results yet. Using Evaluate...ContourPlot3D gives Could not converge... and other errors for 2 variables - I need to work on correcting the code to check whether this method can deal with 2+ variables. $\endgroup$ – Sander Heinsalu Jul 28 '19 at 22:42

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.