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

Posted 9 months ago
1034 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__] := NDSolveFiniteDifferenceDerivative@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
Sort By:
Posted 9 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} = NDSolveFiniteDifferenceDerivative[#, 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 8 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 8 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 9 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 9 months ago
Posted 8 months ago
 Clear[fdd, pdetoode, tooderule, sollst] fdd[{}, grid_, value_, order_] := value; fdd[a__] := NDSolveFiniteDifferenceDerivative@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:
 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}] Is there any way to rectify this? I have also attached my current Mathematica notebook for further reference. Attachments: