Message Boards Message Boards

Get solutions to the gauged GL vortices in (1+1) dimensions?

Posted 6 years ago

I am trying to numerically solve the following differential equations for the profiles of the Gauged GL Vortices:

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])*y[x]^2 == 0;

with the boundary conditions:

inf = 1000;
bc = {y[0] == 0, y[inf] == 1, z[0] == 0, z[inf] == 1};

I have tried a finite difference method but this seems to give large errors even for k = 1 so I am wondering if there is a better numerical method to yield a more accurate solution?

Below is my current function which I employ to solve these coupled differential equations:

Clear[fdd, pdetoode, tooderule, sollst]
fdd[{}, grid_, value_, order_] := value;
fdd[a__] := NDSolve`FiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] := 
  pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], 
   rest];
pdetoode[{func__}[var__], rest__] := 
  pdetoode[Alternatives[func][var], rest];
pdetoode[rest__, grid_?VectorQ, o_Integer] := 
  pdetoode[rest, {grid}, o];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer] := 
  With[{pos = Position[{var}, time][[1, 1]]}, 
   With[{bound = #[[{1, -1}]] & /@ {grid}, 
     pat = Repeated[_, {pos - 1}], 
     spacevar = Alternatives @@ Delete[{var}, pos]}, 
    With[{coordtoindex = 
       Function[coord, 
        MapThread[
         Piecewise[{{1, # === #2[[1]]}, {-1, # === #2[[-1]]}}, 
           All] &, {coord, bound}]]}, 
     tooderule@
      Flatten@{((u : func) | 
            Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat, 
          t_, x2___] :> (Sow@coordtoindex@{x1, x2};

          fdd[{dx1, dx2}, {grid}, 
           Outer[Derivative[dt][u@##]@t &, grid], 
           "DifferenceOrder" -> o]), 
        inde : spacevar :> 
         With[{i = Position[spacevar, inde][[1, 1]]}, 
          Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] := 
  Equal[tooderule[rule][a - b], 0] //. 
   eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@ 
  Reap[expr /. rule]

Clear@pdetoae;
pdetoae[funcvalue_List, rest__] := 
  pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] := 
  pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] := 
 Module[{t}, 
  Function[pde, #[
       pde /. {Derivative[d__][u : func][inde__] :> 
          Derivative[d, 0][u][inde, t], (u : func)[inde__] :> 
          u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
   pdetoode[func[var, t], t, rest]]
POSTED BY: Benjamin Leather
8 Replies

Thank you for response. A further extension to this that I have been attempting is to produce plots of the magnetic field, i.e. plotting the following function based on the numerical solutions: bfield = f'[x]/x;

My attempt so far has yielded interesting results that seem to break down at zero.

bfield = f'[x]/x;
bfieldplot = 
 Plot[Table[bfield, {f, zsol}] // Evaluate, {x, 0, 10}, 
  PlotRange -> {0, 0.6}, 
  PlotStyle -> 
   Table[Blend[{Blue, Cyan}, x], {x, 0, 1, 1/(Length[zsol] - 1)}], 
  GridLines -> Automatic, AxesLabel -> {r, B}]

Plot of magnetic field from numerical solutions

Is there any way to rectify this?

I have also attached my current Mathematica notebook for further reference.

Attachments:
POSTED BY: Benjamin Leather
Attachments:
POSTED BY: Benjamin Leather
Attachments:
POSTED BY: Paritosh Mokhasi
POSTED BY: Benjamin Leather

Benjamin, could you please add the code that you ran using your functions?

POSTED BY: Paritosh Mokhasi

The following could be used as a starting point for solving the system numerically:

Clear["Global`*"];
k = 2.;
inf = 1000.;
xg = N@Range[0, inf, inf/12000];
xg1 = xg + 10^-50.;
{eye, dX, d2X} = NDSolve`FiniteDifferenceDerivative[#, xg, "DifferenceOrder" -> 6]["DifferentiationMatrix"] & /@ {0, 1, 2};
mu = -3*10^2.;
mat = d2X + mu*eye;
mat[[{1, -1}]] = eye[[{1, -1}]];
lfun = LinearSolve[mat];

ysol = zsol = Erf[xg/3];

Do[
 rhs1 = (((k - zsol)^2)*ysol/xg1^2) - ((1 - ysol^2)*ysol)/2 - (dX.ysol)/xg1 + mu*ysol;
 rhs1[[{1, -1}]] = {0., 1.};
 rhs2 = -((k - zsol)*ysol*ysol) + (dX.zsol)/xg1 + mu*zsol;
 rhs2[[{1, -1}]] = {0., 1.};

 ysol1 = lfun[rhs1];
 zsol1 = lfun[rhs2];

 If[Max[Abs[ysol1]] >= 10 || Max[Abs[zsol1]] >= 10, Print["unstable.."]; Break[]];
 err = Max[Norm[ysol1 - ysol], Norm[zsol1 - zsol]];
 If[err <= 10^-7, Break[]];
 If[Mod[iter, 1000] == 0, Print["err ->", err]];
 ysol = ysol1;
 zsol = zsol1;
 , {iter, 1000000}]

The parameter mu will have to be changed as the parameter k is changed. A larger value of k would require a larger value of mu. Also, for larger k-values, the convergence associated with these iterations is slow.

POSTED BY: Paritosh Mokhasi

Dear Paritosh Mokhasi,

Is there any description/explanation of this method of solving the differential equations you could provide? For example a paper or link explaining it.

POSTED BY: Benjamin Leather

The method I am using here is called Successive Over Relaxation (SOR). There are plenty of references that you can find on the web on this. My implementation is a very simple version. There are much more sophisticated methods in SOR that could be used.

POSTED BY: Paritosh Mokhasi
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract