# how to solve ODE with boundary at infinity

y''[x]-x y[x]==0
y[0]==AiryAi[0], y[infinity]==0


the analytic solution to this ODE is the Airy function

y[x]=AiryAi[x]


if I use NDSolve to solve this ODE, I do not know how to express the boundary condition at infinity. I first tried to set y[100]=0, the error message from mathematica is:

NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. The computed solution may match the boundary conditions poorly. >>

NDSolve::berr: There are significant errors {0.,5.04843*10^273} in the boundary value residuals. Returning the best solution found. >>

then I tried to replace x by x=tan[u], with 0<u<Pi/2, now my boundary condition is at finite location, however I still receive error message when I use NDSolve

NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. The computed solution may match the boundary conditions poorly. >>

NDSolve::mxst: Maximum number of 10000 steps reached at the point u == 1.5651256862447316. >>

NDSolve::berr: There are significant errors {0.,1.331263934430340*10^660} in the boundary value residuals. Returning the best solution found. >>

NDSolve::mxst: Maximum number of 10000 steps reached at the point u == 1.5651764393985212. >>

the equation I am posing here can be analytically solved, however my really equations can only be solved by NDSolve with boundary condition at Infinity.

so I am asking for the possibility to solve this equation numerically.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

According to different kindly advice, numerically solve an ODE with unbounded region is impossible. A more practical way is by shooting method. By replacing y[Infinity]==0 to y[100]==10^{-10}, the numerical result agrees very well with the exact solution y[x]=AiryAi[x]

-
@MichaelE2, unfortunately I have to use NDSolve – 3c. Aug 13 '14 at 0:31
NDSolve uses various methods to solve for numerical approximations to DEs over some fixed, bounded interval. Your boundary conditions define the interval -- you cannot use NDSolve for an unbounded interval. (That's not to say it's impossible, there are plenty of hacks to make a function that extends the interval when necessary or that assumes for x>>0, y(x)=0. But it's not clear what you want in this case.) – Kellen Myers Aug 13 '14 at 1:13
@KellenMyers, thanks. I think I need just to use shooting method to approximate the solution. – 3c. Aug 13 '14 at 1:32

You can use ParametricNDSolve to implement a shooting method.

Define a finite version of "infinity".

inf = 5;


Define the differential equation and its initial conditions, parameterised by the initial gradient y'[0] == dy0. For simplicity, I set y[0] == 1.

deqn = {y''[x] - x y[x] == 0, y[0] == 1, y'[0] == dy0};


Compute the numerical solution parameterised by dy0.

ydysol = ParametricNDSolve[deqn, y, {x, 0, inf}, dy0][[1]]


Find the value of the initial gradient dy0 that makes the solution go to zero at "infinity". I choose to minimise y[inf]^2 w.r.t. dy0.

dysol = FindMinimum[((y[dy0] /. ydysol)[inf])^2 // Evaluate, {dy0, -1}]

(* {3.35979*10^-25, {dy0 -> -0.729012}} *)


Sanity check this solution.

AiryAi'[0]/AiryAi[0] // N

(* -0.729011 *)


Plot the solution.

Plot[(y[dy0] /. ydysol /. dysol[[2]])[x] // Evaluate, {x, 0, inf}]


You can use larger values of inf, but the above approach throws warning messages and eventually becomes unstable near x == inf, which is to be expected.

-
I now see that the OP notes that the solution can be obtained directly with NDSolve, which presumably does an implicit shooting method, so there is no need for the explicit ParametricNDSolve shooting method in this case. Oh well, at least it shows how ParametricNDSolve can be used to implement the shooting method if you really want to do it yourself. So much automation ... – Stephen Luttrell Aug 13 '14 at 11:54

The finite element method can be used on this problem if we make a change of variables to convert the domain $[0, \infty)$ to a finite interval. I believe only MachinePrecision is available in FEM. Since AiryAi vanishes so rapidly, it will make a precise result for a large argument difficult to obtain.

Another difficulty in obtaining a precise solution is that all but one solution of the differential equation diverges at infinity. So any numerical error tends to lead to divergence. This makes the finite element method seem a better thing to try than the shooting method. (Of course, as shown in another answer and elsewhere, one can approximate the shooting method by not insisting on carrying the solution all the way to infinity. One of the difficulties is determining an appropriate value for say y[100]. The value y[100] == 10^-10 is way off the true value, in terms of precision. It would be difficult to know that if we did not know the exact solution. See below.)

First we transform the differential equation with the substitutions $x = \tan t$ and $u(t) = y(\tan t)$. (You might want to evaluate the NestList separately if you cannot see what it does.)

eqn = (-x[t] y[x[t]] + D[D[y[x[t]], t]/D[x[t], t], t]/D[x[t], t]) x'[t]^3 == 0 /. x -> Tan;
bvp = {eqn, u[0] == AiryAi[0], u[Pi/2] == 0} /.
Solve[NestList[D[#, t] &, u[t] == y[Tan[t]], 2], {y''[Tan[t]], y'[Tan[t]], y[Tan[t]]}]
(*
{{Sec[t]^6 (-Tan[t] u[t] - 2 Cos[t]^3 Sin[t] u'[t] + Cos[t]^4 u''[t]) == 0,
u[0] == 1/(3^(2/3) Gamma[2/3]), u[π/2] == 0}}
*)


It might seem nice to humans to get rid of the Sec[t]^6 factor, but Mathematica does not seem to care.

A straight forward application produces an OK solution.

airy = u[ArcTan[#]] & /. First@NDSolve[bvp, u, {t, 0, Pi/2}, Method -> {"FiniteElement"}];


To get a better solution, refine the mesh with something like "MeshOptions" -> {MaxCellMeasure -> 0.002}:

airy = u[ArcTan[#]] & /.
First@NDSolve[bvp, u, {t, 0, Pi/2},
Method -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> 0.002}}];


We can compare the precision and accuracy of the two solutions with the following ploting function. (One of the precision plots had a poorly chosen automatic range; hence, the complexity of the command.)

precaccplot[domain : {t1_?NumericQ, t2_?NumericQ} : {0, 10},
precopts : OptionsPattern[LogPlot]] := GraphicsRow[{
LogPlot[Abs[(airy[t] - AiryAi[t])/AiryAi[t]], {t, t1, t2},
PlotLabel -> Precision, precopts],
LogPlot[Abs[airy[t] - AiryAi[t]], {t, t1, t2},
PlotLabel -> Accuracy]}]


The first, less precise solution:

precaccplot[PlotRange -> {10^-8, 10^5}]


The second solution with a finer mesh:

precaccplot[]


Beyond some point, the precision loss is going to be catastrophic, although the accuracy might be acceptable (if, say, you're adding the solution to a fairly large number).

precaccplot[{0, 100}]
AiryAi[100.]
airy[100.]


(*
2.63448*10^-291
4.56296*10^-49
*)

-