Spektralmetoder för linjära elliptiska partiella differentialekvationer i fri rymd

IN

DEGREE PROJECT TECHNOLOGY, FIRST CYCLE, 15 CREDITS

STOCKHOLM SWEDEN 2019,

Spectral methods for linear elliptic partial differential equations in free space

ARIEL ABEDALI CHRISTINA GHAWI

KTH ROYAL INSTITUTE OF TECHNOLOGY

INOM

EXAMENSARBETE TEKNIK, GRUNDNIVÅ, 15 HP

STOCKHOLM SVERIGE 2019,

Spektralmetoder för linjära elliptiska partiella

differentialekvationer i fri rymd

ARIEL ABEDALI CHRISTINA GHAWI

KTH

Spectral methods for linear elliptic partial differential equations in free space

ARIEL ABEDALI CHRISTINA GHAWI

SA114X Degree Project in Engineering Physics Date: May 27, 2019

Supervisors: Fredrik Fryklund and Federico Izzo Examiner: Martin Wiklund

School of Engineering Sciences

Swedish title: Spektralmetoder för linjära elliptiska partiella

Abstract

The purpose of this study is to solve linear elliptic partial differential equa- tions in free space and to achieve spectrally accurate numerical solutions, for smooth and compactly supported inhomogeneities in two and three dimen- sions. This is made by using results from theory of Green’s functions and rewriting the differential operator’s free space Green’s function to a trun- cated spectral representation, by utilizing the inhomogeneity’s compact sup- port; then, using results from Fourier analysis and properties of convolution, calculations are performed using a fast Fourier transform. Although partial differential equations often require non trivial solution methods, this pow- erful approach results in a simple and fast way of achieving highly accurate numerical solutions.

Keywords: Partial differential equation, spectral accuracy, fast Fourier transform, compactly supported smooth inhomogeneity, Green’s function

Sammanfattning

Syftet med den här studien är att lösa linjära elliptiska partiella differentialek- vationer i fri rymd och att uppnå spektralkonvergenta numeriska lösningar, för glatta och kompakt stödda inhomogeniteter i två och tre dimensioner.

Detta utförs genom att använda teori om Greenfunktioner och göra en om- skrivning av differentialoperatorns Greenfunktion i fri rymd till en trunkerad spektralrepresentation, genom att nyttja inhomogenitetens kompakta stöd;

därefter, genom att använda resultat från Fourieranalysen och egenskaper av faltning, beräknas lösningen med hjälp av en snabb Fouriertransform. Trots att partiella differentialekvationer ofta kräver icke-triviala lösningsmetoder, resulterar detta kraftfulla tillvägagångssätt i ett simpelt och snabbt sätt att uppnå spektralkonvergenta numeriska lösningar.

Nyckelord: Partiell differentialekvation, spektral noggrannhet, snabb Fouriertransform, kompakt stödd glatt inhomogenitet, Greenfunktion

Acknowledgements

We take this opportunity to express our gratitude to our supervisors, Fredrik Fryklund and Federico Izzo, for their help, guidance and support. We are deeply thankful for their sincere encouragement and continuous support throughout the project.

Stockholm, Sweden Ariel Abedali and Christina Ghawi May 2019

Contents

1 Introduction 1

1.1 Background . . . 1

1.2 Purpose . . . 1

1.3 Problem formulation . . . 1

2 Preliminaries 2 2.1 Theoretical background . . . 2

2.1.1 Partial differential equation . . . 2

2.1.2 Green’s function . . . 3

2.1.3 Fourier analysis . . . 4

2.1.4 Convolution . . . 5

2.2 Numerical background . . . 6

2.2.1 Discretization . . . 6

2.2.2 Discrete Fourier transform . . . 6

2.2.3 Fast Fourier transform . . . 7

2.2.4 Sampling . . . 7

2.2.5 Convergence theory . . . 7

2.2.6 Spectral accuracy . . . 7

3 Method 9 3.1 Theoretical method . . . 9

3.2 Numerical method . . . 14

4 Results 16 4.1 Convergence analysis . . . 16

5 Conclusion 21

References 22

A Truncated spectral representation of Green’s functions 24

B Code listing 30

Chapter 1 Introduction

1.1 Background

A partial differential equation (PDE) is a multivariable differential equa- tion containing partial derivatives of an unknown function. Partial differen- tial equations are the most common method of modelling physical problems.

Hence, they are used in many applications and can be used to describe a wide variety of phenomena e.g. heat transfer, fluid dynamics, electrodynamics and quantum mechanics. It is often impossible to formulate an explicit expression of the analytic solution to a partial differential equation. Therefore, there are several numerical methods for solving PDEs. Though, methods used for truncated domains are rarely suitable in free space. Furthermore, using e.g. finite element method and perfectly matched layers to mimic function behavior at infinity in a truncated domain is difficult to implement and re- quires a high computational cost. Therefore, it is desirable to solve the PDE numerically using a simple, effective and accurate method in free space.

1.2 Purpose

The purpose of this bachelor thesis is to solve linear elliptic PDEs in two and three dimensions by implementing the method proposed by Vico et al.

(2016) and achieving spectrally accurate numerical solutions.

1.3 Problem formulation

Solving the Poisson equation and the inhomogeneous Helmholtz equation in free space with spectral accuracy in two and three dimensions, for smooth and compactly supported inhomogeneities on an ideal domain.

Chapter 2

Preliminaries

This chapter aims to describe the essential theory needed in order to im- plement the method and understand the results. Both the theoretical and numerical theory are presented.

2.1 Theoretical background

In order to analytically solve linear elliptic partial differential equations, one must first get familiar with the mathematical theory behind it. Since the mathematical theory is technical, only the most essential definitions and theorems are presented. First, a short introduction to PDE is given, followed by Green’s functions theory. Furthermore, theory from Fourier analysis is presented, including the Convolution theorem.

2.1.1 Partial differential equation

Definition 2.1. A partial differential equation of the function u(x) = u(x₁, ...x_n), is an equation on the form

f

x₁, ..., x_n; u, ∂u

∂x₁, ..., ∂u

∂x_n; ∂²u

∂x₁∂x₁, ..., ∂²u

∂x₁∂x_n; ...

= 0, (2.1) where f is a known function and u is an unknown function.

In this thesis, two linear elliptic partial differential equations will be studied in Rⁿ, n = 2, 3: the Helmholtz equation and the Poisson equation.

The inhomogeneous Helmholtz equation represents a time independent form of the wave equation, defined as

Definition 2.2 (Helmholtz’ equation). Helmholtz’ equation is defined as

∆u(x) + k²u(x) = −f (x), (2.2) where ∆ is the Laplace operator, k is the wave number and f : Rⁿ−→ C is a function with compact support, i.e. f is zero outside of a compact set.

The Poisson equation is a generalization of the Laplace equation and used in models containing electric potentials. The Poisson equation can be obtained using k = 0 in (2.2), as

Definition 2.3 (Poisson’s equation). Poisson’s equation is defined as

∆u(x) = −f (x), (2.3)

where f : Rⁿ −→ C is a function with compact support.

2.1.2 Green’s function

Assume that u ∈ C²(Rⁿ) is the solution of the problem

L_xu = −f in Rⁿ, (2.4)

where Lx is a linear differential operator, n = 2, 3. Then, the following theorem can be stated

Theorem 2.1. The solution to (2.4) is u(x) =

Z

Rⁿ

G(x − y)f (y)dy, (2.5)

where G(x) is the free space Green’s function of the problem.

Proof. See Evans (2010) or Sparr & Sparr (2000).

The two and three dimensional free space Green’s functions of the Poisson equation and the Helmholtz equation are listed in table (2.1), where H₀⁽¹⁾ denotes the zeroth order Hankel function of the first kind, r = |r|.

Differential operator 2D 3D

∆ G₀(r) = − log(r)

2π G₀(r) = 1 4πr

∆ + k² G_k(r) = i

4H₀⁽¹⁾(kr) G_k(r) = e^ikr 4πr

Table 2.1: Free space Green’s functions for the Poisson (Laplace) operator and the Helmholtz operator.

For derivations of the free space Green’s functions of the Laplace operator, see Sparr & Sparr (2000).

Remark. The method of Green’s functions (2.5) can also be applied to the linear elliptic differential operators ∆² and ∆(∆ + k²).

2.1.3 Fourier analysis

Definition 2.4. Assume that f ∈ L¹(R). Then the continuous Fourier transform of f is defined as

f (s) =b Z

R

f (x)e^−isxdx , s ∈ R. (2.6) In this report, the Fourier transform of f is denoted bf or F [f].

Remark. The Fourier transform (2.6) can be extended to higher dimension.

In three dimensions, assuming f (r) is a radially symmetric function, the Fourier transform is also a radially symmetric function. Using s = |s| and r = |x|, it can be shown that

f (s) =b









 4π

Z ∞ 0

sin(sr)

sr f (r)r²dr in R³, 2π

Z ∞ 0

J₀(sr)f (r)rdr in R².

(2.7)

Lemma 2.2 (Riemann-Lebesgue lemma). If f satisfies definition (2.4), then lim

|s|→∞

f (s) −b → 0. (2.8)

Proof. See Vretblad (2003).

The Riemann-Lebesgue lemma (2.2) states that the Fourier coefficients of f tend to zero as s tends to infinity. The decay of the Fourier coefficients isb dependent of the smoothness of f , which is stated in the following theorem.

Theorem 2.3. Let bf (s) denote the Fourier transform as in definition (2.4).

If f ∈ C^p−1 for some p ≥ 1, then

f (s) = O(|s|b ^−p−1) as |s| −→ ∞. (2.9) Proof. See Vretblad (2003).

If f is a smooth function, i.e. f ∈ C^∞, then the Fourier coefficients decay at the highest possible rate as p −→ ∞ in (2.9).

Theorem 2.4 (Inversion theorem). Suppose that f ∈ L¹(R) and that f is continuous except for a finite number of finite jumps in any finite interval, and that f (x) = ¹₂ lim_a→x⁺f (a) + lim_a→x⁻f (a)
∀x. Then,

f (x₀) = lim

A→∞

1 2π

Z A

−A

f (s)eb ^isx0ds (2.10) for every x0 where f has left and right derivatives.

Proof. See Vretblad (2003).

If (2.10) holds ∀x₀ ∈ R and bf ∈ L¹(R), then it can be written as f (x) = 1

2π Z ∞

−∞

f (s)eb ^isxds. (2.11) Let F⁻¹[ bf ] denote the inverse Fourier transform of bf . From (2.11), f (x) = F⁻¹[ bf ].

Definition 2.5. Let f be a L-periodic function that is Riemann integrable in [−^L₂;^L₂]. The Fourier series representation of f is defined as

X

m∈Z

c_me^i2πmx/L, (2.12)

where c_m are the Fourier coefficients of f , defined as c_m = 1

L Z L/2

−L/2

f (x)e^−i2πmx/Ldx, m ∈ Z. (2.13) The Fourier series in (2.12) will converge uniformly if the Fourier coeffi- cients tend to zero. Thus, a faster decaying sequence {cm}_m∈Z, obtained by a higher order of differentiability class p, yields a faster convergence rate.

2.1.4 Convolution

Definition 2.6. Let f ∈ L¹(R) and g ∈ L¹(R). The convolution f ∗ g is defined as

(f ∗ g)(x) :=

Z

R

f (x − y)g(y)dy = Z

R

f (y)g(x − y)dy. (2.14) Theorem 2.5 (Convolution theorem). LetF [f] and F [g] denote the Fourier transform of f and g according to (2.6). Then,

F [f ∗ g] = F [f]F [g]. (2.15)

Proof. Using definition (2.4) and (2.6), one obtains F [f ∗ g](s) =Z

R

e^−isx

 Z

R

f (x − y)g(y)dy



dx = (2.16)

= Z Z

R²

e^{−is(x−y+y)}f (x − y)g(y)dxdy = (2.17)

= Z

R

e^−isyg(y)dy Z

R

e^−is(x−y)f (x − y)dx = (2.18)

= Z

R

e^−isyg(y)dy Z

R

e^−isxf (x)dx =F [f](s)F [g](s). (2.19)

The Convolution theorem states that a convolution in a real domain cor- responds to a multiplication in the Fourier domain. Thus, this is a useful tool when evaluating the convolution in equation (2.5). Also, using the Inversion theorem yields

u(x) = Z

Rⁿ

G(x − y)f (y)dy = (G ∗ f )(x) = (2.20)

=F⁻¹[F [G]F [f]] = F⁻¹[ bG · bf ]. (2.21)

2.2 Numerical background

In order to numerically solve partial differential equations and complete a convergence analysis, a computer program is needed, e.g. matlab. matlab contains functions that can be used as tools for solving differential equations.

This chapter presents the relevant matlab implemented functions, along with numerical tools and convergence theory. First, the discrete Fourier transform is stated, along the matlab function Fast Fourier Transform.

Then, convergence theory is presented, including the definition of spectral accuracy.

2.2.1 Discretization

In numerical evaluation, it is needed to use discrete data instead of continu- ous data. Discretization of data can be made in multiple ways and results in a finite number of discrete values, e.g. linearly equispaced data. The process results in a discretization error, which depends on the number of discretiza- tion points N and the distribution of N . In general, a finer linearly spaced mesh will result in a smaller discretization error, in exact arithmetics.

2.2.2 Discrete Fourier transform

Definition 2.7. For N discretization points, assume that f_m = f (m/N ) and {f_m}^{N −1}_m=0 is a finite sequence on R. Then, the discrete Fourier transform of f = {f_m}^{N −1}_m=0 is defined as

fb_ξ =

N −1

X

m=0

f_me^{−i2πξLm/N}, ξ = −N/2, . . . , N/2 − 1. (2.22)

2.2.3 Fast Fourier transform

Fast Fourier transform, fft, is a Matlab implemented algorithm based on the Fastest Fourier Transform in the West. The fft algoritm is used to quickly and accurately compute the discrete Fourier transform (2.22), resulting in O(N log(N )) operations for N discretization points (Trefethen 2000).

2.2.4 Sampling

Sampling, or zero padding, is a process of converting a continuous signal into a discrete signal. The Nyquist-Shannon sampling theorem (Trefethen 2000) establishes a condition on a sample rate, or sampling factor, that is sufficient for reconstructing the original continuous signal, i.e. information loss in the sampling process will be insignificant. The Nyquist-Shannon sampling theorem is used for functions with a Fourier transform that is zero outside a specific frequency region, in order to achieve more accurate amplitude of the sampled signal. This is made by adding zeros to the end of a discretized signal in order to increase its length. Furthermore, if data is discretized using N points in each dimension, af Klinteberg et al. (2017) has shown that a sampling factor of 2.5N in two dimensions and 3N in three dimensions is sufficient for achieving adequate results in the setting used in this thesis.

2.2.5 Convergence theory

In order to complete a convergence analysis, useful definitions are stated below.

Definition 2.8. The `²-norm of a vector x with n elements is defined as

|x| = v u u t

n

X

i=1

|x_i|². (2.23)

Definition 2.9. The maximum absolute row sum of an n × m matrix X is defined as

|X|∞ = max

1≤i≤n m

X

j=1

|x_ij|. (2.24)

2.2.6 Spectral accuracy

Spectral methods are well-used for numerically solving partial differential equations, particularly if the data of the problem are smooth and if one

wants to achieve high accuracy (Trefethen 2000). The idea of spectral meth- ods is to write the numerical solution as a sum of basis functions and identify sum coefficients in order to satisfy the partial differential equation. The Fast Fourier Transform (2.2.3) can be used in such a method, whereas the ba- sis are exponential functions due to the Fourier series representation (2.12).

According to Trefethen (2000), spectral methods can often achieve ten dig- its of accuracy and goes to zero faster than any polynomial, as the typical convergence rate for N discretization points is

O(N^−m) ∀m ≥ 0, for smooth functions. (2.25) This convergence behavior is called spectral. For derivation of convergence rate, see Trefethen (2000).

Chapter 3 Method

In this chapter, the method used to solve linear elliptic partial differential equations in free space, for compactly supported, smooth inhomogeneities, presented by Vico et al. (2016) is described.

The method is based on computing the solution of a PDE as a convolution of the free space Green’s function G and the inhomogeneity f . By utilizing the fact that f is zero outside of a certain domain, G is rewritten to a truncated Green’s function, yielding equivalent solution.

Furthermore, both the theoretical and the numerical method are outlined.

In particular, focus will be on solving the Helmholtz equation and the Poisson equation in Rⁿ for n = 2, 3.

3.1 Theoretical method

Consider a linear elliptic partial differential equation in Rⁿ, n = 2, 3, e.g. the Poisson equation or the Helmholtz equation. According to (2.5), the solution to the PDE can be written as

u(x) = Z

Rⁿ

G(x − y)f (y)dy = (G ∗ f )(x), (3.1) where G(r) is the free space Green’s function of the problem in Rⁿ and f is a smooth, compactly supported inhomogeneity in Ω ⊂ Rⁿ, i.e. f is zero outside of Ω. As seen in the (2.20), by using the Convolution theorem and the Inversion theorem, the solution can be written as

u(x) =F⁻¹[ bG · bf ], (3.2) where F denotes the Fourier transform. Furthermore, since f has compact support, one can rewrite the free space Green’s function G as a truncated Green’s function G^L, whose support contains the support of f as seen in figure (3.1). Define the truncated Green’s function G(r) as

G^L(r) = G(r) · rect r 2L



, (3.3)

Figure 3.1: Illustration of f, G^L, G and respective domain.

where

rect (x) =

(1 |x| ≤ 1/2

0 |x| > 1/2 (3.4)

and L ≥ d, where d is the maximum distance between a source point and a target point in Ω. Due to the compact support of f , the rewriting of G to G^L yields equivalent solution, i.e.

u(x) = Z

Rⁿ

G(x − y)f (y)dy = Z

D

G^L(x − y)f (y)dy. (3.5)

Example 3.1. Consider the inhomogeneous Helmholtz’ equation¹ in R³,

∆u(x) + k²u(x) = −f (x). (3.6) Using G(r) according to table (2.1), the truncated Green’s function for Helmholtz’ equation in three dimensions is defined as

G^L_k(r) = e^ikr 4πrrect

 r 2L



. (3.7)

Furthermore, it is essential to Fourier transform G^L_k, as the solution (3.2) is expressed in terms of the Fourier transforms. As G^L_k is radially symmetric,

1Recall that the Poisson equation is obtained by using k = 0 in Helmholtz’ equation.

Thus, the solution method for the Poisson equation is analogous, only redefining the Green’s function.

using (2.7) results in an evaluation of the Fourier transform as

Gb^L_k(s) =F G^Lk(r)
(s) = 4π Z L

0

sin(sr) sr

e^ikr

4πrr²dr = (3.8)

= Z L

0

sin(sr)

s e^ikrdr = (3.9)

= e^ikr ik

sin(sr) s

L 0

− Z L

0

e^ikr

ik cos(sr)dr = (3.10)

= e^ikL

iks sin(Ls) + e^ikr

k² cos(sr)

L 0

+ Z L

0

e^ikr

k² s sin(sr)dr. (3.11) Evaluating the last integral term results in

Z L 0

sin(sr)

s e^ikrdr = e^ikL

iks sin(Ls) + e^ikL

k² cos(Ls) − 1

k²+ (3.12) +

Z L 0

e^ikr

k² s sin(sr)dr. (3.13)

Rearranging and grouping up terms yield

 1 − s²

k²

 Z L 0

sin(sr)

s e^ikrdr = e^ikL cos(Ls) − i^k_ssin(Ls)
− 1

k² , (3.14)

and thus Z L

0

sin(sr)

s e^ikrdr = e^ikL cos(Ls) − i^k_ssin(Ls)
− 1 k²

 k²− s² k²

⁻¹

. (3.15) Thus, the resulting Fourier transformed truncated Green’s function is

Gb^L_k(s) := −1 + e^ikL cos(Ls) − i^k_ssin(Ls)

k²− s² . (3.16)

Furthermore, there are three values of s for which the function (3.16) becomes singular. When s → 0, the limit is

lims→0Gb^L_k(s) = −1 + e^ikL(1 − ikL)

k² , (3.17)

which is achieved by using the standard limit lims→0

sin(Ls)

s = L. (3.18)

For s → k, using L’Hôpital’s rule yields lim

s→k

Gb^L_k(s) = lim

s→k d

ds −1 + e^ikL cos(Ls) − i^k_ssin(Ls)

d

ds(k²− s²) = (3.19)

= lim

s→k

e^ikL((Ls²− ik) sin(Ls) + ikLs cos(Ls))

2s³ = (3.20)

= e^ikL

2k³[Lk²(sin(Lk) + i cos(Lk))

| {z }

=i exp(−iLk)

−ik sin(Lk)] = (3.21)

= 1

2k³ iLk²− ik sin(Lk)e^ikL
= (3.22)

= 2iLk + 1 − e^2ikL

4k² . (3.23)

The limit s → −k is derived analogously (s → −k in (3.20)). Hence, lim

s→−k

Gb^L_k(s) = 2iLk + 1 − e^2ikL

4k² . (3.24)

The Fourier transformed truncated Green’s functions used in this thesis are listed in tables (3.1- 3.2). The derivation for Helmholtz’ equation in R² and Poisson’s equation in R², R³ is analogous, see Appendix (A).

Let f be a smooth function with compact support in Ω ⊂ R³. After Fourier transforming f , the solution is calculated by using (3.2) and (3.5) as u(x) =F⁻¹[ bG^L_k(s) · bf (s)]. (3.25) In summary, the idea of this method is simply to solve the PDE in free space and utilizing the fact that f has compact support. Thus, it is pos- sible to rewrite the free space Green’s function G to a truncated Green’s function G^L, yielding equivalent solution. The spectral representation of the Green’s function of the Helmholtz operator in three dimensions is presented in figure (3.2). By comparing with the free space kernel, it can be seen that smoothness is obtained at the cost of introducing oscillations.

Differential operator Spectral representation in R²

∆ Gb^L₀(s) = 1 − J₀(Ls)

s² − LJ₁(Ls) log(L) s

∆ + k² Gb^L_k(s) =

1 + iπ

2LsJ1(Ls)H₀⁽¹⁾(Lk) − iπ

2LkJ0(Ls)H₁⁽¹⁾(Lk) s² − k²

Table 3.1: Spectral representation of truncated Green’s functions for the Laplace operator and the Helmholtz operator in R².

Differential operator Spectral representation in R³

∆ Gb^L₀(s) = 2 sin(Ls/2)

s

!2

∆ + k² Gb^L_k(s) = − 1 + e^ikL cos(Ls) − i^k_ssin(Ls)
k²− s²

Table 3.2: Spectral representation of truncated Green’s functions for the Laplace operator and the Helmholtz operator in R³.

Figure 3.2: Spectral representation of free space Helmholtz kernel and truncated Helmholtz kernel in R³ for L = 1.8 and k = 4.

3.2 Numerical method

After following the theoretical approach, deriving the spectral truncated Green’s function and calculating the limit values at singular points, there is a need of implementing a numerical solver. For simplicity, assume that the inhomogeneity f in (3.25) has compact support in a square domain Ω = [−R/2, R/2]ⁿ ⊂ Rⁿ, n = 2, 3. Then, the maximum distance between a source point and a target point in Ω is d = R√

n. Thus, L ≥ R√ n.

Introduce a uniform computational grid x_j =



−R

2, . . . ,R 2 − 1

N



, j = 1, 2, . . . , N, (3.26) where N is even and denotes the number of grid points in each side, that permits f to be discretized on Ω. Furthermore, as seen in figure (3.2), the truncated kernel introduce a slight oscillatory behavior. Thus, in order to sufficiently sample the truncated Green’s function bG^L, there is a need of a

sampling factor. In this thesis, the sufficient sampling factor according to (2.2.4) is

sf :=

(2.5 for n = 2,

3 for n = 3. (3.27)

In order to compute bG^L, introduce N_G = sf ·N and a uniform computational grid

ξ_i = 2π R · sf



−N_G

2 , . . . ,N_G 2 − 1



, i = 1, 2, . . . , N_G. (3.28)

After computing bG^L, replace the singular points in bG^L with the analytically calculated limit values.

As bG^L contains more points than f , f must be zero-padded in order to have the same size N_G as bG^L before computing the discrete Fourier trans- form. Then, compute the discrete Fourier transform of f by a fast Fourier transform. Finally, compute the solution u as in (3.25), i.e.

u =F⁻¹[ bG^L· bf ], (3.29) where the inverse Fourier transform is evaluated by an inverse fast Fourier transform.

Remark. The choice of L is dependent of the inhomogeneity and chosen such that f has numerical compact support in Ω, i.e. f is approximated as zero outside of Ω. For example, if Ω is the unit box, then R = 1 and L ≥ √

3.

Chapter 4 Results

In this chapter, the results gathered from the numerical method described in section (3.2) are presented in a convergence analysis. The inhomogenities are chosen such that an analytic solution is easily retrievable, in order to compare it with the numerical solution. Furthermore, the inhomogenities do not have analytic compact support. Though, choosing L as described in section (3.2), f becomes sufficiently small to be approximated as zero outside of some bounded domain and thus, f achieves numerical compact support.

4.1 Convergence analysis

Consider the following inhomogenities for which the analytic solution is known

f =









 4



r²+ k² 4 − 1



e^−r2 in R² 4



r²+ k² 4 − 3

2



e^−r2 in R³,

(4.1)

where r² = |r|² and k is the wave number in the Helmholtz equation¹. The absolute `²-norm error and the maximum absolute row sum error are presented in figure (4.1) and figure (4.2). As seen in figures (4.1-4.2), the errors rapidly tend to machine precision,  ≈ 10⁻¹⁶ in Matlab, for both the Helmholtz equation and the Poisson equation in two and three dimensions.

Also, the method achieves more than ten digits of accuracy and converges as O(N^−m) for every m ≥ 0. According to section (2.2.6), this behavior is known as spectral convergence.

Furthermore, the decay of discrete Fourier coefficients is presented in figure (4.3) for the two dimensional case, and figure (4.4) for the three di- mensional case. As seen in figures (4.3-4.4), the Fourier coefficients rapidly decay as |ξ| → ∞. According to section (2.1.3), the Fourier series converges uniformly as the Fourier coefficients tend to zero. This confirms the theorem (2.3), since f is a smooth function, i.e. p = ∞, the Fourier coefficients decay exponentially.

1Recall that k = 0 yields the Poisson equation.

Figure 4.1: Absolute `²-norm error for different N , f as in (4.1). For Helmholtz’ equation, k = 4. The line 1/N is used for reference.

Figure 4.2: Maximum absolute error for different N , f as in (4.1). For Helmholtz’ equation, k = 4. The line 1/N is used for reference.

Figure 4.3: Decay of Fourier coefficients in R², f as in (4.1), for N = 256.

Figure 4.4: Decay of Fourier coefficients in R³, f as in (4.1), for N = 128.

Chapter 5 Conclusion

Partial differential equations are often impossible to solve analytically. Thus, it is required to solve them numerically, although achieving highly accurate solutions fast enough can be complicated.

In this thesis, the Poisson equation and the Helmholtz equation were solved in two and three dimensions for smooth and compactly supported in- homogeneities, using the method presented by Vico et al. (2016). A linear elliptic partial differential equation can be solved using a convolution of the Green’s function G and the inhomogeneity f . The fundamental idea of the method is to rewrite the free space Green’s function G of the differential operator, to a truncated Green’s function G^L whose support contains the support of f . Thereafter, computing the solution by utilizing the Convolu- tion theorem and the Inversion theorem.

As seen in the results, the method yields spectrally accurate solutions for linear elliptic PDEs in two and three dimensions. By truncating the Green’s function, it is noted that its spectral representation becomes smooth at the cost of introducing oscillations. The kernel becomes smooth as the limit val- ues at singular points in bG^L is well defined for the different operators. Thus, this method solves the problem of singularities.

The method is easy to implement as there is no need to artificially mimic infinity as e.g. perfectly matched layers. Also, it requires few lines of code and a low computational cost. This is due to only computing relevant points as well as that since the truncated kernel is smooth, few discretization grid points N are needed in each dimension, though still yielding spectral con- vergence. Furthermore, the method is beneficial for solving different linear elliptic PDEs, as only the truncated Green’s function needs to be redefined.

Thus, it is suitable for other linear elliptic differential operators, e.g. ∆² and

∆(∆ + k²), as seen in Vico et al. (2016).

In conclusion, this is a powerful method for solving linear elliptic PDEs in free space with spectral convergence, for smooth and compactly supported inhomogenities.

References

af Klinteberg, L., Shamshirgar, D. S. & Tornberg, A.-K. (2017), ‘Fast ewald summation for free-space stokes potentials’, Research in the Mathematical Sciences 4(1), 1.

URL: https: // doi. org/ 10. 1186/ s40687-016-0092-7

DLMF (n.d.), ‘NIST Digital Library of Mathematical Functions’, http://dlmf.nist.gov/, Release 1.0.22 of 2019-03-15. F. W. J. Olver, A. B.

Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, eds.

URL: http: // dlmf. nist. gov/

Evans, L. C. (2010), Partial Differential Equations, 2nd edn, American Math- ematical Society, Providence, R.I.

Sparr, G. & Sparr, A. (2000), Kontinuerliga System, 2nd edn, Studentlitter- atur AB, Lund.

Trefethen, L. N. (2000), Spectral Methods in MATLAB, SIAM, Philadelphia.

Vico, F., Greengard, L. & Ferrando, M. (2016), ‘Fast convolution with free- space green’s functions’, J. Comput. Phys. 323(C), 191–203.

URL: https: // doi. org/ 10. 1016/ j. jcp. 2016. 07. 028

Vretblad, A. (2003), Fourier Analysis and Its Applications, Springer-Verlag New York, Inc., New York.

Appendix A

Truncated spectral representation of Green’s functions

In this appendix, the spectral representation of the truncated Green’s func- tions and corresponding limit values at singular points are derived, along a comparison to the free space Green’s functions. Throughout this appendix, r = |r| and s = |s| will be used, as well as the Fourier transform (2.7) for radially symmetric functions.

Laplace 3D

The truncated Green’s function for the Laplace operator in 3D is G^L₀(r) = 1

4πrrect r 2L



. (A.1)

A Fourier transform yields

Gb^L₀(s) :=F G^L₀(r)
(s) = 4π Z L

0

sin(sr) sr

1

4πrr²dr = Z L

0

sin(sr)

s dr. (A.2) Evaluating the integral, using the half-angle formulae, results in

Z L 0

sin(sr)

s dr = −cos(Ls) s² + 1

s² = 2 sin(Ls/2) s

2

. (A.3)

It is easy to verify that

lims→0 2 sin(Ls/2) s

2!

= L²

2 , (A.4)

and therefore

Gb^L₀(s) = 2 sin(Ls/2) s

2

. (A.5)

The spectral representation of the Laplace kernel in R³ is presented in figure (A.1).

Figure A.1: Spectral representation of free space Laplace kernel and truncated Laplace kernel in R³ for L = 1.8.

Laplace 2D

The truncated Green’s function for the Laplace operator in 2D is G^L₀(r) = −log(r)

2π rect r 2L



. (A.6)

A Fourier transform results in Gb^L₀(s) = −

Z L 0

J₀(sr) log(r)rdr = (A.7)

= − rJ₁(sr) log(r) s

L 0

+ Z L

0

J₁(sr)

s dr = (A.8)

= −LJ₁(Ls) log(L)

s + lim

r→0⁺

 rJ₁(sr) log(r) s

 +



−J₀(sr) s²

L 0

, (A.9) where the following relations for Bessel functions of the first kind

Z

zJ₀(z)dz = zJ₁(z), (A.10)

Z

J₁(z)dz = −J₀(z), (A.11)

were used (see (DLMF n.d., Eq. 10.22.1)). Furthermore, the middle term that contains the limit as r → 0⁺ can be shown to vanish using L’Hôpital’s rule. Hence, the resulting spectral representation of the Green’s function is¹

Gb^L₀(s) = 1 − J₀(Ls)

s² − LJ₁(Ls) log(L)

s . (A.12)

In the limit s → 0, two applications of L’Hôpital’s rule yields lims→0Gb^L₀(s) = lim

s→0 d

ds(¹₂L²s(−J₀(Ls) + J₂(Ls)) log(L) + J₁(Ls)(L − L log(L)))

d ds(2s)

(A.13)

= lim

s→0

1

4(L²log(L)(J₂(Ls) − J₀(Ls)) + L²s(log(L)(LJ₁(Ls)+

(A.14) + 1

2L(J1(Ls) − J3(Ls))) + L(L − L log(L)(J0(Ls) − J2(Ls))) = (A.15)

= 1

4L²(1 − 2 log(L)) , (A.16)

where the following recurrence relations for Bessel functions were used (see (DLMF n.d., Eq. 10.6.1))

J_ν⁰(z) = J_ν−1(z) − (ν/z)J_ν(z), (A.17) J_ν⁰(z) = −J_ν+1(z) + (ν/z)J_ν(z). (A.18) The spectral representation of the Laplace kernel in R² is presented in figure (A.2).

Helmholtz 2D

The truncated Green’s function for the Helmholtz operator in 2D is G^L_k(r) = i

4H₀⁽¹⁾(kr), (A.19)

where H₀⁽¹⁾ denotes the zeroth order Hankel function of the first kind, defined as

H_ν⁽¹⁾(z) := J_ν(z) + iY_ν(z). (A.20)

1Recall that limz→0J0(z) → 1, limz→0Jν(z) → 0 for ν = 1, 2, . . . (DLMF n.d., Eq. 10.7.1).

Figure A.2: Spectral representation of free space Laplace kernel and truncated Laplace kernel in R² for L = 1.8.

A Fourier transform yields

Gb^L_k(s) = iπ 2

Z L 0

J₀(sr)H₀⁽¹⁾(kr)rdr. (A.21) To evaluate this, the following integral relation for Bessel functions (DLMF n.d., Eq. 10.22.4) can be used

Z

zGν(az)Dν(bz)dz = z (aCν+1(az)Dν(bz) − bGν(az)Dν+1(bz))

a²− b² , (A.22)

where Gν and Dν denotes J_ν, Y_ν, Hν⁽¹⁾ and Hν⁽²⁾. Thus,

Gb^L_k(s) = iπ 2

"

r(sJ₁(sr)H₀⁽¹⁾(kr) − kJ₀(sr)H₁⁽¹⁾(kr)) s² − k²

#L 0

= (A.23)

= iπ 2

L(sJ₁(Ls)H₀⁽¹⁾(Lk) − kJ₀(Ls)H₁⁽¹⁾(Lk))

s²− k² − (A.24)

− iπ 2 lim

r→0⁺

"

r(sJ₁(sr)H₀⁽¹⁾(kr) − kJ₀(sr)H₁⁽¹⁾(kr)) s²− k²

#

. (A.25)

Now, the last term in brackets can be calculated using definition (A.20) and by inserting the power series for Bessel functions (see (DLMF n.d., Sec. 10.8)).

Then, it can be shown that the only non-zero term is

limr→0

r(sJ₁(sr)H₀⁽¹⁾(kr) − kJ₀(sr)H₁⁽¹⁾(kr))

s²− k² = lim

r→0

−rkJ₀(sr)iY₁(kr) s²− k² =

(A.26)

= 2i π

1

s²− k². (A.27) Thus, insertion of (A.27) into (A.25) yields

Gb^L_k(s) =

1 + iπ

2L(sJ₁(Ls)H₀⁽¹⁾(Lk) − kJ₀(Ls)H₁⁽¹⁾(Lk))

s²− k² . (A.28)

For s → k, using L’Hôpital’s rule and the recurrence relation (A.17) yields

lim

s→k

Gb^L_k(s) = lim

s→k

iπLJ1(Ls)



(1 − L)H₀⁽¹⁾(Lk) + LkH₁⁽¹⁾(Lk)



4s + (A.29)

+ lim

s→k

iπLsJ₀(Ls)H₀⁽¹⁾(Lk)

4s = (A.30)

=

iπLJ₁(Lk)

(1 − L)H₀⁽¹⁾(Lk) + LkH₁⁽¹⁾(Lk)

4k + (A.31)

+ iπLJ₀(Lk)H₀⁽¹⁾(Lk)

4 . (A.32)

Analogously for s → −k, using that J₀(s) is an even function and J₁(s) is an odd function yields

s→−klim Gb^L_k(s) =

iπLJ₁(Lk)

(1 − L)H₀⁽¹⁾(Lk) + LkH₁⁽¹⁾(Lk)

4k + (A.33)

+iπLJ₀(Lk)H₀⁽¹⁾(Lk)

4 . (A.34)

The spectral representation of the Helmholtz kernel in R² is presented in figure (A.3).

Figure A.3: Spectral representation of free space Helmholtz kernel and truncated Helmholtz kernel in R² for L = 1.8 and k = 4.

Appendix B Code listing

matlab code, solver in three dimensions

1 % % A l g o r i t h m for s o l v i n g l i n e a r e l l i p t i c PDE 2 % P a r a m e t e r e x p l a n a t i o n :

3 % N is the n u m b e r of grid p o i n t s in each d i m e n s i o n

4 % op is the diff . operator , 1 ( H e l m h o l t z ) or 2 ( P o i s s o n ) 5 % k is the w a v e n u m b e r

6 % f is the i n h o m o g e n e i t y

7 % L is the max d i s t a n c e b e t w e e n a s o u r c e & t a r g e t p o i n t 8 % sf is the s a m p l i n g f a c t o r

9 % G_L is the s p e c t r a l r e p r e s e n t a t i o n of the t r u n c a t e d Green ’ s f u n c t i o n

10 f u n c t i o n u = a l g o r i t m ( N , op , k )

11 R = 12; % C h o o s e R such that f has n u m e r i c a l c o m p a c t s u p p o r t

12 L = R *sqrt(3) ; % M a x i m u m d i s t a n c e when d o m a i n is a box

13 sf = 3;

14 Ng = sf * N ;

15 T = 2*pi/( sf * R ) ;

16 xv = l i n s p a c e( - R /2 , R /2 , N +1) ; 17 xv (end) = [];

18 [ X , Y , Z ] = m e s h g r i d( xv , xv , xv ) ; 19 sv = T *( - Ng /2: Ng /2 -1) ;

20 [ SX , SY , SZ ] = m e s h g r i d( sv , sv , sv ) ; 21 mid = Ng /2 + 1;

22 fpad = z e r o s( Ng , Ng , Ng ) ;

23 f = -4*( X .^2+ Y .^2+ Z .^2+ k ^ 2 / 4 - 3 / 2 ) .*exp( X .^2+ Y .^2+ Z .^2) ;

24 if op == 1

25 G_L = 2*(sin(( L *sqrt( SX .^2+ SY .^2+ SZ .^2) ) /2) ./sqrt ( SX .^2+ SY .^2+ SZ .^2) ) .^2;

26 G_L (i s n a n( G_L ) ) = L ^ 2 / 2 ;

27 end

28 if op == 2

29 G_L = ( -1 + exp(1 i * L * k ) *(cos( k *sqrt( SX .^2+ SY .^2+

SZ .^2) ) - 1 i * k *sin( L *sqrt( SX .^2+ SY .^2+ SZ .^2) ) ./(sqrt( SX .^2+ SY .^2+ SZ .^2) ) ) ) ./( k ^2 -( SX .^2+ SY .^2+ SZ .^2) ) ; 30 % In this case , d i s c r e t i z a t i o n only hits one

s i n g u l a r p o i n t in s =0

31 G_L (i s n a n( G_L ) ) = ( -1 + exp(1 i * k * L ) *(1 -1 i * k * L ) ) / k

^2;

32 end

33 fpad ( mid - N /2: mid + N /2 -1 , mid - N /2: mid + N /2 -1 , mid - N /2:

mid + N /2 -1) = f ;

34 F = i f f t s h i f t ( fftn (f f t s h i f t( fpad ) , [ Ng Ng Ng ]) ) ; 35 U = ( G_L .* F ) ;

36 upad = real(( i f f t s h i f t ( i f f t n (f f t s h i f t( U ) ) ) ) ) ;

37 u = upad ( mid - N /2: mid + N /2 -1 , mid - N /2: mid + N /2 -1 , mid - N /2: mid + N /2 -1) ; % S o l u t i o n u

38 end

Listing B.1: Solver for three dimensional linear elliptic PDE.