# 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. – Ben Leather Nov 29 '17 at 0:50
• The Power::infy error comes from dividing by x == 0. I don't think that's the whole problem, though. – aardvark2012 Nov 29 '17 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>>]. – Ben Leather Nov 29 '17 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. – xzczd Nov 29 '17 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. – Ben Leather Dec 1 '17 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 – xzczd Dec 1 '17 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]];? – Ben Leather Dec 2 '17 at 16:20