I am tring to solve the following ODE with NDsolve


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, 

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!


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}}]
  • $\begingroup$ Hi. Thanks! But what about the issue with the 0.999999999 and not 1? $\endgroup$ – Hamurabi Nov 4 '15 at 12:04
  • $\begingroup$ 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. $\endgroup$ – Lotus Nov 4 '15 at 12:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.