I have a nonlinear set of equations with a boundary condition at infinity. Consequently I have to shoot for the boundary. This question has some good ideas but there is an extra complication in my case; my boundary conditon at infinity is the second derivative going to zero and this is approached asymptotically. I will have more problems the same and this is minimum example I can demonstrate.

The equations and initial conditions at zero are here

eqns = {-Derivative[3][a[0]][η] == (-(1/2))*α*
           (-1 + 2*Derivative[1][a[0]][η]^2 + Derivative[1][a[1]][η]^2 + 
              Derivative[1][a[2]][η]^2 + Derivative[1][b[1]][η]^2 + 
              Derivative[1][b[2]][η]^2 - 
       2*a[0][η]*Derivative[2][a[0]][η] - 
              a[1][η]*Derivative[2][a[1]][η] - a[2][η]*Derivative[2][a[2]][
                  η] - b[1][η]*Derivative[2][b[1]][η] - 
       Derivative[1][b[1]][η] - Derivative[3][a[1]][η] == 
         (1/2)*(-4*α*Derivative[1][a[0]][η]*Derivative[1][a[1]][η] - 
              2*α*Derivative[1][a[1]][η]*Derivative[1][a[2]][η] - 
              2*α*Derivative[1][b[1]][η]*Derivative[1][b[2]][η] + 
              2*α*a[1][η]*Derivative[2][a[0]][η] + 
                Derivative[2][a[1]][η] + α*a[2][η]*
        Derivative[2][a[1]][η] + 
              α*a[1][η]*Derivative[2][a[2]][η] + α*b[2][η]*
                Derivative[2][b[1]][η] + α*b[1][η]*
       1 - Derivative[1][a[1]][η] - Derivative[3][b[1]][η] == 
         (1/2)*(-4*α*Derivative[1][a[0]][η]*Derivative[1][b[1]][η] + 
              2*α*Derivative[1][a[2]][η]*Derivative[1][b[1]][η] - 
              2*α*Derivative[1][a[1]][η]*Derivative[1][b[2]][η] + 
              2*α*b[1][η]*Derivative[2][a[0]][η] + 
              α*b[2][η]*Derivative[2][a[1]][η] - α*b[1][η]*
                Derivative[2][a[2]][η] + 
       2*α*a[0][η]*Derivative[2][b[1]][η] - 
              α*a[2][η]*Derivative[2][b[1]][η] + α*a[1][η]*
                Derivative[2][b[2]][η]), 2*Derivative[1][b[2]][η] - 
           Derivative[3][a[2]][η] == 
         (1/2)*(α - α*Derivative[1][a[1]][η]^2 - 
              4*α*Derivative[1][a[0]][η]*Derivative[1][a[2]][η] + 
              α*Derivative[1][b[1]][η]^2 + 
                  η] + α*a[1][η]*Derivative[2][a[1]][η] + 
              2*α*a[0][η]*Derivative[2][a[2]][η] - 
       -2*Derivative[1][a[2]][η] - Derivative[3][b[2]][η] == 
         (1/2)*(-2*α*Derivative[1][a[1]][η]*Derivative[1][b[1]][η] - 
              4*α*Derivative[1][a[0]][η]*Derivative[1][b[2]][η] + 
              2*α*b[2][η]*Derivative[2][a[0]][η] + 
              α*b[1][η]*Derivative[2][a[1]][η] + α*a[1][η]*
                Derivative[2][b[1]][η] + 2*α*a[0][η]*Derivative[2][b[2]][
ic0 = {a[0][0] == 0, Derivative[1][a[0]][0] == 0, a[1][0] == 0, 
       b[1][0] == 0, Derivative[1][a[1]][0] == 0, Derivative[1][b[1]][0] == 
         0, a[2][0] == 0, b[2][0] == 0, Derivative[1][a[2]][0] == 0, 
       Derivative[1][b[2]][0] == 0}; 

I then have a shooting function that determines the values of the second derivative at a point that is getting towards infinity.

shoot[{(a0_)?NumberQ, a1_, b1_, a2_, b2_}, inf_] := 
   (sol = First[NDSolve[Evaluate[Join[eqns, ic0, 
shoot[{(a0_)?NumberQ, a1_, b1_, a2_, b2_}, inf_] := 
   (sol = First[NDSolve[Evaluate[Join[eqns, ic0, 
                 {Derivative[2][a[0]][0] == a0, Derivative[2][a[1]][0] == a1, 
                   Derivative[2][b[1]][0] == b1, Derivative[2][a[2]][0] == a2, 
                   Derivative[2][b[2]][0] == b2}] /. {α -> 0.1}], 
           {a[0], a[1], b[1], a[2], b[2], Derivative[2][a[0]], 
             Derivative[2][a[1]], Derivative[2][b[1]], Derivative[2][a[2]], 
             Derivative[2][b[2]]}, {η, 0, inf}]]; 
     res = {Derivative[2][a[0]], Derivative[2][a[1]], Derivative[2][b[1]], 
           Derivative[2][a[2]], Derivative[2][b[2]]} /. sol; 
     (#1[inf] & ) /@ res)

By playing around I have found some good starting values and when I use these with infinity set to 9 everything works.

starts = {{a0, 0.054}, {a1, 0.704}, {b1, -0.702}, {a2, 0.021}, {b2, 
inf = 9;
rts = FindRoot[shoot[{a0, a1, b1, a2, b2}, inf] == {0, 0, 0, 0, 0}, 
shoot[{a0, a1, b1, a2, b2} /. rts, inf]

{-6.60792*10^-14, 5.00619*10^-13, 3.11087*10^-13, -3.76011*10^-13, 1.65523*10^-13}

If I plot the second derivative I can see that it goes to zero.

Join[{Plot[Evaluate[{a[0]''[η]} /. sol], {η, 0, inf}, 
   PlotRange -> All]},
  Plot[Evaluate[{a[n]''[η], b[n]''[η]} /. sol], {η, 0, 
    inf}, PlotRange -> All],
  {n, 1, 2}

Mathematica graphics

However my infinity is not big enough. If I increase it to say 10 or more I get

inf = 10;
rts = FindRoot[shoot[{a0, a1, b1, a2, b2}, inf] == {0, 0, 0, 0, 0}, 
shoot[{a0, a1, b1, a2, b2} /. rts, inf]

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

If one runs the integration on from the results with inf = 9 then one can see that the second derivative was made to pass through the point 9 but is not zero subsequently. Here I have run the inf = 9 case to get the roots but now solve and plot using these values for inf = 15.

inf = 15;
sol = First@NDSolve[Evaluate[Join[eqns, ic0,
       {(a[0]'')[0] == 
         a0, (a[1]'')[0] == 
         a1, (b[1]'')[0] == b1,
        (a[2]'')[0] == 
         a2, (b[2]'')[0] == b2}] /. 
      Join[{α -> 0.1}, rts]],
    {a[0], a[1], b[1], a[2], b[2], a[0]'', a[1]'', b[1]'', a[2]'', 
     b[2]''}, {η, 0, inf}];
Join[{Plot[Evaluate[{a[0]''[η]} /. sol], {η, 0, inf}, 
   PlotRange -> All]},
  Plot[Evaluate[{a[n]''[η], b[n]''[η]} /. sol], {η, 0, 
    inf}, PlotRange -> All],
  {n, 1, 2}

Mathematica graphics

So I have failed to get the asymptotic solution. I have tried integrating over an interval at large values of η and then using FindMinimum but this does not give good results.

Do you have any further ideas? Thanks


1 Answer 1


Once again, let's turn to the primitive finite difference method (FDM). I'll use pdetoae for the generation of difference equation:

var = Flatten@{a /@ Range[0, 2], b /@ Range@2}
bc = D[Through@var[t] == 0 // Thread, t, t] /. t -> inf

inf = 15;
points = 200;
difforder = 4;
grid = Array[# &, points, {0, inf}];
ptoafunc = pdetoae[var[η], grid, difforder];
del = #[[3 ;; -2]] &;
ae = del /@ ptoafunc@eqns;
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
aebc = ptoafunc@{ic0, bc};

sollst = Partition[
    With[{guess = 1}, 
      FindRoot[{ae, aebc} /. {α -> 1/10}, 
       Flatten[Table[{var[[i]][η], guess}, {i, 5}, {η, grid}], 1]]][[All, -1]],
     points]; // AbsoluteTiming


The following way to calculate sollst doesn't require one to use del:

fullsys = Flatten@ptoafunc@{eqns, ic0, bc};

lSSolve[obj_List, constr___, x_, opt : OptionsPattern[FindMinimum]] := 
 FindMinimum[{1/2 obj^2 // Total, constr}, x, opt]
lSSolve[obj_, rest__] := lSSolve[{obj}, rest]

sollst = Partition[
    With[{guess = 1}, 
      lSSolve[Subtract @@@ fullsys /. {α -> 1/10}, 
       Flatten[Table[{var[[i]][η], guess}, {i, 5}, {η, grid}], 1]]][[2, 
      All, -1]], points]; // AbsoluteTiming

But notice it's slower and less robut compared to the FindRoot approach.

solfunclst = ListInterpolation[#, grid] & /@ sollst;

Plot[D[Through@solfunclst[x], x, x] // Evaluate, {x, 0, inf}, PlotRange -> All]

Mathematica graphics

FindRoot still spits out lstol, but further check by adjusting points, difforder, etc. shows the result seems to be accurate enough.

  • $\begingroup$ This is wonderful. It is amazing that it does better than 'NDSolve'. I guess that it is the advantage of working on a grid. Many thanks. $\endgroup$
    – Hugh
    Nov 29, 2017 at 13:01
  • $\begingroup$ I have been looking at the output. I am most interested in the derivatives of the solution. I have been interpolating and then taking the derivative. Can the derivative of the solution be extracted from the method directly? The interpolation always adds some numerical noise. $\endgroup$
    – Hugh
    Nov 29, 2017 at 13:03
  • $\begingroup$ @hugh Actually tutorial/NDSolveBVP has mentioned that, "some initial value problems with growing modes are inherently unstable even though the BVP itself may be quite well posed and stable", and I think your equation system happens to belong to this type. As to the derivative, sadly it's not possible to extract them directly AFAIK. (You see, FDM only solves for the function value. ) $\endgroup$
    – xzczd
    Nov 29, 2017 at 13:11
  • $\begingroup$ I have just been back through SE and looked at all the references to pdetoae. It is very impressive. What are the limitations? I guess it can only be used on rectangular grids. However, the ability to solve using 'FindRoot' for nonlinear problems is outstanding. In trying to understand the method I think I see that you put the boundary conditions in as additional equations rather than altering the difference equations. Does this come to the same thing? How do you solve nonlinear problems on non-rectangular grids? Can you get something out of the Finite Element Method? $\endgroup$
    – Hugh
    Dec 1, 2017 at 10:15
  • $\begingroup$ @hugh 1. pdetoae is not limited to rectangular grid, it can handle all the regular domain (e.g. a circle) in principle. 2. I'm not quite sure about the meaning of "altering the difference equations", but the 2 methods should be equivalent if you mean "eliminating some of the variables from the difference equation system using b.c.", I prefer using the discretized b.c. as additional equations because it doesn't disarrange the variables so the rebuilding process will be easier. 3. Sadly I have no tool for nonlinear problem on irregular grid at the moment. $\endgroup$
    – xzczd
    Dec 1, 2017 at 10:56

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.