I am trying to solve the following nonlinear ODE, but Mathematica takes forever to give a result. Any idea why?

eq = f”[r] + r^-1 f'[r] – r^-2 f[r] + (1 – r^2) f[r] == 0

g[r_] = f[r] /. NDSolve[{eq, f[0] == 1, f[1] == 1}, f[r], {r, 0, 1}] //First;

Plot[g[r], {r, 0, 1}]

=================

1

I recommend to use DSolve without boundary conditions. Then try to expand solutions with Series eg. Series[sol,{r,0,3}]. Then it should appear resonable what might go wrong (branch cuts of complex functions). There is also option in NDSolve to avoid symbolic preprocessing.

– Artes

Jan 19 at 17:04

Your code throws some warnings, but runs on my computer in 1 second.

– Chris K

Jan 19 at 19:43

=================

2 Answers

2

=================

The problem here is that your equation has terms singular at r=0, while you fix one of the boundary conditions exactly in this point. The workaround may be in shifting the first boundary condition a bit out of this point:

ClearAll[r];

eq = f”[r] + f'[r]/r – f[r]/r^2 + (1 – r^2) f[r] == 0;

nds = NDSolve[{eq, f[0.1] == 1, f[1] == 1}, f[r], {r, 0, 1}] // First

Plot[f[r] /. nds, {r, 0, 1}]

yielding this:

Alternatively you can regularize the equation by dividing terms by r+0.001 instead of dividing by r:

ClearAll[r];

eq = f”[r] + f'[r]/(r + 0.001) –

f[r]/(r + 0.001)^2 + (1 – r^2) f[r] == 0;

nds = NDSolve[{eq, f[0] == 1, f[1] == 1}, f[r], {r, 0, 1}] // First

Plot[f[r] /. nds, {r, 0, 1}]

returning this:

It is up to you to decide which approach is better.

Have fun!

Because something is wrong with the boundary at x=0x=0. Let’s solve the equation with DSolve. Since DSolve has trouble in dealing with f[0] == 1, we first use a more general b.c. f[a] == b:

eq = f”[r] + r^-1 f'[r] – r^-2 f[r] + (1 – r^2) f[r] == 0

{generalasol} = f[r] /. DSolve[{eq, f[a] == b, f[1] == 1}, f[r], {r}];

and take the limit:

asol = Limit[asolgeneral, a -> 0]

(E^(1/2 – r^2/2) LaguerreL[1/4, -1, r^2])/(r LaguerreL[1/4, -1, 1])

Somewhat surprisingly, we found the solution doesn’t depend on the value of b at all, and the resulting graph is:

Plot[asol, {r, 0, 1}]

It’s actually in agreement with the numeric result given by Alexei. If you’ve tried different shift for the left boundary, you would have found the offset drastically influence the result:

Manipulate[Plot[

f[r] /. NDSolve[{eq, f[e] == 1, f[1] == 1}, f[r], {r, e, 1}] // First // Evaluate, {r,

e, 1}, PlotRange -> {{0, 1}, {0, 1}}], {e, 10^-4, 0.1}]