I am trying to solve a particular Cauchy problem given by

Nonlinear Stiff ODE

I found from a particular paper that the solutions looks like enter image description here

For the auxiliary conditions

enter image description here

For only specific values of $a_{i}$

I found from another paper that the problem can be solved with the Mathematica code


createLargestRegion[h_, d_, mult_] :=

 Block[{ bound = Min[mult*Abs[d]/h, 19], n = 1, result, previous }, 
     Check[result[r_] := 
        First[\[Theta][r] /.

          NDSolve[{\[Theta]''[r] + (1/r ) \[Theta]'[r] -
              (0.5/(r^2)) Sin[
                2 \[Theta][r]] == ((h /2.0) Sin[\[Theta][r]] + 
               (d /r) Sin[\[Theta][r]]^2 ), \[Theta][0.001] == 
             Pi, \[Theta][bound] == 0.1}, \[Theta][r], r, 
           "Method" -> {"Shooting", 
             "StartingInitialConditions" -> { \[Theta][0.001] == 
                Pi, \[Theta]'[0.001] == (-Pi/10)}}]]], False]] || 
    FindMaximum[{result[t], t > 0.1 && t <= bound}, t] > Pi || 
    result [0.5] > result [0.001], bound = bound - 0.1;

hdInterpFuncs = 
 Table[Table[{i*0.1, j*-0.1, 
    Evaluate[createLargestRegion[i*0.1, j*-0.1, 19.0]]}, {j, 3, 18, 
    3}], {i, 3, 18, 3}]

I completeley understand NDSolve and Shooting method, but i am unable to understand the rest of the code. Also when i run the code i am getting

NDSolve::berr: The scaled boundary value residual error of 6.892078251140068`*^7 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.

Any help in helping me understand the code and or figuring out why i am getting the error will be thankful. Also if there are other ways to solve the problem with better code will be thankful too.

(English is not my first language forgive grammar errors.)

Related Paper here


1 Answer 1


The Mathematica code in the question addresses the equation

θ''[r] + θ'[r]/r == Sin[2 θ[r]]/(2 r^2) + h Sin[θ[r]]/2 + d Sin[θ[r]]^2 /r

for 0.3 < h < 1.8 and -1.8 < d < -0.3. In contrast, the LaTeX version of the equation assumes d == 1. I am confident that the figure in the question cannot be reproduced with d == 1. Solving the equation for negative d can be accomplished as follows.

Since the question seeks solutions which asymptotically approach zero at very large r, it makes sense first to linearize the equation about zero and obtain its solution symbolically,

θ''[r] + θ'[r]/r == θ[r]/r^2 + h θ[r]/2

with solutions C[1] BesselI[1, r Sqrt[h/2]] + C[2] BesselK[1, r Sqrt[h/2]]. The second term approaches zero exponentially at large r and is the desired specific solution. However, the presence of the first term causes infinitesimal numerical errors to grow exponentially, unless great care is taken. (This is the typical behavior of the separatrix of a nonlinear ODE.) Matching the numerical solution at large r to C[2] BesselK[1, r Sqrt[h/2]] is an effective boundary condition, and eliminating C[2] yields

θ'[rm] == c θ[rm]
c = N[D[BesselK[1, r Sqrt[h/2]], r]/BesselK[1, r Sqrt[h/2]] /. r -> rm, 30];

with rm the value of r at which the numerical calculation ends.

The equation now can be solved by

Clear[h, d, θp0]; r0 = 10^-3; rm = 20;
sp = ParametricNDSolveValue[{
    c = N[D[BesselK[1, r Sqrt[h/2]], r]/BesselK[1, r Sqrt[h/2]] /. r -> rm, 30]; 
    θ''[r] +  θ'[r]/r  == Sin[2 θ[r]]/(2 r^2) + h  Sin[θ[r]]/2 + d Sin[θ[r]]^2 /r, 
    θ[r0] == Pi + r0 θ'[r0], θ'[rm] == c θ[rm]}, {θ[r], θ'[r]}, {r, r0, rm}, {h, d, θp0}, 
    StartingStepSize -> r0/20, WorkingPrecision -> 30,
    Method -> {"Shooting", 
        "StartingInitialConditions" -> {θ[r0] == Pi + r0 θp0, θ'[r0] == θp0}}];

The challenge, of course, is to find a good initial guess for θp0, the slope of θ at r0. For a single set of parameters h, d, this can be accomplished by scanning the range of θp0 from -0.1 to -0.4 in steps of 0.05, which can be done in a few minutes or less. For instance,

param = {3/10, -1, -28/100}; s = sp @@ param;
Plot[First@s, {r, r0, rm}, PlotRange -> All]        
c = N[D[BesselK[1, r Sqrt[h/2]], r]/BesselK[1, r Sqrt[h/2]] /. r -> rm, 30] /. 
    Thread[{h, d, θp0} -> param];
{(Last@s - c First@s) /. r -> rm, Last@s /. r -> r0}

enter image description here

(* {-1.17413566404716061614*10^-11, -0.32320602499612} *)

Note that the format of this and all subsequent plots is based on

Themes`AddThemeRules["mystyle", AxesStyle -> Directive[Black, Bold, Medium], 
    ImageSize -> Medium, PlotRange -> {0, Pi}, Ticks -> {Automatic, {0, Pi/2, Pi}}];
$PlotTheme = "mystyle";

Now that this solution has been found, and with it the actual value of θ'[r0], -0.32320602499612, nearby solutions in parameter space can be obtained economically by extrapolating from this value of θ'[r0]. For instance, solutions over a range of d can be obtained fairly quickly from

start = Rationalize[{{30/100, -1, Last@sp[30/100, -1, -32/100] /. r -> r0}, 
    {30/100, -99/100, Last@sp[30/100, -99/100, -33/100] /. r -> r0}}, 0];
new = 1; i = 3;
While[new != 0 && start[[i - 1, 2]] < 0, tem = 2 start[[i - 1]] - start[[i - 2]]; 
    new = Rationalize[Last@Check[sp @@ tem, {0, 0}] /. r -> r0, 0]; 
    If[new != 0, AppendTo[start, Join[Most@tem, {new}]]]; i++];
start30 = start;

Show[Plot[sp @@ #, {r, r0, 8.6}, Ticks -> {{0, 5}, {0, Pi/2, Pi}}, 
    PlotStyle -> Directive[Thickness[.005], Hue[-#[[2]]]], 
    AspectRatio -> .65] & /@ start30[[1 ;; 95 ;; 10]]]

enter image description here

The curves, top to bottom, correspond to increasing d from d == -1 to d == - 1/10 in steps of -1/10. Solutions for still larger values of d are exceedingly narrow in r.

With solutions now available for many d, further solutions for ranges of h and each value of d can be computed in the same way. Typical plots similar to that in the question are, for d == -1 and h == {3/10, 41/100, 7/10} (Red, Green, and Blue, respectively)

Show[Plot[sp @@ #, {r, r0, 8.6}, Ticks -> {{0, 5}, {0, Pi/2, Pi}}, 
PlotStyle -> Directive[Thickness[.005], 
    Switch[#[[1]], 3/10, Red, 41/100, Green, 7/10, Blue]]] & 
    /@ {starth[[1]], starth[[12]], starth[[41]]}, AspectRatio -> .65]

enter image description here

and also for d == {-3/4, -1/2}, respectively.

enter image description here

enter image description here

The d == -3/4 curves appear to match the figure in the question best, although differences are evident, expecially near the right of the plots. I would speculate that the calculations used to obtained the figure in the queston employed too small a value of rm.

Returning now to the code in the question, I would observe that it has two problems. First, it uses as a large r boundary condition θ[bound] == 0.1, which is inaccurate. Second, it uses as an initial guess for all h, d the value θ[0.001] == Pi/10, which also is inaccurate. It then iteratively reduces bound until NDSolve converges, often at unrealistically small values of bound. Incidentally, the repeated convergences failures give rise to the error messages mentioned in the question. They can be eliminated, if desired, by placing Quiet@ immediately before NDSolve. Doing so does not, of course, eliminate the inaccuracies just described.

  • $\begingroup$ Could you please kindly give more details on the way you used to replace the boundary condition in infinity. Instead of theta(infinity)==0 I naively would use a boundary condition theta(rm)==0 for a large enough rm. Could you please comment: (1) Do you see some problems with the condition theta(rm)==0? (2) What advantages give your approach of using a linearized equation instead. $\endgroup$ Jan 30, 2018 at 10:28
  • 1
    $\begingroup$ Dear @bbgodfrey I would like to kindly ask you to return to my questions formulated above. I find the approach you offered to be exceptionally important, and (I hope) a more detailed discussion of this question should be of interest for lots of people. Maybe a good idea would be to formulate the problem (of a numerical solution of ODE or PDE with a boundary condition in infinity) as a stand-alone question, which you might both formulate and give the answer with details and examples that you have already collected. $\endgroup$ Feb 7, 2018 at 7:43
  • $\begingroup$ @AlexeiBoulbitch I agree that addressing this question generically makes sense and I have been thinking about how to do so. However, I have not had time to do the problem justice. I certainly intend to do so, but I probably shall not be able to for a week or more . Thanks for your patience. $\endgroup$
    – bbgodfrey
    Feb 7, 2018 at 14:43
  • $\begingroup$ Of course, take your time. $\endgroup$ Feb 7, 2018 at 14:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.