High-order methods for computational wave propagation and scattering

The American Institute of Mathematics

The following compilation of participant contributions is only intended as a lead-in to the AIM workshop “High-order methods for computational wave propagation and scattering.”

This material is not for public distribution.

Corrections and new material are welcomed and can be sent to workshops@aimath.org Version: Tue Aug 28 18:34:49 2007

Table of Contents

A. Participant Contributions . . . . . . . . . . . . . . . . . . . . . . . . 3 1. Beale, J. Thomas 2. Bokil, Vrushali

3. Bruno, Oscar

4. Cakoni, Fioralba

5. Chu, Hanyou

6. Colton, David

7. Dauge, Monique

8. Davies, Penny

9. Ganesh, Mahadevan 10. Gibson, Nathan

11. Hagstrom, Thomas

12. Hiptmair, Ralf

13. Kress, Rainer

14. Lechleiter, Armin

15. Nicholls, David

16. Pieper, Martin

17. Smirnova, Alexandra

18. Strain, John

19. Turc, Catalin

Chapter A: Participant Contributions

A.1 Beale, J. Thomas

I have developed a relatively simple and direct approach for computing singular and nearly singular integrals such as those representing layer potentials on a closed surface in 3D or a closed curve in 2D. At this workshop I hope to see if this approach can be useful in electromagnetic problems. The papers of mine listed below can be found at my web site,

www.math.duke.edu/faculty/beale as well as at the journal sites.

The essence of this method is to replace the singularity with a regularized version, discretize with a standard quadrature, and add corrections based on local analysis near the singularity. Uniform accuracy can be obtained for values on or close to the surface. Typically the error is about O(h^{3}) or better. For a closed surface in 3D the integrals are reduced to sums over regular grid points in overlapping coordinate patches with a partition of unity. No special grid is needed near the singularity. The accuracy depends on coordinate patches without severe distortion. If the coordinate systems are good enough in this sense, only a correction for the regularization is needed, and the quadrature error is smaller. I have concentrated on potentials defining harmonic functions, but a similar approach can be applied to the Helmholtz equation with qualifications. This approach seems to have certain advantages due to its simplicity, but it does not appear straightforward to apply it to cases with corners and edges.

In a recent Ph.D. project with me, Michael Nicholas used such an approach to calculate electromagnetic waves of definite frequency scattered from a doubly periodic material representing a photonic crystal. This work was closely related to work of my colleague Stephanos Venakides and his coworkers, referenced below. They have studied resonances and near-resonances theoretically and computationally in 2D. Nicholas calculated scattered waves and transmission coefficients in some 3D cases with normal incidence and with simple geometry (a slab with one coordinate a function of the periodic coordinates). He found strong dependence on frequency, especially associated with Wood’s anomaly.

Two other current graduate students are working on related projects. Matthew Surles is developing a method for potential integrals over piecewise smooth curves, without using regularization. Jason Wilson is developing an algorithm for introducing coordinates on a given closed surface with low distortion (first fundamental form not far away from the identity).

1. J. T. Beale, A grid-based boundary integral method for elliptic problems in three dimensions, SIAM J. Numer. Anal., 42 (2004), 599-620.

2. J. T. Beale and M.-C. Lai, A method for computing nearly singular integrals, SIAM J. Numer. Anal. 38 (2001), 1902-25.

3. J. T. Beale, A convergent boundary integral method for three-dimensional water waves, Math. Comp. 70 (2001), 977-1029.

4. S. Shipman and S. Venakides, Resonant transmission near nonrobust periodic slab modes, Phys. Rev. E, 71 (2005), pp. 026611(1-10).

5. S. Shipman and S. Venakides, Resonance and bound states in photonic crystal slabs, SIAM J. Appl. Math., 64 (2003), 322-342.

A.2 Bokil, Vrushali

I am interested in developing efficient numerical methods for wave propagation problems in dispersive media, with an emphasis on understanding scattering in complex biological tissues in applications such as tumor detection. Specifically, I intend to investigate higher order finite element methods for time domain wave propagation in general dispersive media, homogenization techniques for complex heterogeneous materials, and related inverse problems for damage detection or material characterization.

