# Error when solving 't Hooft-Polyakov radial equations using NDSolve

I'm trying to solve 2 coupled nonlinear ODEs using NDSolve, but the solution fails when the parameter $\lambda$ increases. The equations are $$r^2 a''(r)= a(r) [a(r)-1] [a(r)-2]- r^2 [1-a(r)] h(r)^2,$$ $$\left(r^2 h'(r)\right)'= 2 [1-a(r)]^2 h(r)+ \lambda\,r^2 [h(r)^2-1]\, h(r).$$ And the boundary conditions are $$a(0)=h(0)=0,$$ $$a(\infty)= h(\infty)=1.$$ I've tried to solve it in the following way:

\[Lambda] = 0.1;
r1 = 10^(-6); r2 = 9.5;
eqn = {r^2*D[D[a[r], r], r] == a[r]*(a[r] - 1)*(a[r] - 2) - r^2*h[r]^2*(1 - a[r]),
D[r^2*D[h[r], r], r] == 2*h[r]*(1 - a[r])^2 + \[Lambda]*r^2*(h[r]^2 - 1)*h[r]};
bc = {a[r1] == 0, a[r2] == 1, h[r1] == 0, h[r2] == 1};
sol = NDSolveValue[{eqn, bc} /. {a -> (#1^2*g[#1] & ), h -> (#1*j[#1] & )}, {g, j}, {r, r1, r2},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5}, WorkingPrecision -> MachinePrecision];
Plot[{1, 1 - r^2*sol[][r], r*sol[][r]}, {r, r1, r2}]


And, indeed, I get a solution. This works for $\lambda=0$ and $\lambda=0.1$. However, I would like to get the solution for a larger range of $\lambda$, for example $\lambda=1$. When I try to use the code above and set $\lambda=1$ Mathematica gives me the following

NDSolveValue: At r == 1.2395178650126792, step size is effectively zero; singularity or stiff system suspected.

NDSolveValue: The scaled boundary value residual error of 62.94578090521123 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.


As a consequence, the plot does not correspond to a solution.

So, my question is: How can I get a solution to these equations when $\lambda=1$, for example? Is it possible to calculate for larger values?

Besides that, how can I use the results of NDSolve in an integral of the form

$$E=4\pi \int_{0}^{\infty}dr\,\left[(a')^2 + \frac{1}{2}\frac{a^2(a-2)^2}{r^2}+ \frac{1}{2}r^2(h')^2 + h^2(a-1)^2 + \frac{\lambda}{4}r^2(h^2-1)^2 \right]?$$

I know this problem has a solution, because the equations are those of the 't Hooft-Polyakov monopole and the last integral is its mass as a function of $\lambda$. But, I'm new in Mathematica and I don't know how to correctly reproduce the solutions which are available at textbooks.

Any hint is valuable to me!

***Update:***Just to make it clear, I would like to stress that I know this problem has been solved. The plots of the solution are available at textbooks on topological solitons, such as [Manton, Sutcliffe] and [Shnir]. What I'm trying to say is that I would like to understand how to obtain this solution in Mathematica. The question marked as a possible duplicate is the solution of the Abelian string (also ANO vortex), which is different. One can see this when looking, for example, into the kinetic part given in those terms with derivatives. Could you, please, give some hints on how I should improve my code?

• The short answer is that it seems to be very sensitive to initial conditions. You should look up the "Shooting" method and adjust "StartingInitialConditions". You can use a nearby solution for them. A lambda that works may be increased a little and a new solution produced. Little by little you might be able to push lambda up to 1. – Michael E2 Feb 28 '18 at 0:56
• I'm inclined to say that this question is morally a duplicate of this question. This set of differential equations was also studied in this paper. – QuantumDot Feb 28 '18 at 13:50
• Possible duplicate of Numerical solution of coupled ODEs with boundary conditions – QuantumDot Feb 28 '18 at 13:51
• Dear @QuantumDot, I know this paper, but there they do not explicitly show how the calculation was done and, besides that, I was trying to reproduce it in Mathematica. This monopole solution is well-known, but books on topological solitons do not discuss in detail how it is done. I edited the end of my question, if you can take a look at it. – MLPhysics Feb 28 '18 at 18:57
• Thank you for the hint, @Michael. As I mentioned in the question, I'm new in Mathematica (just a few days using) and I didn't know about "Shooting" in Wolfram Language. What I have seen was your answer to the question  (mathematica.stackexchange.com/questions/91854/…), but I couldn't generalize the method for my 2 coupled equations. So, about my question: could you, please, just explain me a little more on pushing lambda up? – MLPhysics Feb 28 '18 at 20:11

Once again, compared to "Shooting" method that is the default and currently the only available method for solving nonlinear boundary value problem (BVP) inside NDSolve, finite difference method (FDM) seems to be a better choice for solving this problem. I'll use pdetoae for the generation of difference equations.

r1 = 10^(-6); r2 = 9.58;
eqn = {r^2*D[D[a[r], r], r] == a[r]*(a[r] - 1)*(a[r] - 2) - r^2*h[r]^2*(1 - a[r]),
D[r^2*D[h[r], r], r] == 2*h[r]*(1 - a[r])^2 + λ*r^2*(h[r]^2 - 1)*h[r]};
bc = {a[r1] == 0, a[r2] == 1, h[r1] == 0, h[r2] == 1};
{neweqn, newbc} = {eqn, bc} /. {a -> (#1^2*g[#1] &), h -> (#1*j[#1] &)};
difforder = 4;
domain = {r1, r2};
points = 50;
grid = Array[# &, points, domain];

(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[{g, j}[r], grid, difforder];
delete = #[[2 ;; -2]] &;
ae = delete /@ ptoafunc@neweqn;
aebc = ptoafunc@newbc;

initialfunc[g][r_] = 1/2;
initialfunc[j][r_] = 1/2;

sollst = Partition[
Block[{λ = 1},
FindRoot[{ae, aebc},
Flatten[Table[{var[index], initialfunc[var][index]}, {var, {g, j}}, {index, grid}],
1]]][[All, -1]], points];

{solg, solj} = ListInterpolation[#, domain] & /@ sollst;

Plot[{1, 1 - r^2*solg[r], r*solj[r]} // Evaluate, {r, r1, r2}] As one can see, the initial guesses used here i.e. initialfunc[g][r] and initialfunc[j][r] are rather rough, but we still obtain the desired result.

As to the integral, I think the straightforward solution isn't bad:

Ε =
With[{a = r^2*solg[r], h = r*solj[r], λ = 1},
4 Pi NIntegrate[
D[a, r]^2 + 1/2 a^2 (a - 2)^2/r^2 D[h, r]^2 +
h^2 (a - 1)^2 + λ/4 r^2 (h^2 - 1)^2, {r, r1, r2}]]
(* 9.74871 *)


You may adjust points etc. to make the result preciser.

# Update

It turns out that

1. the transformation $$a(r)=r^2 g(r)$$ and $$h(r)=r j(r)$$

2. approximating $$r_1=0$$ with $$r_1=10^{-6}$$

is not necessary:

r1 = 0; r2 = 95/10;
eqn = {r^2*D[D[a[r], r], r] == a[r]*(a[r] - 1)*(a[r] - 2) - r^2*h[r]^2*(1 - a[r]),
D[r^2*D[h[r], r], r] == 2*h[r]*(1 - a[r])^2 + λ*r^2*(h[r]^2 - 1)*h[r]};
bc = {a[r1] == 0, a[r2] == 1, h[r1] == 0, h[r2] == 1};
{neweqn, newbc} = {eqn, bc};
difforder = 4;
domain = {r1, r2};
points = 50;
grid = Array[# &, points, domain];

(*Definition of pdetoae isn't included in this post,please find it in the link above.*)

ptoafunc = pdetoae[{a, h}[r], grid, difforder];
delete = #[[2 ;; -2]] &;
ae = delete /@ ptoafunc@neweqn;
aebc = ptoafunc@newbc;
initialfunc[a][r_] = 1/2;
initialfunc[h][r_] = 1/2;
sollst = Partition[
Block[{λ = 1},
FindRoot[{ae, aebc},
Flatten[Table[{var[index], initialfunc[var][index]}, {var, {a, h}}, {index, grid}],
1]]][[All, -1]], points];
{sola, solh} = ListInterpolation[#, domain] & /@ sollst;

Plot[{1, 1 - sola[r], solh[r]} // Evaluate, {r, r1, r2}]


The result is the same so I'd like to omit it here.

• Thank you very much @xzczd! I will try to understand it, since it is a sophisticated code for me. – MLPhysics Mar 1 '18 at 17:01
• Have a look at the comment I posted. Do you perhaps have an opinion on that? – user21 Mar 7 '18 at 8:06

This is not an answer but a remark about the boundary conditions.

Let's look at the equation:

r1 = 10^(-6);
r2 = 95/10;
lambda = 1/10;
eqn = {r^2*D[D[a[r], r], r] ==
a[r]*(a[r] - 1)*(a[r] - 2) - r^2*h[r]^2*(1 - a[r]),
D[r^2*D[h[r], r], r] ==
2*h[r]*(1 - a[r])^2 + \[Lambda]*r^2*(h[r]^2 - 1)*h[r]};
bc = {a[r1] == 0, a[r2] == 1, h[r1] == 0, h[r2] == 1};
{teqn, tbc} = {eqn, bc} /. {a -> (#1^2*g[#1] &), h -> (#1*j[#1] &)};
teqn /. {g -> u, j -> v} // Expand;
{gsol, jsol} =
NDSolveValue[{teqn, tbc} /. \[Lambda] -> lambda, {g, j}, {r, r1,
r2}, Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5},
WorkingPrecision -> MachinePrecision];


If we look at the plot:

Plot[{gsol[r], jsol[r]}, {r, r1, r2},
PlotRange -> {{-0.05, 0.15}, All}] We note a jump at the left hand side. The solver has trouble satisfying this BC (and the solver xzczd wrote has a similar problem).

The actual values:

{gsol[r] - 0, jsol[r] - 0} /. r -> r1
{gsol[r] - 1/90.25, jsol[r] - 1/9.5} /. r -> r2

{0.00147074, 1.38778*10^-16}
{-1.23805*10^-10, 2.58372*10^-11}


Do not quite match the requested bc at r1 for gsol. If one looks at:

{{r^2 gsol[r] - 0, r jsol[r] - 0} /.
r -> r1, {r^2 gsol[r] - 1, r jsol[r] - 1} /. r -> r2}

{{1.47074*10^-15, 1.38778*10^-22}, {-1.11734*10^-8, 2.45454*10^-10}}


Things look better. But the issue of the kink remains.

If you were, just for fun (this is not the right thing to do), to replace the bcs in the following way:

{gsol2, jsol2} =
NDSolveValue[{teqn, {g[1/1000000] == gsol[0.001],
90.25 g[9.5] == 1, j[1/1000000] == jsol[0.001],
9.5 j[9.5] == 1}} /. \[Lambda] -> lambda, {g, j}, {r, r1, r2},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5},
WorkingPrecision -> MachinePrecision];


You will note that the convergence is much quicker then with the above BCs. Where did you find these BCs. Are sure they make sense and are correct?

• I'm afraid you've made a mistake on checking, the correct one should be {{r^2 gsol[r] - 0, r jsol[r] - 0} /. r -> r1, {r^2 gsol[r] - 1, r jsol[r] - 1} /. r -> r2} :) – xzczd Mar 7 '18 at 8:27
• @xzczd thanks. Even if one takes r^2 gsol[r] the kink per se remains and indicates trouble that that boundary. – user21 Mar 7 '18 at 8:42
• I think the kink is reasonable, because r1 = 10^(-6) is probably an approximation of r1 = 0. For r1 = 0, a and j can be non-zero. – xzczd Mar 7 '18 at 9:02
• @xzczd, but the problem statement is that the bc a(0)=0`. Too me this looks odd. – user21 Mar 7 '18 at 9:09
• Oh, and I'm sorry to bother you, I know this website is not "MATLAB SE"! – MLPhysics Mar 12 '18 at 22:57