# Normalizing singularities in NDSolve

I've tried to create the following example. Suppose that I have the differential equation:

$$U''(x) = \frac{U'(x) - U(x)}{x},$$

which I know 1 boundary condition and know that this function should behave normal at the origin. Thus, I would assume that to keep the function normal it should be true that:

$$U'(0) = U(0),$$

with say my 1 boundary condition is:

$$U(1) = 10,$$

or something of the sort. I would have then expected I could solve this differential equation with:

soln = NDSolve[{U''[x] == (U'[x] - U[x])/x, U'[0] == U[0],
U[10] == 10},
U,
{x, 0, 1},
Method -> {"Shooting",
"StartingInitialConditions" -> {U[0] == 1}}]


however, this immediately gives an infinite expression. Looking at the trace it seems that it does not consider that the numerator is $$+1-1$$ at all:

(-1. + 1.)/0.
Message[Power::infy, 1/0.]


Any thoughts on how to analyze such a differential equation? Any response is greatly appreciated.

Since the deq is singular at zero, the easiest workaround is to start the solution a little above zero.

Assign an small value

\[Epsilon] = 10^-10

soln = NDSolve[{U''[x] == (U'[x] - U[x])/x, U[\[Epsilon]] == 1,
U[1] == 10}, U[x], {x, \[Epsilon], 1},
Method -> {"Shooting",
"StartingInitialConditions" -> {U'[\[Epsilon]] == 1}},
WorkingPrecision -> 25] // Flatten


I raised the working precision a little to get rid of some warnings about precision.

U[x_] = U[x] /. soln


This particular problem can also be solved analytically:

sol = DSolve[{deq, u'[0] == u[0], u[0] == 1, u[1] == 10},
u[x], {x, 0, 1}] // Flatten

u[x_] = u[x] /. sol


Check

Limit[u'[x] - u[x], x -> 0]
(*  0  *)

u[1]
(*  10  *)

Limit[u[x], x -> 0]
(*  1  *)


Mathematica gets infinity for x = 0, but taking the limit yields correct values.

Plot the numeric and symbolic solutions.

Plot[{U[x], u[x]}, {x, 0, 1}]


And they overlay.

The ODE may be solved symbolically with

uFN = DSolveValue[{U''[x] == (U'[x] - U[x])/x}, U, {x, 0, 2}]

(* Function[{x}, 2 x (BesselJ[2, 2 Sqrt[x]] C[1] - BesselY[2, 2 Sqrt[x]] C[2])]  *)


The singularity enforces the desired "boundary condition" at x == 0 on the whole solution space:

Limit[uFN'[x] - uFN[x], x -> 0, Direction -> "FromAbove"]
(*  0  *)


The BC at x == 10 reduces the dimension of the solution space by 1:

uFN10 = DSolveValue[{U''[x] == (U'[x] - U[x])/x, U[1] == 10},
U, {x, 0, 2}]
(*
Function[{x},
2 x (BesselJ[2, 2 Sqrt[x]] C[1] +
(BesselY[2, 2 Sqrt[x]] (5 - BesselJ[2, 2] C[1]))/BesselY[2, 2])]
*)

Plot[Table[uFN10[x] /. C[1] -> c, {c, 0, 10, 2}] // Evaluate, {x, 0, 1.1}]


• BTW, to see the family of solutions with NDSolve, just use an IVP at U[10] == 10 (or u[1] == 10 as above), with a varying value for U'[10]. This make an nice, suggestive image: ListLinePlot@Table[NDSolveValue[{U''[x] == (U'[x] - U[x])/x, U'[10] == m, U[10] == 10}, U, {x, 0 + \$MachineEpsilon, 25}], {m, -3, 3}] Dec 7 '21 at 2:59

Start NDSolve at a slightly negative x with m as parameter for U'[1] == m , and you have no problems at x == 0. Determine the corresponding m at x==0 for U[m][0]== "anyValue".

As @MichaelE2 says, "The singularity enforces the desired "boundary condition" at x == 0"

sol = Solve[U''[x] == (U'[x] - U[x])/x, U[x]]
sol /. x -> 0
(*   {{U[0] -> Derivative[1][U][0]}}   *)

Usol[m_?NumericQ] :=
U /. First@
NDSolve[{U''[x] == (U'[x] - U[x])/x, U'[1] == m, U[1] == 10},
U, {x, -10^-14, 3/2}]

Plot[Evaluate[Table[Usol[m][x], {m, -50, 50, 11/3}]], {x, 0, 1.5}]


Max@Abs@Table[Usol[m][0] - Usol[m]'[0], {m, -50, 50, 11/3}]
(*   2.30926*10^-13   *)

m1 = m /. FindRoot[Evaluate[Usol[m][0] == 1], {m, 3}]
(*   13.5113   *)

m1000 = m /. FindRoot[Evaluate[Usol[m][0] == 1000], {m, 3}]
(*   -2817.85   *)

Plot[Evaluate[Usol[m1000][x]], {x, 0, 1}]