A.3 Bruno, Oscar

My research efforts relate to the development of numerical methods and related theoretical aspects concerning fast, high-order solution of Partial Differential Equations. The pursuit of such theoretical issues and algorithm development for problems of acoustic and electromagnetic scattering, with an emphasis on applicability to complex engineering geometries and high-frequency configurations, constitute a significant focus of my work in recent years - thus my research interests are closely aligned with the topics to be considered in this workshop.

A.4 Cakoni, Fioralba

Primarily, I am interested in theoretical and computational aspects of inverse scattering problems in acoustic and electromagnetic wave propagation. Im also interested in computational methods in wave propagation based on integral equation approaches. I am expecting to hear a few presentations on the above areas. Particularly important to me in this workshop will be group, panel and in pair extensive discussions about topics and specific problems that interest me.

A.5 Chu, Hanyou

We are interested in efficient treatment of corners, edges, and interface of more than two materials.

A.6 Colton, David

My primary interest is in the numerical solution of the inverse scattering problem for acoustic and electromagnetic waves. I am particularly interested in sampling and probe methods for the reconstruction of the support of the scatterer. Since there are many participants in the workshop with similar interests, I would be interested in a discussion(perhaps a panel discussion?) of the relevant merits and/or difficulties associated with the various sampling and probe methods that are available. Of particular interest here are buried objects, limited aperature data and possible “target signatures”.

A.7 Dauge, Monique

My interest in relation with the workshop concerns harmonic Maxwell equations for possibly inhomogeneous materials, with special focus on asymptotics and approximation issues in presence of singularities or small parameters.

I have contributed, jointly with M. Costabel and other co-authors, several papers on the corner-edge structure of solutions, and on their approximation by standard nodal elements (the weighted regularization).

Presently, with a post-doc and a PhD student, we are investigating the skin effect in presence of highly conducting material, and the amplification of the electromagnetic field around a thin conducting wire exposed to the magnetic field of a MRI device. The latter question arises in medical applications where some unwanted heating of catheters have been observed.

From a theoretical point of view, these questions are related to the computations of complex resonances, and to the diffraction by a thin straight wire. I expect to discuss these topics with participants in the workshop.

References: Most of my papers can be retrieved from my web-page

http://perso.univ-rennes1.fr/monique.dauge

A.8 Davies, Penny

My particular interest is in the numerical approximation of time-dependent boundary integral equations. Especially in developing approximations which are stable and high order in time.

A.9 Ganesh, Mahadevan

Understanding of many physical phenomena and processes in atmospheric science, climatology, and astronomy can be enhanced through simulation of scattered electromagnetic waves from ensembles of particles such as atmospheric aerosols, dust in planetary rings, or interstellar dust. Computer simulations of the radar cross section (RCS) of three dimensional model targets through electromagnetic scattering are a cost effective tool for designing stealth objects. Such simulations are also useful in medical sciences. For example, it is common to use scattering measurements in laboratories to count erythrocytes (red bloods cells), for image-guided neurosurgery using brain tumor models, for detection and classification of tumorous tissue in the breast, for monitoring of the oxygenation level in infant brain tissue, and in functional brain activation studies. The forward and inverse electromagnetic scattering simulations are also required in defense science and efficient materials testing and in detection and identification of visually obscured targets.

Contribution of M. Ganesh in the workshop will be on the development, analysis, and implementation of spectrally accurate high-order algorithms to simulate the interaction of electromagnetic waves by a collection of three dimensional obstacles. The three dimensional, electromagnetic-size and shape, and changes in configurations aspects of the electromagnetic scattering problems require spectrally accurate forward scattering algorithms to efficiently simulate the monostatic RCS of the ensemble and for practical realization of the accurate and reduced data based inverse obstacle scattering iterative Newton-type methods, to reconstruct the shape of unknown scatterers using a limited class of electric far field informational.

A.10 Gibson, Nathan

