Solving an ODE using shooting method

1. Aug 23, 2015

spaghetti3451

1. The problem statement, all variables and given/known data

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.$

2. Relevant equations

My solution is supposed to reproduce the following plots:

3. The attempt at a solution

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

\[Alpha] = 0.99;
\[CapitalPhi]lower = 0;
\[CapitalPhi]upper = 5;
For[counter = 0, counter <= 198, counter++,
\[CapitalPhi]0 = (\[CapitalPhi]lower + \[CapitalPhi]upper)/2;
r0 = 0.00001;
\[CapitalPhi]r0 = \[CapitalPhi]0 + (1/16) (r0^2) (2 (\[CapitalPhi]0) - 3 (\[CapitalPhi]0^2) + \[Alpha] (\[CapitalPhi]0^3));
\[CapitalPhi]pr0 = (1/8) (r0) (2 (\[CapitalPhi]0) - 3 (\[CapitalPhi]0^2) + \[Alpha] (\[CapitalPhi]0^3));
diffeq = {-\[CapitalPhi]''[r] - (3/r) \[CapitalPhi]'[r] + \[CapitalPhi][r] - (3/2) (\[CapitalPhi][r]^2) + (\[Alpha]/2) (\[CapitalPhi][r]^3) == 0, \[CapitalPhi][r0] == \[CapitalPhi]r0, \[CapitalPhi]'[r0] == \[CapitalPhi]pr0};
sol = NDSolve[diffeq, \[CapitalPhi], {r, r0, 200}, Method -> "ExplicitRungeKutta"];
\[CapitalPhi]test = \[CapitalPhi][200] /. sol[[1]];
\[CapitalPhi]upper = If[(\[CapitalPhi]test < 0) || (\[CapitalPhi]test >
1.2), \[CapitalPhi]0, \[CapitalPhi]upper];
\[CapitalPhi]lower = If[(\[CapitalPhi]test < 1.2) && (\[CapitalPhi]test > 0), \[CapitalPhi]0, \[CapitalPhi]lower];
]
Plot[Evaluate[{\[CapitalPhi][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?

Last edited: Aug 23, 2015
2. Aug 24, 2015

soarce

The Mathematica code is not readable, try to reformat it or export it as text of LaTeX.

For $\alpha=0.99$ the decay is very steep and maybe the Runge-Kutta may not handle it very good. Try to use a steep method instead.

3. Aug 24, 2015

spaghetti3451

My edit button is missing, so let me repost my code below:

$$α = 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]$$

4. Aug 24, 2015

spaghetti3451

Might you perhaps provide examples of steep methods or links to available resources? I looked online, but couldn't find resources.

5. Aug 24, 2015

soarce

I mean methods for stiff problems, problems which display steep changes. My mistake.

6. Aug 24, 2015

spaghetti3451

Do you mean that the decay for $\alpha = 0.99$ is not reasonable, perhaps because the solution suddenly takes a plunge downwards after being stable for quite a while?

7. Aug 24, 2015

soarce

Yes, for problems where you have sudden variations some of the numerical methods won't work. However, you function remembers me of the solution of the cubic-quintic nonlinear Schroedinger equation for which we used some numerical relaxation method to get the top-flat profiles.

8. Aug 24, 2015

spaghetti3451

So, do you mean that the steep decay itself is not a problem at all, but rather that some numerical methods such as Runge-Kutta do not give accurate plots in such cases?

9. Aug 24, 2015

soarce

10. Aug 24, 2015

spaghetti3451

Thanks! I'll look into it.

I've also heard of other numerical methods such as finite differences, collocation, multiple shooting, decoupling, and fixed point iteration. Do you think these methods might be any good?

11. Aug 24, 2015

soarce

The shooting method is very simple, but is higly sensitive to the initial condition. You "$\alpha=0.99$" is only an aproximated value because internally the number is stored in floating point convention.
1) try to increase the precision
2) what happens if you increase the total number of iterations ?
3) as your function needs to aproach zero (asimptotically) I would test the value $\Phi(r=130)$ against zero, remove the "1.2 comparison".

I think that shooting method shold work here, let's first check that we do everithing right. I never used Mathematica, but only Matlab and some C code from Numerical Recipes. With finite differences you are locked to a fixed mesh and you will to readjust it for each computation (evaluate the spatial extension of the function, how rapidly changes in space), multiple shooting it's shooting with two or more parameters and you can't apply directly the bisection method and the fixed point iteration depends on the problem you want to solve. For the top-flat profiles the finite diferences and fixed point methods need a good aproximation of your solution in order to converge.

12. Aug 25, 2015