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
7 days 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
7 days 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
2 days ago

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

POSTED BY: Paritosh Mokhasi
Answer
2 days 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
1 day 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
1 day ago

Group Abstract Group Abstract