I have interests in developing and analyzing numerical methods for wave propagation problems in scattering media, with an emphasis on applications of Terahertz frequency interrogation. Specifically, I plan to study higher order finite element methods for time domain wave propagation in powders and foam-like materials. I intend to utilize simulations as forward solves in an inverse problem formulation for damage characterization or substance identification.

A.11 Hagstrom, Thomas

The fundamental technological and scientific importance of wave propagation phenomena is the fact that waves propagate long distances relative to their characteristic dimension, the wavelength, carrying information about the objects which scattered them. It is precisely this fact that makes the numerical simulation of waves so challenging. Even if discretization methods of optimal efficiency can be developed, the need to alot two degrees-of-freedom per wavelength may be prohibitive. We thus require algorithms which allow us to avoid sampling waves away from regions where strong scattering processes are taking place, relying on efficient (unsampled) representations of the solution. In the frequency domain a number of such algorithms exist, for example the fast-multipole algorithm [FMMlow] and equivalent source algorithms [BrunoKunyansky01].

The most limited problem of this type is that of imposing accurate near-field radiation boundary conditions for time-domain scattering problems. We believe that we have a solution to this problem with optimal computational complexity (storage and work/step proportional to ln \frac{cT}{λ}) for isotropic systems which are essentially equivalent to the scalar wave equation. It is based on complete plane wave representations involving both propagating and evanescent modes [complete-theory,complete-Max]. These take the form:

u(x, y, z, t) = ∫0^{\frac{π}{2}} Φ (ct − x cos ϕ, y, z, ϕ) dϕ  (1)
+ ∫_{0}^{∞} e^{−σx}Υ(y, z, t, σ)dσ,  (2)

where u satisfies the wave equation. Applying an appropriate quadrature rule we obtain the finite representation:

u(x, y, z, t) ≈ ∑_{j=0}^{n_{p}−1} w_{j}Φ_{j}(ct − x cos θ_{j}, y, z) + ∑_{j=1}^{n_{e}} d_{j}e^{−σ_{j}x}Υ_{j}(t, y, z).  (3)

Natural questions are:

i.: Can we stably extract the amplitudes in the plane wave representations to build a working algorithm of the type we envision for problems with multiple scatterers? An issue is the fact that each term in the plane wave representation is noncausal. We note that Grote and Kirsch [Grotemult] have constructed an effective multiple scattering algorithm based on multipole expansions, but these do not seem to possess the optimal complexity or geometric flexibility of the plane wave approach. Michielssen’s PWFTD algorithm [PWFTD] is a fast time-domain algorithm based on plane waves. However in this algorithm the evanescent spectrum is avoided via careful space/time decompositions of the sources.

ii.: Is it better to use representations where each term is causal? For example one could use equivalent source representations which are currently being studied by Bruno and Hoch [BrunoHoch]. These take the form:

u(x, y, z, t) ≈ ∑_{j=0}^{n_{s}} \frac{Ψ_{j}(ct − R_{j})}{R_{j}}.  (4)

iii.: Can methods developed for isotropic problems be extended to anisotropic or inhomogeneous systems? In the former case it is possible for the directions of the phase and group velocities to be mismatched, in which case the plane wave representations look strange. (Outgoing waves are represented by apparently incoming plane waves.) Can alternative representations directly involving group velocities be found? Could something like curvelets [CurveletW] prove useful?

iv.: Can we complete the mathematical analysis of any of these methods? For example, can we provide an analysis of the PWFTD algorithm or the complete radiation boundary conditions which seem to work well in prectice?

Bibliography

[FMMlow] L. Greengard, J. Huang, V. Rokhlin, and S. Wandzura, Accelerating fast multipole methods for low frequency scattering, IEEE Comput. Sci. Engng. 5 (1998), 32–38. [BrunoKunyansky01] O. Bruno and L. Kunyansky, A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications, J. Comput. Phys. 169 (2001), 80–110. [complete-theory] T. Hagstrom and T. Warburton, Complete radiation boundary conditions: minimizing the long time error growth of local methods, in preparation. [complete-Max] T. Warburton and T. Hagstrom, Complete radiation boundary conditions for electromagnetic scattering, in preparation. [Grotemult] M. Grote and C. Kirsch, Nonreflecting boundary condition for time dependent multiple scattering, J. Comput. Phys. 221 (2007), 41–62.

[PWFTD] A. Ergin, B. Shanker, and E. Michielssen, Fast evaluation of three-dimensional transient wave fields using diagonal translation operators, J. Comput. Phys. 146 (1998), 157–180. [BrunoHoch] O. Bruno and D. Hoch, Exact nonreflecting boundary conditions based on equivalent sources for time dependent scattering problems, In preparation. [CurveletW] E. Candés and L. Demanet, The curvelet representation of wave propagators is optimal, Comm. Pure and Appl. Math. 58 (2005), 1472–1528.

A.12 Hiptmair, Ralf

(joint work with Ilaria Perugia, Dipartimento di Matematica, Universitá di Pavia, e–mail: ilaria.perugia@unipv.it)

Standard low order Lagrangian finite element discretization of boundary value problems for the Helmholtz equation −Δu − ω^{2}u = f are afflicted with the so-called pollution phenomenon: though for sufficiently small hω an accurate approximation of u is possible, the Galerkin procedure fails to provide it. Attempts to remedy this have focused on incorporating extra information in the form of plane wave functions x ↦ exp(iωd · x), |d| = 1, into the trial spaces. Prominent examples of such methods are the plane wave partition of unity finite element method of Babuska and Melenk [BAM97], and the ultra-weak Galerkin discretization due to Cessenat and Despres [CED98,CED03]. Both perform well in computations, see the articles by Monk and Hutunen [HMK02,HUM07,HMM06] for computational results for the ultra-weak approach.

It turns out that the latter method can be recast as a special so-called discontinuous Galerkin (DG) method employing local trial spaces spanned by a few plane waves. In a sense, this is a special case of a Trefftz-type approximation. This perspective paves the way for marrying plane wave approximation with many of the various DG methods developed for 2nd-order elliptic boundary value problems. Ilaria Perugia and I have pursued this for a generic mixed DG method and a primal DG method, which generalize the ultra-weak scheme.

For these methods we have developed a convergence analysis for the h-version, which achieves convergence through mesh refinement. Key elements are approximation estimates for plane waves and sophisticated duality techniques. The latter entail estimating how well local plane waves can approximate the solution of a dual problem. Unfortunately, we could not help invoking general polynomial estimates in Sobolev spaces for this purpose. This incurs unsatisfactory pollution-affected final error bounds O(ω^{2}h). On the other hand, in 1D a detailed analysis confirms that the plane wave DG method does not suffer from pollution.

However, 1D is a very special case, because there are only two possible wave directions, which can always be resolved. Conversely, in higher dimensions only a few of infinitely many wave directions will be represented in the trial space, usually equidistributed on unit circle/unit sphere. If the solution of a boundary value problem contains plane wave components falling in between these directions, numerical dispersion will be incurred. Then the pollution phenomenon also affect the plane wave method.

In the course of this recent research we have encountered many open problems:

• Can the above reasoning be made rigorous to show the inevitability of pollution also for plane wave based volume discretizations for the Helmholtz equation, c.f. [BAS00] ?

• Can non-linear, e.g., self-adaptive, plane wave approximation provide a remedy for pollution?

• What about properties and convergence of a “p-version” plane wave DG scheme?

Bibliography

[BAM97] I. Babuška and J. Melenk, The partition of unity method, Int. J. Numer. Methods Eng., 40 (1997), pp. 727–758.

[BAS00] I. Babuška and S. Sauter, Is the pollution effect of the FEM avoidable for the Helmholtz equation?, SIAM Review, 42 (2000), pp. 451–484.

[CED98] O. Cessenat and B. Despres, Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz equation, SIAM J. Numer. Anal., 35 (1998), pp. 255–299.

[CED03] , Using plane waves as base functions for solving time harmonic equations with the ultra weak variational formulation, J. Computational Acoustics, 11 (2003), pp. 227–238.

[HMM06] T. Huttunen, M. Malinen, and P. Monk, Solving maxwell’s equations using the ultra weak variational formulation, J. Comp. Phys., (2006).

