The method of finite differences for computational electromagnetics

The partial differential equation of Helmholtz is approximated by a system of linear algebraic equations by using finite differences. It is then easily solved by matrix methods. Although some memory management is involved for 3-d problems, the math is very similar:

∂²φ/∂x²+∂²φ/∂y²+ε_dk²φ=0
φ_x-1,y+φ_x+1,y+φ_x,y-1+φ_x,y+1+(h²ε_dk²-4)φ_x,y=0

where

• h: integration step
• k: wave number (2π/λ)
• λ: (free space) wavelength

Incidentally, the formulation using finite elements is quite similar, but triangular elements are involved in two dimensions, tetrahedrons in three.

The Dirichlet boundary condition φ=0 is pertinent to E waves, while the Von Neumann boundary condition, ∂φ/∂n=0, applies to H fields.

Symmetry should be exploited at all times in field problems: For a 3-d problem, only one eighth of the work may have to be done. A higher resolution is alternatively available for the same amount of work.

As quantisation is involved both in the signal domain and maybe in the spatial domain, if part of the boundary is curved, convergence must be verified: The number of elements is increased until two consecutive results are much the same. Extra elements are only necessary in regions of rapid field variation- but will not help much near singularities.

The solution to the 1-d problem generally really isn't worthy of finite differences, but the simulation error can be calculated easily, since analytical results are available at once. This program provided the results shown in the following table, corresponding to a simple rectangular waveguide of cross section length a in the TE₁₀ mode:

Elements/ (2x/a)	1	.8	.6	.4	.2	0
6	1	.953	.811	.589	.310	0
18	1	.951	.809	.588	.309	0
sin(πx/a)	1	.951	.809	.588	.309	0

The cutoff wavelength for this waveguide is plain to see, but for more sophisticated waveguides, it will have to be calculated: A thirty-liner routine will (usually) provide the eigenvalue of largest absolute value, but this is not useful in terms of waveguides. This should have corresponded to a multiple of the cutoff frequency, but doesn't, because of the inaccuracy introduced by quantisation. On the other hand, the smallest (non zero) eigenvalue does yield the cutoff wave number. Therefore, the matrix is inverted, as well as the resulting eigenvalue.

The next example is a little more interesting: A dielectric loaded, rectangular waveguide: The uniform dielectric is centrally placed along 33% of the cross section: Now this next program gave the results in the following table for λ/(2a):

Elements/ εd	1	3	5	10
6	1.003	1.42	1.75	2.38
12	1.0007	1.46	1.82	2.49
18	1.0003	1.48	1.84	2.53
Extrapolation	1.0002	1.50	1.85	2.55
Probable value	(1)	1.50	1.91	2.64

The results have been improved by using Aitken's δ² process:

x₄=(x₃x₁-x₂x₂)/(x₃-2x₂+x₁)

There are other methods for finding the wave amplitude and cutoff frequency: For example, an iterative method will start with plausible values for the field amplitude and cutoff wavenumber, and alternatively improve upon them, also exploiting the Rayleigh quotient, but this has not been tried.

It doesn't help that for H-modes zero frequency is also a solution (corresponding to constant field within the waveguide). The trick is to employ an educated guess for the eigen wavenumber in the system matrix: The double inversion is again in order, but now the 'eigenvalue' found is really the difference between the assumed and actual values.

This so called inverse method is taken advantage of in the ridge waveguide program by which the data in the following table was procured (ridge dimensions d x d/2):

Elements	kd
35	2.32
59	2.29
89	2.28
Extrapolation	2.27
Reference	2.25

The field values for the ridge waveguide dominant H-mode are shown on a 9x9 grid (for 1/2 of the cross section), but have actually been calculated on a 69x69 grid. Otherwise, there is no difference between the program used and the program shown.

797	784	745	681	592	481	355	217	73
809	796	758	693	603	489	360	219	73
832	820	783	720	627	506	368	223	74
864	854	822	765	670	528	377	227	75
901	893	870	830	777	573	392	231	76
938	933	918	898	882
969	966	958	947	941
990	988	983	977	974
1000	996	991	987	984

The cutoff wavelength and field amplitude for the circular waveguide ( diameter d) have been obtained in much the same way: The eigenCirc program verified the reference value for the fundamental frequency, while the circular program created the figures in the next table: Again, one quarter of a cross section is shown on a 9x9 mesh, but the program calculations for the principal TE mode were performed on a 67x67 mesh (for the entire circle). In every other respect, the program supplied is exactly like the program used:

Elements	kd
51	3.48
68	3.54
96	3.56
Extrapolation	3.57
Reference	3.682

And the wave amplitude:

						491	999	1000
				881	926	960	981	988
			791	842	882	913	931	937
		686	731	772	806	831	847	853
	563	599	637	670	698	719	732	736
	452	485	515	541	562	578	588	592
287	326	350	371	389	404	416	423	425
165	185	199	211	222	230	237	240	242
35	38	41	43	45	47	48	49	49

When an antenna is not long in terms of the wavelength, further techniques include the finite difference time domain, transmission line matrix (TLM) and Moment method. The method of moments seems to be in class of its own, since it does not have to model the space between or around the antenna elements, and is therefore appropriate for Yagi type antennae. It does, however, involve exacting numeric integration, as well as matrix inversion.

The TLM method does not use either, while FDTD and finite elements or finite differences do not need heavy duty numeric integration. Finite elements or finite differences involve matrix inversion, though. More sections for computational electromagnetics are in preparation.

Reference: J. B. Davies and C. A. Muilwyk: 'Numerical solution of uniform hollow waveguides with boundaries of arbitrary shape', Proc. IEE, Vol. 113, No. 2, pp. 277-283

• Site map
• Finite elements for electrostatics
• Moment method
• Long file splitter

Present random quote:

Ignorance is not bliss

Up to this point: Valid XHTML 1.0!