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$):

enter image description here

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":

enter image description here

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?


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?

  • $\begingroup$ I suppose that your Partial Differential Equation (PDE) and its approximation is standard to be sure, because some of your notations are nonstandard (you might mean $\Delta$ or $\nabla^2$ instead of $\nabla$, the discrete (approximative) values for $E$ are mostly denoted with always two indices, e.g. $E_{i j}$, or always with one index, e.g. $E_k$). I think the community could help you if you can specify exactly how the boundary conditions are defined and approximated for all 4 edges of domain. I personally can not recognize it only from the picture of nonzero entries in matrix. $\endgroup$ Jul 5, 2016 at 11:06
  • $\begingroup$ Thank you @PeterFrolkovič - you are right. I have made some corrections and added some more description to my boundary definition. is this enough? If not, I can also post the Matlab code $\endgroup$
    – Kosha Misa
    Jul 5, 2016 at 12:01
  • $\begingroup$ Is E the electric field as you seem to imply or the electrical potential? If the former, why aren't there two equations at each mesh point? $\endgroup$ Jul 6, 2016 at 11:43
  • $\begingroup$ E is the Electric field. Do you mean 2x equations for each spatial direction (x and y) ? If so, then only one direction should be enough to calculate the right cuttoff frequencies $\endgroup$
    – Kosha Misa
    Jul 6, 2016 at 17:10

1 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.

  • $\begingroup$ Thanks a lot, it seems that I really might have misunderstood how to implement the PEC (Dirichlet). To calculate eigenvalues and vectors I am using a MATLAB method which takes the main matrix as an input. So I need to modify it. Please take a look at my example in the main description (I made an update at the end of the last update) $\endgroup$
    – Kosha Misa
    Jul 6, 2016 at 18:31
  • $\begingroup$ Im new at this and wanted to grand you a +1 but I only can check your post as answer $\endgroup$
    – Kosha Misa
    Jul 6, 2016 at 19:28
  • $\begingroup$ Your considerations for the node 7 are not correct. If there is no left neighbor it does not mean you simply skip it, you can make it if that left point is a node where you know the value of $E$ is zero. But then you suppose that $E$ is zero not only at the boundary, but also outside of domain (and that is not the case for eigenfunction of your test example). So take simply A(7,7)=1. The influence of the fact that $E_7=0$ will occur in the discrete equation for the node 8, but you need not to modify the row for node 8 because of that. $\endgroup$ Jul 7, 2016 at 6:07
  • $\begingroup$ Again taking my example from above: If I don't set the matrix entry (the coefficient) of the left neighbor to 0 the equation for this row will be: $(−4E_7+E_1+E_8+E_{13}+E_6)=Δx^2λE_7$ And due to my understanding this $E_6$ should not be there. Now if I set $A(7,7)=1$ the equations becomes $(E_7+E_1+E_8+E_{13}+E_6)=Δx^2λE_7$ but the $E_6$ value and coefficient stays there. Please help me understand where is my mistake - should I just ignore the fact that $E_6$ is popping up in this equation and if so why? $\endgroup$
    – Kosha Misa
    Jul 7, 2016 at 7:40
  • $\begingroup$ As you proposed I've modified my matrix. I've also included $E_x$ & $E_y$ in my matrix $A=[E_x, sparse(np,np) ....sparse(np,np), E_y]$ (so my matrix is 2*np) and changed the diagonal coefficients to $1$ for $E_x$ on the top and bottom and for $E_y$ for the left and right boundaries and now I do get my lowest cut off frequency (at about $1.49Mhz$) and even the E-Field distribution looks almost right (it looks rotated by $90°$) but It still don't get why it works. See my question from the prev. comment $\endgroup$
    – Kosha Misa
    Jul 7, 2016 at 7:47

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.