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:
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.
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?