Solving an ODE using shooting method

[spaghetti3451]

Homework Statement

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

Homework Equations

My solution is supposed to reproduce the following plots:

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));
differential equation = {-\[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?

 

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

 

[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));
differential equation = {-Φ''[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]
$$

 

[spaghetti3451]

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

 

[soarce]

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

 

[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?

 

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

 

[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?

 

[soarce]

Indeed, the numericla method is not able to keep track of the steep variations. There are special numerical methods, but I think that Mathematica must have support for stiff problems. Try to investigate a little more here http://reference.wolfram.com/language/tutorial/NDSolveStiffnessSwitching.html

 

[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?

 

[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 approach 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.

 

[BvU]

[edit] misread - after many readings !