The Nyström method
In this post, we will shortly describe the Nyström method [1], a valid alternative to the Method of Moments (MoM) [2] for numerically solving integral equations and we will use it to solve an electromagnetic scattering problem from an indefinitely long, perfectly conducting cylinder.
The Nyström method evaluates integrals using quadrature rules so that the unknown is not represented by the expansion coefficients of the induced current density over a certain basis, but rather by the current density at the quadrature nodes.
The Nyström Method should be applied with care to cases when the integral equation exhibits a singular kernel as in the case of the Electric Field Integral Equation [1,2] dealt with in this post.
This post is divided into three parts: theory, validation formulas and numerical implementation with Python. The validation formulas refer to the scattering by a circular cylinder. The Python implementation is available at this link and comprises both the Nyström Method code and the validation code.
Theory
The Electric Field Integral Equation (EFIE) for an indefinitely long, perfectly conducting cylinder embedded in free space under Transverse Magnetic (TM) incidence is the following:
In such an equation, C represents the contour of the cylinder’s cross section, ω is the angular frequency, μ0 is the free-space permeability, Jz is the only (z) component of the current density induced over the cylinder’s surface, x and y represent generic observation and source points over C, respectively, R=|x-y|, G is the Green’s function
provided by an Hankel function of zero-th order and second kind, k0 is the free space wavenumber and the left hand side of the EFIE is the only (z) component of the impinging electric field.
After having introduced the following parameterization of the contour C,
the EFIE can be rewritten as
where
An accurate application of the Nyström method requires to extract the logarithmic singularity of M(t,t’). For this reason, such a function is rewritten as
where
is the coefficient to give to the logarithmic singularity and
is the regular part of the kernel M(t,t’) which can be expressed as
and C=0.5772156649.
The EFIE can be then rewritten as:
The first integral refers to the irregular part of the kernel, while the second integral refers to the regular part.
As far as the second integral is concerned, it can be evaluated by a simple quadrature rule involving 2N points. In other words:
In order to properly deal with the logarithmic singularity, the first integral is evaluated by approximating Jz using a global trigonometric polynomial which is integrated exactly, see [1]:
where
By a point matching approach in t, the following system of linear equations is finally defined:
The above system can be recast in matrix form as:
In such a system of linear equations, the data are represented by the impinging field while the unknowns are the samples of the current density distribution.
Validation
A validation test of the method can be set up in the case of a circular cylinder with radius a. In this case, indeed, both the scattered field and the induced current density can be represented by cylindrical harmonics for which the expansion coefficients can be given in closed form [2].
Under the incidence of a plane wave propagating along the x1 axis
the z-component of the scattered field can be expanded as:
where Jn is the Bessel function of n-th order and
Finally, the z-component of the induced current density can be represented as
Python code
The first step of the Python code is to define the electromagnetic constants:
Later on, the parameters defining the object’s contour are defined. In the available code, the contour is assumed to be elliptical so that the lengths of the two half-axes need to be defined:
Once defined the contour, we need to select the number of quadrature nodes 2N. The number N is decided thanks to an iterative approach. The surface discretization level is fixed and the number of quadrature points is increased until their distance is larger than the fixed discretization level:
After having decided the number of quadrature nodes 2N, their positions on the object’s contour as well as the derivatives of the contour parameterization at the quadrature nodes are computed:
In the above snippet, the quadrature nodes are uniformly spaced in the independent variable of the contour representation parameter.
It is now time to fill in the matrix of the linear system representing the Nyström method. This is done by the following snippet implementing the equations detailed above:
In the above code snippet, the trigInterp function represents the trigonometric interpolation function used by the Nyström method and is below reported:
The last steps to achieve the unknown current distribution are represented by defining the impinging field and by the solution of the linear system:
Following the unknown determination, the series for the calculation of the reference current is summed up in order to check for the reconstruction accuracy:
Finally, starting from the calculation of the unknown current density, the radiated field is computed by the relation:
Such a computation is worked out in the far-field, namely when
where D is the diameter of the smaller circle enclosing the scatterer and λ is the wavelength:
The reference scattered field is calculated by summing up the above indicated cylindrical waves series
References
[1] D. Colton, R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Springer, Fourth Edition.
[2] C.A. Balanis, Advanced Engineering Electromagnetics, J. Wiley & Sons.
