On the accuracy of the Complex-Step-Finite-Difference method

Accepted Manuscript On the accuracy of the Complex-Step-Finite-Difference method Rafael Abreu, Zeming Su, Jochen Kamm, Jinghuai Gao

PII: DOI: Reference:

S0377-0427(18)30131-6 https://doi.org/10.1016/j.cam.2018.03.005 CAM 11552

To appear in:

Journal of Computational and Applied Mathematics

Received date : 23 August 2017 Revised date : 19 January 2018 Please cite this article as: R. Abreu, Z. Su, J. Kamm, J. Gao, On the accuracy of the Complex-Step-Finite-Difference method, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.005 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

On the accuracy of the Complex-Step-Finite-Difference method Rafael Abreu 1

∗,1,2

, Zeming Su3 , Jochen Kamm1 and Jinghuai Gao3

Institut f¨ ur Geophysik, Westf¨ alische Wilhelms-Universit¨ at M¨ unster, Corrensstraße 24, D-48149 M¨ unster, Germany 2

Instituto Andaluz de Geof´ısica, Universidad de Granada, Campus de Cartuja s/n, E-18071 Granada, Spain 3

Xi’an Jiaotong University, Xi’an, Shaanxi, China

Abstract The Complex-Step-Finite-Difference method (CSFDM) is a very simple methodology that can be implemented in well known numerical techniques helping to improve, for instance in the wave propagation problem, time and/or space derivative based wavefields. We clarify differences between the CSFDM and previous implementations of the Complex-Step (CS) derivative approximation in well known numerical techniques. We study dispersion properties for the one-way and two-way wave equations using the Finite-Difference method (FDM), the Pseudospectral method (PSM), the Finite-Element method (FEM) and the CSFDM, under the influence of a plane wave and Ricker source time functions. We show the gain in numerical accuracy offered by the methodology of the CSFDM over the FDM, PSM and FEM. We finally discuss consequences of the CSFDM in future scenarios and propose directions of study in this area. Key words: numerical solution, wave propagation, Finite-Difference method, Pseudospectral method, Finite-Element method, Complex-Step method, Complex-Step-Finite-Difference method.

∗ Corresponding

author: [email protected]

1

1

Introduction

The concept of the derivative can be considered as one of the pillars of modern human society. Since its invention in the times of Isaac Newton and Gottfried Leibniz, a countless number of different theories have been developed, e.g., continuum and quantum mechanics, that have helped to develop modern technology as we know it today. It thus becomes imperative to evaluate derivatives for many different reasons, for instance, the evaluation of a building’s safety, the design of microchips, the calculation of barrier options in financial engineering, the propagation of seismic waves through the entire Earth, inverse problems (data assimilation), and design optimization (simulation parameter choice) and many more. The numerical evaluation of derivatives was firstly done using the Finite-Difference method (FDM). Despite its limitations, the FDM is the simplest, easiest and most intuitive numerical technique. The numerical implementation of the FDM can be problematic since the accuracy of the results will depend on the size of the perturbation chosen. This perturbation (∆x) should be selected to avoid numerical cancellation errors and ensuring suitable results [4]. Different methods have been developed to improve accuracy, reduce human error in the implementation and computational cost, for instance, automatic differentiation (AD) [47] and the Complex-Step method (CSM) [38]. The fundamental differences between the FDM, CSM and AD can be intuitively understood using the Taylor’s series. A Finite-Difference (FD), Complex-Step (CS) and AD approximations for the first-order derivative f 0 (x) can be respectively obtained using the Taylor series expansion of f (x + ∆x), f (x + i∆x) and f (x + ∆x) as follows f (x + ∆x) = f (x) + ∆xf 0 (x) + f (x + i∆x) = f (x) + i∆xf 0 (x) − f (x + ∆x) = f (x) + ∆xf 0 (x),

∞ X ∆xn n ∆xn n ∆x2 00 f (x) + ... + f (x) = f (x), 2 n! n! n=0

∞ X ∆x2 00 (i∆x)n n (i∆x)n n f (x) + ... + f (x) = f (x), 2! n! n! n=0

(1)

where i such that i2 = −1 is the imaginary unit and  such that 2 = 0 with  6= 0 is the dual unit. The FDM uses the first expression in (1) to obtain an approximation for f 0 (x). The CSM uses the imaginary part Im of the second expression in (1) to obtain an approximation for f 0 (x) and the method of AD uses the dual part Du of the last expression in (1) to obtain the exact value of f 0 (x). The method of AD and the CSM are based on the same principle: making a perturbation of the function in a direction different from the real axis. In the case of AD, the perturbation is made using a dual number  (2 = 0 with  6= 0) and in the case of the CSM, the perturbation is made using an imaginary number i (i2 = −1). Despite AD gives exact derivative values, its numerical implementation can be far from obvious [36]. Instead, the CSM is a simpler technique, straightforward to implement in a computer and it has been recently generalized and proposed to solve the wave propagation problem [5]. The Complex-Step-Finite-Difference method (CSFDM) introduced by [5] addresses the problem of finding improved numerical techniques for solving the wave propagation problem while trying to reduce or maintain the computational requirements. The CSFDM represents a simple methodology which uses any other efficient, accurate and stable numerical scheme for the acoustic (or elastic) wave equation, like for instance the FiniteDifference method (FDM), the Spectral-Element method (SEM) and/or the Discontinuous-Galerkin method (DGM) and propagates the imaginary part of an imaginary perturbation of the displacement u in time— Im(ut+i∆t )—or space—Im(utx+i∆x ). This simple approach allows us to obtain velocity and acceleration fields x up to sixth-order of accuracy within the use of three time levels of the temporal discretization. In this study, we aim to explore dispersion errors when propagating the imaginary part of the imaginary perturbation (in time and/or space) of the displacement field u in synthetic 1D models where the analytical solution of the wave propagation is known. To this end, we begin with a general introduction establishing and clarifying previous implementations of the CS derivative approximations and the difference with the work introduced with the CSFDM, bringing to light the main differences presented by the CSFDM. We next analyze differences between analytical and numerical solutions of the first-order and second-order wave equations using the FDM, the PSM, the FEM and the CSFDM (Section 4). We show that CS derivative approximations are always more accurate than conventional FD approximations in presence of numerical dispersion errors. In Section 6, we draw the main conclusions of this study and propose future directions of research. Due to the lack 2

of a reference, in Appendix A we compute the analytical solution for the second-order wave equation under the influence of a Ricker source.

2

The Complex-Step method, the Complex-Step-Finite-Difference method and further implementations

The CSM for computing first-order derivatives was formally introduced by [38] in a very simple way using the Taylor series expansion of an imaginary perturbation of an analytic function f (x + i∆x) and taking the imaginary part Im, leading to the following expression f 0 (x) =


Im(f (x + i∆x)) + O ∆x2 , ∆x

(2)

where i is the imaginary unit (i2 = −1), the operator Im refers to the imaginary part, the symbol O refers to the order of the error made in the approximation and ∆x, with ∆x → 0, refers to the differential step (spacing) of the approximation. To establish future comparisons, we now recall the most basic FD approximations to the first-order derivative given by the following expressions f 0 (x) =


f (x + ∆x) − f (x) + O ∆x ∆x

f 0 (x) =

and


f (x + ∆x) − f (x − ∆x) + O ∆x2 . 2∆x

(3)

Approximation ((3) left) is called forward because only a differential step in the positive direction (+∆x) is taken, while approximation ((3) right) is called centered because it takes differential steps in the positive (+∆x) and negative (−∆x) directions. Two main differences are noticeable between the CS approximation (2) and the FD approximations (3): (i) both FD approximations in (3) show a subtraction in the numerator (hence its name finite difference), while the CS approximation (2) shows a single term in the numerator. This avoids numerical cancellation errors inherent
to all classical FD approximations1 and (ii) the CS approximation (2) is second-order accurate O ∆x2 at the cost of a single imaginary step (i∆x) unlike the second-order FD approximation ((3) right) at the cost of two steps in the positive (+∆x) and negative (−∆x) directions. The CS derivative approximation (2) has proved to be more efficient than conventional FD techniques in many applications (e.g. [7, 6, 43, 25, 26, 10, 18, 21, 8, 46, 35]). Several generalizations to the initial CS derivative approximation (2) have been made (e.g. [23, 9, 2, 4, 16, 15, 24]). A simple generalization made by [4] consisted into taking a complex step (h + iv) in a strict sense in the original CS approximation (2) as follows f 0 (x) =


Im(f (x + h + iv)) + O h, v 2 , v

(4)

where h, v ∈ R and h, v → 0. Although discretization (4) is less accurate compared to the original CS discretization (2), it allows a generalization of the original CS derivative concept to a wider range of possible higher-order accurate approximations using the imaginary part of the function, like for instance (for a full list see [5, 4, 3])  2    Im(f (x + h + iv)) + Im(f (x − h + iv)) 3h − v 2 f 0 (x) = − f 000 (x) + O 5h4 − 10h2 v 2 + v 4 (5) 2v 6 and f 0 (x) =



 i 3h2 − v 2 v h Im(f (x + iv)) + Im(f (x + h + iv)) + Im(f (x − h + iv)) 3h2 v 6h2  2   4 2  2 2
v (5h − 7v ) 7h v − 35h2 v 4 + 18v 6 v − f (x) + f vii (x) + O h8 , v 8 . 5!3 7!3

(6)

1 the numerical cancellation error refers to when the machine precision is reached and the computer can no longer distinguish between two small numbers, so the differences in expression (3) become zero.

3

Note that in approximations ((5)-(6)), we can choose the differential steps h, v and the approximations will remain second- and fourth-order accurate respectively. If we choose a certain relation between h, v like: 3h2 −
2 4 v = 0, the error becomes fourth-order O h in approximation (5) and choosing 5h2 − 7v 2 = 0, the error
6 becomes of sixth-order O h in approximation (6), within compact three levels stencils2 . It is here where the generalization of the CSM made by [4] gains its importance. At a first glance, it becomes very attractive to apply a higher-order accurate discretizations like (4) to solve certain PDE. Although this may look very promising, it was studied in [3] for the wave propagation problem, without significant practical results. However, it was found by [5] that a practical application for approximations like ((5)-(6)) happens to be extremely simple. The Complex-Step-Finite-Difference method (CSFDM) introduced by [5], combines generalized approximations ((5)-(6)) with any well known numerical technique applied to solve the wave propagation, and improves the accuracy of the solution of the wavefields. The use of the original CS approximation (4) in conjunction with other numerical techniques like the FEM is already well known (see [39, 40, 43, 44, 45, 20, 30, 11, 32, 31, 14] among many others). We emphasize that the main difference in all of these previous implementations and the CSFDM is that the CSFDM does not modify the structure of any numerical technique in order to gain more accuracy in the approximation of the partial derivatives. Instead, the CSFDM uses any known efficient numerical technique and modifies the motion condition only (initial conditions in time and/or space). In the case of computational wave propagation, an initial plane wave and/or a source time function located at a single point or at different spatial locations. The main advantage and simplicity of the CSFDM is that it does not involve tedious, complex and time consuming modifications of pre-existing numerical codes. In most cases, these modifications increase the computational requirements of the numerical simulation and also change its dispersion properties. However, despite of this fact, in many papers, the study of the dispersion properties of the numerical technique is completely omitted. In this sense, it is proper to perform an analytical study of the dispersion properties of the modified scheme in order to really state whether the introduction of the CS derivative calculation improves numerical results. It is already well known that the increase in accuracy is not the only factor that makes a numerical technique better, but the reduction of errors due to numerical dispersion plays also a fundamental role [27, 5].

3

Finite-difference and complex-step-finite-difference approximations

As mentioned before, it has been shown by [5] that to obtain higher-order accurate solutions of the velocity and/or acceleration fields in the numerical solution of the seismic wave propagation, we can use any discretization given by the FDM, PSM, FEM, SEM and DGM and propagate the imaginary part of the imaginary perturbation in time and/or space. This simple and novel methodology, called by [5] the CSFDM, is valid from a mathematical and computational point of view but introduces a different physical perspective of the PDE under study. In this section, we illustrate this by studying differences when using the FDM and the CSFDM for two simple PDEs: the 1D first-order and second-order wave equations.

3.1

The first-order wave equation and the FDM and CSFDM

To illustrate differences in the physical interpretations of the FDM and CSFDM, we first consider the 1D one-way wave equation given by the following mathematical expression ∂u ∂u =c , ∂t ∂x

(7)

where u = u(x, t) is the displacement field, x is the space variable, t is the time variable and c the propagation velocity. The most basic FD discretizations to the one-way wave equation (7) are called (i) upwind, resulting from applying the forward approximation ((3) left) in time and space to (7) and (ii) the leapfrog approximation, resulting from applying the centered approximation ((3) right) in time and space to (7). The upwind and the leapfrog approximations are respectively given by the following expressions

ut+∆t = S utx+∆x + 1 − S utx + O ∆t, ∆x , (8) x 2 here we mean three-level stencil by considering the complex differential steps as h + iv, iv and −h + iv. Note that the locations change in the real value h but not in the imaginary value iv.

4

c∆t ∆x



ut+∆t = S utx+∆x − utx−∆x + ut−∆t + O ∆t2 , ∆x2 , x x

(9)

is called the Courant number and/or stability factor. The parameter ∆x is referred as the grid where S = spacing and ∆t as the time step. Expressions ((8)-(9)) can be respectively written, without including the error term, in a general operator form as follows t + t ∆+ t ux = S ∆x ux

and

with ∆+ x f (x) = f (x + ∆x) − f (x)

2t utx = S 2x utx ,



2nx f (x) = f x + n∆x − f x − n∆x .

and

(10) (11)

The CSFDM simply takes any efficient numerical discretization given by any numerical technique and propagates the imaginary part of the imaginary perturbation in time (Im(ut+i∆t )) or space (Im(utx+i∆x )). Following x this, we can write (without including the error terms) two different complex-step-finite-difference (CSFD) versions of discretizations (10), for imaginary perturbations in time and space as follows



+ t+i∆t + t ∆+ =0 and ∆+ (12) t − S ∆x Im ux t − S ∆x Im ux+i∆x = 0,

2t − S 2x Im uxt+i∆t = 0

and



2t − S 2x Im utx+i∆x = 0.

(13)

Although, we are able to compute more accurate velocity (and acceleration) fields using discretizations ((12)(13)) (see [5]), the PDE that they mathematically represent may be different from the initial PDE used in the FD discretization (7). After some algebra it can be shown that eqs. (12) lead to the following PDEs respectively ∂2u ∂2u ∂2u ∂2u = c and = c . (14) ∂t2 ∂x∂t ∂t∂x ∂x2 Equations (14) are different from the initial PDE (7), simply representing the time and spatial derivatives of the original equation (7). It is clear that the PDEs in (14) can be reduced to the original one (7) only in the case of homogeneous media.

3.2

The second-order wave equation and the FDM and CSFDM

We apply the FDM to the 1D second-order wave equation given by the following mathematical expression ∂2u ∂2u = c2 2 , (15) 2 ∂t ∂x where u = u(x, t) is the displacement field, x is the space variable, t is the time variable and c the propagation velocity. We recall the simplest FD approximation for the second-order derivative f 00 (x) given by the following expression
f (x + ∆x) − 2f (x) + f (x − ∆x) f 00 (x) = + O ∆x2 . (16) 2 ∆x Applying the FD approximation (16) in time and space to the second-order wave equation (15), we find the leapfrog FD approximation


ut+∆t = S2 utx+∆x + utx−∆x + 2 1 − S2 utx − ut−∆t + O ∆t2 , ∆x2 . (17) x x Equation (17) can be written, without including the error term, in an operator form as follows
t
t − 2 + − 1/2 2 1/2 1/2 ∆+ or 21/2 ux = 0, t ◦ ∆t − S ∆x ◦ ∆x ux = 0, x ◦ 2x − S 2x ◦ 2x

(18)

where ◦ is the convolution operator and ∆, 2 are given in (11). Depending whether an imaginary perturbation in time (i∆t) or space (i∆x) is used, two different CSFD versions of discretization (18) are obtained



− − 2 + − t+i∆t 2 + − t ∆+ = 0, ∆+ (19) t ◦ ∆t − S ∆x ◦ ∆x Im ux t ◦ ∆t − S ∆x ◦ ∆x Im ux+i∆x = 0.

The previous CSFD schemes in fact respectively solve the following PDEs

∂3u ∂3u ∂3u ∂3u = c2 2 and = c2 3 . (20) 3 2 ∂t ∂x ∂t ∂t ∂x ∂x It is noticeable that the corresponding PDEs (20) are not the same as the initial two-way wave equation (15) and they can be equivalent only in the case of homogeneous media. 5

4

Accuracy of the Complex-Step-Finite-Difference approximations

The main purpose of every numerical technique is to give an approximation to the exact solution of a certain mathematical equation. From a strict point of view, the most appropriate numerical method, for a certain physical configuration, is the one that reproduces the analytical solution within the lowest possible error band and with the smallest requirement of computational memory. In the case of computational wave propagation, numerical errors are related to spatial and temporal dispersion-dissipation properties of the numerical technique (for a detailed study of dispersion-dissipation properties of the FDM refer to [27]). In this study, we simply refer to both dispersion-dissipation properties as numerical dispersion. In practice, there are two main different elements that may induce numerical dispersion: physical and numerical heterogeneities. In the case of computational wave propagation, physical heterogeneities refer to differences in the model material parameters (density, velocities, elastic parameters, etc), which in fact induce physical dispersion in the wave propagation. As users of the numerical technique, we are interested into reproducing (up to a certain level of accuracy) those physical phenomena with the numerical solutions. The second one, called in this study numerical heterogeneities, appear when users do not set the correct (or most appropriate) parameter values previous running the numerical simulation. Setting incorrect values will certainly induce dispersion into the numerical solution without any physical meaning. Errors induced from numerical heterogeneities come by the wrong selection of the grid spacing (∆x) and/or the wrong time step (∆t). These parameter values are connected by the Courant number S (S = c ∆t/∆x). This means that the best parameter values that can be chose, are given by the optimum Courant number S that the numerical technique can handle. In the next sections, we first illustrate dispersion properties of the FDM applied on the one-way and two-way wave equations. Later, we perform different numerical experiments using the CSFDM.

4.1

Dispersion properties of the upwind and leapfrog FD schemes

To gain insight into the wave propagation phenomena using the upwind (8) and leapfrog (17) schemes, we assume a plane wave solution of the form u(x, t) = Aei(kx+ωt) ,

(21)

where k is the wavenumber, ω the frequency and A the amplitude. Substituting the plane wave (21) into the upwind (8) and leapfrog (17) discretizations, leads to the following expressions respectively eiω∆t − 1 = S(eik∆x − 1)

cos ω∆t − 1 = S2 (cos k∆x − 1),

and

(22)

∆t . The first step of the dispersion analysis is to determine the maximum value of the Courant with S = c ∆x number S that the numerical technique can handle. After some algebra we obtain 0 ≤ S ≤ 1 for both schemes. After having obtained the values that the Courant number may adopt, we proceed to analyze the numerical errors for different values of the Courant number. This is physically translated into how good the numerical method performs in presence of spacial heterogeneities (c = c(x)). The reader should note that by varying the value of the velocity field c(x) over the spatial coordinate, it is mathematically translated into varying the value of the Courant number S = c(x)∆t/∆x. At this point we use the concept of the ratio R between the numerical cnum and analytical c phase velocities, given by the following expression

R=

cnum cnum ∆t ω∆t = = . c S∆x S k∆x

(23)

For finding values of the ratio R, we assume the Courant number and k∆x values, and using expressions (22) we obtain ω∆t values needed to compute the ratio R using (23). Figure 1-a shows values of the ratio R (23) as a function of the Courant number S, ω∆t and k∆x for the upwind FD scheme (8). The real part of R—Re(R)—is equal to one for only a single value of the Courant number (S = 1). Having R = 1 means that the numerical and the analytical solutions are the same. The imaginary part of R—Im(R)—is not zero. This is mathematically translated into having imaginary values of 6

a)

Dispersion plot − Real part

b)

Dispersion plot − Imaginary Part

1.03

Courant Courant Courant Courant Courant

0.5

1.02 1.01

0.25 0.5 0.7 0.9 1

Dispersion plot − Real part

Dispersion plot − Imaginary Part

1.005

0.05

1

0.04 0.03

0.995

0.4

Courant Courant Courant Courant Courant

1

0.25 0.5 0.7 0.9 1

0.02

0.98 0.97

0.3

0.2

imag(Ratio Factor)

0.99

real(Ratio Factor)

imag(Ratio Factor)

real(Ratio Factor)

0.99 0.985 0.98

0.01 0 −0.01

0.975 0.96

−0.02 0.1

0.97

0.95 0.94 0.93 0

c)

0.2

0.4

k∆x

0.6

0.8

1

0

0.2

0.4

k∆x

0.6

0.8

0.96 0

1

d)

Dispersion plot − Imaginary Part

Dispersion plot − Real part

−0.03

0.965

0

1.5

Courant Courant Courant Courant Courant

0.5 1.4

0.25 0.5 0.7 0.9 1

−0.04

0.2

0.4

k∆x

0.6

0.8

−0.05 0

1

Dispersion plot − Real part

0.2

0.4

k∆x

0.6

0.8

1

Dispersion plot − Imaginary Part 0.1

Courant Courant Courant Courant Courant

1

0.98

0.25 0.5 0.7 0.9 1

0.4 0.96

0.3

0.2

imag(Ratio Factor)

1.2

real(Ratio Factor)

0.05

imag(Ratio Factor)

real(Ratio Factor)

1.3 0.94

0.92

0.9

1.1 0.1

0

0.88

1 0.86 0 0.9 0

0.2

0.4

k∆x

0.6

0.8

1

0

Dispersion plot − Real part

e)

0.2

0.4

k∆x

0.6

0.8

0.84 0

1

Dispersion plot − Imaginary Part 0.6

Courant Courant Courant Courant Courant

1 0.5

0.98 0.96

0.4

k∆x

0.6

0.8

−0.05 0

1

Dispersion plot − Real part

f)

0.25 0.5 0.7 0.9 1

0.2

0.2

0.4

k∆x

0.6

0.8

1

Dispersion plot − Imaginary Part

1.06

0.05

Courant Courant Courant Courant Courant

0.04 1.05 0.03

0.25 0.5 0.7 0.9 1

0.4 1.04

0.02

0.9 0.88

0.3

0.2

0.1

0.86

imag(Ratio Factor)

0.92

real(Ratio Factor)

imag(Ratio Factor)

real(Ratio Factor)

0.94 1.03

1.02

0.01 0 −0.01 −0.02

1.01 0

−0.03

0.84

0.8 0

1

−0.1

0.82

0.2

0.4

k∆x

0.6

0.8

1

−0.2 0

−0.04

0.2

0.4

k∆x

0.6

0.8

0.99 0

1

0.2

0.4

k∆x

0.6

0.8

1

−0.05 0

0.2

0.4

k∆x

0.6

0.8

1

Figure 1: Dispersion plots for the (a) upwind (8), (b) leapfrog (17), (c) Lax-Friedrichs (35), (d) Lax-Wendroff (37), (e) FEM (39) and (f) PSM (41) discretizations. On the left the real part and on the right the imaginary part of the ratio R (23).

the frequency ω and physically translated into numerical dissipation for different values of the Courant number not equal to one S 6= 1. In general, we can conclude that the upwind scheme reproduces the analytical solution (Re(R) = 1 and Im(R) = 0) for S = 1 only, and for smaller values of the Courant number (S < 1) the upwind scheme presents significant numerical dispersion. Figure 1-b shows values of the ratio R (23) as a function of the Courant number S, ω∆t and k∆x for the leapfrog scheme (17). Like in the upwind case, the real part of R—Re(R)—is equal to one for only a single value of the Courant number (S = 1) and unlike the upwind case, the imaginary part of R—Im(R)—is zero. This is physically translated into having numerical dispersion but not dissipation for different values of the Courant number not equal to one (S 6= 1). This means that the leapfrog scheme (17) reproduces the analytical solution for only S = 1 and for smaller values of the Courant number (S < 1), the numerical errors due to dispersion become more significant.

7

4.2

Accuracy of Complex-Step-Finite-Difference approximations for the first-order wave equation

To analyze the accuracy of the CSFDM applied to the one-way wave equation, we consider the upwind approximation (8). We have seen that the stability condition of the upwind approximation is given by satisfying the following condition for the Courant number: 0 ≤ S ≤ 1 and, for the value of S = 1, the upwind approximation (8) gives the exact solution. We analyze dispersion errors for the one-way wave equation using the FDM and CSFDM, considering the following initial value problem ∂u ∂u =c , ∂t ∂x

with

u(x, 0) = sin πx

and

0 ≤ x ≤ 3π.

(24)

The analytical displacement and velocity solutions for the initial value problem (24) are given by the following expressions respectively u(x, t) = sin (πx + πct),

u(x, ˙ t) = πc cos (πx + πct),

with

0 ≤ x ≤ 3π

and t > 0.

(25)

As pointed out before, the CSFDM simply takes the numerical scheme (upwind discretization (8) in this case) and propagates the imaginary part of the imaginary perturbation (in space in this case since we have an initial condition over the space instead of a time dependent function), that is, u(x, 0) = Im(sin (π(x + i∆x)))

with

0 ≤ x ≤ 3π.

The upwind FD and CSFD discretizations may be then respectively written as follows

+ + ∆+ and ∆+ t − S ∆x sin πx = 0 t − S ∆x Im(sin (π(x + i∆x))) = 0,

(26)

(27)

where, we again emphasize that the only difference is in the value to be propagated. At this point we can pose the following argument: if with the FDM we are able to propagate the initial displacement condition (24) up to the time level of interest and, with those displacement values we are able to compute velocity and acceleration fields (up to second-order accuracy within three time levels) and, analogously the CSFDM tells that we are able to propagate the imaginary part of the imaginary perturbation (in time or space) to get higher-order accurate velocity or acceleration fields. Then, we should be able to propagate the value of our preference, despite the physical interpretation of the numerical discretization. Following the previous conjecture, we can generalize the upwind numerical scheme (8) to the following expressions
t
t + + ∆+ ˙ x = 0, and/or ∆+ (28) t − S ∆x u t − S ∆x ux,x = 0,

where the u˙ tx stands for the displacement derivative with respect to time t and utx,x for the derivative with respect to space x. We can write the corresponding velocity/gradient propagation schemes (28) using the initial motion condition given in (24), which leads to the following expression
+ ∆+ (29) t − S ∆x π cos πx = 0.

The natural question that follows our conjecture is as follows: which is the most appropriate discretization when we are interested in computing velocity and/or acceleration fields? To answer this, we have to take into account dispersion errors intrinsic to the numerical discretization. To this end, we consider the case when we are in the presence of numerical heterogeneities, that is, the numerical artifacts related to errors in the selection of the simulation parameters by the user of the numerical technique. We recall, numerical heterogeneities as defined in this study. We emphasize that considering numerical heterogeneities is synonymous to considering heterogeneous models (c = c(x)) in the simulation, which is simply translated into a variable Courant number S. To induce numerical heterogeneities, we select the wrong (or not optimum) time step (∆t) at the beginning of the numerical simulation, which again, is directly reflected in the value of the Courant number S. To explore 8

differences in the FD and CSFD numerical schemes, we compute errors in amplitudes between the analytical and numerical waveforms. Errors in amplitudes EA are computed using the following expression X Aanalytical − Anumerical i i EA =

i

X Aanalytical i

× 100%,

(30)

i

and Anumerical are the analytical and numerical wave amplitudes respectively evaluated at the where Aanalytical i i location i. Figure 2-a shows the obtained numerical results. A complex behavior of the error function can be observed: on one hand, when the Courant number S is close to the optimum value (S ≈ 1), all approximations behave as expected; we get the characteristic solution propagating the velocity value in the initial motion condition (eq. (29), see Fig. 2-a left). On the other hand, for smaller values of the Courant number (S → 0), all CSFD approximations give smaller errors (Fig. 2-a right) and also, the solution with the lowest error amplitudes is the fourth-order CSFD approximation (5). The FD and CSFD approximations used are given by eqs. (3)-(6) and assuming the most appropriate relationship between h, v.

4.3

Accuracy of Complex-Step-Finite-Difference approximations for the secondorder wave equation

We explore the accuracy of CSFD approximations for the second-order wave equation (15). We use a Ricker function as time dependent source given by the following expression 2

2

f (xs , t) = δ(x − xs ) (1 − 2f02 (t − t0 )2 ) e−f0 (t−t0 ) ,

(31)

where δ is the Dirac delta function, f0 is the dominant frequency, t0 the time delay and xs the location of the source. The analytical solution is computed in Appendix A. For the numerical implementation of the Dirac delta function δ, we use δ(x − xs ) =

1 + O(∆x2 ) ∆x

at

x = xs ,

(32)

which is a second-order approximation [37]. This means that the numerical solution for the initial value problem (A) cannot be more accurate than second-order due to the numerical implementation of the Dirac delta function in the FD scheme. Figure 2-b shows errors in amplitude for the FDM and the CSFDM using the leapfrog scheme (17). Like in the case of the first-order wave equation studied in the previous section, the approximation that gives smaller errors in presence of numerical heterogeneities (S 6= 1), is the fourth-order CS approximation (5). We emphasize that we compare only higher-order approximations against the analytical solution since we use a second-order approximation for the location of the Dirac delta function δ, for which second-order approximations are evidently the most accurate ones. The FD and CSFD approximations used are given by eqs. (3)-(6) and employing the most appropriate relationship between h, v. As a second illustration, we solve the second-order wave equation (15) subject to the boundary condition π

u(0, t) = −u0 ei(ωt+ 2 ) ,

(33)

where u0 is an amplitude. The analytical solution is given by the following expression π

u(x, t) = −u0 ei(kx+ωt+ 2 ) .

(34)

Figure 2-c shows the obtained results. We observe that like in previous illustrations, the CSFDM gives the most accurate results. Once again, the fourth-order CSFD approximation (5) is more accurate than the sixth-order CSFD (6) for small values of the Courant number.

9

velocity amplitude errors

2.5

velocity amplitude errors

b)

0.236

0.2335

2.3

2

1.5

1

velocity amplitude errors (zoom)

velocity propagation amplitude errors 4−order CSFD amplitude errors 6−order CSFD amplitude errors

0.234

amplitude error percentage

amplitude error percentage

velocity amplitude errors (zoom) 2.4

2−order FD amplitude erros velocity propagation amplitude errors 2−order CSFD amplitude errors 4−order CSFD amplitude errors 6−order CSFD amplitude errors

0.232

2.2

2.1

0.23

amplitude error percentage

3

amplitude error percentage

a)

0.228

0.226

0.224

0.222

0.2335

0.2335

0.2335

2 0.22

0.2335

0.5 0.218

1.9 0 0.2

0.3

c)

0.4

0.5

0.6

0.7

Courant number S

0.8

0.9

1

0.365

velocity amplitude errors 0.35

0.375

0.38

0.385

0.39

0.395

Courant number S

0.4

0.216 0.2

0.405

d)

velocity amplitude errors (zoom)

0.2335 0.3

0.4

0.5

0.6

0.7

0.8

Courant number S

0.9

1

0.22730.22730.22730.22730.22730.22740.22740.22740.2274

Courant number S

velocity amplitude errors 18

2−order FD amplitude erros velocity propagation amplitude errors 2−order CSFD amplitude errors 4−order CSFD amplitude errors 6−order CSFD amplitude errors

0.3

0.37

0.3

16

velocity amplitude errors (zoom)

2−order FD amplitude erros velocity propagation amplitude errors 2−order CSFD amplitude errors 4−order CSFD amplitude errors 6−order CSFD amplitude errors

12.95

12.9

14

0.15

0.25

0.2

amplitude error percentage

0.2

amplitude error percentage

amplitude error percentage

amplitude error percentage

0.25

12 10 8 6

0.1

12.85

12.8

12.75

4 0.05

12.7

2 0.15 0.4

0.5

0.6

0.7

Courant number S

0.8

0.9

1

0.3

0.5

f)

5

0.04

0.03

0.02

0.04

0.035

0.4

0.6

Courant number S

0.8

1

0.2385

0.239

Courant number S

velocity amplitude errors 6

0.045

amplitude error percentage

amplitude error percentage

0.45

0.05

2−order FD amplitude erros velocity propagation amplitude errors 2−order CSFD amplitude errors 4−order CSFD amplitude errors 6−order CSFD amplitude errors

0.05

0.4

Courant number S

velocity amplitude errors (zoom)

velocity amplitude errors 0.06

0.35

0.2395

velocity amplitude errors (zoom)

2−order FD amplitude erros velocity propagation amplitude errors 2−order CSFD amplitude errors 4−order CSFD amplitude errors 6−order CSFD amplitude errors

3.8 3.78 3.76

amplitude error percentage

e)

0.3

0 0.2

amplitude error percentage

0 0.2

4

3

2

3.74 3.72 3.7 3.68 3.66 3.64

0.01

3.62

1

0.03

3.6 0 0.2

0.4

0.6

Courant number S

0.8

1

0.35

velocity amplitude errors

g)

0.07

0.4

0.45

0.5

Courant number S

0 0.2

0.4

0.6

Courant number S

0.8

1

0.2505 0.251 0.2515 0.252 0.2525 0.253 0.2535

Courant number S

velocity amplitude errors (zoom)

velocity propagation amplitude errors 4−order CSFD amplitude errors 6−order CSFD amplitude errors

0.026

0.06

0.026

0.0259

amplitude error percentage

amplitude error percentage

0.05

0.04

0.03

0.0259

0.0259

0.0259

0.0259

0.02 0.0258

0.0258

0.01

0.0258 0 0.2

0.3

0.4

0.5

Courant number S

0.6

0.7

0.3728

0.3729

0.3729

0.373

Courant number S

0.373

Figure 2: Errors in amplitude for the one-way wave equation solved with the (a) upwind (8), (d) Lax-Friedrichs (35), (e) Lax-Wendroff (37) and (f) FEM (39) under an initial condition of motion given by a plane wave. Errors in amplitude for the two-way wave equation solved with the FDM and the CSFDM using the leapfrog scheme (17) under the presence of a Ricker source time dependent function (b) and an harmonic boundary condition (c). Errors in amplitude for the two-way wave equation solved with the (g) PSM (41).

10

4.4

Further numerical approximations for the first-order wave equation

There are several techniques for solving the first-order wave equation. We explore three more different popular schemes: Lax–Friedrichs, Lax-Wendroff and the FEM and the solutions obtained using the CSFDM. The Lax–Friedrichs discretization [41] solves the first-order wave equation (7) in the following way ut+∆t = x

 S  1 t ux+∆x + utx−∆x + utx+∆x − utx−∆x . 2 2

(35)

Inserting the plane wave solution (21) into (35), we can find the dispersion relation for the Lax–Friedrichs discretization given by the following expression eiω∆t = cos k∆x + iS sin k∆x.

(36)

Figure 1-c shows values of the ratio R (23) as a function of the Courant number S, ω∆t and k∆x for the Lax–Friedrichs scheme (35). The real part of R—Re(R)—is equal to one for only a single value of the Courant number (S = 1) and the imaginary part of R—Im(R)—is not zero. This means that the Lax-Friedrichs scheme presents numerical dissipation errors in heterogeneous media. Solving the initial value problem (24) leads to similar behavior compared to the upwind scheme (see Fig. 2-d). The Lax-Wendroff discretization [41] solves the first-order wave equation (7) in the following way ut+∆t = utx + x

 S2   S t ux+∆x − utx−∆x + utx+∆x − 2utx + utx−∆x . 2 2

(37)

Inserting the plane wave solution (21) into (37), we can find the dispersion relation for the Lax–Wendroff discretization given by the following expression
eiω∆t = 1 + iS sin k∆x + S2 cos k∆x − 1 . (38)

Figure 1-d shows values of the ratio R (23) as a function of the Courant number S, ω∆t and k∆x for the Lax–Wendroff scheme (37). The real part of R—Re(R)—is equal to one for only a single value of the Courant number (S = 1) and the imaginary part of R—Im(R)—is not zero. Compared to the upwind and Lax-Friedrichs schemes, the Lax-Wendroff scheme presents considerably smaller numerical dissipation and dispersion errors. This means that the Lax–Wendroff scheme is a more appropriate choice to solve the wave propagation in heterogeneous media, compared to the upwind and Lax-Friedrichs schemes. Solving the initial value problem (24), we observe that the CSFDM shows smaller numerical errors (see Fig. 2-e). We finally evaluate one (out of many) Finite-Element discretization given by the following expression (see [19] page 109)   1 t 2 t 1 t 1 t 1 t t+∆t ux = ux−∆x + ux + ux+∆x + S u − u . (39) 6 3 6 2 x+∆x 2 x−∆x Inserting the plane wave solution (21) into (39), we can find the dispersion relation given by the following expression eiω∆t =

 1 cos k∆x + 2 + iS sin k∆x. 3

(40)

Figure 1-e shows values of the ratio R (23) as a function of the Courant number S, ω∆t and k∆x for the FEM (39). The real part of R—Re(R)—is never equal to one and the imaginary part of R—Im(R)—is not zero. The scheme presents similar dissipation errors compared to the Lax-Friedrichs, despite that the real part of R presents small numerical dispersion errors. Solving the initial value problem (24), we observe that the FEM shows smaller numerical errors not for the highest value of the Courant number but for the value of S = 0.6 (see Fig. 2-f). This is of particular interest when modeling wave propagation in highly heterogeneous media. Like in the previous simulations, the CSFDM gives the most accurate velocity results.

11

4.5

Further numerical approximations for the second-order wave equation

The pseudospectral discretization solves the second-order wave equation (15) in the following way (see [17])
ut+∆t − 2 utx + ut−∆t x x = c2 PS2(utx ) + O ∆t2 ∆t2

(41)

where PS2 refers to the second-order Fourier derivatives (see [13, 17]). Inserting the plane wave solution (21) into (41), we can find the dispersion relation for the PSM given by the following expression cos ω∆t = 1 −

c2 k 2 ∆t2 S2 k 2 ∆x2 =1− . 2 2

(42)

After some algebra, it can be shown that 0 < S ≤ 0.64 for discretization (41) (see [17]). Figure 1-f shows values of the ratio R (23) as a function of the Courant number S, ω∆t and k∆x for the PSM (41). The real part of R—Re(R)—shows considerable small errors and the imaginary part of R—Im(R)—is zero. Solving the initial value problem (24), we observe that the CSFDM shows smaller numerical errors (see Fig. 2-g).

5

Computing the imaginary part of an imaginary perturbation for discrete signals

There is still the open problem of how to apply the CSM when only real and discrete values of the functional f (x) in study are known. For the application of the CSFDM, the imaginary part of imaginary perturbation in time or space for signals have to be known. Since no clear analytic expression is available when only discrete values of a certain function are given, we can resort to the description of the imaginary part of the imaginary perturbation written in terms of the Taylor’s series as follows Im(f (x + i∆x)) = ∆xf 0 (x) −

∞ X (−1)n+1 ∆x2n−1 2n−1 ∆x3 000 f (x) + ... = f (x), 3! (2n − 1)! n=1

(43)

where Im(f (x)) = 0 because x is set to be a real number and f (x) is real on the entire real axis. To obtain the best approximation to the imaginary perturbation of the available discrete data, we may employ an accurate numerical technique to compute a derivative of a function: the Fourier derivative. Despite its limitations to homogeneous spacing ∆x, the Fourier derivative has shown to be very accurate in the seismic wave propagation problem [12, 13]. The derivative of a function f (x) is calculated using the differentiation theorem of the Fourier transform dn f (x) = IF T [(−i k)n F (k)], dx

(44)

where IF T holds for the inverse Fourier transform, k the wavenumber and F (k) is the Fourier transform of f (x). Using the discrete Fourier transform of functions, we can obtain exact derivatives up to the Nyquist wavenumber kN = π/∆x. Since the Fourier series is based on sine and cosine functions, it is implied that f (x) has to be periodically continued outside of the domain of interest. We evaluate a 2π-periodic Gauss function in the interval x ∈ (0, 2π) given by the following expression f (x) = e−1/σ

2

(x−x0 )2

,

(45)

with x0 = π and the analytical derivative f 0 (x) = −2

x − x0 −1/σ2 (x−x0 )2 e . σ2

(46)

In our numerical illustration, we select a domain x ∈ [−1, 1], with the number of space discretizations of N = 150. The Fourier derivative is computed using IPython [34] and the NumPy library [42]. Fig. 3 shows the Gauss waveform and the comparisons between the analytical Im(f (x + i∆x)) given by the following expression 2
2 Im(f (x + i∆x)) = Im e−1/σ (x−x0 +i∆x) , (47) 12

Figure 3: Left: Gauss function. Right: analytical and numerical imaginary part of the imaginary perturbation of the Gauss function.

and the approximation using eq. (43) truncated at the 9th derivative. It is clear that this approximation holds for a discretization enough to describe the Gauss waveform. After having computed the approximate values of the imaginary perturbation, the interpolation error may act as noise. In realistic applications, the noise level might be unknown, however, there may be cases where the order of magnitude of the noise is known. In cases when there is no information about the noise magnitude, procedures to deal with the unknown noise can be first applied. The work presented by [33, 22] may be combined with noise estimation techniques to evaluate CS derivatives when only information of discrete real values of the functional f (x) are given.

6

Discussion and conclusions

We have performed numerical experiments comparing popular numerical discretizations for the first- and secondorder wave equations and the CSFDM. Although it is well known that simulations of wave propagation depend on the number of points per wavelength assumed and the total distance traveled by the wave [17], our results show that the CSFDM shows always more accurate fields in comparison with respect to those obtained by the FDM on well established numerical discretizations of the wave equations. We have discussed the difference of the CSFDM with respect to previous implementations of CS derivative approximations in well known numerical techniques like the FEM, emphasizing that the CSFDM only modifies the motion conditions, unlike previous implementations which modify the mathematical structure of the numerical technique, leading to different dispersion properties that may induce important differences in the numerical schemes. The methodology presented in this study can improve numerical results in presence of heterogeneities, which we have linked to differences in the values of the Courant number S. This suggests that the CSFDM improves velocity results in presence of structural heterogeneities, which has been one of the main challenges in the numerical solution of wave propagation. The increase in accuracy (temporal and/or spatial) does not imply that the solution offered by the numerical technique is more adequate to solve the wave propagation, emphasizing that the main attribute to focus our attention is to reduce numerical dispersion/dissipation errors. We have proposed a methodology to compute the imaginary part of an imaginary perturbation when only discrete real values of the function is known. This however, requires further studies to quantify the implications of the numerical noise present in the approximation made and further numerical investigations need to be done. More detailed studies are required and will be performed in future contributions, where a complete dispersion analysis of popular numerical techniques like the SEM and/or the DGM, often employed for numerical simulations in highly heterogeneous media, will be performed. This will allow us to quantitatively establish the worst case Courant number S value (or values) that each numerical technique, used for realistic applications, can handle. These future studies will give to the user of the numerical technique a value (or values) of the most appropriate time step ∆t that should be set in the numerical simulation depending on the model in study. In the case of computational geophysics, alternative techniques for numerical computation of derivatives like automatic differentiation tools are being used [36, 29, 28, 1]. The CSM still remain not widely spread in the 13

geoscientific community. There is hope that results presented in this work help to clarify potential uses of the CS derivative based approximations.

7

Acknowledgments

We would like to thank two anonymous reviewers for constructive comments. Rafael Abreu was partly supported by the Spanish Ministry of Economy and Competitiveness (project CGL2015-67130-C2-2-R).

References [1] A. A. Abokhodair. Numerical tools for geoscience computations: semiautomatic differentiation—SD. Computational Geosciences, 11(4):283–296, 2007. [2] A. A. Abokhodair. Complex differentiation tools for geophysical inversion. Geophysics, 74(2):H1–H11, 2009. [3] R. Abreu. Complex Steps Finite Differences with applications to seismic problems. PhD thesis, Universidad de Granada, 2014. [4] R. Abreu, D. Stich, and J. Morales. On the generalization of the complex step method. Journal of Computational and Applied Mathematics, 241(0):84 – 102, 2013. [5] R. Abreu, D. Stich, and J. Morales. The complex-step-finite-difference method. Geophysical Journal International, 202(1):72–93, 2015. [6] A. H. Al-Mohy and N. J. Higham. The complex step approximation to the fr´echet derivative of a matrix function. Numerical Algorithms, 53(1):133–148, 2010. [7] W. K. Anderson, J. C. Newman, D. L. Whitfield, and E. J. Nielsen. Sensitivity analysis for navierstokes equations on unstructured meshes using complex variables. American Institute of Aeronautics and Astronautics, 39(1):56–63, 2001. [8] C. Burg and J. Newman. Computationally efficient, numerically exact design space derivatives via the complex Taylor’s series expansion method. Computers and Fluids, 32(3):373–383, 2003. [9] L. I. Cervi˜ no and T. R. Bewley. On the extension of the complex-step derivative technique to pseudospectral algorithms. Journal of Compututational Physics, 187:544–549, 2003. [10] G. Dziatkiewicz. Complex variable step method for sensitivity analysis of effective properties in multi-field micromechanics. Acta Mechanica, 227(1):11–28, 2016. [11] R. Fielder, A. Montoya, H. Millwater, and P. Golden. Residual stress sensitivity analysis using a complex variable finite element method. International Journal of Mechanical Sciences, 2017. [12] B. Fornberg. The pseudospectral method: Comparisons with finite differences for the elastic wave equation. Geophysics, 52(4):483–501, 1987. [13] T. Furumura, B. L. N. Kennett, and M. Furumura. Seismic wavefield calculation for laterally heterogeneous whole earth models using the pseudospectral method. Geophysical Journal International, 135(3):845–860, 1998. [14] A. Gomez-Farias, A. Montoya, and H. Millwater. Complex finite element sensitivity method for creep analysis. International Journal of Pressure Vessels and Piping, 132:27–42, 2015. [15] A. H¨ urkamp, M. Tanaka, and M. Kaliske. Complex step derivative approximation of consistent tangent operators for viscoelasticity based on fractional calculus. Computational Mechanics, 56(6):1055–1071, 2015. [16] R. W. Ibrahim and H. A. Jalab. The fractional complex step method. Discrete Dynamics in Nature and Society, 2013, 2013. 14

[17] H. Igel. Computational Seismology: A Practical Introduction. Cambridge University Press, 2016. [18] W. Jin, B. Dennis, and B. Wang. Improved sensitivity analysis using a complex variable semi-analytical method. Structural and Multidisciplinary Optimization, 41:433–439, 2010. [19] M. Kawahara. Finite Element Methods in Incompressible, Adiabatic, and Compressible Flows: From Fundamental Concepts to Applications. Springer, 2016. [20] H. Kim and M. Cho. Study on the design sensitivity analysis based on complex variable in eigenvalue problem. Finite Elements in Analysis and Design, 45(12):892 – 900, 2009. [21] J. Kim, D. G. Bates, and I. Postlethwaite. Nonlinear robust performance analysis using complex-step gradient approximation. Automatica, 42(1):177 – 182, 2006. [22] N. Kreji´c, Z. Luˇzanin, F. Nikolovski, and I. Stojkovska. A nonmonotone line search method for noisy minimization. Optimization Letters, 9(7):1371–1391, 2015. [23] K.-L. Lai and J. Crassidis. Extensions of the first and second complex step derivative approximations. Journal of Computational and Applied Mathematics, 219:276–293, 2008. [24] G. Lantoine, R. P. Russell, and T. Dargent. Using multicomplex variables for automatic computation of high-order derivatives. ACM Transactions on Mathematical Software, 38(3):16:1–16:21, 2012. [25] J. Martins, P. Sturdza, and J. J. Alonso. The connection between the complex-step derivative approximation and algorithmic differentiation. American Institute of Aeronautics and Astronautics, 921:2001, 2001. [26] J. R. Martins, P. Sturdza, and J. J. Alonso. The complex-step derivative approximation. ACM Transactions on Mathematical Software (TOMS), 29(3):245–262, 2003. [27] P. Moczo, J. Kristek, and M. G´alis. The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures. Cambridge University Press, 2014. [28] C. Molkenthin, F. Scherbaum, A. Griewank, N. Kuehn, and P. Stafford. A study of the sensitivity of response spectral amplitudes on seismological parameters using algorithmic differentiation. Bulletin of the Seismological Society of America, 104(5):2240–2252, 2014. [29] C. Molkenthin, F. Scherbaum, A. Griewank, H. Leovey, S. Kucherenko, and F. Cotton. Derivative-based global sensitivity analysis: upper bounding of sensitivities in seismic-hazard assessment using automatic differentiation. Bulletin of the Seismological Society of America, 2017. [30] J. F. Monsalvo, M. J. Garc´ıa, H. Millwater, and Y. Feng. Sensitivity analysis for radiofrequency induced thermal therapies using the complex finite element method. Finite Elements in Analysis and Design, 135:11–21, 2017. [31] A. Montoya, R. Fielder, A. Gomez-Farias, and H. Millwater. Finite-element sensitivity for plasticity using complex variable methods. Journal of Engineering Mechanics, 141(2):04014118, 2014. [32] A. Montoya and H. Millwater. Sensitivity analysis in thermoelastic problems using the complex finite element method. Journal of Thermal Stresses, 40(3):302–321, 2017. [33] F. Nikolovski and I. Stojkovska. Complex-step derivative approximation in noisy environment. Journal of Computational and Applied Mathematics, 327:64 – 78, 2018. [34] F. P´erez and B. E. Granger. Ipython: a system for interactive scientific computing. Computing in Science & Engineering, 9(3):21–29, 2007. [35] M. S. Ridout. Statistical applications of the complex-step method of numerical differentiation. The American Statistician, 63(1):66–74, 2009. [36] M. Sambridge, P. Rickwood, N. Rawlinson, and S. Sommacal. Automatic differentiation in geophysical inverse problems. Geophysical Journal International, 170(1):1–8, 2007. 15

[37] P. Smereka. The numerical approximation of a delta function with application to level set methods. Journal of Computational Physics, 211(1):77 – 90, 2006. [38] W. Squire and G. Trapp. Using complex variables to estimate derivatives of real functions. SIAM Journals Online, 40(1):110–112, 1998. [39] M. Tanaka, M. Fujikawa, D. Balzani, and J. Schr¨ oder. Complex-step derivative approximation schemes for the robust calculation of numerical constitutive tangent moduli. Proceedings of Applied Mathematics and Mechanics, 1:167–168, 2013. [40] M. Tanaka, M. Fujikawa, D. Balzani, and J. Schr¨ oder. Robust numerical calculation of tangent moduli at finite strains based on complex-step derivative approximation and its application to localization analysis. Computer Methods in Applied Mechanics and Engineering, 269:454–470, 2014. [41] J. W. Thomas. Numerical Partial Differential Equations: Finite Difference Methods. Springer, 1995. [42] S. van der Walt, S. C. Colbert, and G. Varoquaux. The numpy array: a structure for efficient numerical computation. Computing in Science & Engineering, 13(2):22–30, 2011. [43] A. Voorhees, H. Millwater, and R. Bagley. Complex variable methods for shape sensitivity of finite element models. Finite Elements in Analysis and Design, 47(10):1146 – 1156, 2011. [44] A. Voorhees, H. Millwater, R. Bagley, and P. Golden. Fatigue sensitivity analysis using complex variable methods. International Journal of Fatigue, 40:61 – 73, 2012. [45] D. Wagner and H. Millwater. 2D weight function development using a complex Taylor series expansion method. Engineering Fracture Mechanics, 86:23 – 37, 2012. [46] B. P. Wang and A. P. Apte. Complex variable method for eigensolution sensitivity analysis. American Institute of Aeronautics and Astronautics, 44(12):2958–2961, 2006. [47] R. E. Wengert. A simple automatic derivative evaluation program. Communications of the ACM, 7(8):463– 464, 1964.

A

Ricker analytical solution

Due to the lack of an accessible reference, we present the derivation of the analytical solution for the second-order wave equation in presence of a Ricker source time dependent function. The initial value problem is given by the following expression  ∂2u ∂2u    − a2 2 = f (x, t), −∞ < x < ∞, t > 0,  2  ∂t ∂x (48) u(x, 0) = ϕ(x), −∞ < x < ∞,     ∂u(x, 0) = ψ(x),  −∞ < x < ∞. ∂t According to d’Alembert’s formula, the solution of (48) can be expressed as follows " # # Z x+at Z x+at Z t" Z x+a(t−τ ) ∂ 1 1 1 u(x, t) = ϕ(ξ) dξ + ψ(ξ) dξ + f (ξ, τ )dξ dτ ∂t 2a x−at 2a x−at 2a x−a(t−τ ) 0 Z x+at Z t Z x+a(t−τ ) 1 1 1 = [ϕ(x + at) + ϕ(x − at)] + ψ(ξ) dξ + dτ f (ξ, τ )dξ. 2 2a x−at 2a 0 x−a(t−τ )

If we set the initial motion conditions as    ϕ(x) = 0, ψ(x) = 0,   f (x, t) = δ(x − x0 )R(t − t0 ),

(49)

(50) with

2 t2 ) exp(−π 2 f 2 t2 ), R(t) = (1 − 2π 2 fM M

where fM is the dominant frequency, t0 is the time delay, δ is Dirac delta function, x0 the source location and R(t) represents a Ricker wavelet. Herein, we can write the general solution (49) as follows Z t Z x+a(t−τ ) Z t Z x+a(t−τ ) 1 1 (51) u(x, t) = dτ f (ξ, τ )dξ = R(τ − t0 ) dτ δ(ξ − x0 )dξ. 2a 0 2a 0 x−a(t−τ ) x−a(t−τ )

16

Before continuing, we recall from the properties of the Dirac delta function δ, we admit that ∀ε > 0, Z x0 +ε δ(x − x0 )f (x) dx = f (x0 ).

(52)

x0 −ε

To simplify our derivation of the analytical solution, we first consider the case where x0 < x, in other words, we only consider the integral respect to ξ, that is, within the limits [x − a(t − τ ), x + (at − τ )]. Our primary goal is to verify whether x0 is in the range (x − a(t − τ ), x + (at − τ )). If x0 is not in the range, then the second integration is Z x+a(t+τ ) δ(ξ − x0 )dξ = 0. (53) x−a(t−τ )

x − x0 = τ0 , we can always find ε that 0 < ε < x0 − (x − a(t − τ )). The range We know that if x0 > x − a(t − τ ), that is, τ < t − a of integration will transform from [x − a(t − τ ), x + a(t − τ )] to [x − a(t − τ ), x0 − ε] ∪ (x0 − ε, x0 + ε) ∪ [x0 + ε, x + a(t − τ )]. We can write the integral corresponding to the variable ξ in eq. 51 as follows Z x+a(t−τ ) Z x0 −ε Z x0 +ε Z x+a(t−τ ) δ(ξ − x0 )dξ = δ(ξ − x0 )dξ + δ(ξ − x0 )dξ + δ(ξ − x0 )dξ x−a(t−τ )

x0 −ε

x−a(t−τ )

x+ε

(54)

=0+1+0 =1

if

τ
x − x0 . a

We write Z

x+a(t−τ )

x−a(t−τ )

δ(ξ − x0 )dξ =

  0  1

if if

x − x0 , a x − x0 τ t−

(55)

x0 − x . If x0 ∈ [x − at, x + at], that is, x0 − at < x < x0 + at, a |x − x0 | there must exist τ0 ∈ [0, t] subject to x + a(t − τ0 ) = x0 or x − a(t − τ0 ) = x0 , then τ0 = t − and we can write the following a Z t Z x+a(t−τ ) 1 u(x, t) = dτ δ(ξ − x0 )R(τ − t0 )dξ 2a 0 x−a(t−τ ) Z τ0 Z t Z x+a(t−τ ) Z x+a(t−τ ) 1 1 = dτ dτ δ(ξ − x0 )R(τ − t0 )dξ + δ(ξ − x0 )R(τ − t0 )dξ 2a 0 2a τ0 x−a(t−τ ) x−a(t−τ ) Z τ0 Z τ0 1 1 2 2 R(τ − t0 ) dτ = (1 − 2π 2 fM (τ − t0 )2 ) exp(−π 2 fM (τ − t0 )2 ) dτ = 2a 0 2a 0 2 π 2 (t − τ )2 ) τ0 2 2 2 2 2 2 (τ − t0 ) exp(−fM |x − x0 | 0 = (τ0 − t0 ) exp(−fM π (t0 − τ0 ) ) + t0 exp(−fM π t0 ) , = with τ0 = t − . 2a 2a a We have derived the solution for x < x0 , which involves τ0 = t −

0

(56) This solution also indicates the physical meaning of wave equation: τ0 is the moment at which the wave reached the point x and priori to this moment, u(x, t) = 0 and if x ≤ x0 − at or x ≥ x0 + at, the displacement is zero u(x, t) = 0. Finally we can write the solution (51) as follows  2 π 2 (τ − t )2 ) + t exp(−f 2 π 2 t2 )  0 0 0  (τ0 − t0 ) exp(−fM 0 M , x ∈ B(x0 , at), u(x, t) = (57) 2a  0, otherwise, where B(x0 , at) ⇔ x0 − at < x < x0 + at and τ0 = t −

|x − x0 | . a

17