# Getting error NDSolve::ndsz from my ODE

I'm trying to figure out why I am consistently getting singularity error when I'm trying to solve an ODE using NDSolve. Any help will be appreciated.

EEr = 0.5*((1 + u'[#1])^2/(1 + om*con[#1])^(2/3) - 1) &;
EEt = 0.5*((1 + u[#1])^2/(1 + om*con[#1])^(2/3) - 1) &;

Pr = (1 + om*con[#1])^(1/3)*Ev*((1 + v)/v*EEr[#1] + EEt[#1])*(1 +u'[#1]) &;
Pt = (1 + om*con[#1])^(1/3)*Ev*(EEr[#1] + (1 + v)/v*EEt[#1])*(1 + u[#1]/#1) &;

EE = 80*10^9;
v = 0.25;
om = 5*10^-5;
Ev = EE*v/(1 + v)/(1 - 2 v);
R = 1;
con = (1 +1000 #1) &;
rs = $MachineEpsilon; equ = {Pr'[r] + (Pr[r] - Pt[r])/r == 0, u[rs] == rs, Pr[1] == 0}; sol = NDSolve[equ, {u}, {r, rs, 1}, Method -> {"Shooting", "StartingInitialConditions" -> {u[rs] == rs, u'[rs] == rs}}];  Evaluating the above code produces a host of errors of the form NDSolve::ndsz: At r == 2.220446049*10^-16 , step size is effectively zero; singularity or stiff system suspected To be clear, the equation itself is solvable, all I need is to avoid the singular point and any tips will be very appreciated! • Check what you posted here because I am getting a NDSolve::dvnoarg: The function u appears with no arguments. error. – Jose Enrique Calderon Nov 1 '18 at 9:20 • Thank you for noticing that, I've changed the code and hope this won't disturb you – Chengjun xu Nov 1 '18 at 10:00 • Is this Prandtl number Pr? – Alex Trounev Nov 1 '18 at 10:34 • The changes rs = Sqrt@$MachineEpsilon and Pr[1]/Ev == 0 get rid of the error. But maybe the solution is not good enough? – Michael E2 Nov 1 '18 at 11:11
• From the approximate solution, it looks like u'[r] goes to infinity at r == 0. A large error is somewhat expected (and even the error computed by NDSolve is approximate and probably subject to large errors near r == 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 than u[rs] == rs and possibly use a higher WorkingPrecision to reduce rounding error. (2) Find a locally exact solution at r == 0 and use NDSolve to integrate the difference. [Site tip: Use @MichaelE2 to make sure I'm notified of your reply.] – Michael E2 Nov 1 '18 at 12:51

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

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_] :=
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]}]  • To be honest, I am just in the first grade for my doctor degree and your code is a bit complicated to me. However, I will continue to work hard, and your professional attitude motivates me to discover a broader science world. Your words mean a lot to me and thanks again. @Michael E2 – Chengjun xu Nov 7 '18 at 3:08 Michael explained almost everything. I will add a few words. In this problem, the main source of errors is a large parameter Ev = 3.2 * 10 ^ 10 by which all parts of the equation are multiplied, but on which the solution does not depend. We reduce this parameter, and we solve the equation with high precision, for which we present all parameters in a rational form EEr = 1/2*((1 + u'[#1])^2/(1 + om*con[#1])^(2/3) - 1) &; EEt = 1/2*((1 + u[#1])^2/(1 + om*con[#1])^(2/3) - 1) &; Pr = (1 + om*con[#1])^(1/3)*((1 + v)/v*EEr[#1] + EEt[#1])*(1 + u'[#1]) &; Pt = (1 + om*con[#1])^(1/3)*(EEr[#1] + (1 + v)/v*EEt[#1])*(1 + u[#1]/#1) &; EE = 80*10^9; v = 1/4; om = 5*10^-5; Ev = EE*v/(1 + v)/(1 - 2 v); R = 1; con = (1 + 1000 #1) &; rs = Rationalize[Sqrt[$MachineEpsilon]];

equ = {Pr'[r] + (Pr[r] - Pt[r])/r == 0, u[rs] == rs, Pr[1] == 0};

sol = NDSolve[equ, {u}, {r, rs, 1},
Method -> {"Shooting",
"StartingInitialConditions" -> {u[rs] == rs,
u'[rs] == Rationalize[13.3709]}}, WorkingPrecision -> 30];
{Plot[Evaluate[u[r] /. sol], {r, rs, 1}, AxesLabel -> {"r", "u"}],
LogLogPlot[Evaluate[u'[r] /. sol], {r, rs, 1}, PlotRange -> All,
AxesLabel -> {"r", "u'"}]}


Calculate

Pr[1] /. sol

Out[]= {-4.9329723767604719*10^-14}


At first glance, the solution is obtained with high accuracy, but if we calculate

Pr[1]*Ev /. sol

Out[]= {-0.00157855116056335102}


then we will see that zero is still very far. Raise WorkingPrecision -> 50, then we calculate

Pr[1]*Ev /. sol

Out[]= {-5.370961405513687831138142980*10^-13}


This is close to zero and can be considered a solution to the problem.

• With the increased WorkingPrecision you can return to the OP's BC rs = Rationalize[$MachineEpsilon, 0]. Basically I expect that rs can be set to be about$10^{-wp/2}$. (I think the desired BC is that$u \rightarrow 0$as$r \rightarrow 0\$.) – Michael E2 Nov 1 '18 at 23:57
• The explanation is very clear and thoughful. l learn a lot from you and I’m vrey grateful. Finally, really thank you very much and hope this not didn't spend your time too much @Alex Trounev – Chengjun xu Nov 2 '18 at 4:10