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.


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.

  • 3
    $\begingroup$ Have you tried to code anything? If so you could share that. $\endgroup$
    – user21
    Feb 16, 2018 at 7:46
  • $\begingroup$ 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. $\endgroup$ Feb 17, 2018 at 13:04
  • $\begingroup$ @MatthiasBernien Does not your Dirichlet condition impose $E_x = E_y = 0$ on the whole domain? $\endgroup$
    – anderstood
    Feb 18, 2018 at 15:38
  • $\begingroup$ @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. $\endgroup$ Feb 18, 2018 at 19:42
  • 1
    $\begingroup$ @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. $\endgroup$ Feb 22, 2018 at 7:04


Your Answer

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

Browse other questions tagged or ask your own question.