1
$\begingroup$

I have two coupled ODEs that I am trying to solve numerically. It appears that there is a singularity in the solution to the equations which I am unsure how to get past. Both functions $\alpha$ and $\phi$ should asymptotically converge to 1 at infinity, however I can't set my boundary conditions to be more than 3.4; I get the error:

NDSolve::ndsz: At r == 0.5165576868139004`, step size is effectively zero; singularity or stiff system suspected.

I have seen suggestions that I should start away from a singularity (i.e. r = 0) or change boundary conditions etc, but my boundary conditions are fixed.

How do I get around this issue?

Code below:

(* Set Constants *)
n = 1;
λ = 1;

(* Set functions *)
p[r_] = 1/r^2;
q[r_] = 1/r;

(* ODE's, which for some reason I called pde1/2 *)
pde1 = D[ϕ[r], {r, 2}] + q[r]*D[ϕ[r], {r, 1}] - p[r]*(n - α[r])^2 ϕ[r] + λ/2 (1 - ϕ[r]^2) ϕ[r] == 0
pde2 = D[α[r], {r, 2}] - q[r]*D[α[r], {r, 1}] + (n - α[r]) ϕ[r]^2 == 0

(* Min/max numbers *)
min = 0.000001;
max = 3.4; (* Define a "finite infinity" - I would like to set this higher *)

(* Boundary Conditions *)
bc = {ϕ[min] == 0, α[min] == 0, ϕ[max] == 1 - 1/max Exp[-max], α[max] == 1 - 5/max Exp[-max]};

sol = NDSolve[
       {pde1, pde2, bc}, {ϕ, α}, {r, min, max}, 
       Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}}, 
       AccuracyGoal -> 15, PrecisionGoal -> 15, MaxSteps -> 20000
      ]

Plot[Evaluate[{ϕ[r], α[r]} /. sol], {r, 0, 4}]
$\endgroup$
5
$\begingroup$

In this BVP, Mathematica uses the shooting method, and you're shooting at an unstable solution. You can iteratively approach the solution, extending the interval of integration at each step, but at machine precision, you can extend it only so far.

I'll set up bc to depend on symbolic min and max, so we can fiddle with them at each iteration:

(*Boundary Conditions*)
Clear[min, max];
bc = {ϕ[min] == 0, α[min] == 0, ϕ[max] == 1 - 1/max Exp[-max],
  α[max] == 1 - 5/max Exp[-max]};

Then with max = 3.4, we get an initial solution, just as the OP does.

Block[{min = 0.000001, max = 3.4},
 sol1 = NDSolve[{pde1, pde2, bc}, {ϕ, α}, {r, min, max}]
 ]

From this point, the basic iterative step uses the previous step to determine the starting initial conditions for the shooting method. So thus:

Block[{sol0 = First@sol1, min = 0.000001, max = 7.5}, 
 sol2 = NDSolve[{pde1, pde2, bc}, {ϕ, α}, {r, min, max}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" ->
       {ϕ[min] == (ϕ[min] /. sol0), ϕ'[min] == (ϕ'[min] /. sol0),
        α[min] == (α[min] /. sol0), α'[min] == (α'[min] /. sol0)}}]
 ]

Mathematica graphics

The setting max = 7.5 was determined by running NDSolve first out to 10; but it developed a singularity at 7.7, so I stepped it back a little. (Push too close to the limit and the shooting method might have problems.)

Iterate again:

Block[{sol0 = First@sol2, min = 0.000001, max = 15}, 
 sol3 = NDSolve[{pde1, pde2, bc}, {ϕ, α}, {r, min, max}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" ->
       {ϕ[min] == (ϕ[min] /. sol0), ϕ'[min] == (ϕ'[min] /. sol0),
        α[min] == (α[min] /. sol0), α'[min] == (α'[min] /. sol0)}}]
 ]

Mathematica graphics

Hmm, no error but something looks wrong.

Plot[Evaluate[{ϕ[r], α[r]} /. sol3], {r, 0, 15}, GridLines -> {None, {1}}]

Mathematica graphics

Well, that might be as good as you can do at a WorkingPrecision of MachinePrecision.

If you use WorkingPrecision -> 20, you get the following and probably could push the integration farther than r == 15.:

Mathematica graphics

An alternative approach is to implement one's own shooting method, and adapt the initial conditions to integrate as far as possible. This was one of the methods used in my answer to Nonlinear differential equation: numerical solution.

$\endgroup$
  • $\begingroup$ Thank you, this answer is fantastically explained. $\endgroup$ – Akoben May 6 '16 at 8:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.