[HUM07] T. Huttunen and P. Monk, The use of plane waves to approximate wave propagation in anisotropic media, J. Computational Mathematics, 25 (2007), pp. 350–367.

[HMK02] T. Huttunen, P. Monk, and J. Kaipio, Computational aspects of the ultraweak variational formulation, J. Comp. Phys., 182 (2002), pp. 27–46.

A.13 Kress, Rainer

My main research interests are in inverse obstacle scattering for time-harmonic acoustic and electromagnetic waves. Form the workshop I expect a fruitful interaction on efficient methods for solving boundary integral equations for forward obstacle scattering problems and their application in inverse algorithms.

For an overview on my research I refer to the monograph

Inverse Acoustic and Electromagnetic Scattering Theory. Applied Mathematical Sciences 93, Springer-Verlag, Heidelberg, 2nd Edition (1998) (with D. Colton). and the two recent surveys

Using fundamental solutions in inverse scattering. Inverse Problems 22, R49-R66 (2006) (with D. Colton).

Uniqueness and numerical methods in inverse obstacle scattering. In: Inverse Problems in Applied Sciences–towards Breakthrough (Cheng et al, eds). J. Physics Conference Series 73, 012003 (2007).

A.14 Lechleiter, Armin

I’m interested in direct and inverse scattering problems in the frequency domain. Concerning inverse problems, a lot of my work deals with Factorization methods, for instance contamination detection in rough surface scattering or reconstruction of periodic penetrable surfaces. Both of the latter problems are motivated by physicists working on photonic crystals. Through this physicists group I also got interested in waveguide scattering and one goal would is the analysis of scattering in a periodic waveguide.

A.15 Nicholls, David

To date my research has focused on the numerical and analytical aspects of linear and nonlinear wave motion in fluid mechanics (water waves), linear acoustics, and electromagnetics. In particular I have worked on “Boundary Perturbation” methods which rely upon analyticity properties of field quantities to justify highly accurate numerical algorithms for the simulation of these physical problems. However, these methods, as applied to acoustic and electromagnetic problems, are not specially designed for high-frequency radiation in their original formulation. Recently, in collaboration with F. Reitich (Minnesota), we have shown how these algorithms can be extended to high-frequency configurations using a phase extraction approach similar in spirit to the recent work of Bruno, Geuzaine, Munro, and Reitich. Thus far we have developed the algorithm for shallow diffraction gratings which precludes the possibility of multiple scattering, and are now investigating how to build an iterative algorithm (e.g. based upon Geometrical Optics) to account for multiple reflections. Given the expertise and interests of the conference participants this seems to be an ideal opportunity to make significant progress on this project.

A.16 Pieper, Martin

My work in the last years was focused on spectral boundary integral methods for electromagnetic scattering problems. Especially I’m interested in methods, which use so-called vector spherical harmonics. These special vectorfields can be seen as the vectorial analogon to the common scalar spherical harmonics (see Freeden, Gervens, Schreiner “Constructive Approximation on the Sphere”). Some participants of the workshop have published papers on such methods, so it could be interesting to compare the different ideas and results.

To get convergence results for these spectral methods, one needs a sharp norm estimate of the used hyperinterpolation operator. Up to now an estimate, which is sharp enough, is not yet available and it is questionable, if it is actually possible to get such an estimate.

There are some recent results on this topic, but there are also still a lot of open questions, which are worthwhile to discuss.

References:

1) “Vector hyperinterpolation on the sphere”, submitted to Journal of Approximation Theory, preprint available on http://www.num.math.uni-goettingen.de/bibi/preprintserie/2007.html

A.17 Smirnova, Alexandra

Assume that a nonlinear operator F acts on a pair of Hilbert spaces (H, H_{1}), that is F : D(F) ⊂ H → H_{1}, and that F is Fréchet differentiable in its domain D(F) without such structural assumptions as monotonicity, or invertibility of F^{′∗}( ·)F′( ·), etc. We consider a general problem of minimizing the functional

J(q) := ||F (q) − g_{δ}||_{H_{1}}^{2},  (5)

