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, 2019 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, 2019 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, 2019 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, 2019 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, 2021 at 5:49

1 Answer 1


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, 2019 at 7:08
  • $\begingroup$ What version are you using? $\endgroup$
    – user21
    Feb 20, 2019 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, 2019 at 16:30
  • 1
    $\begingroup$ @AsafMiron I updated the code, corrected Norm[]->Norm[]^2 and added fig. 2. $\endgroup$ Feb 20, 2019 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, 2019 at 16:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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