Message Boards Message Boards

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

GROUPS:

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
Answer
6 months ago

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
Answer
6 months ago

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
Answer
5 months ago

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
Answer
5 months ago

Unfortunately, the numerical solutions associated with your method are actually worse than my original finite difference method which I have shown above, rather than better.

POSTED BY: Benjamin Leather
Answer
5 months ago

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

POSTED BY: Paritosh Mokhasi
Answer
5 months ago
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]]

Clear@k;
Clear@y;
Clear@z;
Clear@\[Phi];
inf = 1000;
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;
eqns = Join[{eqn1, eqn2}];
bc = {y[0] == 0, y[inf] == 1, z[0] == 0, z[inf] == k};
set = Join[eqns, bc];

points = 1000;
domain = {0, inf};
grid = Array[# &, points, domain];
difforder = 4;
(*Definition of pdetoode is above.*)

ptoafunc = pdetoae[{y[x], z[x]}, grid, difforder];

delete = #[[2 ;; -2]] &;
ae = delete /@ ptoafunc[Map[x^2 # &, eqns, {2}] // Simplify];
WorkingPrecision -> 30;

sollst = Table[
  ListInterpolation[#, grid] & /@ Partition[#, points] &@
   With[{initial = 1}, 
     FindRoot[{ae, bc}, 
      Flatten[{{y@#, initial} & /@ grid, {z@#, initial} & /@ grid}, 
       1]]][[All, -1]], {k, 3}]
Attachments:
POSTED BY: Benjamin Leather
Answer
5 months ago

The boundary conditions that you had presented for z[inf] was 1, but in your notebook I see that it is dependent on k. That might have been one reason for the discrepancy. The attached notebook shows the same tests that you performed and looks at the residual errors, and it seems to be as expected.

Attachments:
POSTED BY: Paritosh Mokhasi
Answer
5 months ago

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
Answer
5 months ago

Group Abstract Group Abstract