where g_{δ} approximates the exact data g with the accuracy δ, i.e.,

||g − g_{δ}|| ≤ δ.  (6)

(5) is a typical mathematical model of various inverse problems. One of the most important examples is the problem of a numerical reconstruction of the potential by scattering data. The well-known approaches to theoretical investigation of inverse scattering problem are Gelfand-Levitan, Marchenko and Krein methods. However, since the problem is ill-posed, the development of the corresponding numerical procedures, which must be combined with certain regularization techniques, is still a very complicated task.

My major goal for the upcoming workshop is the development of fast, accurate and admitting parallel implementation, iteratively regularized procedures for nonlinear inverse problem (5). In particular, I am extremely interested in studying theoretical convergence results for an Iteratively Regularized Gauss Newton (IRGN) algorithm with a Tikhonov regularization term using a seminorm generated by a linear operator, as well as other stabilizing preconditioned numerical methods combined with proper stopping rules, line search routines and parameter selection algorithms. I plan to illustrate theoretical results by numerical simulations for the inverse scattering problem, and I am hopeful that the new methods will result in much more accurate and efficient reconstructions.

A.18 Strain, John

Classical potential theory converts linear constant-coefficient elliptic problems in complex domains into integral equations on interfaces, and generates robust, efficient numerical methods. The conversion is usually carried out for a particular situation such as the Poisson equation in dimension 2, and the efficiency of the resulting methods then depends on detailed analysis of the appropriate special functions.

I have recently developed a general conversion scheme which leads naturally to a fast algorithm: arbitrary elliptic problems in arbitrary dimension are converted to first-order systems, a periodic fundamental solution is mollified for convergence, and the mollification is locally corrected via Ewald summation. Local linear algebra and the elementary theory of distributions yield a simple boundary integral equation. With the aid of a new nonequidistant fast Fourier transform for piecewise polynomial functions, the resulting numerical methods provide highly accurate solutions to general elliptic systems in complex domains.

A natural extension of this approach would be to nominally elliptic systems such as the Maxwell equations at fixed frequency, describing electromagnetic wave propagation. Two obvious issues are the complexity of Ewald summation at high frequency, and the extension from Fourier series to Fourier transform required by infinite domains. However, possible payoffs include a simple general derivation for robust boundary integral equations of combined-field type, a natural treatment of complex obstacles, and fast wave propagation algorithms.

A.19 Turc, Catalin

My main interests in this upcoming workshop concern the development of novel wellposed integral equation formulations for acoustic and electromagnetic scattering problems and the design of fast, high-order algorithms for the numerical solution of such problems. In addition, these algorithms should have the capability of producing accurate solutions while only requiring small numbers of Krylov-subspace iterations regardless of the geometric complexity of the scatterers. In this regard, our efforts were directed towards the design of certain uniquely solvable Regularized Combined Field Integral Equations (RCFIE) formulations for the cases of sound-hard acoustic configurations and perfectly conducting electromagnetic applications. The aforementioned RCFIE rely on the use of certain regularizing operators that mollify the unfavorable effects of the pseudodifferential operators of order one that enter the classical CFIE, so that the eigenvalues of the integral operators in RCFIE are asymptotically clustered. One possible methodology for the numerical solution of the surface integral equations underlying the RCFIE formulations uses a high-order Nyström approach based on partitions of unity and high-order quadrature methods which can be accelerated via high-order “two-face” equivalent source approximations. A variety of numerical results demonstrate that for a given accuracy, the new equations can give rise to order-of-magnitude reductions in computational costs over those resulting from previous approaches especially in the case of large scatterers that exhibit geometrically complex features. Part of our recent efforts are focused on extending the high-order methodology outlined above to treat problems of scattering from open surfaces and from objects that present geometric singularities such as edges and corners. Our main idea here is to derive and use the leading order asymptotic behavior of the solutions to the boundary integral equations in the neighborhood of the geometric singularities. Another area of interest is the design of integral equation formulations for acoustic transmission problems and scattering from dielectric problems that would share the desirable featured previously discussed and which could furthermore address the issue of “near resonances” that occurs in these cases.