Message Boards Message Boards

GROUPS:

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

Posted 1 year ago
1300 Views
|
8 Replies
|
5 Total Likes
|

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]]
8 Replies

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.

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.

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.

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.

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

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:

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:

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:
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