# Nonlinear 2nd order ODE with regular singularities

I am tring to solve the following ODE with NDsolve

$2x~(1-x)~f''(x)+(3-4x)~f'(x)+a~f(x)+b~f^n(x)=0;~~a,b\in\mathbb{R},~n\in\mathbb{N}$.

The mathematica "code" is:

sol=NDSolve[{2*x (1 - x) *f''[x] + (3 - 4 x) f'[x] +a*f[x] -
b*(f[x])^n == 0, y[0.999999999] == 0,
y'[0.999999999] == 1}, {f[x]}, {x, 0, 0.999999999},
MaxSteps -> Infinity, Method -> "ExplicitRungeKutta",
Method -> {"StiffnessSwitching", "NonstiffTest" -> False}]

Plot[Evaluate[f[x] /. sol], {x, 0, 0.999999999}]


For values of $a,b,n$ and starting from an $x$ close to $1$ one can find solutions while solving towards $x=0$. I am wondering how correct/trustworthy the results are, considering the fact that at $x=0$ one has a regular singularity. Does anybody have an experience or intuition for such a problem?

I can't really understand why NDsolve does only work if I set the boundary to $x=0.999999999$ or closer to $1$ but not if I set $x=1$.

Finally, I would like to compute the $L^2$ norm of the solution with respect to a certain measure using NIntegrate leading to

NIntegrate[(Abs[f[x] /. sol])^2*Sqrt[x]/Sqrt[1 - x], {x, 0,
0.999999999}]


What I see is that depending on the precision set in the NDsolve the result changes significantly and varies over many orders of magnitude. I assume that this is because of the singularity so I am actually seeking for a very stable solution method which works nicely there.

I would be happy with any suggestion. Thanks!

• – Michael E2 Nov 2 '15 at 14:04

In my experience results from NDSolve are generally quite reliable. When there is doubt I would try to understand what the issue is by solving a smaller but similar problem, symbolically if possible. For the example you posed, here is what DSolve gives for specific values of a, b and n.

In:= dsol =
With[{a = 1, b = 1, n = 1, int = Rationalize[0.99999999999, 0]},
DSolve[{2*x (1 - x) f''[x] + (3 - 4 x) f'[x] + a f[x] -
b*(f[x])^n == 0, f[int] == 0, f'[int] == 1}, f[x], x]]

Out= {{f[x] -> -((
199999983448 (2 Sqrt Sqrt[1 - x] - Sqrt[x]))/(
9999998345000068475625 Sqrt[x]))}}


We see there is trouble when x is close to 0.

NDSolve also gives the same thing, only numerically. Try with increased precision and plot the results from DSolve and NDSolve to compare:

sol = With[{a = 1, b = 1, n = 1, int = Rationalize[0.99999999999, 0]},
NDSolve[{2*x (1 - x) f''[x] + (3 - 4 x) f'[x] + a f[x] -
b*(f[x])^n == 0, f[int] == 0, f'[int] == 1}, f[x], {x, 0, int},
MaxSteps -> Infinity, Method -> "ExplicitRungeKutta",
Method -> {"StiffnessSwitching", "NonstiffTest" -> False},
WorkingPrecision -> 20]]

Plot[Evaluate[f[x] /. {sol, dsol}], {x, 0, 0.999999999},
PlotStyle -> {{Red, Dashed}, {Blue, Dotted}}]

• Hi. Thanks! But what about the issue with the 0.999999999 and not 1? – Hamurabi Nov 4 '15 at 12:04
• I believe it is because of the (1-x) term in the coefficient for f''[x], which is the highest derivative term. Try changing it to (2-x) and see what happens. You may have to change the NDSolve limits also to say {x,0,3 int}. You should see messages pertaining to x = 0 and x = 2. – Lotus Nov 4 '15 at 12:48