1
$\begingroup$

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

enter image description here

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

enter image description here

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

$\endgroup$
  • 2
    $\begingroup$ The ODE is singular at y[x] == 1, as can be seen from examining the coefficient of y''[x]. $\endgroup$ – bbgodfrey Nov 23 '17 at 5:32
  • $\begingroup$ Does this mean that there is no 'true' solution for this ODE, or simply that numerical methods won't be able to find it? $\endgroup$ – bonifaz Nov 23 '17 at 20:20
  • 1
    $\begingroup$ 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. $\endgroup$ – bbgodfrey Nov 23 '17 at 20:52
  • $\begingroup$ 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). $\endgroup$ – bonifaz Nov 23 '17 at 22:06
  • $\begingroup$ What boundary conditions would you impose on other solutions? $\endgroup$ – bbgodfrey Nov 23 '17 at 22:59
4
$\begingroup$

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]

enter image description here

The higher the WorkingPrecision, the smaller the values.

$\endgroup$
  • $\begingroup$ It's true that the resulting curve is smooth, but several orders of magnitude off the results obtained from setting $y(1.2)=.99999$. $\endgroup$ – bonifaz Nov 23 '17 at 20:25
  • 1
    $\begingroup$ But if you keep increasing the number of 9's, the numbers keep getting smaller and smaller also $\endgroup$ – Bill Watts Nov 23 '17 at 20:58
  • 1
    $\begingroup$ @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. $\endgroup$ – bbgodfrey Nov 23 '17 at 23:50
  • $\begingroup$ So does this mean that this problem doesn't have a 'true' solution? Or just none that can be computed properly? $\endgroup$ – bonifaz Nov 24 '17 at 7:35
  • 1
    $\begingroup$ @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. $\endgroup$ – bbgodfrey Nov 24 '17 at 19:29
4
$\begingroup$

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]

enter image description here

enter image description here

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]

enter image description here

enter image description here

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]

enter image description here

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.

$\endgroup$

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.