# NDSolve with removable singularity

I am trying to solve the following ODE with NDSolve $$xy''(x)+2y'(x)-xy(x)=x^2$$ with boundary conditions $$y(-1)=3 - 2e$$ and $$y(-0.5)=9/2-4\sqrt{e}$$. I want so solve it on the interval $$[-2,2]$$. The exact solution is given by $$y(x)=\frac{-2-x^2+2e^{-x}}{x}$$. The solution is not defined at $$x=0$$ where it has a removable singularity. I get the error (I guess because of that singularity):

NDSolve::ndsz: At x == -2.27*10^-322, step size is effectively zero; singularity or stiff system suspected.


Here is my code:

NDSolve[{x*D[f[x],{x,2}]+2*D[f[x],x]-x*f[x]==x^2, f[-0.5]==9/2-4*Sqrt[E], f[-1]==3-2*E},f, {x,-2,2}]


Is there a way to avoid that error such that NDSolve just "skips" the singularity and continues beyond?

• I get some extra errors. Are you aware of those?
– bmf
Mar 28, 2022 at 22:54
• @bmf Only warning on V13 due to numerical issue. But I rewrote his original code before trying it to break it down so I can see it more easily. Mar 28, 2022 at 23:03
• @Nasser many thanks. The re-writing you did is what I usually like seeing when solving a d.e
– bmf
Mar 28, 2022 at 23:04

Maybe this?:

ode = x*f''[x] + 2*f'[x] - x*f[x] == x^2;
ic = {f[-1/2] == 9/2 - 4*Sqrt[E], f[-1] == 3 - 2*E};
ysol = NDSolveValue[{ode, ic} /. f -> Function[x, u[x]/x],
Piecewise[{{u[x]/x, x != 0}}, u'[0]], {x, -1, 1}];

Plot[ysol, {x, -1, 1}]


Another way to avoid a singularity is with Piecewise, provided you can fill in the hole. In this case there seems to be numerical instability around the singularity. Using the automatic discontinuity processing resulted in large discontinuities in the solution at x == 0. Automatic adaptive step size couldn't integrate past zero. So here's an adaptive fixed-step approach, that results in a very small bump at x == 0.

(*ode2=f''[x]==Piecewise[{{x+f[x]-(2 f'[x])/x,x!=0}},1/3 (x+f[x])];*)
ode2 = f''[x] ==
PiecewiseExpand[
ode /. f'[x] :> Piecewise[{{f'[x], x != 0}}, x f''[x]] //
Inactivate[#, Equal] &] /.
eq : Inactive[Equal][__] :> First@SolveValues[Activate[eq], f''[x]];

run[step_] :=
NDSolveValue[{ode2, ic}, f, {x, -1, 1},
Method -> {"FixedStep", Method -> "ExplicitRungeKutta",
"DiscontinuityProcessing" -> False}, StartingStepSize -> step,
WorkingPrecision -> 32, MaxSteps -> Ceiling[1 + 2/step]];

{stepsize, ysol2, ysol3} = NestWhile[
{First[#]/2, Last[#], run[First@#/2]} &,
{1/512, run[1/256], run[1/512]},
NDSolveScaledVectorNorm[
2, {10^-8, 10^-8}][#[[2]][
"ValuesOnGrid"] - #[[3]]["ValuesOnGrid"][[;; ;; 2]], #[[2]][
"ValuesOnGrid"]] > 1 &,
1,
5
];

stepsize
NDSolveScaledVectorNorm[2, {10^-8, 10^-8}][
ysol3["ValuesOnGrid"][[;; ;; 2]] - ysol2["ValuesOnGrid"],
ysol2["ValuesOnGrid"]]
(*
1/512
0.0056963837440
*)

sol = -((E^-x (-2 + 2 E^x + E^x x^2))/x);
Plot[ysol3[x] - sol // RealExponent, {x, -1, 1},
PlotRange -> {-16.5, -8}, PlotLabel -> "Log abs. error"]


Update:

Actually, with high enough working precision, the default method of NDSolve succeeds with the piecewise ODE:

ysol4 = NDSolveValue[{ode2, ic}, f, {x, -1, 1},
Method -> {"DiscontinuityProcessing" -> False},
WorkingPrecision -> 50];
Plot[ysol4[x] - sol // RealExponent,
{x, -1, 1},
PlotRange -> {-26.5, 8}, WorkingPrecision -> 50,
PlotLabel -> "Log abs. error"]


Update 2

I've been wondering about (1) the exact boundary conditions and (2) the need for a numerical solution. Well, (2) is not that important to me. However, (1) keeps raising a question in my mind. The ODE is linear, so it's easy to understand first that at $$x=0$$, we must have $$f'(0)=0$$, if $$y=f(x)$$ does not blow up; second, that in the neighborhood of a initial condition $$f(0)=y_0,\, f'(0)=0$$ at $$x=0$$, almost all solutions go to infinity. So choosing a BVP not at $$x=0$$ and which stays finite as $$x$$ passes through $$0$$ is a delicate business. Any error at a step of a solver and the numerical solution almost certainly steps onto a solution that blows up. No wonder it seems nearly impossible to solve numerically. The bumps in the error we see above are from the solution starting to turn to infinity when $$x$$ switches from $$x<0$$ to $$x>0$$. Now the question is, why not solve the problem as an IVP with an initial condition at $$x=0$$? If we have the precise knowledge to set up the BVP, why can't we set up the corresponding IVP $$f(0)=-2,\,f'(0)=0$$? And even if we don't know the IVP, we can shoot for the boundary conditions from an IVP at $$x=0$$, although we have to do this by hand.

It's easy to figure out the starting point for the shooting method. Pick a not-stiff solver and solve the BVP, which will stop just before the singularity when the solution suddenly heads for infinity, near the desired IVP:

ic = {f[-(1/2)] == 9/2 - 4 Sqrt[E], f[-1] == 3 - 2 E};
NDSolveValue[{ode, ic},
f[-1.*^-4],  (* take a value to the left of x == 0 *)
{x, -1, 1},
Method -> "ExplicitRungeKutta"]


NDSolveValue::ndstf: At x == -0.0160629, system appears to be stiff. Methods Automatic, BDF, or StiffnessSwitching may be more appropriate.

InterpolatingFunction::dmval: Input value {-0.0001} lies outside the range of data in the interpolating function. Extrapolation will be used.

(*  -2.  <-- can verify with a plot that
it's a good starting value  *)


Now we can solve the BVP using the piecewise ODE ode2 from above. It works well at machine precision, but also at higher precision. (I ran it both at machine precision and WorkingPrecision -> 32). The result shows that the warning from FindMinimium is negligible. (There is no warning at WorkingPrecision -> 32.)

bcres = ic /. Equal -> Subtract;

obj // ClearAll;
obj[yy_?NumericQ] := Total[bcres /. f -> psol[yy]]^2;
FindMinimum[
obj[y], {y,
Quiet@NDSolveValue[{ode, ic}, f[-1.*^-4], {x, -1, 1},
Method -> "ExplicitRungeKutta"]}]
ysolMP = psol[y] /. Last[%];


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

(*  {1.3113*10^-17, {y -> -2.}}  <-- good residual despite warning *)

Plot[{ysolMP[x], ysol32[x]} - sol // RealExponent // Evaluate, {x, -1,
1}, PlotRange -> {-18.5, -6}, PlotLabel -> "Log abs. error",
PlotLegends -> {HoldForm[WP -> "MP"], HoldForm[WP -> "32"]},
WorkingPrecision -> 32]


• This is why numerical solvers are "messy" :) But they are needed for real life applications ofcourse. Mar 29, 2022 at 2:16

Is there a way to avoid that error such that NDSolve just "skips" the singularity and continues beyond?

Why not just solve it analytically?

ode = x*f''[x] + 2*f'[x] - x*f[x] == x^2
ic = {f[-1/2] == 9/2 - 4*Sqrt[E], f[-1] == 3 - 2*E}
sol= DSolveValue[{ode, ic}, f[x], x]


gives

And

Plot[sol, {x, -2, 2}]


Compared to

 sol = NDSolveValue[{ode, ic}, f[x], {x, -2, 2}]


 Plot[sol, {x, -2, 2}]


• Thanks, but as I wrote in my question I am aware of the easy analytical solution. I still would like to solve it numerically in order to tackle a similar problem where the analytical solution is more complicated. Mar 28, 2022 at 23:06
• @esner1994 perhaps StiffnessSwitching?
– bmf
Mar 28, 2022 at 23:07
• @esner1994 I tried first all the singularity avoidance tricks I know about but non-worked I am afraid. There are lots of questions about how to avoid singularity in numerical solving of ode on this site. some options work and some do not. Mar 28, 2022 at 23:15
• ok thanks. Worst case I divide the intervall and interpolate in between... Mar 28, 2022 at 23:19
• @esner1994 You can try WhenEvent and see if you can jump over the singularity. This might or might not work. I have not tried it. Mar 28, 2022 at 23:25