# Spectral problem for differential vector operator (calculation of EM field in a cavity)

I know that mathematica has a DEigensystem and NDEigensystem which allow one to find eigenfunction and eigenvalue for scalar differential operator. But I want to find eigenfunction and eigenvalue for vector operator with additional relation between component of that wector field.

I need to solve Maxwell's equation in arbitrary shape cavity. $$\Delta \bf{E}=-\omega^2 \bf{E}\quad \bf{\nabla}\cdot \bf{E}=0$$ With boundary conditions $\bf{E}_{\tau}=0$ in a boundary of a cavity, where $\bf{E}_{\tau}$ is a tangent vector of a cavity boundary.

I don't now how to correctly take into account the equation $\bf{\nabla}\cdot \bf{E}=0$. And how two express the boundary condition correctly?

I can express one component of vector field $\bf E$ through another two, but in this case the boundary condition take a more complex form. This boundary conditions can not be expressed via DirichletCondition and NeumannValue.

Edit:

I can express one component of the vector field through another two. For example $$E_z=-\int dz(\partial_x E_x+\partial_y E_y)+g(x,y)$$ Function $g(x,y)$ should obey the same equation $\Delta g=-\omega^2 g$.

The boundary condition can be rewrited in the following form $\bf{E}_\tau=\bf{E}-(\bf{E}\cdot \bf{n})\,n=0$ where $\bf{n}$ is a unit vector perpendicular to the boundary.

It is possible to solve it for rectangle. But it does not make sense because in such case there is a analytical solution. I don't know how to realize this boundary condition in the mathematica.

I agree, that the additional conditions $\bf{\nabla}\cdot \bf{E}=0$ it is a not big problem here. But the boundary conditions it is a big problem, at least for me.

Edit2 For simplisity I try to solve the same 2 dimentional problem. If the boundary can be set in the following form $y=g(x)$, one can calcuete the tangent vector analyticaly $\bf{\tau}=\bf{e}_x+g'(x)\bf{e}_y$. I try the realize it in the Mathematica:

gg[x_] := 1 - x^2/2
region = ImplicitRegion[x > -1 && x < 1 && y < gg[x] && y > -gg[x], {x, y}];
NDEigensystem[{-Laplacian[{ux[x, y], uy[x, y]}, {x, y}],
DirichletCondition[uy[x, y] == 0, x == -1],
DirichletCondition[uy[x, y] == 0, x == 1],
DirichletCondition[{ux[x, y], uy[x, y]}.{1, D[gg[x], x]} == 0,
y > 0 && x > -1 && x < 1],
DirichletCondition[{ux[x, y], uy[x, y]}.{1, -D[gg[x], x]} == 0,
y > 0 && x > -1 && x < 1]}, {ux[x, y],
uy[x, y]}, {x, y} \[Element] region, 6]


But it does not work.

Also for solution inside 3D sphere

region = ImplicitRegion[x^2 + y^2 + z^2 <= 1, {x, y, z}];
g1 = (x* Ex[x, y, z] + y *Ey[x, y, z] + z* Ez[x, y, z])/(x^2 + y^2 + z^2);
NDEigensystem[{-Laplacian[{Ex[x, y, z], Ey[x, y, z], Ez[x, y, z]}, {x,
y, z}], DirichletCondition[{Ex[x, y, z] - x*g1 == 0,
Ey[x, y, z] - y*g1 == 0, Ez[x, y, z] - z*g1 == 0}, True]}, {Ex[x,  y, z], Ey[x, y, z], Ez[x, y, z]}, {x, y, z} \[Element] region, 6]


It also does not work.

• Have you tried to code anything? If so you could share that. – user21 Feb 16 '18 at 7:46
• I give up. This is how far I got.op = Laplacian[{Ex[x, y], Ey[x, y]}, {x, y}] + {Ex[x, y], Ey[x, y]}; bc = DirichletCondition[{Ex[x, y] == 0, Ey[x, y] == 0}, True]; {ev, ef} = NDEigensystem[{op, bc}, {Ex[x, y], Ey[x, y]}, {x, y} \[Element] Disk[], 9]; VectorPlot[#, {x, y} \[Element] Disk[]] & /@ ef I tried to impose the boundary condition by defining a normal vector and using the skalar product. But the DirichletCondition is only allowed to be linear. These lecture notes may be helpful. – Matthias Bernien Feb 17 '18 at 13:04
• @MatthiasBernien Does not your Dirichlet condition impose $E_x = E_y = 0$ on the whole domain? – anderstood Feb 18 '18 at 15:38
• @anderstood Yes. What I wanted to do was nvec={x,y}; DirichletCondition[{Ex[x, y], Ey[x, y]}.nvec/(Norm[nvec]Norm[{Ex[x, y], Ey[x, y]}])==1, True]. But that did not work. – Matthias Bernien Feb 18 '18 at 19:42
• @Peter Besides adding the second component of the field everywhere you need to split the eigenvectors into two halfs:evIFx = ElementMeshInterpolation[{mesh}, #[[1 ;; 539]]] & /@ eigenVectors; evIFy = ElementMeshInterpolation[{mesh}, #[[540 ;; -1]]] & /@ eigenVectors;. Then use VectorPlot. – Matthias Bernien Feb 22 '18 at 7:04