*2.3. Numerical Implementation*

Introducing the finite differences and derivatives can be summarized as follows:

$$\begin{aligned} \frac{\partial u\_x}{\partial x} &\approx \frac{1}{2\Delta x} \left[ u\_x^{i+1,j} - u\_x^{i-1,j} \right], \\\\ \frac{\partial u\_x}{\partial x} &\approx \frac{1}{2\Delta x} \left[ u\_x^{i,j+1} - u\_x^{i,j-1} \right], \\\\ \frac{\partial^2 u\_x}{\partial x^l} &\approx \frac{1}{\Delta x^l} \left[ u\_x^{i+1,j} - 2u\_x^{i,j} + u\_x^{i-1,j} \right], \\\\ \frac{\partial^2 u\_y}{\partial x^2} &\approx \frac{1}{\Delta x^2} \left[ u\_x^{i,j+1} - 2u\_x^{i,j} + u\_x^{i,j-1} \right], \\\\ \frac{\partial^2 u\_y}{\partial x \partial z} &\approx \frac{1}{4\Delta x \Delta z} \left[ u\_x^{i+1,j+1} - u\_x^{i+1,j-1} - u\_x^{i-1,j+1} + u\_x^{i-1,j-1} \right], \end{aligned} \tag{11}$$

where Δ*x* and Δ*z* are the horizontal and vertical grid sizes, respectively; and the value of grid sizes must be small enough so that there are enough grid points per wavelength [10]; *i* and *j* are the indices of the horizontal and vertical grids. The finite differences of *uz* are similar to the way of the horizontal one.

Using the finite-difference star stencil [10]:

*∂*2*ux*

$$
\begin{array}{ccccc}
\mathbf{M}\_1 & \mathbf{M}\_4 & \mathbf{M}\_7 \\
\mathbf{M}\_2 & \mathbf{M}\_5 & \mathbf{M}\_8 \\
\mathbf{M}\_3 & \mathbf{M}\_6 & \mathbf{M}\_9
\end{array}.\tag{12}
$$

Substituting the finite difference of derivatives and each grid point of Equation (12) according to Figure 2 can be collected as follows:

After using the finite difference method in the frequency domain, the final equation can be written as Equation (14). Thus, for each frequency component, the matrix equation can be given as:

$$\mathbf{M}(\omega)\mathbf{U}(\omega) = \mathbf{F}(\omega). \tag{14}$$

The impedance matrix **M**(*ω*) can be obtained by Equation (13), **U**(*ω*) is the monochromatic wavefield, and **F**(*ω*) is the source term. The frequency-domain Ricker wavelet is used in the following numerical examples, and LU decomposition is used to solve the matrix Equation (14) to get the monochromatic wavefield.

**<sup>M</sup>**<sup>1</sup> <sup>=</sup> (*C*13+*C*44)*i*,*<sup>j</sup>* 4Δ*x*Δ*z*(*exez* )*i*,*<sup>j</sup>* 0 1 1 0 , **M**<sup>2</sup> = ⎡ ⎢ ⎢ ⎢ ⎣ (*C*11)*i*,*<sup>j</sup>* Δ*x*(*ex*)*i*,*<sup>j</sup>* <sup>2</sup> <sup>+</sup> (*C*11)*i*,*<sup>j</sup>* 2Δ*x*(*ex*) 3 *i*,*j* · *∂*(*ex*)*i*,*<sup>j</sup> ∂x* (*C*44)*i*,*<sup>j</sup>* 2Δ*x*(*ez e*<sup>2</sup> *<sup>x</sup>* )*i*,*<sup>j</sup>* · *∂*(*ex*)*i*,*<sup>j</sup> ∂z* (*C*13)*i*,*<sup>j</sup>* 2Δ*x*(*ez e*<sup>2</sup> *<sup>x</sup>* )*i*,*<sup>j</sup>* · *∂*(*ex*)*i*,*<sup>j</sup> ∂z* (*C*44)*i*,*<sup>j</sup>* (*ex*)*i*,*<sup>j</sup>* Δ*x* <sup>2</sup> <sup>+</sup> (*C*44)*i*,*<sup>j</sup>* 2Δ*x*(*ex*) 3 *i*,*j* · *∂*(*ex*)*i*,*<sup>j</sup> ∂x* ⎤ ⎥ ⎥ ⎥ ⎦ , **<sup>M</sup>**<sup>3</sup> <sup>=</sup> <sup>−</sup>(*C*13+*C*44)*i*,*<sup>j</sup>* 4Δ*x*Δ*z*(*exez* )*i*,*<sup>j</sup>* 0 1 1 0 , **M**<sup>4</sup> = ⎡ ⎢ ⎢ ⎢ ⎣ (*C*44)*i*,*<sup>j</sup>* Δ*z*(*ez* )*i*,*<sup>j</sup>* <sup>2</sup> <sup>+</sup> (*C*44)*i*,*<sup>j</sup>* 2Δ*z*(*ez* ) 3 *i*,*j* · *∂*(*ez* )*i*,*<sup>j</sup> ∂z* (*C*13)*i*,*<sup>j</sup>* 2Δ*z*(*exe*<sup>2</sup> *<sup>z</sup>* )*i*,*<sup>j</sup>* · *∂*(*ez* )*i*,*<sup>j</sup> ∂x* (*C*44)*i*,*<sup>j</sup>* 2Δ*z*(*exe*<sup>2</sup> *<sup>z</sup>* )*i*,*<sup>j</sup>* · *∂*(*ez* )*i*,*<sup>j</sup> ∂x* (*C*33)*i*,*<sup>j</sup>* (*ez* )*i*,*<sup>j</sup>* Δ*z* <sup>2</sup> <sup>+</sup> (*C*33)*i*,*<sup>j</sup>* 2Δ*z*(*ez* ) 3 *i*,*j* · *∂*(*ez* )*i*,*<sup>j</sup> ∂z* ⎤ ⎥ ⎥ ⎥ ⎦ , **M**<sup>5</sup> = ⎡ ⎢ ⎢ ⎢ ⎣ *<sup>ω</sup>*2*ρi*,*<sup>j</sup>* <sup>−</sup> <sup>2</sup>(*C*11)*i*,*<sup>j</sup>* Δ*x*(*ex*)*i*,*<sup>j</sup>* <sup>2</sup> <sup>−</sup> <sup>2</sup>(*C*44)*i*,*<sup>j</sup>* Δ*z*(*ez* )*i*,*<sup>j</sup>* <sup>2</sup> 0 <sup>0</sup> *<sup>ω</sup>*2*ρi*,*<sup>j</sup>* <sup>−</sup> <sup>2</sup>(*C*44)*i*,*<sup>j</sup>* Δ*x*(*ex*)*i*,*<sup>j</sup>* <sup>2</sup> <sup>−</sup> <sup>2</sup>(*C*33)*i*,*<sup>j</sup>* Δ*z*(*ez* )*i*,*<sup>j</sup>* 2 ⎤ ⎥ ⎥ ⎥ ⎦ , **M**<sup>6</sup> = ⎡ ⎢ ⎢ ⎢ ⎣ (*C*44)*i*,*<sup>j</sup>* Δ*z*(*ez* )*i*,*<sup>j</sup>* <sup>2</sup> <sup>−</sup> (*C*44)*i*,*<sup>j</sup>* 2Δ*z*(*ez* ) 3 *i*,*j* · *∂*(*ez* )*i*,*<sup>j</sup> ∂z* −(*C*13)*i*,*<sup>j</sup>* 2Δ*z*(*exe*<sup>2</sup> *<sup>z</sup>* )*i*,*<sup>j</sup>* · *∂*(*ez* )*i*,*<sup>j</sup> ∂x* −(*C*44)*i*,*<sup>j</sup>* 2Δ*z*(*exe*<sup>2</sup> *<sup>z</sup>* )*i*,*<sup>j</sup>* · *∂*(*ez* )*i*,*<sup>j</sup> ∂x* (*C*33)*i*,*<sup>j</sup>* (*ez* )*i*,*<sup>j</sup>* Δ*z* <sup>2</sup> <sup>−</sup> (*C*33)*i*,*<sup>j</sup>* 2Δ*z*(*ez* ) 3 *i*,*j* · *∂*(*ez* )*i*,*<sup>j</sup> ∂z* ⎤ ⎥ ⎥ ⎥ ⎦ , **<sup>M</sup>**<sup>7</sup> <sup>=</sup> <sup>−</sup>(*C*13+*C*44)*i*,*<sup>j</sup>* 4Δ*x*Δ*z*(*exez* )*i*,*<sup>j</sup>* 0 1 1 0 , **M**<sup>8</sup> = ⎡ ⎢ ⎢ ⎢ ⎣ (*C*11)*i*,*<sup>j</sup>* Δ*x*(*ex*)*i*,*<sup>j</sup>* <sup>2</sup> <sup>−</sup> (*C*11)*i*,*<sup>j</sup>* 2Δ*x*(*ex*) 3 *i*,*j* · *∂*(*ex*)*i*,*<sup>j</sup> ∂x* −(*C*44)*i*,*<sup>j</sup>* 2Δ*x*(*ez e*<sup>2</sup> *<sup>x</sup>* )*i*,*<sup>j</sup>* · *∂*(*ex*)*i*,*<sup>j</sup> ∂z* −(*C*13)*i*,*<sup>j</sup>* 2Δ*x*(*ez e*<sup>2</sup> *<sup>x</sup>* )*i*,*<sup>j</sup>* · *∂*(*ex*)*i*,*<sup>j</sup> ∂z* (*C*44)*i*,*<sup>j</sup>* (*ex*)*i*,*<sup>j</sup>* Δ*x* <sup>2</sup> <sup>−</sup> (*C*44)*i*,*<sup>j</sup>* 2Δ*x*(*ex*) 3 *i*,*j* · *∂*(*ex*)*i*,*<sup>j</sup> ∂x* ⎤ ⎥ ⎥ ⎥ ⎦ , **<sup>M</sup>**<sup>9</sup> <sup>=</sup> (*C*13+*C*44)*i*,*<sup>j</sup>* 4Δ*x*Δ*z*(*exez* )*i*,*<sup>j</sup>* 0 1 1 0 . (13)
