Integrating second order differential equation

Questions
, ,
1

Hi everyone, I’m a new user of Mathematica, and i’m having difficulties in solving a particular problem. My aim is to integrate this differential equation:

\frac{d^{2 }U}{dx^{2}} + \frac{2}{x} \frac{dU}{dx} = - \frac{I(U)}{I(U_{0})} - \lambda \frac{I(\mu U)}{I(\mu U_{0})}

with boundary conditions:

U(0)=U_{0} \; \frac{dU}{dx}=0 \; \mathrm{at} \;x=0

where I(U) is defined as:

I(U)= \frac{3}{4} \sqrt{\pi} -\frac{2}{3} e^{U} \Gamma(\frac{5}{2}, U)

where \Gamma(a,U) is the incomplete gamma function.
In Mathematica I tried this code, with given values for the parameters above:

f[x_]:= 3/4 Sqrt[Pi] - 2/3 Exp[x]* Gamma[5/2, x];
g[x_, p_, l_, m_]:= (f[x]/f[p]) + l* (f[m*x]/f[m*p]);

DSolve[{y''[x] + (2/x) y'[x] + g[y[x], 2, 2, 3]==0, y[0]==2, y'[0]==0}, y, {x,0, 1}]

but i get the following error: the function y appears with no arguments. Could you please give me any help? Please, also consider that my final aim is then to plot g(U(x)) as defined at different values of U_{0}, \lambda and \mu. Thank you.

2

Hello,
I believe that you may have previous definitions of f,g, or y.
When I try your code from a fresh kernel, I get the DSolve back, indicating that Mathematica failed to find a closed form solution.

Try:

Clear[f,g,y];
DSolve[{y''[x] + (2/x) y'[x] + g[y[x], 2, 2, 3]==0, y[0]==2, y'[0]==0}, y, {x,0, 1}]

I believe you may have to resort to NSolve:

Clear[y, f, g];
f[x_] := 3/4 Sqrt[Pi] - 2/3 Exp[x]*Gamma[5/2, x];
g[x_, p_?NumericQ, l_?NumericQ, m_?NumericQ] := (f[x]/f[p]) + 
   l*(f[m*x]/f[m*p]);
nsol = y /. 
  NDSolve[{y''[x] + (2/x) y'[x] + g[y[x], 2, 2, 3] == 0, y[1] == 2, 
     y'[1] == 0}, y, {x, 1, 2}][[1]]
Plot[nsol[t], {t, 1, 2}]

I changed the domain of x because there is a singularity at x=0.
You could generalize nsol to take arguments for numerical values of p,l,and m.

Note the use of the ?NumericQ pattern restriction is because NSolve only plays nice when it calls functions that are defined for numerical values.

3

Thank you very much for you help, I tried your suggestions. However, my specific problem requests conditions for x=0, is there a way to do this avoiding the singularity?

4

Hello,
There is probably a clever way to do this by looking at the limits as t_{init}\rightarrow 0, but this will give you a visual clue about the behavior:

Clear[nsol]

nsol[init_ ] := 
 y /. NDSolve[{y''[x] + (2/x) y'[x] + g[y[x], 2, 2, 3] == 0, 
     y[init] == 2, y'[init] == 0}, y, {x, init, 2}][[1]]

Manipulate[
 Plot[nsol[init][t], {t, init, 2}, PlotRange -> {{0, 2}, All}], {init,
   0.0001, .1}]

or,

nsol2[eps_ ] := 
 y /. NDSolve[{ y''[x] + (2/(x + eps)) y'[x] + g[y[x], 2, 2, 3] == 0, 
     y[0] == 2, y'[0] == 0}, y, {x, 0, 2}][[1]]

Manipulate[
 Plot[nsol2[eps][t], {t, 0, 2}, 
  PlotRange -> {{0, 2}, All}], {{eps, 0.01}, .0001, .1}]