1
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Nov 2 '15 at 14:02
  • $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this this meta Q&A helpful $\endgroup$ – Michael E2 Nov 2 '15 at 14:03
  • $\begingroup$ Related: mathematica.stackexchange.com/questions/96799/…, and perhaps mathematica.stackexchange.com/questions/91854/… $\endgroup$ – Michael E2 Nov 2 '15 at 14:04
2
$\begingroup$

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}}]
$\endgroup$
  • $\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.