Open access

Simulation and Analysis of Transient Processes in Open Axially-symmetrical Structures: Method of Exact Absorbing Boundary Conditions

Written By

Olena Shafalyuk, Yuriy Sirenko and Paul Smith

Submitted: 27 October 2010 Published: 21 June 2011

DOI: 10.5772/18329

From the Edited Volume

Electromagnetic Waves

Edited by Vitaliy Zhurbenko

Chapter metrics overview

2,527 Chapter Downloads

View Full Metrics

1. Introduction

Present day methodologies for mathematical simulation and computational experiment are generally implemented in electromagnetics through the solution of boundary-value (frequency domain) problems and initial boundary-value (time domain) problems for Maxwell’s equations. Most of the results of this theory concerning open resonators have been obtained by the frequency-domain methods. At the same time, a rich variety of applied problems (analysis of complex electrodynamic structures for the devices of vacuum and solid-state electronics, model synthesis of open dispersive structures for resonant quasi-optics, antenna engineering, and high-power electronics, etc.) can be efficiently solved with the help of more universal time-domain algorithms.

The fact that frequency domain approaches are somewhat limited in such problems is the motivation for this study. Moreover, presently known remedies to the various theoretical difficulties in the theory of non-stationary electromagnetic fields are not always satisfactory for practitioners. Such remedies affect the quality of some model problems and limit the capability of time-domain methods for studying transient and stationary processes. One such difficulty is the appropriate and efficient truncation of the computational domain in so-called open problems, i.e. problems where the computational domain is infinite along one or more spatial coordinates. Also, a number of questions occur when solving far-field problems, and problems involving extended sources or sources located in the far-zone.

In the present work, we address these difficulties for the case of TE0n- and TM0n-waves in axially-symmetrical open compact resonators with waveguide feed lines. Sections 2 and 3 are devoted to problem definition. In Sections 4 and 5, we derive exact absorbing conditions for outgoing pulsed waves that enable the replacement of an open problem with an equivalent closed one. In Section 6, we obtain the analytical representation for operators that link the near- and far-field impulsive fields for compact axially-symmetrical structures and consider solutions that allow the use of extended or distant sources. In Section 7, we place some accessory results required for numerical implementation of the approach under consideration. All analytical results are presented in a form that is suitable for using in the finite-difference method on a finite-sized grid and thus is amenable for software implementation. We develop here the approach initiated in the works by Maikov et al. (1986) and Sirenko et al. (2007) and based on the construction of the exact conditions allowing one to reduce an open problem to an equivalent closed one with a bounded domain of analysis. The derived closed problem can then be solved numerically using the standard finite-difference method (Taflove & Hagness, 2000).

In contrast to other well-known approximate methods involving truncation of the computational domain (using, for example, Absorbing Boundary Conditions or Perfectly Matched Layers), our constructed solution is exact, and may be computationally implemented in a way that avoids the problem of unpredictable behavior of computational errors for large observation times. The impact of this approach is most significant in cases of resonant wave scattering, where it results in reliable numerical data.


2. Formulation of the initial boundary-value problem

In Fig. 1, the cross-section of a model for an open axially-symmetrical (/ϕ0) resonant structure is shown, where {ρ,ϕ,z} are cylindrical and {ρ,ϑ,ϕ} are spherical coordinates. By Σ=Σϕ×[0,2π] we denote perfectly conducting surfaces obtained by rotating the curve Σϕ about the z-axis; Σε,σ=Σϕε,σ×[0,2π]is a similarly defined surface across which the relative

Figure 1.

Geometry of the problem in the half-planeϕ=π/2.

permittivity ε(g) and specific conductivity σ0(g)=η01σ(g) change step-wise; these quantities are piecewise constant inside Ωint and take free space values outside. Here,g={ρ,z}; η0=(μ0/ε0)1/2is the impedance of free space;ε0, and μ0 are the electric and magnetic constants of vacuum.

The two-dimensional initial boundary-value problem describing the pulsed axially-symmetrical TE0n- (Eρ=Ez=Hϕ0) and TM0n- (Hρ=Hz=Eϕ0) wave distribution in open structures of this kind is given by

{[ε(g)2t2σ(g)t+2z2+ρ(1ρρρ)]U(g,t)=F(g,t),t>0,gΩU(g,t)|t=0=φ(g),tU(g,t)|t=0=ψ(g),g={ρ,z}Ω¯Etg(p,t)|p={ρ,ϕ,z}Σ=0,t0Etg(p,t)andHtg(p,t)are continuous when crossingΣε,σU(0,z,t)=0,|z|<,t0D1[U(g,t)Ui(1)(g,t)]|gΓ1=0,D2[U(g,t)]|gΓ2=0,t0,E1

where E={Eρ,Eφ,Ez} and H={Hρ,Hφ,Hz} are the electric and magnetic field vectors; U(g,t)=Eϕ(g,t)for TE0n-waves and U(g,t)=Hϕ(g,t) for TM0n-waves (Sirenko et al., 2007). The SI system of units is used. The variable t which being the product of the real time by the velocity of light in free space has the dimension of length. The operatorsD1, D2will be described in Section 2 and provide an ideal model for fields emitted and absorbed by the waveguides.

The domain of analysis Ω is the part of the half-plane ϕ=π/2 bounded by the contours Σϕ together with the artificial boundaries Γj (input and output ports) in the virtual waveguidesΩj,j=1,2. The regions Ωint={g={r,ϑ}Ω:r<L} and Ωext (free space), such thatΩ=ΩintΩextΓ, are separated by the virtual boundaryΓ={g={r,ϑ}Ω:r=L}.

The functionsF(g,t), φ(g), ψ(g), σ(g), and ε(g)1 which are finite in the closure Ω¯ of Ω are supposed to satisfy the hypotheses of the theorem on the unique solvability of problem (1) in the Sobolev spaceW21(ΩT), ΩT=Ω×(0;T)where T< is the observation time (Ladyzhenskaya, 1985). The ‘current’ and ‘instantaneous’ sources given by the functions F(g,t) andφ(g), ψ(g)as well as all scattering elements given by the functionsε(g), σ(g)and by the contours Σϕ and Σϕε,σ are located in the regionΩint. In axially-symmetrical problems, at points g such thatρ=0, only Hz or Ez fields components are nonzero. Hence it follows thatU(0,z,t)=0;|z|<, t0in (1).


3. Exact absorbing conditions for virtual boundaries in input-output waveguides



in (1) give the exact absorbing conditions for the outgoing pulsed waves Us(1)(g,t)=U(g,t)Ui(1)(g,t) and Us(2)(g,t)=U(g,t) traveling into the virtual waveguides Ω1 andΩ2, respectively (Sirenko et al., 2007). Ui(1)(g,t)is the pulsed wave that excites the axially-symmetrical structure from the circular or coaxial circular waveguideΩ1. It is assumed that by the time t=0 this wave has not yet reached the boundaryΓ1.

By using conditions (2), we simplify substantially the model simulating an actual electrodynamic structure: the Ωj-domains are excluded from consideration while the operators Dj describe wave transformation on the boundaries Γj that separate regular feeding waveguides from the radiating unit. The operators Dj are constructed such that a wave incident on Γj from the region Ωint passes into the virtual domain Ωj as if into a regular waveguide – without deformations or reflections. In other words, it is absorbed completely by the boundaryΓj. Therefore, we call the boundary conditions (2) as well as the other conditions of this kind ‘exact absorbing conditions’.

In the book (Sirenko et al., 2007), one can find six possible versions of the operators Dj for virtual boundaries in the cross-sections of circular or coaxial-circular waveguides. We pick out two of them (one for the nonlocal conditions and one for the local conditions) and, taking into consideration the location of the boundaries Γj in our problem (in the plane z=L1 for the boundary Γ1 and in the plane z=L2 forΓ2) as well as the traveling direction for the waves outgoing through these boundaries (towards z= for Γ1 and towards z= forΓ2), write (2) in the form:


(nonlocal absorbing conditions) and


(local absorbing conditions). The initial boundary-value problems involved in (5) and (6) with respect to the auxiliary functions Wj(ρ,t,φ) must be supplemented with the following boundary conditions for all timest0:


(on the boundaries ρ=0 and ρ=aj of the region Ωj for a circular waveguide) and


(on the boundaries ρ=bj and ρ=aj of the region Ωj for a coaxial waveguide).

In (3) to (8) the following designations are used: J0(x)is the Bessel function, ajand bj are the radii of the waveguide Ωj and of its inner conductor respectively (evidently, bj0if only Ωj is a coaxial waveguide), {μnj(ρ)}and {λnj} are the sets of transverse functions and transverse eigenvalues for the waveguideΩj.

Analytical representations for μnj(ρ) and λnj are well-known and for TE0n-waves take the form:

{μnj(ρ)=J1(λnjρ)2[ajJ0(λnjaj)]1,n=1,2,...,ρ<aj(circular waveguide)λnj>0are the roots of the equationJ1(λaj)=0,E9
{μnj(ρ)=G1(λnj,ρ)2[aj2G02(λnj,aj)bj2G02(λnj,bj)]1/2,n=1,2,...,bj<ρ<aj(coaxial waveguide)λnj>0are the roots of the equationG1(λ,aj)=0Gq(λ,ρ)=Jq(λρ)N1(λbj)Nq(λρ)J1(λbj).E10

For TM0n-waves we have:

{μnj(ρ)=J1(λnjρ)2[ajJ1(λnjaj)]1,n=1,2,...,ρ<aj(circular waveguide)λnj>0 are the roots of the equationJ0(λaj)=0,E11
{μnj(ρ)=G˜1(λnj,ρ)2[aj2G˜12(λnj,aj)bj2G˜12(λnj,bj)]1/2,n=1,2,...μ0(ρ)=[ρln(aj/bj)]1,bj<ρ<aj(coaxial waveguide)λnj>0(n=1,2,...)are the roots of the equationG˜0(λ,bj)=0,λ0j=0G˜q(λ,ρ)=Jq(λρ)N0(λaj)Nq(λρ)J0(λaj).E12

Here Nq(x) are the Neumann functions. The basis functions μnj(ρ) satisfy boundary conditions at the ends of the appropriate intervals (ρ<ajorbj<ρ<aj) and the following equalities hold


in the circular or coaxial waveguide, respectively.


4. Exact radiation conditions for outgoing spherical waves and exact absorbing conditions for the artificial boundary in free space

When constructing the exact absorbing condition for the wave U(g,t) crossing the artificial spherical boundaryΓ, we will follow the sequence of transformations widely used in the theory of hyperbolic equations (e.g., Borisov, 1996) – incomplete separation of variables in initial boundary-value problems for telegraph or wave equations, integral transformations in the problems for one-dimensional Klein-Gordon equations, solution of the auxiliary boundary-value problems for ordinary differential equations, and inverse integral transforms.

In the domainΩext=Ω\(ΩintΓ), where the field U(g,t) propagates freely up to infinity ast, the 2-D initial boundary-value problem (1) in spherical coordinates takes the form


Let us represent the solution U(r,ϑ,t) asU(r,ϑ,t)=u(r,t)μ(ϑ). Separation of variables in (14) results in a homogeneous Sturm-Liouville problem with respect to the function μ˜(cosϑ)=μ(ϑ)


and the following initial boundary-value problem foru(r,t):


Let us solve the Sturm-Liouville problem (15) with respect to μ˜(cosϑ) andλ. Change of variablesx=cosϑ, μ˜(x)=μ˜(cosϑ)yields the following boundary-value problem forμ˜(x):


With λ2=λn2=n(n+1) for each n=1,2,3, equation (17) has two nontrivial linearly independent solutions in the form of the associated Legendre functions Pn1(x) andQn1(x). Taking into account the behavior of these functions in the vicinity of their singular points x=±1 (Bateman & Erdelyi, 1953), we obtain


Here {μ˜n(cosϑ)}n=1,2, is a complete orthonormal (with weight functionsinϑ) system of functions in the space L2[(0<ϑ<π)] and provides nontrivial solutions to (15). Therefore, the solution of initial boundary-value problem (14) can be represented as


where the space-time amplitudes un(r,t) are the solutions to problems (16) forλ2=λn2.

Our goal now is to derive the exact radiation conditions for space-time amplitudes un(r,t) of the outgoing wave (19). By defining wn(r,t)=run(r,t) and taking into account thatλn2=n(n+1), we rewrite equation (16) as


Now subject it to the integral transform


where the kernel Zγ(ω,r)=ra[α(ω)Jγ(ωr)+β(ω)Nγ(ωr)] satisfies the equation (Korn & Korn, 1961)


Here α(ω),β(ω) are arbitrary functions independent ofr, and a is a fixed real constant. Applying to (20) the transform (21) with a=1/2 andγ=n+1/2, we arrive at


Since the ‘signal’ wn(r,t) propagates with a finite velocity, for any t we can always point a distance r such that the signal has not yet reached it, that is, for these t and r we havewn(r,t)0. Then we can rewrite equation (23) in the form


From (24) the simple differential equation for the transforms w˜n(ω,t) of the functions wn(r,t) follows:


In this equation, the values α(ω) and β(ω) entering into Zγ(ω,r) are not defined yet. With α(ω)=1 andβ(ω)=0, we have


The last integral is the Hankel transform (Korn & Korn, 1961), which is inverse to itself, and

By χ we denote the Heaviside step-function
χ(r)={0 for r<01 for r0.E29

Taking into account (26), equation (25) can be rewritten in the form




and the symbol ‘’ denotes derivatives with respect to the whole argumentωL.

If G is a fundamental solution of the operator B[G] (i.e.,B[G(t)]=δ(t) , where δ(t) is the Dirac delta function), then the solution to the equation B[U(t)]=g(t) can be written as a convolution U=(Gg) (Vladimirov, 1971). For [2/t2+ω2]G(t)=δ(t) we haveG(t)=χ(t)ω1sinωt, and then

Applying the inverse transform (28) to equation (32), we can write


Let us denote


Then from (Gradshteyn & Ryzhik, 2000) we have for r>L>0


where Pγ(x) and Qγ(x) are the Legendre functions of the first and second kind, respectively. Forγ=n+1/2, we can rewrite this formula as


where q=[r2+L2(tτ)2]/2rL and Pn(q) denotes a Legendre polynomial. Considering that


(Janke et al., 1960), anddχ(x)/dx=δ(x), we can derive


Finally, taking into account the relationwn(r,t)=run(r,t), we have from (33)


By using (19), we arrive at the desired radiation condition:


By passing to the limit rL in (40), we obtain


Formula (41) represents the exact absorbing condition on the artificial boundaryΓ. This condition is spoken of as exact because any outgoing wave described by the initial problem (1) satisfies this condition. Every outgoing wave U(g,t) passes through the boundary Γ without distortions, as if it is absorbed by the domain Ωext or its boundaryΓ. That is why this condition is said to be absorbing.


5. On the equivalence of the initial problem and the problem with a bounded domain of analysis

We have constructed the following closed initial boundary-value problem

{[ε(g)2t2σ(g)t+2z2+ρ(1ρρρ)]U(g,t)=F(g,t),t>0,gΩintU(g,t)|t=0=φ(g),tU(g,t)|t=0=ψ(g),gΩ¯intEtg(p,t)|p={ρ,ϕ,z}Σ=0,U(0,z,t)=0,|z|L,t0Etg(p,t)andHtg(p,t)are continuous when crossingΣε,σD1[U(g,t)Ui(1)(g,t)]|gΓ1=0,D2[U(g,t)]|gΓ2=0,t0D[U(g,t)]|gΓ=0,E42

where the operator D is given by (41). It is equivalent to the open initial problem (1). This statement can be proved by following the technique developed in (Ladyzhenskaya, 1985). The initial and the modified problems are equivalent if and only if any solution of the initial problem is a solution to problem (42) and at the same time, any solution of the modified problem is the solution to problem (1). (In the Ωext-domain, the solution to the modified problem is constructed with the help of (40).) The solution of the initial problem is unique and it is evidently the solution to the modified problem according construction. In this case, if the solution of (42) is unique, it will be a solution to (1). Assume that problem (42) has two different solutions U1(g,t) andU2(g,t). Then the function u(g,t)=U1(g,t)U2(g,t) is also the generalized solution to (42) forF(g,t)=Ui(1)(g,t)=φ(g)=ψ(g)0. This means that for any function γ(g,t)W21(ΩT) that is zero att=T, the following equality holds:


Here, ΩintT=Ωint×(0,T)and ΣintT are the space-time cylinder over the domain Ωint and its lateral surface; cos(n,ρ)and cos(n,z) are the cosines of the angles between the outer normal n to the surface ΣintT and ρ- and z-axes, respectively; the element dg of the end surface of the cylinder ΩintT equalsρdρdz.

By making the following suitable choice of function,


it is possible to show that every term in (43) is nonnegative (Mikhailov, 1976) and therefore u(g,t) is equal to zero for all gΩint and0<t<T, which means that the solution to the problem (42) is unique. This proves the equivalency of the two problems.


6. Far-field zone problem. Extended and remote sources

As we have already mentioned, in contrast to approximate methods based on the use of the Absorbing Boundary Conditions or Perfectly Matched Layers, our approach to the effective truncation of the computational domain is rigorous, which is to say that the original open problem and the modified closed problem are equivalent. This allows one, in particular, to monitor a computational error and obtain reliable information about resonant wave scattering. It is noteworthy that within the limits of this rigorous approach we also obtain, without any additional effort, the solution to the far-field zone problem, namely, of finding the field U(g,t) at arbitrary point in Ωext from the magnitudes of U(g,t)on any arcr=ML, 0ϑπ, lying entirely in Ωint and retaining all characteristics of the arcΓ. Thus in the case considered here, equation (39) defines the diagonal operator SLr(t) such that it operates on the space of amplitudes u(r,t)={un(r,t)} of the outgoing wave (19) according the rule


and allows one to follow all variations of these amplitudes in an arbitrary region ofΩext. The operator


given by (40), in turn, enables the variations of the fieldU(g,t), gΩext, to be followed.

It is obvious that the efficiency of the numerical algorithm based on (42) reduces if the support of the function F(g,t) and/or the functions φ(g) and ψ(g) is extended substantially or removed far from the region where the scatterers are located. The arising problem (the far-field zone problem or the problem of extended and remote sources) can be resolved by the following straightforward way.

Let us consider the problem

{[ε(g)2t2σ(g)t+2z2+ρ(1ρρρ)]U(g,t)=F(g,t)+F˜(g,t),t>0,gΩU(g,t)|t=0=φ(g)+φ˜(g),tU(g,t)|t=0=ψ(g)+ψ˜(g),g={ρ,z}Ω¯Etg(p,t)|p={ρ,ϕ,z}Σ=0,t0Etg(p,t)andHtg(p,t)are continuous when crossingΣε,σU(0,z,t)=0,|z|<,t0D1[U(g,t)Ui(1)(g,t)]|gΓ1=0,D2[U(g,t)]|gΓ2=0,t0,E47

which differs from the problem (1) only in that the sources F˜(g,t)andφ˜(g), ψ˜(g)are located out of the domain Ωint enveloping all the scatterers (Fig. 1). The supports of the functionsF˜(g,t), φ˜(g), and ψ˜(g) can be arbitrary large (and even unbounded) and are located in Ωext at any finite distance from the domainΩint.

Let the relevant sources generate a field Ui(g,t) in the half-planeΩ0={g:ρ>0,|z|<}. In other words, let the function Ui(g,t) be a solution of the following Cauchy problem:


It follows from (47), (48) that in the domain Ωext the function Us(g,t)=U(g,t)Ui(g,t) satisfies the equations


and determines there the pulsed electromagnetic wave crossing the artificial boundary Γ in one direction only, namely, from Ωint intoΩext.

The problems (49) and (14) are qualitatively the same. Therefore, by repeating the transformations of Section 4, we obtain


or, in the operator notations, D[U(g,t)Ui(g,t)]|gΓ=0, – the exact absorbing condition allowing one to replace open problem (47) with the equivalent closed problem

{[ε(g)2t2σ(g)t+2z2+ρ(1ρρρ)]U(g,t)=F(g,t),t>0,gΩintU(g,t)|t=0=φ(g),tU(g,t)|t=0=ψ(g),gΩ¯intEtg(p,t)|p={ρ,ϕ,z}Σ=0,U(0,z,t)=0,|z|L,t0Etg(p,t)andHtg(p,t)are continuous when crossingΣε,σD1[U(g,t)Ui(1)(g,t)]|gΓ1=0,D2[U(g,t)]|gΓ2=0,t0D[U(g,t)Ui(g,t)]|gΓ=0.E51

7. Determination of the incident fields

To implement the algorithms based on the solution of the closed problems (42), (51), the values of the functions Ui(1)(g,t) and Ui(g,t) as well as their normal derivatives on the boundaries Γ1 and Γ are required (see formulas (3), (5), (50)). Let us start from the functionUi(1)(g,t). In the feeding waveguideΩ1, the field Ui(1)(g,t) incoming on the boundary Γ1 can be represented (Sirenko et al., 2007) as


Here (see also Section 3), n=0,1,2,...only in the case of TM0n-waves and only for a coaxial waveguideΩ1. In all other casesn=1,2,3,.... On the boundaryΓ1, the wave Ui(1)(g,t) can be given by a set of its amplitudes{vn1(L1,t)}n. The choice of the functionsvn1(L1,t), which are nonzero on the finite interval0<T1tT2<T, is arbitrary to a large degree and depends generally upon the conditions of a numerical experiment. As for the set{vn1(z,t)/z|z=L1}n, which determines the derivative of the functiоn Ui(1)(g,t) onΓ1, it should be selected with consideration for the causality principle. Each pair Vn1(ρ,t)={vn1(L1,t)μn1(ρ);(vn1(z,t)/z|z=L1)μn1(ρ)}is determined by the pulsed eigenmode Uni(1)(g,t) propagating in the waveguide Ω1 in the sense of increasingz. This condition is met if the functions vn1(L1,t) and vn1(z,t)/z|z=L1 are related by the following equation (Sirenko et al., 2007):


The function Ui(g,t) generated by the sourcesF˜(g,t), φ˜(g), and ψ˜(g) is the solution to the Cauchy problem (48). Let us separate the transverse variable ρ in this problem and represent its solution in the form (Korn & Korn, 1961):


In order to find the functionsvλ(z,t), one has to invert the following Cauchy problems for one-dimensional Klein-Gordon equations:


Here, Fλ(z,t), φλ(z)и ψλ(z,t) are the amplitude coefficients in the integral presentations (54) for the functionsF˜(g,t), φ˜(g), andψ˜(g).

Now, by extending the functions Fλ(z,t) and vλ(z,t) with zero on the intervalt<0, we pass on to a generalized version of problems (56) (Vladimirov, 1971)

(δ(1)(t)is the generalized derivative of the functionδ(t)). Their solutions can be written by using the fundamental solution G(z,t,λ)=(1/2)χ(t|z|)J0[λ(t2z2)1/2] of the operator B(λ) as follows:

Equations (54) and (58) completely determine the desired functionUi(g,t).


8. Conclusion

In this paper, a problem of efficient truncation of the computational domain in finite-difference methods is discussed for axially-symmetrical open electrodynamic structures. The original problem describing electromagnetic wave scattering on a compact axially-symmetric structure with feeding waveguides is an initial boundary-value problem formulated in an unbounded domain. The exact absorbing conditions have been derived for a spherical artificial boundary enveloping all sources and scatterers in order to truncate the computational domain and replace the original open problem by an equivalent closed one. The constructed solution has been generalized to the case of extended and remote field sources. The analytical representation for the operators converting the near-zone fields into the far-zone fields has been also derived.

We would like to make the following observation about our approach.

In our description, the waveguide Ω1 serves as a feeding waveguide. However, both of the waveguides can be feeding or serve to withdraw the energy; also both of them may be absent in the structure.

The choice of the parameters α(ω) and β(ω) determining Zγ(ω,r) (see Section 4) affects substantially the final analytical expression for the exact absorbing condition on the spherical boundaryΓ. When constructing boundary conditions (41), (50), we assumed that α(ω)=1 andβ(ω)=0. In (Sirenko et al., 2007), for a similar situation, the exact absorbing conditions for outgoing pulsed waves were constructed with the assumption that α(ω)=Nγ(ωL) andβ(ω)=Jγ(ωL). With such α(ω) andβ(ω), equation (21) is the Weber-Orr transform (Bateman & Erdelyi, 1953). However, the final formulas corresponding to (39), (40) for this case turn into identities asrL, which present a considerable challenge for using them as absorbing conditions. In addition, the analytical expressions with the use of Weber-Orr transform are rather complicated to implement numerically.

The function Ui(g,t) (see Section 7) can be found in spherical coordinates as well. In this situation, we arrive (see Section 4) at the expansions like (19) with the amplitude coefficients vn(r,t) determined by the Cauchy problems


whereFn(r,t), φn(r), and ψn(r) are the amplitude coefficients for the functionsF˜(g,t), φ˜(g), andψ˜(g).

The standard discretization of the closed problems (42), (51) by the finite difference method using a uniform rectangular mesh attached to coordinates g={ρ,z} leads to explicit computational schemes with uniquely defined mesh functionsU(j,k,m)=U(ρj,zk,tm). The approximation error isO(h¯2), where h¯ is the mesh width in spatial coordinates, l¯=h¯/2for θ=maxgΩint[ε(g)]<2 or l¯<h¯/2 for θ2 is the mesh width in time variablet;ρj=jh¯, zk=kh¯, andtm=ml¯. The range of the integersj=0,1,...,J, k=0,1,...K, and m=0,1,...M depends both on the size of the Ωint domains and on the length of the interval [0,T] of the observation timet. The condition providing uniform boundedness of the approximate solutions U(j,k,m) with decreasing h¯ and l¯ is met (see, for example, formula (1,50) in (Sirenko et al., 2007)). Hence the finite-difference computational schemes are stable, and the mesh functions U(j,k,m) converge to the solutions U(ρj,zk,tm) of the original problems (42), (51).

As opposed to the well-known approximate boundary conditions standardly utilized by finite-difference methods, the conditions derived in this paper are exact by construction and do not introduce an additional error into the finite-difference algorithm. This advantage is especially valuable in resonant situations, where numerical simulation requires large running time and the computational errors may grow unpredictably if an open problem is replaced by an insufficiently accurate closed problem.


  1. 1. BatemanH.ErdelyiA. 1953 Higher transcendental functions, McGraw-Hill, ISBN 048644614X, New York
  2. 2. BorisovV. V. 1996 Electromagnetic fields of transient currents, St. Petersburg Univ. Press, 5-28801-256-3 Petersburg
  3. 3. GradshteynI. S.RyzhikI. M. 2000 Table of Integrals, Series, and Products, 0-12294-757-6 Press, San Diego, London
  4. 4. JankeE.EmdeF.LöschF. 1960 Tables of Higher Functions, McGraw-Hill, 0-07032-245-7 York
  5. 5. KornG. A.KornT. M. 1961 Mathematical Handbook for Scientists and Engineers, Mc Graw-Hill, 0-07035-370-0 York
  6. 6. LadyzhenskayaO. A. 1985 The boundary value problems of mathematical physics, Springer-Verlag, 3-54090-989-3 York
  7. 7. MaikovА. R.SveshnikovA. G.YakuninS. A. 1986 Difference scheme for the Maxwell transient equations in waveguide systems. Journal of Computational Mathematics and Mathematical Physics, 26 851863 , 0965-5425
  8. 8. MikhailovV. P. 1976 Partial Differential Equations. Nauka, 20203-032-13-76
  9. 9. SirenkoY. K.StromS.YashinaN. P. 2007 Modeling and Analysis of Transient Processes in Open Resonant Structures. New Methods and Techniques, Springer, 0-387-30878-4New York
  10. 10. TafloveA.HagnessS. C. 2000 Computational Electrodynamics: the Finite-Difference Time-Domain Method, Artech House, 0-89006-792-9
  11. 11. VladimirovV. S. 1971 Equations of mathematical physics, Dekker, 978-0-82471-713-1 New York

Written By

Olena Shafalyuk, Yuriy Sirenko and Paul Smith

Submitted: 27 October 2010 Published: 21 June 2011