finite difference frequency domain Eigenvalue matrix to get eigenmodes

Question

I am using the FDFD method to calculate eigenmodes for an empty wave guide. It's a 1x1meter structure with a PEC boundary. (Here I have 6x6 points to make it simple). Sounds simple, but I can't get it right.

At each mesh point my source free E-Field is allocated and calculated using the Helmholtz-equation: $(∇^2-\lambda)E(x,y) = 0$

And since I can use a central differential quotient $∇E_{ij}(x,y) = (-4E_{ij} +E_{i+1,j}+E_{i-1,j}+E_{i,j+1}+E_{i,j-1})/\Delta x^2$

I am getting a matrix like this one (on the central diagonal the values are $4$ and the rest is of the diagonals are $-1$):

So far so good. But now I need to apply 2 changes:

1. eliminate non existing points (on the boundary each point has only 3x dots of the central differential quotient and not 4, so one of them must be set to 0). And on the edges there are only 2x dots instead of 4.

2. Apply PEC for which $ E_{tan}(boundary)=0$ which means that on the boundary all tangetial fields must be equal 0.

Doing so I get the "1. A-Matrix":

But here comes my problem:

Eigenvalues (for me here Eigenmodes) that I get are not the ones, that I calculate analytically. My guess is that the 1st or 2nd adjustment that I am applying is somehow wrong. Could you please tell me whether I am apply my matrix adjustments and PEC implementation correctly?

UPDATE

Thanks to Peter Frolkovič, I noticed that there were some mistakes in my notations. I corrected the $∇$ to $∇^2$ and the $E_k$ to $E_{ij}$. Also regarding the boundary conditions:

Since it's PEC I'm after, the E-Field is 0 inside PEC. So naturally all E-Field components on all boundary points should be 0. But if I do so, I really have to set all boundary points to zero in the Matrix A and this gives me the "2. A-Matrix" (see the 2nd image above)

But solving the eigenvalue problem with this matrix results in some very weird Eigenvalues. So that's why I thought that not all E-Field components on the boundary should be 0 but only the tangetial ones. So I only set $ E_{tan}(boundary)=0$. But as described this didn't help much.

Here an Example with my Grid (see first image above):

In the above image you can see how my grid is defined (blue numbers are the node numbers). As described, at each node I solve the Helmholtz PDE from above. Lets take node number $7$ as an example. The corresponding equation would be:

$(∇^2-\lambda)E_7 = 0$

with $∇^2E_7 = (-4E_7+E_1+E_8+E_{13})/ \Delta x^2$ (here are only 3 nodes of the central differential quotient, since node number 7 doesn't have a left node) $(∇^2-\lambda)E_7 = 0$

$(-4E_7+E_1+E_8+E_{13}) = \Delta x^2 \lambda E_7 $

Now for Dirichlet $E_7 =0$ which means that the equation simplifies to $(E_1+E_8+E_{13}) = 0 $

Do I understand it right? If so, then the only way I know how to implement this in an eigevalue equation is by editing the main matrix - so by setting the corresponding matrix value (node) to 0... so in MATLAB $A(7,7)=0$.

As Peter Frolkovič pointed out, matrix entries are only the coefficients of the unknowns - I agree, but how else can I influence the node $E_7$ in this case, if not by doing as I described?

Answer 1

Looking to your "2. A-Matrix" my first guess is that you might misunderstand the implementation of so called Dirichlet boundary conditions when one prescribes given values for the unknown solution at boundary. For instance the zero Dirichlet boundary conditions mean $$ E_{i,j}=0 $$ for some indices $(i,j)$ for which the corresponding grid points lie on boundary.

One of the standard treatment is that you view the above condition as a linear algebraic equation (although a trivial one), so you include it into your system of linear algebraic equations, so a corresponding row and column must be present in the matrix. This row has a trivial form having $1$ on the diagonal and zero elsewhere. The other rows are not influenced by this boundary condition! Especially it does not mean that you replace in the matrix the entries corresponding to boundary nodes by zeroes, the matrix entries are defined by coefficients before the unknowns, not by the unknowns themselves.

So my suggestion is that you check this issue in your code. In "2. A-Matrix" it seems that you have zero rows that can not happen.

The other standard approach to treat Dirichlet boundary conditions is that you eliminate these values, it means you do not include them into unknowns, so no rows and columns for them exist in the matrix. If you have the zero Dirichlet boundary conditions, it means simply that you delete the rows and columns that would correspond to these values. In general you modify all linear algebraic equations for $i,j$ when at least one neighbor lies on the boundary (but not $i,j$-node itself, it is inner node). You might follow this approach in the "1. A-Matrix", but I am not sure if completely.