# Numerically solving two coupled nonlinear PDEs

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

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.

• Have you tried anything yourself, at least to reproduce your equation in Mathematica code? – MarcoB Feb 18 at 23:34
• How did this problem arise? It is necessary to indicate the beginning of the discussion on mathematica.stackexchange.com/questions/191456/… – Alex Trounev Feb 18 at 23:55
• Thanks @AlexTrounev, you're absolutely right. I've clarified the background in an edit to my original post. – Asaf Miron Feb 19 at 8:23
• The boundary conditions are erroneous. Apparently, it is necessary to set the angle in local coordinates on each circle. – Alex Trounev Feb 19 at 20:41

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 =
ImplicitRegion[
x^2 + y^2 <= R^2 && (x - d)^2 + y^2 >= r^2 && (x + d)^2 + y^2 >=
r^2, {x, y}];
DiscretizeRegion[A]

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}] +
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,
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],
DirichletCondition[
u[x, y] == ArcTan[x - d, y], (x - d)^2 + y^2 == r^2],
DirichletCondition[
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,
k}]


{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]"}]}


• What are you using "InterpolationOrder" -> {u -> 2, v -> 2} for? – user21 Feb 20 at 7:08
• What version are you using? – user21 Feb 20 at 7:21
• @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} – Alex Trounev Feb 20 at 16:30
• @AsafMiron I updated the code, corrected Norm[]->Norm[]^2 and added fig. 2. – Alex Trounev Feb 20 at 17:10
• @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/… – Alex Trounev Feb 21 at 16:08