• Support PF! Buy your school textbooks, materials and every day products Here!

Solving an ODE using shooting method

  • #1
1,344
32

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:

Capture.jpg




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:

phi_0_50.jpg


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

phi_0_99.jpg

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:

Answers and Replies

  • #2
155
24
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
1,344
32
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
1,344
32
Might you perhaps provide examples of steep methods or links to available resources? I looked online, but couldn't find resources.
 
  • #5
155
24
I mean methods for stiff problems, problems which display steep changes. My mistake.
 
  • #6
1,344
32
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
155
24
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
1,344
32
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?
 
  • #10
1,344
32
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
155
24
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
BvU
Science Advisor
Homework Helper
2019 Award
13,361
3,165
[edit] misread - after many readings !
 
Last edited:

Related Threads on Solving an ODE using shooting method

Replies
4
Views
8K
Replies
0
Views
6K
Replies
2
Views
5K
Replies
2
Views
7K
Replies
4
Views
2K
Replies
1
Views
657
  • Last Post
Replies
0
Views
884
Replies
1
Views
1K
Replies
0
Views
3K
Top