Message Boards Message Boards

Plot a hamiltonian density for gauged vortices?

Posted 6 years ago

I have numerically solved two coupled differential equations for Gauged Vortices in (2+1) dimensions using the following code:

SolveGaugedGLVortices[k_Integer, inf_?NumericQ, tol_?NumericQ, 
   damping_?NumericQ] := 
  Block[{xg, xg1, eye, dX, d2X, mu, mat, lfun, ysol, zsol, rhs1, rhs2,
     ysol1, zsol1, err, i, residual1, residual2, yif, zif},
   xg = N@Range[0, inf, inf/12000];
   xg1 = xg + 10^-50.;
   {eye, dX, d2X} = 
    NDSolve`FiniteDifferenceDerivative[#, xg, "DifferenceOrder" -> 6][
       "DifferentiationMatrix"] & /@ {0, 1, 2};
   mu = -Abs[N[damping]];
   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., k};
    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 <= Abs[tol], Break[]];
    If[Mod[iter, 10000] == 0, dPrint["err ->", err]];
    ysol = ysol1;
    zsol = zsol1;
    i = iter;
    , {iter, 1000000}
    ];
   residual1 = 
    xg*(d2X.ysol) + dX.ysol - (ysol*((k - zsol)^2))/xg1 + 
     xg*(1 - ysol^2)*ysol/2;
   residual1[[{1, -1}]] = 0.;
   residual2 = xg*(d2X.zsol) - dX.zsol + xg*(k - zsol)*ysol^2;
   residual2[[{1, -1}]] = 0.;
   {yif, zif} = 
    Interpolation[Transpose[{xg, #}], Method -> "Spline"] & /@ {ysol, 
      zsol};
   {{yif, zif}, {Mean[Abs[residual1]], Mean[Abs[residual2]]}, mu, i}
   ];

sol = Table[SolveGaugedGLVortices[k, 1000, 10^-7, 3*(10^k)]
   , {k, 3}];

{ysol, zsol} = Transpose[sol[[All, 1]]];

to yield solutions for k=1,2,3 for the two variables y and z in the lists ysol and zsol respectively.

Now I want to plot the following function, the hamiltonian density:

hdens = (f'[x])^{2}/x^2 + (u'[x])^{2} + (k - f[x])^{2}*
u[x]/x^2 + (1 - u[x]^2)^{2};

for each value of k using the respective numerical solutions such that f = zsol and u = ysol; i.e. for k = 1, f = ysol[[1]][[1]] and y = zsol[[1]][[1]].

I would like some way of iterating over each list and plotting the hamiltonian density for each k value on the same plot. I've also attached my Mathematica Notebook for further reference. Any help on this matter would be greatly appreciated.

Attachments:
POSTED BY: Benjamin Leather
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