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.

Improve this question
$\endgroup$
2
  • $\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, 2017 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, 2017 at 1:51

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
12
  • $\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, 2017 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, 2017 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, 2017 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, 2017 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, 2017 at 16:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.