# An Improved Lax-Wendroff Scheme for Two-Dimensional Transient Thermal Simulation
## 2 Methods
### Difference Operators and Coefficient Formulas
#### 2.3.2 HONFD Operator and Coefficient Formulas  
[PERSON] and [PERSON] [30] give fourth-order-accuracy FD operators to approximate the second-order spatial derivative in the boundary region. Inspired by the reference [31], we develop a formula for the coefficients of the HONFD operator, and use a strategy of combining the HONFD and HOCFD approximation to calculate the spatial derivatives of the boundary region. The HONFD operator of the fourth-order spatial derivative along the \(x\) coordinate with \((N-3)\)-order accuracy, which can be written as  
\[D_{x}^{4}[T^{n}(i,j)]=\frac{1}{\Delta x^{4}}\sum\limits_{i=1}^{N}p_{i}T^{n}(i +m_{i},j) \tag{11}\]  
where \(p_{i}\) represents the HONFD coefficient, which can be calculated through Equation (12) (see Appendix B for detailed derivation)  
\[p_{i}=4!\left(\frac{\prod\limits_{ \begin{subarray}{c}n-5\cdots>n_{1}=2\\ n_{N-5}\cdots>n_{1}=2\end{subarray} \prod\limits_{j_{1},\cdots,j_{N-5}\
eq i}q_{ r_{1}}q_{r_{2}}\cdots q_{r_{N-5}}}}{\prod\limits_{k=1}^{N}(q_{k}-q_{i})}\right) \tag{12}\]  
where the stencil vector \(q=[q_{1},q_{2},\cdots,q_{N}]\), \(N\) represents the length of the operator, \(q_{i}=m_{i}\), \(m=[-L_{p},-L_{p}+1,\cdots,0,1,2,\cdots,L_{q}]\), and \(L_{p}+L_{q}=N-1\), \(L_{p}\), \(L_{q}\) are all positive integers and \(m\) is also a vector. \(\Delta x\) represents the spatial interval in the \(x\) direction. \(m_{i}=0\) represents the derivative at the mesh point in the boundary area, and also represents the zero point position of the operator vector. The left boundary area of the two-dimensional model in the \(x\) direction, \(L_{p}\), represents the number of meshes from the derivative mesh point or operator zero point to the left boundary point. \(L_{q}\) represents the length to the right side of the zero value point of the operator, \(L_{q}=(N-1)-L_{p}\). The right boundary area, \(L_{q}\), represents the number of meshes from the zero position of the operator to the right boundary point, and the \(L_{p}\) represents the length to the left side of the zero position of the operator, \(L_{p}=(N-1)-L_{q}\). The derivatives in the \(y\) direction in the boundary area are calculated the same as in the \(x\) direction. \(i\) represents the sequence number of the difference operator vector element from left to right \(i=1,2,\cdots,N\). The sixth-order HONFD operators and coefficients for the fourth-order spatial derivative are shown in Table 1. The operators' lengths are \(N=7\). The \(L_{p}\) and \(L_{q}\) of different operators with serial numbers 1-5 are 1, 2, 5, 4, 3 and 5, 4, 1, 2, 3, respectively. Table 1 only shows the sixth-order HONFD operators and their difference coefficients for the left boundary zone (serial numbers 1, 2) and right boundary zone (serial numbers 3, 4) as well as the internal zone (serial number 5) of the two-dimensional model in the \(x\) direction. The sixth-order HONFD operator in the top and bottom boundary zone and its difference coefficients in the \(y\) direction are the same as in the \(x\) direction. As can be seen from Table 1, the HONFD coefficients of the right boundary zone are consistent with the difference coefficients of the same operator (such as serial numbers 1 and 3) in the left boundary zone arranged in reverse order. Therefore, the two-dimensional model only needs to calculate the HONFD coefficients for either boundary region without an additional computational cost. Moreover, Table 1 demonstrates that, when the HONFD operator becomes an HOCFD operator (see Table 1 serial number 5), The coefficients calculated using the coefficient formula (see Equation (12)) of the HONFD operator are consistent with the results calculated using the coefficient formula (see Equation (9b)) of the HOCFD operator.