1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solving an ODE using shooting method

  1. Aug 23, 2015 #1
    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:

    Capture.jpg



    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:

    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: Aug 23, 2015
  2. jcsd
  3. Aug 24, 2015 #2
    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.
     
  4. Aug 24, 2015 #3
    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]
    $$
     
  5. Aug 24, 2015 #4
    Might you perhaps provide examples of steep methods or links to available resources? I looked online, but couldn't find resources.
     
  6. Aug 24, 2015 #5
    I mean methods for stiff problems, problems which display steep changes. My mistake.
     
  7. Aug 24, 2015 #6
    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?
     
  8. Aug 24, 2015 #7
    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.
     
  9. Aug 24, 2015 #8
    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. Aug 24, 2015 #9
  11. Aug 24, 2015 #10
    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?
     
  12. Aug 24, 2015 #11
    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.
     
  13. Aug 25, 2015 #12

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    [edit] misread - after many readings !
     
    Last edited: Aug 25, 2015
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted