Main Nomenclature

g(x) Unknown function over the interval .

f(x) Known function over the interval .

K(x,u) Kernel of the integral equation.

g_Ω (x) Truncated function of over the interval x∈R .

f^∞ (x) Extended function of over the interval x∈R.

g_Ω (x_n) Instantaneous value of at sample x_n∈R.

g ̂ _Ω (x_n) Instantaneous approximation of g_Ω (x_n ).

Ψ(x) Nonideal sampling function over the interval x∈R.

ϑ Sampling rate.

p_n (x;Δ) Generalized interpolation function over the interval x∈R .

h ̃(v) Fourier transform of h(x).

g ̃^Ω (v)Fourier transform of g_Ω (x).

BW₂ ^h Essential bandwidth of ∈L² .

L^1∧2 Means L¹∩L² .

Mathematical definitions and other notations can be found in appendicesA.1 and A.2.

Introduction

Integral equations (IEs) allow the modeling of a large number of complex problems in several areas of engineering and physics such as heat and mass transfer, oscillation theory and electromagnetic fields. An IE is an equation in which an unknown function related to an integral needs to be found under a known excitation . Some examples include one-dimensional IEs with forms ∃g(∙):∫_a ^b K(x,u)g(u)du=f(x) and ∃g(∙):∫_a ^xK(x,u)g(u)du=f(x), which are respectively known as the Fredholm and Volterra equations of the first kind, in which is usually known as the kernel.

Several techniques have been developed over the years to solve this kind of problem using exact, asymptotic and numerical approaches. A comprehensive summary of these techniques for different problems can be found in [1]–[3]. In particular, a well-known numerical technique to solve electromagnetic problems using IEs is called the method of moments (MoM) [4], [5]. This approach allows the development of a general framework including different types of modern techniques, known as the method of weighted residuals (MWR). For example, the finite element method (FEM) and the finite difference method (FDM) can be represented using the MWR [6], [7]. In recent times, the MoM has been used to solve multidimensional open boundary problems in electrodynamics, and it is included in many CAD platforms such as HFSS, IE3D, FEKO, and NEC.

In [8], we proposed an original approach to solve integral equations, based on the nonideal instantaneous sampling theory, which involved developing a novel framework for solving linear operators using Fourier analysis. The new approach was based on the fact that a Fredholm equation of the first kind with a difference kernel K(x,u)=K(x-u) can be expressed as an extended convolutional problem with form K(x)*g_Ω (x)=f^∞ (x) , where where is a truncated function of , is an extended function of , , and is the standard convolution operator. As a result, although the problem must be solved via matrix techniques using the original domain, the convolution properties and Fourier analysis can be used to design an efficient solution. This new property, known as the convolutional equivalence, allows the establishment of design rules that cannot be determined when using the MoM. These results are detailed in this paper through a new framework termed as the method of sampling for one dimension (MoS!-1D).

Most previous works aimed at solving Fredholm equations in the context of electrodynamics were focused mainly on the MoM. The different approaches established benchmarks for computational techniques (e.g., the comparison between MoM, FEM and FD-TD) and experimental techniques (e.g., the wire antenna) and identified particular computational issues (e.g., extremely ill-conditioned matrices), and they have been the cornerstone for generalizations in two and three dimensions. For example, [4] and [9] solved the problem of wire antennas with arbitrary shapes; [10] and [11] used telegraphist’s equations to evaluate transmission lines by using the Green function; [12] proposed the use of semiorthogonal compactly supported spline wavelets to evaluate an infinitely long metallic cylinder illuminated by a TM plane wave; [13] found the unknown current of a linear dipole with MFIE using high-order basis functions; and [14] and [15] solved the scattering for an infinite curved smooth strip considering the TM and TE cases by using the Nystrom method. However, few discussions exist regarding the relationship between the system theory based on the convolution operator and the Fredholm equation in the context of electrodynamics. Notably, [16] defined an effective bandwidth of a scattered field using the bandlimited behavior of the far field. Subsequently, [17] and [18] used this property to define basis functions with better convergence based on the entire domain and quasi-localized bandlimited interpolation functions. Furthermore, [19] used this behavior to simplify the Green function by removing the unexpected spectral components. This approach was supported by the similarity of the problem with a convolutional formulation, which was not very well defined.

The main contributions of this paper are new design paradigms and strategies for solving a Fredholm equation of the first kind, which covers the MoM methodology. For instance, we found that the functions used to determine a linear equation system can be designed as generalized interpolation functions without the assumption of continuity or smoothness, the weighting functions are sampled linear filters with particular specifications, and the new matrix coefficients Smn may contain information concerning the supradomain (x ∈ Ω^c). In the same manner, we proved that the generalized interpolation functions can be different from the functions used to approximate the unknown continuous function, and it is possible to reduce the computational cost based on the properties of linear and invariant systems.

The remaining sections of the article are organized as follows: in method of sampling section describes the justification and development of the novel approach using the nonideal instantaneous sampling theory. In analysis and evaluation section, several aspects of the matrix approach are discussed, and two examples are numerically evaluated. In conclusions, we present the conclusions of this research. Finally, appendices A and B detail several definitions, theorems and properties needed to make this paper self-contained.

Method of Sampling

Convolutional Equivalence

Assume a convolutional integral equation (CIE) with the form ∃g(∙):∫_ab K(x,u)g(u)du=f(x) with a<b ∀x∈Ω≜[a,b]⊂R, and K(x,u)=K(x-u). . Let K∈L² and g_Ω ∈ L^1∧2 with the essential bandwidths being and (see definition 8).

Consider the extended equivalent problem ∃g_Ω (∙):∫_-∞ ^∞ K(x,u)g(u)du=f^∞ (x)=f_Ω ^∞ (x)+f_Ωc ^∞ (x) , with f_Ω ^∞ (x)=f_Ω (x) and Ω^c=R-Ω [8]. From this formulation, the extended function f^∞ (x) has two components: The first (f_Ω ^∞) is the known function because f_Ω (x)={f(x),∀x∈Ω and zero elsewhere ; the second (f_Ωc ^∞) is the extended function f^∞ (x) evaluated in the supradomain (x∈Ω^c),which is a degree of freedom with mathematical sense allowing the hypothesis g_Ω K→ f^∞.For a major discussion concerning this equivalence,see [8].

The analysis of this new problem can be performed in the original domain (x-domain) or the Fourier domain (v-domain) by using

where f^∞∈L² with BW₂ ^f∞ . The critical bandwidth of this new problem is denoted by BW=min⁡{BW₂ ^f∞,BW₂ ^fΩ∞} .

The solution of this problem must be in the domain because the degree of freedom f(_Ω ^c)^∞ (x) is unknown, and even though this function makes mathematical sense, it is not related to the physics of the original problem.

Asymptotic Reduction

Let xn ∈ {...,x_-1,x₀,x₁,…} be sampling points such that S≜{x₁,…,x_N}⊆Ω and S^c≜{…,x_-1,x₀,x_N+1,x_N+2,…}⊈Ω . Let Ψ(x) be a nonideal sampling function defined by

where g ̂_Ω (x_n ) is an instantaneous approximation of g_Ω (x_n ), and p_n (x;Δ) is a generalized interpolation function [8].

In contrast with other computational methods developed so far, the nonideal sampling function Ψ(x) is not necessarily a satisfactory representation of g_Ω(x) in the conventional sense (e.g., variational sense). This fact means that Ψ(x):∆ >0 does not necessarily minimize a continuous error, such as min_∀x∈Ωmax|Ψ(x)-g_Ω (x)| , min⁡∫_Ω|Ψ(x)-g_Ω (x)|dx or min⁡∫_Ω(Ψ(x)-g_Ω (x))² dx , as can be inferred from figure 1. We use functions such that g ̂_Ω(x) and ϵ_Δ (x) are asymptotically (Δ→0⁺) separable through the v-domain. As a result, Ψ(x) is not necessarily conformed by a combination of orthogonal or smooth functions with domain Ω (see figure 1).

If we use uniform generalized interpolation functions p_n (x;Δ)=p(x-x_n;Δ)=p(x;Δ)*δ(x-x_n) with domain Ω_pn ⊆R and uniform samples x_n=x_n-1+Δ=x₀+nΔ, the spectral error is (see theorem 11)

Therefore, the spectral error function ϵ ̃_Δ (v) has three spectral components: The first component is the distortion of {∆^-1∙p ̃(v;Δ)-1} in |v|≤BW. The second component is the interference of in ρ ̃_∆ (v) in |v|>BW. From (7), this interference may be centered at the harmonics ±kϑ , and it may be located in the intervals |v±kϑ|≤BW . The third component is the aliasing of ρ ̃_∆ (v) in |v|≤BW . Consequently, if

is a nascent Delta-Dirac function, then the distortion and aliasing are asymptotically removed via ϵ ̃_Δ (v) for g_Ω ∈ L^1∧2 using the Nyquist criterion ϑ≜Δ^-1>2BW . If the Nyquist criterion is not satisfied, the aliasing cannot be controlled, and therefore, it is not possible for convergence to occur in the problem.

v-domain to be defined as follows:

The novel approach is based on the solution of this new problem assuming the asymptotic reduction of the spectral error, summarized by

As a result, g ̂_Ω (x)→g_Ω (x) and g ̂_Ω (x_n )→g_Ω (x_n ) almost everywhere. Therefore, the relations among the interpolation function p(x;Δ) , kernel K(x) and sampling rate (ϑ≜Δ^-1) to achieve this particular behavior are not arbitrary.

Remark 1.Because we do not know g ̂_Ω (x) a priori, our methodology only obtains the instantaneous values g ̂_Ω (x_n ) [8]. Consequently, a smooth interpolation algorithm, denoted by J{∙} , must be used to find the unknown continuous approximation given by g ̂_Ω (x)=I{g ̂_Ω (x_n )}, under the hypothesis of compliance with the Nyquist criterion.

In contrast with other methodologies, the proposed approach allows the determination of several design strategies and interpretations to solve linear operators using Fourier analysis. For clarity and space reasons, we discuss only two of these strategies below.

Bandwidth Matching Procedure (BM-P)

We propose a general methodology to solve a CIE assuming that the combination of the interpolation functions and the kernel allows the undervaluing of the influence of the spectral error in the problem, for example, by assuming BW₂ ^gΩ ≫BW₂ ^K or ‖(K ̃∙g ̃^Ω )(v)‖₂≫‖(K ̃∙ϵ ̃_Δ )(v)‖₂ .We name this technique the bandwidth matching procedure (BM-P).

x=x_m∈S i Without loss of generality, if we take samples at i in (9),

Rewriting (12) in integral form, we obtain

The linear equation system using BM-P is SΨ=f with Ψ=[g ̂_Ω (x_n )], f=[f(x_m)] and S=[S_mn ], where g ̂_Ω (x_n ) is an instantaneous approximation of g_Ω (x_n ), f(x_m) is the instantaneous value printed by the output, and S_mn is the instantaneous equivalent system between the samples, defined as

Corollary 2. If we design p_n (∙;Δ) as the basis functions and Ω_pn ⊆Ω , the BM-P obtains the same equations as those used by the MoM when using the point matching procedure (PM-P).

If the Fourier transform of the functions are well-defined and holds F{K(x)*p(x;Δ)}=K ̃(v)p ̃(v;Δ) with K ̃∙p ̃∈L¹ , another way to evaluate the system coefficients [8] is by using

Regularization Technique by Convolution (RT)

We propose another general methodology to solve a CIE assuming some spectral anomaly, for instance BW₂ ^gΩ≪B_₂ ^K or ‖(K ̃∙g ̃^Ω )(v)‖₂≪‖(K ̃∙ϵ ̃_Δ )(v)‖₂ . In this case, we propose the regulation of (1) using R(x;Δ) with the convolution operator (literally, filtering the problem) as

such that the new system can be modeled as a well-defined low-pass problem with a new critical bandwidth (BW) ̅<BW. We name this technique as the regularization technique by convolution (RT).

Therefore, to have asymptotic convergence toward the original problem, we design the filter R(x;Δ) as a nascent Delta-Dirac function with nulls specifically located to reduce the spectral anomalies in the problem. In particular, we design these nulls to remove asymptotically the last component of the spectral error, that is, the harmonic interference.

Without loss of generality, if we take samples x = x_m∈ S at in the new equivalent instantaneous problem,

Rewriting (18) in integral form, we obtain

where the function R_m (ξ;Δ) = R(x_m-ξ;Δ) is literally the sampled impulse response of a standard linear filter.

The linear equation system using RT is S ̅Ψ = f ̅ with Ψ = [g ̂Ω (xn )], f ̅=[f ̅∞ (xm)], and S ̅=[S ̅mn ], where g ̂Ω(xn ), is an instantaneous approximation of gΩ (xn) considering the filter R(x;Δ), f ̅∞ (xm) , is the instantaneous value printed by the regularized output, and S ̅_mn is the new instantaneous equivalent system between the samples given by

If we design R_m (ξ;Δ) with subdomains (Ω_Rm⊂Ω) or entire domains (Ω_Rm =Ω) to avoid unknown information of f^∞ (x), the resulting respective techniques are called the subdomain regularization technique (S-RT) and the entire domain regularization technique (E-RT).

Corollary 3. If we design p_n (∙;Δ) as the basis functions, R_m (∙;Δ) as the weighting functions related to an interior product, and Ω_pn,Ω_Rm ⊆Ω, the RT obtains the same equations as those used when applying the general approach of the MoM (and in general, those employed when using the MWR).

If it is possible to use the associativity property R*{K*Ψ}={R*K}*Ψ=K*{R*Ψ} (e.g., see Fubini’s theorem and Hölder’s inequality), the formulation can be simplified and computational cost can be reduced by using ∑_n ^N=₁g_Ω (x _n) K(x)*P_n (x;∆)= f ̅^∞ (x) because P_n (x;∆)≜(R*_pn )(x;∆) and Ω_Pn ≜Ω_R*Ω_pn may have closed forms. The new function P_n (x;∆) is known as an ultrafunction, and the simplified formulation as termed as the ultra-RT. In this case, the system coefficients [8] can be obtained using

Finally, if the Fourier transforms of the functions are well defined and F{R(x;Δ)*K(x)*p(x;Δ)}=R ̃(v;Δ) K ̃(v)p ̃(v;Δ) holds with R ̃∙K ̃∙p ̃∈L¹,

From (10) and (17), if we choose R(x;Δ)=∆^-1∙p(x;Δ), the harmonic interference can be reduced by using a quadratic factor. In the same way, if the related functions are designed as in the MoM, this technique can be considered an alternative explanation of the Galerkin method, and it can easily explain under what conditions the solution can be successfully applied [8].

Design Considerations

The asymptotic reduction is obtained and justified using the sample theory in the following manner. The distortion is asymptotically removed by choosing the interpolation function as a nascent Delta-Dirac function because lim┬(Δ→0⁺)⁡δ ̃_Δ (v)→1:∀|v|≤BW. The interference is reduced by designing the interpolation functions (and the filter) with spectral nulls at the harmonic interference, given by |δ ̃_Δ (v)|≤δ_I≪1,∀|v±kϑ|≤BW, k ∈ N. This interference is asymptotically removed by increasing the sampling rate because the kernel has a finite bandwidth. The spectral nulls should be designed to exactly match with this interference to reduce the possibility of occurrence of ambiguities caused by the nonexistence of f(_Ωc ^∞ (x). Finally, the aliasing is asymptotically removed by increasing the sampling rate because ∀g ̂_Ω ∈ L^1∧2 has the property lim┬(_v→∞)g ̂ ̃^Ω (v)→0. As a result, lim┬(Δ→0)g ̂ ̃^Ω (v±kϑ)〗→0:∀|v|≤BW,k∈N and lim┬(Δ→0⁺)⁡ρ ̃_∆ (v)〗→0:∀|v|≤BW. For details and examples of these considerations, see [8].

Analysis and Evaluation

Novel Approach to Solve Convolutional Integral Equations

Because the approach of the MoS! has a different way to obtain convergence, to the best of the author’s knowledge, our proposal is an innovative framework for solving Fredholm equations of the first kind. The MWR solves ∃g(∙):L{g}=f by using g ̂=∑_n ^N=₁α_n b_n, such that the residual error r=L{ĝ}-f is reduced by means of the interior product ⟨r,R_m ⟩=0 [4]. The novel approach includes the same interpretation because the same equations are determined when Ω_pn,Ω_Rm ⊆Ω, albeit with new paradigms and degrees of freedom. We highlight the following aspects of our methodology.

Before taking the samples x_m in (12) and (18), the functions p(x;Δ) and R(x;Δ) can be designed to reduce the error in both domains because ‖ϵ ̃_Δ (v)‖₂=‖ϵ_Δ (x)‖₂. A more comprehensive approach to design R(x;Δ) (and the problem) using the communication theory may be considered because that function is not necessarily related to a reduction in the residual error in the sense of a geometric projection. For instance, the filter R(∙;Δ) may be designed as a standard receptor system with several roles related to reducing the noise caused by the rounding errors (e. g., pre-emphasis and de-emphasis techniques and matched filter) and transforming a dense matrix in a sparse matrix (e. g., intersymbol interference techniques).

After taking the samples x_m in (12) and (18), the obtained matrices SΨ=f and S ̅Ψ=f ̅ are the formulations of a standard discrete-time system (formally a discrete x-domain system), in which it is now evident that all the information between the samples x_m and x_m+1 is lost because, in general, f^∞ (x) and f ̅^∞ (x) are not bandlimited functions. As a result, Ψ(x) is not necessarily a good continuous approximation of the unknown function g_Ω (x).

Remark 4. The only methodological difference between BM-P and RT is that the original problem is regularized by a filter to obtain a slower problem. As a consequence, the MoM (and, in general, the MWR) is an instantaneous technique, or “point matching procedure”, independent of the weighting functions.

Discussion Concerning Interpolation Functions

A standard interpolation function is defined as a basis function in the context of being used by standard numerical methods such as the MoM. For instance, the well-known rectangular and triangular basis functions are p_n ^Π (x;Δ)=Π_Δ (x-x_n) and p_n ^Λ (x;Δ)=Λ_Δ (x-x_n). Our approach allows the analysis of these functions because their Fourier transforms are well defined by p ̃^Π (v;Δ)=ϑ^-1 sinc(πv/ϑ) and p ̃^Λ (v;Δ)=ϑ^-1sinc² (πv/ϑ).

Consequently, to maintain the distortion in 0.95≤|Δ^-1∙p ̃(v;Δ)|≤1,∀|v|≤BW, it is required that BW≤0.175ϑ and BW≤0.124ϑ. As a result, the minimum sampling rate with that distortion for p^Π (∙;Δ) is ϑ= 0.175^-1 BW≥6∙BW and that for p^Λ(∙;Δ) is ϑ= 0.124^-1 BW≥9∙BW.

These results prove that if the kernel has a narrow bandwidth such that it reduces the interference of the sampling process, the rectangular function should unexpectedly obtain a better accuracy than the triangular function does because the distortion is lower at the same sample rate.

In contrast, because the reduction in the first harmonic interference is given by |Δ^-1∙p ̃(v;Δ)|<δ_I ^±1,∀|v±ϑ|<BW with δ_I ^±1.Π=0.2016 and δ_I ^±1.Λ=0.0192, the triangular function should demonstrate superior performance in the kernels with wide bandwidths (e. g., kernels with some singularity or quasi-singularity) because the main harmonic interference of the sampling process is one decade lower at the same sample rate. From our approach, it is now evident that the rectangular function does not necessarily minimize ‖Ψ(x)-g_Ω (x)‖, although the instantaneous value g ̂_Ω (x_n) may be more accurate.

A generalized interpolation function is defined in a different direction using the asymptotic separability in the domain between the unknown approximation and the error. As a result, Ψ(x) may not be a standard representation of g_Ω (x), (see figure 1). Nevertheless, these functions require careful design to avoid the numerical instability caused by the spectral mismatch (e.g., new spectral nulls in the problem). An uncomplicated method to formulate a generalized interpolation function is

where p̃(v;Δ)=∏^M _i=1p ̃ⁱ (v;Δ) and Δ^-1∙p(x;Δ) represent a nascent Delta-Dirac function with nulls located strictly in the harmonic interference of the sampling process. To demonstrate this methodology, the authors propose p^e (x;Δ)=(p¹⁾*p²⁾)(x;Δ) as an example, such that p¹⁾ (x;Δ)=Π_Δ (x)-c/2 Π_Δ (x+Δ/2)-c/2 Π_Δ (x-Δ/2) and p²⁾ (x;Δ)= t₁∙Π_Δ (x), as shown in figure 2a.

The formulation of this unusual combination in the v-domain is p ̃^e (v;Δ)=(p ̃¹)∙p ̃²) )(v;Δ)=t₁ ϑ^-2 sinc² (πv/ϑ){1-c∙cos⁡〖(πv/ϑ)〗}; therefore, to achieve the nascent Delta-Dirac behavior, the constant t₁=ϑ/(1-c). Figure 2b shows that if c=1/3, p ̃^e (v;Δ) increases the linearity regarding p ̃^Π (v;Δ), and the spectral nulls have a similar attenuation to p ̃^Λ (v;Δ). As a result, using low-pass kernels, this interpolation function ought to exhibit better performance and low computational cost. The function defined in this example is similar to those defined by the wavelet theory, although it does not necessarily imply a good continuous representation of the original function or a multiresolution signal decomposition.

Remark 5. Because the rectangular and triangular functions are nascent Delta-Dirac functions with nulls at ±kϑ, that basis functions are generalized interpolation functions, as well. The main difference between MoM and MoS! in these cases it is that the instantaneous values g ̂_Ω (x_n ) may allow a better continuous approximation than Ψ(x) does.

Finally, if |K ̃(v)| is strongly decreasing toward zero at high frequencies, it is possible to define

because (46) is proved from (45). As a result, some low-pass kernels may have a convergent closed-form matrix using high sampling rates , which allows the study of an analytic or asymptotic inverse for the problem, using the digital signal processing framework.

Example 1: ANarrow Band Kernel

The problem with is the CIE of the differential equation , where is not an arbitrary function. This function must have some properties [1] such as appropriate boundaries , and it must be a bounded twice-differentiable function with a bounded first derivative. As a theoretical and numerical benchmark, we propose and

Using our approach, the extended convolutional model is

where Ω=[0,1], K(x)=|x|, f_Ω ^∞ (x)=f_Ω (x), and f_(Ωc) ^∞ (x) is a degree of freedom with mathematical sense. Because d/dx |x|=sgn(x),∀x≠0, K ̃(v)=-1/(2π² v² ),∀v≠0.

In this problem, it not possible to use the v-domain to calculate the matrix coefficients because K ̃(v) has a distribution not calculated at v=0. However, the Fourier domain interpretation is reasonably valid because the Fourier transform is well defined using the distribution theory, and the convolution operator is well established in this problem for a finite interval Ω. As a result, by using the spectral interpretation for ∀v≠0, it is valid to assume that this problem is a narrowband kernel for BW≥π^-1 2^-1/2.

The matrix formulation simplifying Ω=(0,1) for p^Π (x;Δ) is (I+4B)Ψ=4ϑ² f, where B=[B_mn ]=[|m-n|] and I is the identity matrix; for p^Λ (x;Δ) is (I+3B)Ψ=3ϑ² f, and for p^Π (x;Δ) and R^Π (x;Δ)=Δ^-1 Π_Δ (x) is (I+3B)Ψ=3ϑ² f ̅, where f ̅_m=π-¹-2π^-2 sin⁡(πx_m)sinc(1/2 π∆). From the MoM approach, the first two formulations pertain to the PM-P, and the last one pertains to the Galerkin method.

The matrix formulation using the MoS!-1D with p^e (x;Δ) is (T+B)Ψ=ϑ² f, where T is an uniform tridiagonal matrix with coefficients T_1,1=12^-1 (4-7c)(1-c)^-1 and T_1,2=24^-1 (24-25c) (1-c)-¹-1. We emphasize that this solution is different from the MoM approach because the integrals have different domains.

Figure 3 compares the different estimations through the relative residual error

for Θ≜[10,1-10Δ]; in this case, several new paradigms are evident. Because a low-pass kernel is more relevant for the distortion rather than the interference, the function p^Π (∙) finds a better outcome than p^Λ (∙) does, and the new proposal p^e (∙) has a significantly better performance for Θ, with a different slope (logarithm) of convergence. For instance, the number of unknown variables for an instantaneous relative residual error of 10^-7 is N≈50 for p^e (∙) using the MoS! and N≈2000 for the classical low order functions using the MoM.

Nevertheless, we note that p^Λ (∙) still exhibits better elimination of the first harmonic interference in figure 2b. This little energy is accumulated at the boundaries at x_n∈(0,9Δ]∪[1-9Δ,1). As a result, this location is the only one at which p^Λ (∙) exhibits better performance in this problem.

Considering special interpolation functions, we found an inverse with linear complexity using p^δ (x;Δ). The solution in this case is Ψ=∆^-2 B^-1 f, where B^-1 is detailed in appendix B. This matrix solution is similar to the matrix used by FDM for g(x)=1/2 d²/dx² f(x). This solution exhibits the same performance as p^Λ (∙), thereby proving conclusively that the interference is not dominant in this problem. The theoretical matrix B^-1 allows the solution of the general case (I+βB)Ψ=βϑ² f by applying the Neuman series by means of Ψ≈{β^-1 B^-1+⋯+(-1)^k+1 β^-k B^-k}βϑ² f for N≥3 and ρ^{+ (β^(-1) B^(-1) )}<2β^-1<1. As shown in figure 3, this algorithm converges with k=3 in the same way as Gaussian elimination for the rectangular (β=4) and triangular (β=3) cases.

This problem exhibits ill conditioning at high sampling rates (or breakdown segmentation) because this kernel has extremely low spectral information at extremely high frequencies. From figure 3, the breakdown for the classical basis functions can be noted at and that for p^e (∙) can be noted at N≳400

In summary, several of the new paradigms were distinguished in this example. First, we demonstrated that Ψ(x) is not necessarily the best way to represent the unknown function because the rectangular function exhibits a better instantaneous performance than the triangular function does. Second, we demonstrated that the Galerkin method does not increase the accuracy because this problem is not anomalous at high frequencies. Consistent with the ultraformulation, the same system matrix is found using the triangular function with PM-P and the rectangular function with the Galerkin method. As a consequence, any difference between both results can only be interpreted by filtering of the source. Third, considering the example of p^e (x;Δ), the use of generalized interpolation functions with supradomain information can significantly increase the accuracy. Finally, we demonstrated that the formulation based on p^δ (x;Δ) can be used to design asymptotic techniques to solve this narrowband problem.

Example 2: AWideband Kernel

To discuss the novel approach using a standard electromagnetic benchmark, the authors solved the unknown current of a linear dipole with radio 0<a≤10^-2 λ and length L≫a using the reduced or approximate Pocklington equation under a particular printed source E(x). Without loss of generality, we use the Richmond simplification [20] to use a standard CIE. The extended convolutional problem is

where x∈R, Ω=(-L/2,L/2), K₀=j∙(2πfϵ)^-1, j=√(-1), c=λf, k=2π/λ, f [Hz] is the frequency, c [m/s] is the speed and λ [m] is the wavelength in a medium with permittivity ϵ [F/m], E_Ω ^∞ (x)=E_Ω (x) is the printed source in Ω, and E_Ωc ^∞ (x) is a degree of freedom with mathematical sense.

This quasi-singular kernel can be understood from the Fourier transform K ̃_ap (v)=F{K_ap (x)} approximated through the DFT with samples x_n=a/150 n. For instance, figure 4c details log₁₀ |K ̃_ap (v_n )| using a=0.05λ and λ=1 m, in which several behaviors are evident.

As expected, this kernel has a large bandwidth caused by its quasi-singularity. We find a local minimum at v=±1 (and, in general, at v=±λ^-1) related to the main spatial harmonic of the unknown current, and a local maximum at v≈±50 (and, in general, at v=±1/4 a^-1). Moreover, the transfer function has an exponential decay for |v|≳200, whose value tends rapidly toward zero. Although the bounded oscillations in |v|≳600 have a pseudorandom appearance, they are caused by an infinite x-domain function with a fast decay toward zero at high frequencies, which is approximated by a truncated function.

To compare the MoM and MoS!, we find the input impedance (Z_in) of a center-fed dipole using the magnetic frill source simplified by

for all x∈Ω, b=2.3a and V_S=1. Because the complex power provided by the source is S=1/2 ∫_-L/2 ^L/2 E(x)I* (x)dx, and it can also be defined using S=1/2 |I(0)|² Z_in, we estimate the input impedance with the MoM using I_Ω (x)~Ψ(x) as

In the MoS!, any continuous operation with the unknown current requires the definition of an interpolator because Ψ(x) loses information between discrete samples. Without loss of generality, we find I ̂_Ω (x) through the piecewise cubic Hermite interpolating polynomials [21] using I ̂_Ω (x_n) and the boundaries I ̂(±L/2)=0. Consequently, we estimate the input impedance with the MoS! by using

Furthermore, the difference between the MoS! and MoM is particularly accentuated when we use bandlimited interpolation functions. For instance, and only to differentiate between the methods, if p(x;Δ)=sinc(πx∆^-1) with x_n=-L/2+∆∙n and ∆=L/(N+1), PM-P obtains an impedance matrix (Z) with coefficients

Therefore, O(N²) is required to ensure computational storage and O(N³) is required to solve the linear equation system using Gaussian elimination because the impedance matrix does not have useful properties. In contrast, the system matrix (S) obtained by (15) is a symmetric Toeplitz matrix, which only requires O(N) for computational storage, and several efficient algorithms are available to solve the linear equation system [22]. We compare Z_in using PM-P and BM-P with L=0.47λ in table 1a.

Another application of the MoS! is to use the system matrix as a preconditioner for the MoM because S^-1 E may lead to a low computational cost and both formulations are closely related at x_m,x_n∈S. As an example, figure 5 shows the relative residual error in solving ZΨ=E using sinc(∙) and the conjugate gradient squared method (CGS) with and without the preconditioner S. As shown, the solution with the preconditioner S has fast convergence, even for a small number of unknown variables (N).

Although the unknown current has a well-identified fundamental harmonic at v=±1, a considerable bandwidth (BW> 1/2Δ|_N=101= N+1/2L |_N=101=108.5) is necessary to obtain an absolute relative error less than 3.5 % in the imaginary part. The analysis of the magnetic frill source shows that its bandwidth is the cause of this situation. If we assume that Z_in is mainly obtained from band base information, filtering the source with the normalized filter R^Λ (x;Δ)=Δ^-1 Λ_Δ (x) should increase the convergence rate, as shown in table 1b.

We compare the results with those obtained using the Galerkin method in table 1c. According to our approach, the Galerkin method obtains an oscillatory convergence rate because filters conformed by truncated sinc functions do not have well located spectral nulls to ensure the efficient elimination of the interference.

The Z_in evaluated using the ultraformulation (23)-(25) and using the v-domain (16), (26) exhibit relative residual errors that are less than 1×10^-5. In particular, using the symmetric of (32), we simplify (16) and (26) using RΛ (∙) by using

Although (40)-(41) have some disadvantages caused by the Fourier transform estimation and the evaluation of an integrand with an oscillatory behavior, they have other advantages owing to the prior theoretical knowledge of R ̃(∙) and p ̃(∙), the finite interval of integration, the use of the FFT for a fast calculation and windowing to reduce the truncation effect over K ̃_ap (v). The system coefficients calculated by v-domain are a relevant option for other more convenient bandlimited interpolation functions, such as the raised cosine functions.

In conclusion, using potentially the most critical bandlimited interpolation function, we demonstrate that the MoS! is a viable alternative with low computational cost because the new matrix has useful properties based on symmetry, and it can be calculated using low cost alternative expressions. Although the MoM has a better convergence rate in this example, the solution of S ̅_N+18 is mostly equivalent to Z ̅_N. Moreover, we show that the system matrix can be used as a preconditioner for the MoM using a well-known Krylov subspace iterative method with very fast convergence.

Furthermore, the MoS! can simplify the solution using subdomain functions because the ultraformulation is well defined for K∈L². For such a case, we emphasize that if Ω_pn,Ω_Rm)⊆Ω, both matrices are equal (S=Z). However, the input impedance is different because both methods employ different approaches for the unknown continuous approximation. In this example, this difference is smaller because the input impedance is calculated using an integral operator (a low-pass operator). To discuss the ultraformulation, we first evaluated Z_in in tables 2a and b using the rectangular and triangular functions with PM-P (and BM-P in the parenthesis) using (15), and MoM (and S-RT in the parenthesis) with the filter R^Λ (∙) using (21). As expected, the rectangular function with BM-P exhibited an oscillatory behavior without convergence because the reduced kernel has low attenuation for its interference (see figure 4c in |v|≤50).

Because (R^Λ*p^Π)(x;∆) and (R^Λ*p^Λ)(x;∆) have closed forms, the ultraformulation has two advantages. The first advantage is the reduction of the computational cost owing to the transformation of the 2D integral to 1D. The second advantage is the reduction of the nontrivial effects of rounding errors. In particular, (24) reduces the jitter produced between the interpolation functions and filters, and (25) reduces the jitter produced by the kernel. For instance, we note that the input impedance in table 2b using p^Π (∙) and R^Λ (∙) must lead to better performance because the kernel and the source have been filtered. Table 2c summarizes the new estimation using (24), in which we remove the lack of convergence without changing the numerical methods. We highlight this nontrivial issue because κ(S_N)≤600 for N≤113.

Continuing the discussion on regularization, we note that filters with subdomains are inefficient in reducing the source’s bandwidth due to the time-frequency relationship (formally, the relationship between the x and v domains). To enhance the convergence, we propose the sampled filter R_Ω ^g (x_m-x) inspired by the standard Gaussian filter

where σ=c∙∆ to obtain a nascent Delta-Dirac function (c≥1).

Table 2d shows the recalculation of Z_in using the entire domain filter with c=1.5. As shown, it is possible to improve the convergence significantly, even for low order interpolation functions (e. g., tables 1 and 2 show that the triangular interpolation function with the quasi-Gaussian filter using N=17 has the same approximate performance as that of the triangular filter using N=41, without filter using N=77, and with the entire domain formulations using the sinc function and N>53). Nevertheless, some disadvantages were observed when filters were used in this problem. First, the conditioning number and the fake oscillations [23] and [24] may be increased because the filters reduce the bandwidth artificially. Second, the use of any filter includes the establishment of a particular model for the discontinuities at low sampling rates, which implies that a specific model may be used for the current at the feed point and the edges of the wire. Lastly, monotone convergence using the quasi-Gaussian filter when c≳1.5 does not occur because spectral mismatches exist in this particular problem (e.g., a relevant transience exists for N<17 and c=1.5).

In summary, we solved a classical electromagnetic problem with several new paradigms. It was noted that the limiting factor for fast convergence is the bandwidth of the printed source. Therefore, the regularization filter, and not the interpolation function, is the design challenge for a fast and numerically stable solution. Moreover, it was demonstrated that the instantaneous values obtained from the system matrix allow the solution of the problem, and less notably, Ψ(x) was considered to be an interpolator with an unnecessary high computational cost. According to our approach, the use of Ψ(x) to calculate other performance parameters based on wideband operators (e.g., vector operators) is not desirable because Ψ(x) may possess information that is not valid at high frequencies. In particular, the information obtained using the reduced kernel at high frequencies by employing any method is strictly incorrect because the real part of this kernel is an almost ideal low-pass filter with essential bandwidth BW₂ ^Re{K_ap}≈λ^-1. For these reasons, we changed the terminology of “basis functions” to “interpolation functions” because these functions do not necessarily provide basis for the unknown function following the standard algebraic approach.

Conclusions

The authors presented a novel approach to solve integral equations with difference kernels, based on the generalization of a Fredholm equation of the first kind, as an extended convolutional problem. The results demonstrated that the method of moments (and in general, the method of weighted residuals) is a particular case, in which new interpretations, degrees of freedoms and design rules can be found using standard matrix techniques.

The new approach can be used to address major methodological questions for the method of moments; these questions include those pertaining to the design of all the functions related with the computational method (e.g., to the best of the author’s knowledge, a general approach for the design of weighting functions has not been reported so far), the establishment of preconditioners in a systematic manner (e.g., in problems using bandlimited basis functions), the reduction of computational cost based on several new properties (e.g., using ultraformulation and functions with an approach that does not necessarily imply orthogonality, differentiability or smooth continuity), the use of the Fourier transform as a fundamental design tool although the solution must be performed in the original domain (e.g., via matrix coefficients calculated by the Fourier domain), and better understanding the formulation based on the interior product as a particular case of regularization using convolution.

References

Notes

* Research article This article is derived from the doctoral research of the first coauthor at the Universidad de los Andes, Colombia. That study was partially financed by the Plan de Formación Permanente del Profesor Javeriano (PFPP) from Pontificia Universidad Javeriana, Colombia.

How to cite this article: C.I.Páez Rued a, andR. Bustamante Miller,“Novel computational approach to solve convolutional integral equations: Method of sampling for one dimension,” Ing. Univ.,vol. 23, no. 2,2019. https://doi.org/10.11144/Javeriana.iyu23-2.ncas

Author notes

^a Autor de correspondencia. Correo electrónico paez.carlos@javeriana.edu.co

Descargar