I am trying to use NDEigensystem to find eigenvalues and eigenmodes of the following Laplacian with boundary conditions:

$$\vec{\nabla}^2 \vec{A}(x,z) = vals \; \vec{A},$$

$$A_{x}(x,z=\pm \frac{d}{2}) = 0,$$

$$A_{y}(x,z=\pm \frac{d}{2}) = A_{y}(x = \pm r,z) = 0,$$

$$A_{z}(x = \pm r,z) = 0,$$

$$\frac{\partial A_{x}(x,z = \pm d/2)}{\partial x} = 0,$$

$$\frac{\partial A_{z}(x = \pm r,z)}{\partial z} = 0,$$

where $\vec{A} = \{A_x,A_y,A_z\}$ and $vals$ is an eigenvalue found by NDEigensystem. The domain is a rectangle centred at $(x,z) = (0,0)$ and with basis $2r$ and height $d$. Importantly, to get these equations I used the fact that

$$\nabla \cdot \vec{A} = 0,$$

and as such I would like the solutions to respect this constraint. However, for (some) of the solutions found by NDEigensystem, $\nabla \cdot \vec{A} \neq 0$ (as you can see in the last plot if you ran the code below). Luckily, these solutions come in groups with the same eigenvalue, and a superposition of them actually can give null divergence. I am confident that I am using NeumannValue in the wrong way. I also tried to impose the null divergence directly in the operator by setting $\partial_x A_x = -\partial_z A_z$, but in that case NDEigensystem was giving a completely crazy solution. In particular, the boundary conditions involving the derivatives were not satisfied anymore (this further suggests that there is a problem with NeumannValue - or better with me trying to use NumannValue).

My question is the following: How can I impose null divergence of the solution? How does NeumannValue works, exactly, in case of coupled differential equations? I have tried to read documentation about NeumannValue, but unfortunately I am a lame physicists and have really not understood if and how it can be used to impose null divergence of a field in my case or in general.

As a small comment, I am aware that my problem can be solved analytically in a relatively simple manner. The reason for which I want to have a numerical solution is that I would like to increase the complexity once I am sure that I understood how NDEigensystem and NeumannValue work.

I have a "bonus question" here:

  • There is a way to change the normalization of the eigenfunctions? I have read that NDEigensystem spits out eigenfunctions $\vec{A}_i$ s.t.: $\int \vec{A}_i^* \cdot \vec{A}_j dxdz = \delta_{i,j}$. I would like to change that to $\int \epsilon_r(x,z) \vec{A}_i^* \cdot \vec{A}_j dxdz = \delta_{i,j}$, where $\epsilon_r(x,z)$ is a certain function. I can do that by hand afterwards, but it seems that there is a method option in NDEigensystem that does that (called VectorNormalization), but I did not understand how it works...

In the following, the block of the code that I am currently using:

d = 3;
r = 1;
\[CapitalOmega] = Rectangle[{-r, -(d/2)}, {r, d/2}];
Graphics[{Opacity[.3], Blue, \[CapitalOmega]}, Axes -> True, 
 AxesLabel -> {"x", "z"}]
(*Differential operator:*)
EqSt[x_, y_, z_] := 
 Laplacian[{Ax[x, z], Ay[x, z], Az[x, z]}, {x, y, z}, "Cartesian"]
(*Dirichlet and Neumann boundary conditions:*)
BndCnd = {DirichletCondition[Ax[x, z] == 0, z^2 == (d/2)^2], 
   DirichletCondition[Ay[x, z] == 0, Or[x^2 == r^2, z^2 == (d/2)^2]], 
   DirichletCondition[Az[x, z] == 0, x^2 == r^2]};
NeumBnd = {NeumannValue[0, x^2 == r^2], 0, 
   NeumannValue[0, z^2 == (d/2)^2]};
(*Getting solutions:*)
{valsZ3, funsZ3} = 
   Flatten[{EqSt[x, y, z] - NeumBnd, BndCnd}], {Ax, Ay, 
    Az}, {x, z} \[Element] \[CapitalOmega], 16];
(*Checling the solutions:*)
ii = 3;
{Plot3D[funsZ3[[ii, 1]][x, z], {x, -r, r}, {z, -(d/2), d/2}, 
  PlotRange -> All, AxesLabel -> {"x", "z"}, PlotLabel -> "Ax[x,z]"],
 Plot3D[funsZ3[[ii, 2]][x, z], {x, -r, r}, {z, -(d/2), d/2}, 
  PlotRange -> All, AxesLabel -> {"x", "z"}, PlotLabel -> "Ay[x,z]"],
 Plot3D[funsZ3[[ii, 3]][x, z], {x, -r, r}, {z, -(d/2), d/2}, 
  PlotRange -> All, AxesLabel -> {"x", "z"}, PlotLabel -> "Az[x,z]"],
   Div[{funsZ3[[ii, 1]][x, z], funsZ3[[ii, 2]][x, z], 
     funsZ3[[ii, 3]][x, z]}, {x, y, z}]], {x, -r, r}, {z, -(d/2), d/
   2}, PlotRange -> All, AxesLabel -> {"x", "z"}, 
  PlotLabel -> 
   "\!\(\*SubscriptBox[\(\[Del]\), \({x, y, \
  ColorFunction -> Function[{x, y, z}, Hue[.65 (1 - z)]]]}

Your Answer

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

Browse other questions tagged or ask your own question.