2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Ben Leather Nov 29 '17 at 0:50
  • $\begingroup$ The Power::infy error comes from dividing by x == 0. I don't think that's the whole problem, though. $\endgroup$ – aardvark2012 Nov 29 '17 at 1:51
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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>>]. $\endgroup$ – Ben Leather Nov 29 '17 at 9:23
  • $\begingroup$ @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. $\endgroup$ – xzczd Nov 29 '17 at 11:16
  • $\begingroup$ @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. $\endgroup$ – Ben Leather Dec 1 '17 at 14:44
  • $\begingroup$ @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 $\endgroup$ – xzczd Dec 1 '17 at 15:02
  • $\begingroup$ 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 form ae = ptoafunc[x^2 # & /@ eqns // Simplify][[2 ;; -2]];`? $\endgroup$ – Ben Leather Dec 2 '17 at 16:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.