2

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] == 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] == 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. Profile solutions plot

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.

xzczd
  • 65,995
  • 9
  • 163
  • 468
Ben Leather
  • 31
  • 4
  • 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 - https://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 00: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 01:51

1 Answers1

1

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] == 0, y[inf] == 1};

points = 100;
domain = {0, inf};
grid = Array[# &, points, domain];
difforder = 4;
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
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}]

Mathematica graphics

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • 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 09: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 01 '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: https://mathematica.stackexchange.com/a/25616/1871 – xzczd Dec 01 '17 at 15:02
  • If eqns = Join[{eqn1, eqn2}]; and one specifies the boundary conditions such that bc = {y[0] == 0, y[inf] == 1, z[0] == 0, z[inf] == 1};, then wouldn't the subsequent code take the formae = ptoafunc[x^2 # & /@ eqns // Simplify][[2 ;; -2]];`? – Ben Leather Dec 02 '17 at 16:20
  • @BenLeather You need delete = #[[2 ;; -2]] &; ae = delete /@ ptoafunc[Map[x^2 # & , eqns,{2}] // Simplify] – xzczd Dec 02 '17 at 16:31
  • This gives me a ‘Part:’ error, saying it cannot take ‘solutions -2 through 2’? Am I missing something? – Ben Leather Dec 03 '17 at 23:03
  • @BenLeather Please notice it's 2 ;; -2, not -2 ;; 2. – xzczd Dec 04 '17 at 05:28
  • I have rectified this but the same error still occurs? – Ben Leather Dec 04 '17 at 08:02
  • @BenLeather Please add the specific code to your question. – xzczd Dec 04 '17 at 08:07
  • delete = #[[2 ;; -2]] &; ae = delete /@ ptoafunc[Map[x^2 # &, eqns, {2}] // Simplify]; - gives the error Part::take: Cannot take positions 2 through -2 in eqns. – Ben Leather Dec 04 '17 at 08:48
  • Let us continue this discussion in chat. – xzczd Dec 04 '17 at 08:57