Solving a second-order nonlinear differential equation

Question

I am trying to solve a particular Cauchy problem given by

I found from a particular paper that the solutions looks like

For the auxiliary conditions

For only specific values of $a_{i}$

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

Off[General::stop]

createLargestRegion[h_, d_, mult_] :=

 Block[{ bound = Min[mult*Abs[d]/h, 19], n = 1, result, previous }, 
  While[Not[
     Check[result[r_] := 
       Evaluate[
        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;
   n++;];
  result[r]]

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

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}

(* {-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]]]

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]

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

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.