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]][η] -
b[2][η]*Derivative[2][b[2]][η]),
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]][η] +
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][η]*
Derivative[2][b[2]][η]),
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 +
2*α*a[2][η]*Derivative[2][a[0]][
η] + α*a[1][η]*Derivative[2][a[1]][η] +
2*α*a[0][η]*Derivative[2][a[2]][η] -
α*b[1][η]*Derivative[2][b[1]][η]),
-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.
ClearAll[shoot];
shoot[{(a0_)?NumberQ, a1_, b1_, a2_, b2_}, inf_] :=
(sol = First[NDSolve[Evaluate[Join[eqns, ic0,
{DerivatClearAll[shoot];
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,
0.020}};
inf = 9;
rts = FindRoot[shoot[{a0, a1, b1, a2, b2}, inf] == {0, 0, 0, 0, 0},
starts];
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]},
Table[
Plot[Evaluate[{a[n]''[η], b[n]''[η]} /. sol], {η, 0,
inf}, PlotRange -> All],
{n, 1, 2}
]
]
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},
starts];
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]},
Table[
Plot[Evaluate[{a[n]''[η], b[n]''[η]} /. sol], {η, 0,
inf}, PlotRange -> All],
{n, 1, 2}
]
]
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
