I'm interested in numerically studying vortex solutions to PDEs. By this I mean solutions for PDEs for a field $\theta$ which satisfy $\oint\nabla\theta\cdot d\vec{\ell}=2\pi m$ where $m=...,-2,-1,0,1,2,...$ and $d\vec{\ell}=rd\phi$ is a vector along a curve made of a circle of fixed radius $r$ where $\phi\in(0,2\pi)$. A solution to this is $\theta=m\phi$. As a first step in this direction, I asked about the vortex solution to the Laplace equation for $\theta$ here: Numerically solving the Laplace equation in a 2d cylinder.

The next step I'm hoping to learn is how to code two coupled, nonlinear PDEs for two fields where one of them has two opposite vortex solutions. Hopefully, after understanding, this I could generalize the code and my understanding to other similar problems.

Here I've chosen, as a toy model, the stationary equations for two fields: $\theta$, whose vortex solutions are of interest, and an additional field $\rho$ that's coupled to $\theta$. The two nonlinear PDEs are $$0=\nabla\rho\cdot\nabla\theta+\frac{\rho}{2}\nabla^{2}\theta$$ and $$0=\nabla^{2}\rho(1-(\nabla\theta)^{2})-2\rho(\nabla^{2}\theta)^{2}.$$

These are derived from the Hamiltonian $\mathcal{H}=\int dr^{2}[\frac{1}{2}\rho^{2}((\nabla\theta)^{2}-1)]$.

To avoid singularities, I want to solve them on a domain that looks like this enter image description here

where the outer radius is $R$, the inner radii are the same and equal to $R_0$, the origin is dead center between the two inner circles and the distance from the origin to the center of either of the inner circles is $d$.

The boundary conditions I'm interested in are: $$\theta(x,y)=\arctan(\frac{y}{x})$$ for $(x-d)^2+y^2=R_{0}^2$ and $$\theta(x,y)=-\arctan(\frac{y}{x})$$ for $(x+d)^2+y^2=R_{0}^2$ and $$\rho(x,y)=\rho_{0}$$ for $x^2+y^2=R^2$ where $\rho_0$ is a constant (which can be set to 1). If additional conditions are required on $\rho$ on the two inner circles then these should be two different constants.

  • $\begingroup$ Have you tried anything yourself, at least to reproduce your equation in Mathematica code? $\endgroup$
    – MarcoB
    Feb 18 '19 at 23:34
  • $\begingroup$ How did this problem arise? It is necessary to indicate the beginning of the discussion on mathematica.stackexchange.com/questions/191456/… $\endgroup$ Feb 18 '19 at 23:55
  • $\begingroup$ Thanks @AlexTrounev, you're absolutely right. I've clarified the background in an edit to my original post. $\endgroup$
    – Asaf Miron
    Feb 19 '19 at 8:23
  • $\begingroup$ The boundary conditions are erroneous. Apparently, it is necessary to set the angle in local coordinates on each circle. $\endgroup$ Feb 19 '19 at 20:41
  • $\begingroup$ I wonder why people on the forum do not have the habit of voting for the questions and answers they read. The people at TeX.stackexchange.com are much more welcoming. $\endgroup$ Apr 15 '21 at 5:49

To solve a system of nonlinear equations with FEM we need to build a converging iterative process. I use the method of the false transient.The solution of the problem converges quickly, but at the first step instability arises.

R = 4; r = 1/2; d = 1; A = 
  x^2 + y^2 <= R^2 && (x - d)^2 + y^2 >= r^2 && (x + d)^2 + y^2 >= 
    r^2, {x, y}];

rho[0][x_, y_] := 1;
theta[0][x_, y_] := 0;
t0 = 1/5; k = 5;
Do[{theta[i], rho[i]} = 
NDSolveValue[{Laplacian[u[x, y], {x, y}] + 
    2*Grad[rho[i - 1][x, y], {x, y}].Grad[
        theta[i - 1][x, y], {x, y}]/rho[i - 1][x, y] == (u[x, y] -
       theta[i - 1][x, y])/t0, 
  Laplacian[v[x, y], {x, y}] - 
    2*rho[i - 1][x, 
      y]*(2*Grad[rho[i - 1][x, y], {x, y}].Grad[
            theta[i - 1][x, y], {x, y}]/rho[i - 1][x, y])^2/(1 - 
        Norm[Grad[theta[i - 1][x, y], {x, y}]]^2) == (v[x, y] - 
      rho[i - 1][x, y])/t0, 
  DirichletCondition[v[x, y] == 1, x^2 + y^2 == R^2], 
   u[x, y] == ArcTan[x - d, y], (x - d)^2 + y^2 == r^2], 
   u[x, y] == -ArcTan[x + d, y], (x + d)^2 + y^2 == r^2]}, {u, 
  v}, {x, y} \[Element] A, 
 Method -> {"FiniteElement", 
   "InterpolationOrder" -> {u -> 2, v -> 2}, 
   "MeshOptions" -> {"MaxCellMeasure" -> 0.001}}], {i, 1, 
k}]; // Quiet

{ContourPlot[theta[k][x, y], {x, y} \[Element] A, Contours -> 20, 
  ColorFunction -> Hue, PlotRange -> All, PlotLegends -> Automatic], 
 ContourPlot[rho[k][x, y], {x, y} \[Element] A, Contours -> 20, 
  ColorFunction -> Hue, PlotRange -> All, PlotPoints -> 50, 
  PlotLegends -> Automatic]}

Table[DensityPlot[theta[i][x, y], {x, y} \[Element] A, 
  ColorFunction -> Hue, PlotRange -> All], {i, 1, k}]
Table[DensityPlot[rho[i][x, y], {x, y} \[Element] A, 
  ColorFunction -> Hue, PlotRange -> All, PlotPoints -> 50], {i, 1, 


{Plot[theta[k][-d, y], {y, r, R}, PlotRange -> All, 
  AxesLabel -> {"y", "\[Theta]"}], 
 Plot[rho[k][-d, y], {y, r, R}, PlotRange -> All, 
  AxesLabel -> {"y", "\[Rho]"}], 
 Plot[theta[k][d, y], {y, r, R}, PlotRange -> All, 
  AxesLabel -> {"y", "\[Theta]"}], 
 Plot[rho[k][d, y], {y, r, R}, PlotRange -> All, 
  AxesLabel -> {"y", "\[Rho]"}]}


  • $\begingroup$ What are you using "InterpolationOrder" -> {u -> 2, v -> 2} for? $\endgroup$
    – user21
    Feb 20 '19 at 7:08
  • $\begingroup$ What version are you using? $\endgroup$
    – user21
    Feb 20 '19 at 7:21
  • 1
    $\begingroup$ @user21 In[1]:= $Version Out[1]= "11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)". My algorithm needs an option "InterpolationOrder" -> {u -> 2, v -> 2} $\endgroup$ Feb 20 '19 at 16:30
  • 1
    $\begingroup$ @AsafMiron I updated the code, corrected Norm[]->Norm[]^2 and added fig. 2. $\endgroup$ Feb 20 '19 at 17:10
  • 1
    $\begingroup$ @AsafMiron 1) I chose the parameters based on the experiences. 2) There are methods for accelerating convergence. In this forum, I posted several examples of the use of the method of the false transient, see mathematica.stackexchange.com/… $\endgroup$ Feb 21 '19 at 16:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.