# Errors when numerically solving differential equations for global vortex profiles

I am currently to numerically solve the following differential equation for the profile of global vortices with a simple complex scalar field:

y''[x] + y'[x]/x  - y[x]/x^2 + (1 - y[x]^2) y[x] == 0,
y == 0, y[Infinity] == 1.


This is my following code in Mathematica to try and obtain the solution.

inf = 100;
eqn = y''[x] + y'[x]/x  - y[x]/x^2 + (1 - y[x]^2) y[x] == 0;
NDSolve[{eqn, y == 0, y[inf] == 1}, y, {x, 0, 10}]


Unfortunately, Mathematica's NDSolve function gives the following errors due to infinities being obtained:

Power::infy: Infinite expression 1/0.^2 encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

Power::infy: Infinite expression 1/0. encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

Power::infy: Infinite expression 1/0.^2 encountered.

General::stop: Further output of Power::infy will be suppressed during this calculation.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

General::stop: Further output of Infinity::indet will be suppressed during this calculation.

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0..


I have tried to specify the method of NDSolve to the "Shooting" method but this yields similar problems. I know this form of differential equation can be solved numerically to give profiles such as the ones given in the figure below. In this plot y[x] is simply labelled f and x is equivalent to r, k is the co-efficient of the third term, -y[x]/x^2, which I have taken to be equal to 1 for convenience. The singularity at x = 0 can be avoided by an appropriate solution y[x].

I am therefore wondering how Mathematica can be used to obtain similar results. Any help on this matter would be greatly appreciated.

• Sorry, it is simply differing notation, I didn't think it was that difficult to understand. The notation is also explained in the question. The plot is taken from a paper - arxiv.org/pdf/0903.1528.pdf - and shown as proof that the differential equation can be solved using numerical methods despite the singularity. Nov 29, 2017 at 0:50
• The Power::infy error comes from dividing by x == 0. I don't think that's the whole problem, though. Nov 29, 2017 at 1:51

A step backward i.e. make use of finite difference method (FDM) seems not to be a bad idea in this case. I'll use pdetoae for the generation of difference equation:

Clear@k;
inf = 100;
eqn = y''[x] + y'[x]/x - k^2 y[x]/x^2 + (1 - y[x]^2) y[x] == 0;
bc = {y == 0, y[inf] == 1};

points = 100;
domain = {0, inf};
grid = Array[# &, points, domain];
difforder = 4;
(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[y[x], grid, difforder];

ae = ptoafunc[x^2 # & /@ eqn // Simplify][[2 ;; -2]];

sollst = Table[
ListInterpolation[
With[{initial = 1}, FindRoot[{ae, bc}, {y@#, initial} & /@ grid]][[All, -1]],
grid], {k, 10}]

(* Alternatively: *)
(*
lSSolve[obj_List, constr___, x_, opt : OptionsPattern[FindMinimum]] :=
FindMinimum[{1/2 obj^2 // Total, constr}, x, opt]
lSSolve[obj_, rest__] := lSSolve[{obj}, rest]

fullae = ptoafunc[x^2 # & /@ eqn // Simplify];
sollst = Table[
ListInterpolation[
With[{initial = 1},
lSSolve[Subtract @@@ Flatten@{fullae, bc}, {y@#, initial} & /@ grid]][[2, All, -1]],
grid], {k, 10}]
*)
Plot[sollst[x] // Through // Evaluate, {x, 0, 50}, PlotRange -> All,
PlotStyle -> Table[Blend[{Blue, Magenta}, x], {x, 0, 1, 1/9}], GridLines -> Automatic,
AxesLabel -> {r, f}] • When I try to use this method within Mathematica I get the following error when trying to output ae: Part::take: Cannot take positions 2 through -2 in pdetoae[y[x],{0,100/99,200/99,100/33,400/99,500/99,200/33,700/99,800/99,100/11,1000/99,100/9,400/33,1300/99,1400/99,500/33,1600/99,1700/99,200/11,1900/99,2000/99,700/33,200/9,2300/99,<<4>>,2800/99,2900/99,1000/33,3100/99,3200/99,100/3,3400/99,3500/99,400/11,3700/99,3800/99,1300/33,4000/99,4100/99,1400/33,4300/99,400/9,500/11,4600/99,4700/99,1600/33,4900/99,<<50>>},4][<<1>>==x <<1>>]. Nov 29, 2017 at 9:23
• @BenLeather As mentioned in the anotation, definition of pdetoae is not included in this post, please find it in the link at the beginning of the post. Nov 29, 2017 at 11:16
• @xzcxd Can you apply the same finite difference method for a system of differential equations such that: eqn1 = y''[x] + y'[x]/x - (k - z[x])^2 y[x]/x^2 + (1 - y[x]^2) y[x]/2 == 0; eqn2 = z''[x] - z[x]/x + (k - z[x]) z[x]^2 == 0. Here, y[x] and z[x] are the two variables we are solving for. Dec 1, 2017 at 14:44
• @BenLeather Use pdetoae[{y[x], z[x]}, …, the subsequent code also needs to be slightly modified depending on your b.c. of course. Well, if you're having difficulty in understanding those /@, @, etc., consider starting from here: mathematica.stackexchange.com/a/25616/1871 Dec 1, 2017 at 15:02
• If eqns = Join[{eqn1, eqn2}]; and one specifies the boundary conditions such that bc = {y == 0, y[inf] == 1, z == 0, z[inf] == 1};, then wouldn't the subsequent code take the form ae = ptoafunc[x^2 # & /@ eqns // Simplify][[2 ;; -2]];? Dec 2, 2017 at 16:20