About two weeks ago I've posted the following question on MathOverflow:

Solving a boundary value problem numerically, with high precision.

That is, the ODE is \begin{equation} y''=y^2-t \tag{1} \end{equation} and we're looking for the two slopes $\alpha_1,\alpha_2$ such that the solutions $y_i$ with the initial conditions \begin{align} y_i(0)&=0 \\y_i'(0)&=\alpha_i \quad i=1,2 \end{align} are asymptotic to $\sqrt{t}$ as $t \to +\infty$.

I've been trying to attack it myself for quite some time, and here's what I came up with:

Using the quadratic Taylor approximation based at a point $t_i$ we get \begin{align} y(t_{i+1}) & \approx y(t_i)+y'(t_i)(t_{i+1}-t_i)+\frac{y''(t_i)}{2}(t_{i+1}-t_i)^2\\ &=y(t_i)+y'(t_i)(t_{i+1}-t_i)+\frac{y^2(t_i)-t_i}{2}(t_{i+1}-t_i)^2 \end{align} where I have used the ODE $(1)$ to reduce the final coefficient. Similarly one gets an approximation for the slope: \begin{align} y'(t_{i+1})& \approx y'(t_i)+y''(t_i)(t_{i+1}-t_i)\\ &=y'(t_i)+\left(y^2(t_i)-t_i \right)(t_{i+1}-t_i) \end{align}

In Mathematica, I've implemented this "time step" as follows (The Clip command is used to avoid overflow)

QuadraticTaylor[{y_, p_, t_, dt_}] := {Clip[y + p dt + 1/2 (y^2 - t) dt^2, {-4, 4}],
p + (y^2 - t) dt, t + dt, dt};

Then after trying out different slopes I using the method in the interval $[0,10]$ with $N=10^5$ subintervals and $\alpha=0.92437$ results in a blowup. That is, the code

