# An Improved Lax-Wendroff Scheme for Two-Dimensional Transient Thermal Simulation
## 2 Methods
### Difference Operators and Coefficient Formulas  
#### 2.3.1 HOCFD Operator and Coefficient Formula  
For the two-dimensional transient HCEs, the 2\(M\)-order HOCFD operator of the second-order spatial derivative along the \(x\) and \(y\) coordinates can be expressed as [26]  
\[ \begin{split} D_{x}^{2}[T^{n}(i,j)]&=\frac{1}{ \Delta x^{2}}(\sum\limits_{m=1}^{M}a_{m}(T^{n}(i+m,j)+T^{n}(i-m,j)-2T^{n}(i,j) ))\\ D_{y}^{2}[T^{n}(i,j)]&=\frac{1}{\Delta y^{2}}( \sum\limits_{m=1}^{M}a_{m}(T^{n}(i,j+m)+T^{n}(i,j-m)-2T^{n}(i,j)))\end{split}  \tag{6}\]  
and the 2\(M\)-order HOCFD operator of the fourth-order spatial derivative along the \(x\) and \(y\) coordinates can be can be expressed as  
\[ \begin{split} D_{x}^{4}[T^{n}(i,j)]&=\frac{1}{ \Delta x^{4}}(\sum\limits_{m=1}^{M+1}b_{m}(T^{n}(i+m,j)+T^{n}(i-m,j)-2T^{n}(i, j)))\\ D_{y}^{4}[T^{n}(i,j)]&=\frac{1}{\Delta y^{4}}( \sum\limits_{m=1}^{M+1}b_{m}(T^{n}(i,j+m)+T^{n}(i,j-m)-2T^{n}(i,j)))\end{split}  \tag{7}\]  
as well as the 2\(M\)-order HOCFD operator of the fourth-order spatial mixed derivative  
\[c_{m,l}(T^{n}(i+m,j+l)+T^{n}(i-m,j+l)+T^{n}(i+m,j-l)+T^{n}(i-m,j-l)) \tag{8}\]  
where \(\Delta x\) and \(\Delta y\) represent the mesh interval along the \(x\) and \(y\) coordinates, respectively. \(M\) represents half the length of the HOCFD operator. There are other similarities along the \(y\) coordinate. The HOCFD operator coefficient \(a_{m}\) can be determined by the analytic formula of Equation (9a) (see [26] for a detailed derivation). \[a_{m}=\left\{ \begin{array}{cc}1&m=1,M=1\\ \frac{(-1)^{n+1}\prod\limits_{i=l,j=m}^{M}i^{2}}{m^{2}\prod\limits_{i=1}^{M} (m^{2}-i^{2})\prod\limits_{i=m+1}^{M}(i^{2}-m^{2})}&m=1,2,\cdots,M,M\geq 2 \end{array} \right.\] (9a) and \[b_{m}\] can be determined by the analytic formula of Equation ( 9b ) (see Appendix A for a detailed derivation). \[ \begin{split}&\text{(1)}\quad M=1,\,b_{1}=-4,\,b_{2}=1\\ &\text{(2)}\quad M=2,\,b_{1}=-13/2,\,b_{2}=2.0,\,b_{3}=1/6\\ &\text{(3)}\quad M\geq 3\end{split} \]\[b_{m}=\left\{ \begin{array}{c}\frac{(-1)^{n}4!}{2}\frac{\prod\limits_{i=3}^{M+1} \beta^{2}+4\prod\limits_{i=-3}^{M+1}\beta^{M+1}\sum\limits_{i=3}^{M+1}1/^{2}}{ \prod\limits_{i=1}^{M}(m^{2}-i^{2})},m=1\\ \frac{(-1)^{n}4!}{2m^{2}}\frac{\prod\limits_{i=1}^{M+1}\beta^{2}+\prod\limits_ {i=-3}^{M+1}\beta^{2}\sum\limits_{i=-3}^{M+1}1/^{2}}{\prod\limits_{i=1}^{M}(m^ {2}-i^{2})\prod\limits_{i=+1}^{M+1}(\beta^{2}-m^{2})},m=2,3,\cdots,M+1\end{array} \right. \tag{9b}\]  
For the fourth-order spatial mixed derivatives \(\partial^{4}T^{n}(i,j)/\partial x^{2}\partial y^{2}\), from \(\partial^{4}/\partial x^{2}\partial y^{2}=(\partial^{2}/\partial x^{2})( \partial^{2}/\partial y^{2})\), using the product of the \(2M\)-order HOFD coefficients of the second-order spatial derivatives along the x and y coordinates (see Equation (9a)), respectively, one can calculate the \(2M\)-order HOFD coefficients matrix of the fourth-order spatial mixed derivative by means of a simple method that can easily treat the boundary (see Equation (10)). \[C_{1/4}= \begin{bmatrix}a_{0}^{2}&a_{0}a_{1}&\cdots&a_{0}a_{M}\\ a_{0}a_{1}&a_{1}^{2}&\cdots&a_{1}a_{M}\\ \vdots&\vdots&\vdots&\vdots\\ a_{0}a_{M}&a_{1}a_{M}&\cdots&a_{M}^{2}\end{bmatrix}  \tag{10}\]  
where \(C_{1/4}\) represents the \(1/4\) block matrix of the coefficient matrix \(C\). The \(2M\)-order HOFD coefficients matrix \(C\) of the fourth-order spatial mixed derivatives is a symmetric matrix, i.e., \(c_{m,l}=c_{l,m}=a_{m}a_{l}\), \(m=1,2,\cdots,M\), \(l=1,2,\cdots,M\), \(a_{m}\), \(a_{l}\) represent the \(2M\)-order HOFD coefficients of the second-order spatial derivatives along the \(x\) and \(y\) coordinates, respectively (see Equation (9a)), and \(a_{0}=-2\sum\limits_{m=1}^{M}a_{m}\).