# Optimizing Finite-Difference Operator in Seismic Wave Numerical Modeling
## 2 Methods  
### Finite-Difference Method  
We start with a brief introduction of the FD method. Based on studies of [PERSON] (2017) and [PERSON] et al. (2018), we describe the derivation process of FD operators as follows. First, the FD operators can be derived by the sinc interpolator. According to [PERSON]'s theorem, a continuous signal \(f(x)\) with limited bandwidth can be reconstructed by interpolating the normalized sinc function. Roughly speaking, the FD method represents waveforms with discrete samples. \[sinc(x)= \begin{cases}1,\text{ if }x=0\\ \frac{sin(\pi x)}{\pi x},\text{ otherwise}\end{cases} \qquad f(x)=\sum_{n=- \infty}^{+\infty}sinc\bigg{(}\frac{x-n\Delta x}{\Delta x}\bigg{)}f(n\Delta x ), \tag{1}\]  
where \(\Delta x\) is the sampling interval. Thus, \(1/\Delta x\) is the sampling frequency, and \(f(n\Delta x)\) is the infinite uniformly spaced discrete sequence of \(f(x)\). The first and second derivatives of \(f\) are derived from Equation (1). The derivation can be simplified by letting \(x=0\) (see Appendix A.1). \[\left.\frac{\partial f}{\partial x}\right|_{x=0}=\frac{1}{\Delta x}\sum_{n=- \infty}^{+\infty}\bigg{(}-\frac{1}{n}\cos(n\pi)\bigg{)}f(n\Delta x), \tag{2}\]\[\left.\frac{\partial^{2}f}{\partial x^{2}}\right|_{x=0}=\frac{1}{\Delta x^{2}} \sum_{n=-\infty}^{+\infty}\left(-\frac{2}{n^{2}}\cos(n\pi)\right)f(n\Delta x). \tag{3}\]  
We use a window function \(\omega(n)\) with length \(N+1\) to truncate the derivative function Equations (2) and (3). \(N\) is an even number, which is related to the accuracy. Windowing is a weighting process that truncates the function of infinite length to \([-N/2,N/2]\). The FD operators can be calculated as:  
\[\left.\frac{\partial f}{\partial x}\right|_{x=0}\approx\frac{1}{\Delta x} \sum_{n=-N/2}^{N/2}\omega(n)\left(-\frac{1}{n}\cos(n\pi)\right)f(n\Delta x), \tag{4}\]  
\[\left.\frac{\partial^{2}f}{\partial x^{2}}\right|_{x=0}\approx\frac{1}{ \Delta x^{2}}\sum_{n=-N/2}^{N/2}\omega(n)\left(-\frac{2}{n^{2}}\cos(n\pi) \right)f(n\Delta x). \tag{5}\]  
[PERSON] and [PERSON] [22] pointed out that \(w(n)\) in the conventional FD operator is a binomial constant coefficient, and they referred to the conventional \(w(n)\) as a binomial window. One method to obtain a better FD operator is to design a more suitable \(w(n)\) function, which requires more control over the parameters. For the optimization algorithms, we need to set an objective function. Consider the coefficients of the function Equations (4) and (5) as a whole:  
\[\left.\frac{\partial f}{\partial x}\right|_{x=0}\approx\frac{1}{\Delta x} \sum_{n=-N/2}^{N/2}b_{n}f(n\Delta x),\qquad\left.\frac{\partial^{2}f}{\partial x ^{2}}\right|_{x=0}\approx\frac{1}{\Delta x^{2}}\sum_{n=-N/2}^{N/2}c_{n}f(n \Delta x), \tag{6}\]  
where  
\[b_{n}=-\frac{1}{n}\cos(n\pi)\omega(n),\qquad c_{n}=-\frac{2}{n^{2}}\cos(n\pi) \omega(n),\]  
\(b_{n}\) and \(c_{n}\) are the finite-difference coefficients for the first and second derivative, respectively. They are related by \(c_{n}=2b_{n}/n\). The FD weights can also be obtained through the Taylor series expansions (Thongyoy and Boonyasiriwat [24]). It reflects the relation between finite-difference operators and the pseudospectral method. The finite-difference operator is formulated as (see Appendix A.2):  
\[\frac{\partial f}{\partial x}\approx\frac{1}{\Delta x}\sum_{n=1}^{N/2}b_{n}(f _{n}-f_{-n}), \tag{7}\]  
\[\frac{\partial^{2}f}{\partial x^{2}}\approx\frac{1}{\Delta x^{2}}\Bigg{[}c_{ 0}f_{0}+\sum_{n=1}^{N/2}c_{n}(f_{n}+f_{-n})\Bigg{]}, \tag{8}\]  
where \(f_{n}=f(x+n\Delta x)\), \(b_{n}\) and \(c_{n}\) are defined similarly as before. \(\sum_{n=-N/2}^{N/2}b_{n}=\sum_{n=-N/2}^{N/2}c_{n}=0\), \(b_{-n}=-b_{n}\), \(c_{-n}=c_{n}\). Hence, for \(n=0\),  
\[b_{0}=0,\qquad c_{0}=-2\sum_{n=1}^{N/2}c_{n}. \tag{9}\]  
Performing the Fourier transform on Equations (7) and (8), the data can be transformed from the conventional spatial and time domains to the frequency-wavenumber domain, which is also known as the F-k domain (see Appendix A.3).