# Weird (?) NDSolve behavior at boundary

I'm trying to numerical solve the following differential equation (also described here):

$$c_1 (y(x)-1)y'(x) - (c_2 x + x^3 (c_3y(x)-c_4))y'(x)^2+c_5 x^3y'(x)^3 + xy''(x)(y(x)-1)(c_6 + x^2(c_7y(x) - c_8))=0$$

which should admit a cdf $y(x)$, with non-negative density $y'(x)$ on it's entire support.

eq = 2.4 (-1 + y[x]) Derivative[1][y][
x] - (1.2 x - 2.28 x^3 + 2.1 x^3 y[x]) Derivative[1][y][x]^2 +
0.24 x^3 Derivative[1][y][x]^3 +
x (-1 + y[x]) (1.2 - 0.88 x^2 + 0.7 x^2 y[x]) Derivative[2][y][x] == 0;
s = NDSolveValue[{eq, y[1.2] == 1, y[.7] == 0}, y, {x, .7, 1.2}];
Plot[s'[x], {x, 1.1, 1.2}, PlotRange -> All]


which yields

Conversely,

s2 = NDSolveValue[{eq, y[1.2] == .999, y[.7] == 0}, y, {x, .4, 1.2}];
Plot[s2'[x], {x, 1.1, 1.2}, PlotRange -> All]


yields the much smoother

What goes wrong if I set the 'true' boundary condition $y(1.2)=1$?

• The ODE is singular at y[x] == 1, as can be seen from examining the coefficient of y''[x]. – bbgodfrey Nov 23 '17 at 5:32
• Does this mean that there is no 'true' solution for this ODE, or simply that numerical methods won't be able to find it? – bonifaz Nov 23 '17 at 20:20
• Nonlinear ODEs can have many solutions, or no solutions, depending on boundary conditions. For instance, y[x] == 1 is a solution of the general equation, but not with your particular boundary conditions. I shall give this more this more thought later today. – bbgodfrey Nov 23 '17 at 20:52
• If my ODE would have multiple solutions, can I also find the other ones using ndsolve? The one currently obtained doesn't really make sense ($(y(x)=1$ for $x$ large enough, although the derivative is non-zero). – bonifaz Nov 23 '17 at 22:06
• What boundary conditions would you impose on other solutions? – bbgodfrey Nov 23 '17 at 22:59

Rationalize the equation and increase the working precision, and you will get a smooth curve, but it also shows that the y is going to 0 for your range of x.

eq = 2.4 (-1 + y[x]) Derivative[1][y][
x] - (1.2 x - 2.28 x^3 + 2.1 x^3 y[x]) Derivative[1][y][x]^2 +
0.24 x^3 Derivative[1][y][x]^3 +
x (-1 + y[x]) (1.2 - 0.88 x^2 + 0.7 x^2 y[x]) Derivative[2][y][
x] == 0 // Rationalize;

s = NDSolveValue[{eq, y[1.2] == 1, y[.7] == 0}, y, {x, .7, 1.2},WorkingPrecision -> 50];

Plot[s'[x], {x, 1.1, 1.2}, PlotRange -> All]


The higher the WorkingPrecision, the smaller the values.

• It's true that the resulting curve is smooth, but several orders of magnitude off the results obtained from setting $y(1.2)=.99999$. – bonifaz Nov 23 '17 at 20:25
• But if you keep increasing the number of 9's, the numbers keep getting smaller and smaller also – Bill Watts Nov 23 '17 at 20:58
• @bonifaz It appears that there is no solution that goes exactly to 1 at x == 1.2, although one can get arbitrarily close. The solution in the answer above loses precision at about x == 1.15; y == 10^-26. – bbgodfrey Nov 23 '17 at 23:50
• So does this mean that this problem doesn't have a 'true' solution? Or just none that can be computed properly? – bonifaz Nov 24 '17 at 7:35
• @bonifaz Your question is difficult to answer without more details. With respect to your earlier comment, I do not believe the your original ODE has a second solution exactly satisfying the two boundary conditions. – bbgodfrey Nov 24 '17 at 19:29

This is, to me at least, an interesting and somewhat unusual ODE.

Higher Precision Solution

The ODE is singular at y[x] == 1, which acts as an attractor when approached from below. Consequently, there appears to be no solution that goes exactly to 1 at x == 1.2, although some can be arbitrarily close. For instance, in the solution given by Bill Watts,

s = NDSolveValue[{eq, y[7/10] == 0, y[12/10] == 1}, y, {x, 7/10, 12/10},
MaxSteps -> 10^6, WorkingPrecision -> 50];
Plot[s[x], {x, .7, 1.2}, PlotRange -> All, AxesLabel -> {x, y},
LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]
LogPlot[s'[x], {x, .7, 1.2}, PlotRange -> All, AxesLabel -> {x, y'},
LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]


y'[x] decreases approximately exponentially as x approaches 1.2, losing precision at about y[1.15] == 10^-27. Note that the slope ofy[x]atx == 0.7 is

s'[7/10]
(* 4.93285063801117625190985639916133706 *)


A likely course of action to obtain a solution approaching y[1.2] == 1 more rapidly is to increase the initial slope and use a larger WorkingPrecision. The built-in "Shooting" method often does not work well with such nonlinear threshold problems, but trial-and-error often does. For instance, with an initial slope of 5.04,

s = NDSolveValue[{eq, y[7/10] == 0, y'[7/10] == 504/100}, y, {x, 7/10, 12/10},
MaxSteps -> 10^6, WorkingPrecision -> 180, Method -> "StiffnessSwitching"];
Plot[s[x], {x, .7, 1.2}, PlotRange -> All, AxesLabel -> {x, y},
LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]
LogPlot[s'[x], {x, .7, 1.2}, PlotRange -> All, AxesLabel -> {x, y'},
LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]


Although the solution may appear to have a discontinuous slope at about x == .9014, expanding the plots near there shows that the curves vary smoothly, and that y'[x] decays approximately exponentially, finally losing precision at a value of about 10^-110, which probably is sufficiently close to zero. Note that Method -> "StiffnessSwitching" is used to avoid the somewhat mysterious NDSolveValue::nderr error. Also, the corresponding computation with an initial slope of 5.05 or larger fails, because the y'[x] becomes infinite before the curve reaches y[x] == 1.

Search for Alternative Solutions

Because the ODE is nonlinear, it is possible that other solutions approaching y[1.2] == 1 at large angles may exist. Before conducting a search, it is useful to observe that the ODE has a second curve on which the coefficient of y''[x] vanishes.

ss = Flatten@Solve[(1.2 - 0.88 x^2 + 0.7 x^2 y[x]) == 0 // Rationalize, y[x]] // Values
(* {(4 (-15 + 11 x^2))/(35 x^2)} *)


This curve crosses y[x] == 1 at x == 2 Sqrt[5/3]. Now perform a sweep of initial slopes.

sp = ParametricNDSolveValue[{eq, y[7/10] == 0, y'[7/10] == yp}, y, {x, 7/10, 30/10}, {yp},
WorkingPrecision -> 60, MaxSteps -> Infinity, Method -> "StiffnessSwitching"];
Show[Plot[Evaluate@Table[sp[yp][x], {yp, -8, 5, 1/5}], {x, 7/10, 30/10},
PlotRange -> {-1.2, 1.2}, AxesLabel -> {x, y},
LabelStyle -> Directive[Bold, Black, Medium]],
Plot[{ss, 1}, {x, 7/10, 30/10}, PlotRange -> All,
PlotStyle -> Directive[Black, Dashed]], ImageSize -> Large]


The ODE can be integrated through the singular curve (shown as a dashed line), because y'[x] == 0 on the singular curve for the solutions computed. Incidentally, curves integrated from, for instance, y[.7] == - 1/10 also pass smoothly through the singular curve, again with y'[x] == 0. Hence, defining initial conditions {y[x], y'[x]} on the singular surface is not sufficient to determine unique solutions originating there. Returning to the original issue, we see that there is no other family of solutions approaching y[1.2] == 1`.