# Instability of FindRoot solution

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][], sol[c][]}, {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][], sol[c][]}, {c, 0, 1}]


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}] 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}] 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][], sol[c][]}, {c, 0, 1},
PlotLegends -> Placed["Expressions", {0.75, 0.6}]] 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][], sol[c][]}, {c, 0, 1},
PlotLegends -> Placed["Expressions", {0.75, 0.6}]] 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 == x0}, x, {a, 0, 1}];


Visualization:

Plot[sol[t], {t, 0, 1}] 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 == p10, p2 == p20}, {p1, p2}, {c, 0, 1}];


Visualization:

Plot[Through @ sol[t], {t, 0, 1}] 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}] The solutions can be evaluated using NDSolve

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

Show[{cp, ListPlot[Table[{a, x /. sol[a]}, {a, 0, 2, .1}]]}] solF[a_] := FindRoot[{ gl[x, a] == 0 }, {x, .5, 0, 1} ] give the same results!

• 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. – Sander Heinsalu Jul 28 '19 at 22:42