# NDSolve: Infinite expression 1/0 encountered

I am trying to solve the following differential equation

EOM=(-Piecewise[{{1.157793814828612*^106, r >= 0 && r <= 1.0408060341965232*^-18}, {(9.734399999999999*^36*(2.7912782245116803*^-35 - 3*r^2))/
(2.7912782245116803*^-35 + r^2)^3, r > 1.0408060341965232*^-18}}, 0])*f[r] + (7*Derivative[1][f][r])/r + Derivative[2][f][r] == 0


with initial conditions

f[0] == 1
Derivative[1][f][0] == 0


numerically using NDSolve

NDSolve[{EOM, {f[0] == 1, Derivative[1][f][0] == 0}}, {f[r]}, {r, 0, 100}]


But I get

Power::infy: Infinite expression 1/0 encountered.

Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered.

NDSolve::ndnum: Encountered non-numerical value for a derivative at r == 0.

I searched other related posts but could not solve the problem. I will be thankful if you help me with this.

• Your EOM contains a term (7*f'[r])/r which results in Infinite expression 1/0 encountered at the initial conditions (r = 0). Jun 24 at 14:00
• Thank you for your reply, I actually tried multiplying the whole equation by r but that did not help and I still get the same error. @Bob Hanlon
– amy
Jun 24 at 15:14
• A little besides the point, but the numbers you have seem a bit wacky - they span 141 orders of magnitude which will likely create some numerical issues for you. Are you sure they're correct? Jun 24 at 16:41
• Change the initial conditions to {f[epsilon] == 1, f'[epsilon] == 0} for a sufficiently small epsilon Jun 24 at 17:18
• Yes, they are correct. The function pieces are of the same order of magnitude and also all the terms added up are of the same order of magnitude. So, in principle the equation is Okay. @bRost03
– amy
Jun 24 at 17:23

This may be an answer to an approximation to your question. I am making assumptions.

In addition to the problem at r=0, the derivative of your piecewise function is zero when r=0. The this gives you a solution:

DSolve[{(7 Derivative[1][f][r])/r + (f^\[Prime]\[Prime])[r] == 0,
f'[0] == 0, f[0] == 1}, f, r]


and so the integrator isn't going to behave well near zero.

It look to me like you are approximating a potential that looks like this:

pot = a/r^5 - b/r^3


which gives you a force that looks like this.

force = -Simplify[D[pot, r]]


which looks like your piecewise function. Please forgive me for making assumptions about what you are doing.

Continuing..

This problem is begging to be rescaled to non-dimensional units. Let's suppose the mininum potential (-eMin) is located at rMin

solNorm =
Solve[{rMin == SolveValues[force == 0, r][[2]],
eMin == pot /. r -> rMin}, {a, b}][[1]]
eNorm = Simplify[((pot /. r -> rho rMin)/eMin) /. solNorm]


That gives a normalized force term (to get units, multiply by Emin/rMin)

fNorm = - D[eNorm, rho]


With the assumption, your equation of motion becomes

eom = phi[rho] fNorm + 7 D[phi[rho], rho]/rho + D[phi[rho], {rho, 2}]


Note that starting the equation of motion at pho=0 is going to make phi fly off to infinity. Using @Bob Hanlon's advice, here is a function to get a numerical solution for different starting points.

phiSol[sp_] :=
phiSol[sp] =
NDSolveValue[{eom == 0, phi'[sp] == 0, phi[sp] == 1},
phi, {rho, sp, 2}]


For example:

phiSol[.1]


We will need to extract the domain of the InterpolationFunction

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];


The behavior can be investigated for different initial points:

Manipulate[
With[{interpolator = phiSol[sp],
domain = Flatten[InterpolatingFunctionDomain[phiSol[sp]]]},
Plot[interpolator[t], {t, domain[[1]], domain[[2]]}, Frame -> True,
FrameLabel -> {"r/rMin", "phi"}]
],
{{sp, .5}, .01, 1}
]


Note that the length unit is rMin which is very small in your original problem.