Nonlinear 2nd order ODE with regular singularities

Question

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!

Answer 1

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[35]:= 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[35]= {{f[x] -> -((
    199999983448 (2 Sqrt[24999997931] 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}}]