Here are the issues, including ones implied in my comment:
- Increasing round-off error as $r \rightarrow 0$.
- Increasing stiffness as $r \rightarrow 0$.
- Poor scaling due to large coefficients.
- Starting integration at the stiff boundary $r \approx 0$.
In my comments I suggested not getting too close to r == 0
, rescaling the BC Pr[1]/Ev == 0
, and possible increasing WorkingPrecision
would help solve the BVP. Here I'll show another approach by starting the integration at the stabler boundary condition.
Generally, one should expect it to be easier to start the integration at the boundary where there is no singularity (or numerics issues arising from one). The trouble is finding good starting initial conditions. One approach to BVP is to use to use FindRoot
to find the Chebyshev series of a function that satisfies the BVP at the Chebyshev collocation points. The BCs are easy to implement, even the one at the singularity r == 0
. One the one hand, this is not a particularly good example for the approach because of the singularity and numeric issue; other the other hand, with just 33 points, we quickly get a decent approximation that yields good starting points for u[1]
and u'[1]
. To implement this approach I modified J.M.'s chebdeval[] to return, the function, derivative, and second derivative values. (It also returns the input x
to make the code in residual[]
simpler. A bit of a kludge, I guess.)
Another technique I use is to use Solve
to rearrange an equation, so that it is better scaled. For instance, Solve[ODE, u''[r]]
yields a result that can be put in the form u''[r] == RHS
such that the error u''[r] - RHS
is scaled such that the difference can be compared to u''[r]
. which is the value NDSolve
is computing from the ODE
. Since NDSolve
does this internally, it is not important for it, but it turns out to be helpful for FindRoot
. Similarly we can rescale the BC at r == 1
with Solve[Pr[1] == 0, u[1]]
, although because of the nonlinear nature of Pr[1]
, it is more convenient here to solve for (1 + u[1])^2
. NDSolve
does not do this, so we have to do it in the call to NDSolve[]
.
Utility definitions and functions:
Clear[equ];
equ[] := Pr'[r] + (Pr[r] - Pt[r])/r == 0; (* Basic ODE *)
equ[rs_, re_: 1] := { (* BVP *)
equ[],
u[rs] == Sqrt@rs, (* BC near r == 0 *)
(1 + u[1])^2 == C[0] /. First@Solve[ (* Effectively rescales BC at r == 1 *)
Rationalize[Rationalize@Pr[1] == 0, 0] /. (1 + u[1])^2 -> C[0],
C[0]]};
reduced = (Subtract @@ First@First@Solve[equ[], u''[r]]); (* Rescales ODE *)
error = reduced/(1*^-8 + Abs[u''[r]]); (* err rel to u''[r] when u''[r] > 10^-8 *)
(* computes {F, F', F'', x} of a Chebyshev series for F over {a,b} *)
Clear[chebdeval];
chebdeval[cof_List, x_, a_: 0, b_: 1] :=
Module[{f0, f1, f2, j, p0, p1, p2, q0, q1, q2, t, dtdx},
t = (2 x - a - b)/(-a + b); dtdx = 2/(b - a);
f2 = f1 = 0; (* f = F *)
p2 = p1 = 0; (* p = F' *)
q2 = q1 = 0; (* q = F'' *)
Do[
q0 = 2 (t q1 + 2 p1) - q2;
q2 = q1; q1 = q0;
p0 = 2 (t p1 + f1) - p2;
p2 = p1; p1 = p0;
f0 = 2 t f1 - f2 + cof[[j]];
f2 = f1; f1 = f0
, {j, Length[cof], 2, -1}];
{t f1 - f2 + First[cof], dtdx (t p1 + f1 - p2), dtdx^2 (t q1 + 2 p1 - q2), x}
];
(* adds BC u[0]==0 to coefficient list add prepending the zero-order
Chebyshev coefficient c[0] to coefficient array c[1], c[2],... *)
addbc = Prepend[#, ((-1)^Range[0, Length@# - 1]).#] &;
(* computes residual errors of the ode at collocation nodes *)
Clear[residual];
residual[c_?(VectorQ[#, NumericQ] &), nodes_] :=
With[{coeff = addbc[c]}, (* add BC at r == 0 *)
Append[
reduced /.
Thread[{u[r], u'[r], u''[r], r} -> chebdeval[coeff, nodes[[2 ;; -2]]]],
Pr[1]/Ev /. (* add BC at r == 1 *)
Thread[{u[1], u'[1], u''[1], r} -> chebdeval[coeff, nodes[[-1]]]]
]
];
Carrying out the solution:
(* construct Chebyshev solution *)
Block[{n = 32},
nodes = Sin[Pi/2 Range[-N@n, n, 2]/n];
cc32 = coeff /. FindRoot[ (* Find Chebyshev series that satisfies ODE at nodes *)
residual[coeff, (nodes + 1)/2] == ConstantArray[0, Length@nodes - 1],
{coeff, ConstantArray[0, Length@nodes - 1]},
PrecisionGoal -> 10];
]
(* Use Chebyshev sol. for "StartingInitialConditions" at r == 1 for NDSolve *)
PrintTemporary@Dynamic@{Clock[Infinity]};
Block[{rs = $MachineEpsilon, wp = 32},
rs = SetPrecision[rs, wp];
sol = First@NDSolve[
Rationalize[Rationalize@equ[rs], 0], {u}, {r, rs, 1},
Method -> {"Shooting",
"StartingInitialConditions" ->
{{u[1], u'[1]} == Rationalize[chebdeval[addbc@cc32, 1][[;; 2]], 0]},
Method -> Automatic},
WorkingPrecision -> wp, PrecisionGoal -> 10,
InterpolationOrder -> All]
] // AbsoluteTiming
Error and visualizations:
NIntegrate @@ {
SetPrecision[error, 24]^2 /. sol,
Flatten@
{r, Sort[Join[
u@"Grid" /. sol,
{{r}} /. FindRoot[u''[r] /. sol, {r, 0.22}, WorkingPrecision -> 24]]]},
MaxRecursion -> 20, PrecisionGoal -> 8, AccuracyGoal -> 16,
WorkingPrecision -> 24} // Sqrt
(* 9.43592806141403019272207*10^-7 *)
Plot @@ {error /. sol // RealExponent, Flatten@{r, u@"Domain" /. sol},
PlotPoints -> 200, AspectRatio -> 1/4}

Plot[{First@chebdeval[addbc@cc32, x], u[x] /. sol}, {x, 0, 1},
PlotStyle -> {AbsoluteThickness[6], AbsoluteThickness[3]}]

rs = Sqrt@$MachineEpsilon
andPr[1]/Ev == 0
get rid of the error. But maybe the solution is not good enough? $\endgroup$u'[r]
goes to infinity atr == 0
. A large error is somewhat expected (and even the error computed byNDSolve
is approximate and probably subject to large errors nearr == 0
as well). I don't have time to investigate, but there are two avenues that occur to me: (1) Come up with a better BC thanu[rs] == rs
and possibly use a higherWorkingPrecision
to reduce rounding error. (2) Find a locally exact solution atr == 0
and useNDSolve
to integrate the difference. [Site tip: Use @MichaelE2 to make sure I'm notified of your reply.] $\endgroup$