3

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
  • 219
  • 1
  • 4
  • 11
nightmarish
  • 141
  • 4
  • 3
    The oscillations are likely due to the solution at large $r$ being extremely sensitive to the initial conditions at $r_0$, meaning that very small errors at in $\Phi(r_0), \Phi'(r_0)$ propagate into very large errors in $\Phi(r)$ as $r\to\infty$. The shooting method is not very suitable for such unstable problems. Try printing plots of solutions for intermediate IC guesses to see how much the solution changes during the binary search, it may be helpful. – Kirill Aug 23 '15 at 20:01
  • Also, how does the bisection work with the test $\Phi_{\mathrm{test}}\in(0,1.2)$? – Kirill Aug 23 '15 at 22:41
  • I've seen that, as $r \rightarrow \infty$, the solution either tends to $\Phi \rightarrow \infty, -\infty,$ or $(3 + \sqrt{9 - 8 \alpha})/(2 \alpha)$ (which is within 0.8 to 1.2 for the values of $\alpha$ I used). The last value of $\Phi$ was found by solving $\Phi - \frac{3}{2}\Phi^{2}+\frac{\alpha}{2}\Phi^{3}=0$, gotten from the differential equation itself. Since my required solution should tend to $0$ as $r \rightarrow \infty$, I know that my required solution is between the solution which shows oscillations about $\Phi \approx 1.2$ and the solution which goes to $-\infty$. – nightmarish Aug 24 '15 at 10:02
  • In fact, if we solve the equation $\Phi - \frac{3}{2}\Phi^{2}+\frac{\alpha}{2}\Phi^{3}=0$, we get $\Phi = 0$, $\frac{1.5 \pm \sqrt{2.25 - 2\alpha}}{\alpha}$. If we use $\Phi(0)$ to be either one of these values, then the solution is an equilibrium solution - a horizontal line. If we choose $\Phi(0)<0$, then the solution diverges to $-\infty$. If we choose $\Phi(0)>\frac{1.5 + \sqrt{2.25 - 2\alpha}}{\alpha}$, then the solution diverges to $\infty$. For all other choices of $\Phi(0)$, the solution either oscillates about $\frac{1.5 - \sqrt{2.25 - 2\alpha}}{\alpha}$, or diverges to $-\infty$. – nightmarish Aug 24 '15 at 10:16
  • In the meantime, I'll do multiple shooting and get back to you. – nightmarish Aug 24 '15 at 10:23

0 Answers0