2

I have been trying to solve the following nonlinear ordinary differential equation:

$$-\Phi''-\frac{3}{r}\Phi'+\Phi-\frac{3}{2}\Phi^{2}+\frac{\alpha}{2}\Phi^{3}=0$$

with boundary conditions

$$\Phi'(0)=0,\Phi(\infty)=0.$$

My solution is supposed to reproduce the following plots:

Now, to produce the plots given above, I wrote the following Mathematica code:

α = 0.99;
Φlower = 0;
Φupper = 5;
For[counter = 0, counter <= 198, counter++,
    Φ0 = (Φlower + Φupper)/2;
    r0 = 0.00001;
    Φr0 = Φ0 + (1/16) (r0^2) (2 (Φ0) - 3 (Φ0^2) + α (Φ0^3));
    Φpr0 = (1/8) (r0) (2 (Φ0) - 3 (Φ0^2) + α (Φ0^3));
    diffeq = {-Φ''[r] - (3/r) Φ'[r] + Φ[r] - (3/2) (Φ[r]^2) + (α/2) (Φ[r]^3) == 0, Φ[r0] == Φr0, Φ'[r0] == Φpr0};
    sol = NDSolve[diffeq, Φ, {r, r0, 200}, Method -> "ExplicitRungeKutta"];
    Φtest = Φ[200] /. sol[[1]];
    Φupper = If[(Φtest < 0) ||   (Φtest > 1.2), Φ0, Φupper];
    Φlower = If[(Φtest < 1.2) && (Φtest > 0),   Φ0, Φlower];
]
Plot[Evaluate[{Φ[r]} /. sol[[1]]], {r, 0, 200}, 
     PlotRange -> All, PlotStyle -> Automatic]

In the code, I used Taylor expansion at $r=0$ due to the $-\frac{3}{r}\Phi'$ term. Moreover, I used shooting method and continually bisected an initial interval from $\Phi_{\text{upper}}=5$ to $\Phi_{\text{lower}}=0$ to obtain more and more precise values of $\Phi(0)$.

With the code above, I was able to produce the plots for $\alpha = 0.50, 0.90, 0.95,0.96,0.97$. For example, my plot for $\alpha = 0.50$ is as follows:

However, my plot for $\alpha = 0.99$ does not converge to the required plot:

Can you suggest how I might tackle this problem for $\alpha = 0.99$? Also, is there an explanation for the plots shooting upwards and oscillating after a prolonged asymptotic trend towards the positive $r$-axis?

Glorfindel
  • 547
  • 1
  • 8
  • 14
nightmarish
  • 283
  • 1
  • 7
  • Thank you for your comment! Now that you mention it, I think my question is not really about the technical details of Mathematica, but rather about the implementation technique of an ODE problem. So, would you suggest me to delete it? – nightmarish Aug 23 '15 at 19:43
  • In any event: what happens if you use a different method, like "StiffnessSwitching"? – J. M.'s missing motivation Aug 24 '15 at 01:38
  • I would expect Φ[0] to be (3 + Sqrt[9 - 8 α])/(2 α). – bbgodfrey Aug 24 '15 at 04:57
  • @ Guess who it is: I tried "Stiffness Switching", bit it does not work for some reason. – nightmarish Aug 24 '15 at 09:50
  • @ bbgodfrey: I understand that you calculated $\Phi[0]=(3 + \sqrt{9 - 8 \alpha})/(2 \alpha)$ using the cubic equation in the differential equation. But, that value of $\Phi[0]$ actually leads to an equilibrium solution - the solution does not tend to $0$ as $r$ tends to $\infty$. – nightmarish Aug 24 '15 at 09:52
  • 5
    The identical equation is solved in question 84131 – bbgodfrey Aug 24 '15 at 18:26

0 Answers0