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: