Return to Answer

Fixed typos

Source Link

edited Sep 9, 2014 at 16:42

244.7k
18
351
774

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.)

A straight forward application produceproduces an OK solution.

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

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

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

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 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. See below.)

A straight forward application produce an OK solution.

To get a better solution, improve the mesh:

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

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

A straight forward application produces an OK solution.

To get a better solution, refine the mesh with something like "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.)

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).

Source Link

created Sep 9, 2014 at 14:37

Michael E2

244.7k
18
351
774

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 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. 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 produce an OK solution.

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

To get a better solution, improve the mesh:

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

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}]

Mathematica graphics

The second solution with a finer mesh:

precaccplot[]

Mathematica graphics

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

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

Mathematica graphics

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