With[{Parabola = ContourPlot[x == y^2, {x, 0, 10}, {y, -5, 5},ContourStyle->Black]}, 
 Manipulate[T = NestList[QuadraticTaylor, {0, \[Alpha], 0, 10/N}, N]; 
   ListPlot[Map[{#[[3]], #[[1]]} &, T]]], {{\[Alpha], 0.92437`}, -10, 
   10}, {{dt, 0.1}, 0, 1}, {{N, 10000}, 1, 1000000, 1}]]

produced the following plot

enter image description here

while using $\alpha=0.92436$ resulted in an undershoot.

enter image description here

Thus, I thought that one of the critical slopes lies in the interval $[0.92436,0.92437]$. However, using a higher degree Taylor approximation (which I could only do with fewer subintervals, since the evaluations are computationally more expensive), it looks like the interval is actually $$[0.92437548,0.92437549] .$$

Here I'm at a point in which I'm not sure what is the optimal tradeoff between the degree of the Taylor approximation, and the number of subintervals $N$. Moreover, these calculations take a pretty long time (over a minute) for $N>10^6$.

I'd love to hear your thoughts on my attempt, also, if you know a better way to do this please let me know! I'd like to get at least 16 decimal digits on both of the $\alpha$s (my attempt only discusses the positive $\alpha$). Thank you!

P.S. I've also tried using NDSolve to do this, but I'm not sure on how to tweak the settings to best fit my problem. Using default settings seems to get the 7th digit wrong.

  • 4
    $\begingroup$ "numerically, accurately and efficiently" - this feels like one of those situations where you can only pick two out of three... in the meantime, you might be interested in this. $\endgroup$ Commented Feb 2, 2017 at 4:36
  • 1
    $\begingroup$ why do you call this BVP? You have initial conditions both at zero? $\endgroup$
    – Nasser
    Commented Feb 2, 2017 at 5:47
  • $\begingroup$ @Nasser the problem is fundamentally a BVP, relating the behaviour of the function $y(t)$ at $t=0$ and $t=\infty$. I tried solving it using a shooting method which reduces it to an IVP. $\endgroup$
    – user1337
    Commented Feb 2, 2017 at 5:50
  • $\begingroup$ But you haven't transform the ode, why only the bcs? $\endgroup$
    – zhk
    Commented Feb 2, 2017 at 7:00

2 Answers 2


It's somewhat surprising that

  1. brute-force attack turns out to be a good enough solution:

    {inf = 20, pre = 32};
    pa = ParametricNDSolveValue[{y''[x] == y@x^2 - x, y[0] == 0, y'[0] == alpha}, 
       y, {x, 0, inf}, alpha, MaxSteps -> Infinity, WorkingPrecision -> pre];
    {i = 5, init = 92437/10^5};
    (alphavalue = 
       NestWhile[If[pa[# + 10^(-i)][inf] < Sqrt[inf], # + 10^(-i), i++; #] &, init, 
     i <= Floor[pre/2] &] // Quiet) // AbsoluteTiming    
    (* {36.075637, 9243754874375653/10000000000000000} *)    
    p1 = Plot[#, {x, 0, 20}] &@{x^(1/2), pa[alphavalue][x], 
          pa[alphavalue + 10^-(pre/2)][x]} // Quiet

    Mathematica graphics

  2. FindRoot seems not to work well on this problem.

  3. though it doesn't help much in finding an accurate enough $\alpha$, finite difference method (FDM) does a good job in solving the BVP. I'll use pdetoae for the generation of difference equations:

    points = 25;
    grid = Array[# &, points, {0, inf}];
    difforder = 4;
    (* Definition of pdetoae isn't included in this code piece,
       please find it in the link above. *)
    ptoa = pdetoae[y[x], grid, difforder];
    set = {y''[x] == y@x^2 - x, y[0] == 0, y[inf] == Sqrt[inf]};
    ae = MapAt[#[[2 ;; -2]] &, ptoa@set, 1];
    sollst = 
      FindRoot[ae, {y/@ grid, Sqrt[grid]}\[Transpose], WorkingPrecision -> 32][[All, -1]];
    sol = ListInterpolation[sollst, grid, InterpolationOrder -> difforder];
    (* 1.019443775765894394411167742141 <- Not accurate at all. 
       This can be improved by using a larger points and difforder,
       but not quite economic. *)
    (* However, the solution i.e. value of y is accurate enough: *)
    p2 = Plot[sol[x], {x, 0, inf}, PlotStyle -> {Red, Thick, Dashed}];
    Show[p1, p2]

    Mathematica graphics

  • $\begingroup$ So what you are trying to say is that, there will be blow up/down using MMA? $\endgroup$
    – zhk
    Commented Feb 2, 2017 at 13:23
  • 3
    $\begingroup$ @mmm To be precise, the solution will blow up/undershoots using any software, if we're solving the equation with initial condition $y(0)=0,\ y'(0)\approx0.92436$. The solution is extremely sensitive to initial value $y'(0)$. ("$α=0.92437$ results in a blowup, $α=0.92436$ resulted in an undershoot. " ) And if I understand the question correctly, an accurate enough $\alpha$ (rather than a stable solution) is exactly what OP needs. $\endgroup$
    – xzczd
    Commented Feb 2, 2017 at 13:41
  • $\begingroup$ In my comment, I have argued that, Maple is able to solve the BVP without any problem, "numerically, accurately and efficiently". $\endgroup$
    – zhk
    Commented Feb 2, 2017 at 13:47
  • $\begingroup$ @MMM Yeah, I saw it, but can a accurate enough $y'(0)$ be calculated from the solution given by Maple? $\endgroup$
    – xzczd
    Commented Feb 2, 2017 at 14:07
  • $\begingroup$ Maple gives `alpha=0.924375469638557'. $\endgroup$
    – zhk
    Commented Feb 2, 2017 at 14:09


Solving IVP

I have converted your BVP in to a IVP, like this,

$$y1^{'}(t)=y2(t),\,\, y2^{'}(t)=y1^{2}(t)-t,\,\, y1(0)=y10, \,\,y2(0)=y20$$

Now, I will try to solve the above system for a set {y10,y20} of initial conditions,

soln[y10_?NumericQ, y20_?NumericQ] :=First@NDSolve[{y1'[t]==y2[t],y2'[t]-y1[t]^2 + t == 0, 
y1[0] == y10, y2[0] == y20}, {y1, y2}, {t, 0, 5}];
Plot[Evaluate[{y1[t], y2[t]} /. soln[#, #] & /@Range[-1.5, 0.4, 0.1]], {t, 0, 5}, 
PlotRange -> All,PlotStyle -> {Red, Green}, AxesLabel -> {"t", "y"},
PlotLegends -> {"y1", "y2"}]

ParametricPlot[Evaluate[{y1[t], y2[t]} /. soln[#, #] & /@Range[-1.5, 0.4, 0.1]], {t, 0, 3},
PlotRange -> All,PlotStyle -> {Red, Green}, AxesLabel -> {"y1", "y2"}]

Doing different experimentation with initial conditions, it is evident that the solution blowup when {y10,y20}>0.4, i.e., Range[0.4, 2, 0.1].

Solving BVP

In maple, we can solve the bvp "numerically, accurately and efficiently" with out any problem,

odeplot(A1, [[t,y(t), color=red],[t,sqrt(t), color=green,linestyle=3]],0..N,axes=boxed)

enter image description here

Here is the value of alpha


enter image description here

I am observing a very strange behavior on Mathematica,

eq = y''[x] == y[x]^2 - x;
ibcs = {y[0] == 0, y[N1] == Sqrt[N1]};
sol = First@NDSolve[Join[{eq}, ibcs], y[x], {x, 0, N1}]
Plot[{y[x]} /. sol, {x, 0, N1}]

For N1=4, I get this

enter image description here

and for N1=5

enter image description here

Is it possible, that there are multiple solutions? or there is something else going on?

I am thankful to @halirutan for teaching me how to input the bc y[x]=Sqrt[x].

  • $\begingroup$ Thank you for your answer. There are indeed two such solutions: one with a positive initial slope, and one with a negative one. $\endgroup$
    – user1337
    Commented Feb 2, 2017 at 14:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.