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}]
Improve this question
$\endgroup$

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you, this answer is fantastically explained. $\endgroup$ Commented May 6, 2016 at 8:45

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.