## 3 Iterative Algorithm of Two-Dimensional Central Difference Galerkin Spectral Method  
Assume that \(\Delta t\) is expressed as a time step, and finite sequences \(\{x_{n}\}\) and \(\{y_{m}\}\) represent nodes in the x and y directions, respectively. Here, Chebyshev-Gauss-Lobatto type configuration points are selected as space nodes. \[x_{n}=-cos\frac{n\pi}{H_{x}},\ n=0,1,2,\cdot\cdot,H_{x},\ y_{m}=-cos\frac{m \pi}{H_{y}},\ m=0,1,2,\cdot\cdot\cdot,H_{y} \tag{4}\]  
By combining the theory of the Galerkin spectrum method and the temperature field function at the time node \(u(x,y,t)\in L^{2}(\Omega;[0,T])\), \(i\Delta t\) can be constructed:  
\[u(x,y,it)\approx u_{H}^{i}(x,y)=\sum_{f=0}^{H}\sum_{g=0}^{H}a_{f,g}^{i}F_{f}( x)F_{g}(y) \tag{5}\]  
where both \(\left\{F_{f}(x)\right\}_{f=0}^{H}\) and \(\left\{F_{g}(y)\right\}_{g=0}^{H}\) are Chebyshev orthogonal polynomials (see, e.g., [25]), and \(\left\{a_{f,g}^{i}\right\}_{f,g=0}^{H}\) is the spectral coefficient, when \(x=x_{n},y=y_{m},t=i\Delta t\), \(u(x,y,t)\)can replace \(u(x_{n},y_{m},i\Delta t)\), the effective approximate discrete form of temperature field \(u(x,y,t)\) on the space node is obtained. That is, the first discrete form:  
\[u(x_{n},y_{m},i\Delta t)\approx u_{H}^{i}(x_{n},y_{m})=\sum_{f=0}^{H}\sum_{g=0}^ {H}a_{f,g}^{i}F_{f}(x_{n})F_{g}(y_{m}) \tag{6}\]  
In addition, from Chebyshev's forward discrete transformation, the following equation can be obtained:  
\[a_{j}=\frac{1}{b_{j}H}\sum_{k}^{H}\frac{1}{b_{k}}u_{H}(x_{n})F_{j}(x_{n}) \tag{7}\]  
Thus, the second discrete form of the effective approximation function can be obtained:  
\[u_{H}^{i}(x_{n},y_{m})=\sum_{l=0}^{H}W_{j}^{i}F(x_{n})\ 0\leq n\leq H \tag{8}\]  
where  
\[W_{j}^{i}=\frac{2}{b_{j}H}\sum_{n=0}^{H}\frac{1}{b_{n}}u_{H}(x_{n},y_{m},i \Delta t)\cos\frac{nj\pi}{H}\ 0\leq j\leq H \tag{9}\]  
as well as,  
\[u_{H}^{i}(x_{n},y_{m})=\sum_{l=0}^{H}V_{j}^{i}F(y_{m})\ 0\leq m\leq H \tag{10}\]  
\[V_{j}^{i}=\frac{2}{b_{j}H}\sum_{n=0}^{H}\frac{1}{b_{m}}u_{H}(x_{n},y_{m},i \Delta t)\cos\frac{mj\pi}{H}\ 0\leq j\leq H \tag{11}\]  
where  
\[b_{j}= \begin{cases}1&1\leq j\leq H-1\\ 2&n=0,H\end{cases}  \tag{12}\]  
Furthermore, the recurrence relation of Chebyshev polynomials is  
\[F_{l+1}(x)=2 xF_{l}(x)-F_{l-1}(x)\ l\geq 1 \tag{13}\]  
From (13), the following expressions for the first derivative and the second derivative of Chebyshev polynomials can be obtained. \[F_{l}^{\prime}(x)=\sum_{ \begin{subarray}{c}i=0\\ i+l=odd\end{subarray} }^{l-1}\frac{2 lF_{l}(x)}{b_{i}}\,F_{l}^{\prime\prime}(x)= \sum_{ \begin{subarray}{c}i=0\\ i+l=even\end{subarray} }^{l-2}\frac{lF_{i}(x)}{b_{i}}(l^{2}-i^{2}) \tag{14}\]  
Obviously, from (11) and (14), it is easy to obtain the approximation function of the second partial derivative of the temperature field in the \(x\) direction and \(y\) direction on the Chebyshev-Gauss-Lobatto configuration point. \[ \begin{array}{c}\frac{\partial^{2}u_{Hx}^{i}(x_{n},y_{m})}{\partial x^{2}} =\sum_{j}^{H}Y_{j}^{i}F_{j}^{\prime\prime}(x_{n})=\sum_{j}^{H}W_{j}^{j-2} \sum_{\begin{subarray}{c}i=0\\ i+l=even\end{subarray}}^{l-2}\frac{lF_{i}(x_{n})}{b_{i}}(j^{2}-i^{2})\\ \\ \frac{\partial^{2}u_{Hy}^{i}(x_{n},y_{m})}{\partial y^{2}}=\sum_{j}^{H}X_{j}^{ i}F_{j}^{\prime\prime}(y_{m})=\sum_{j}^{H}V_{j}^{j-2}\sum_{\begin{subarray}{c}i=0 \\ i+l=even\end{subarray}}^{l-2}\frac{lF_{i}(y_{m})}{b_{i}}(j^{2}-i^{2})\end{array}  \tag{15}\]  
In addition, the first-order partial derivative of the temperature field with respect to the time direction can also be obtained by the central difference approximation method. \[\frac{\partial u_{H}^{i}(x_{n},y_{m})}{\partial t}=\frac{u_{H}^{i+1}(x_{n},y_{m})- u_{H}^{i-1}(x_{n},y_{m})}{2\Delta t}+O(\Delta t^{2})\]  
Assume that \(u_{k}{}^{i}=u^{i}(x_{n},y_{m})\), \(u_{H_{k}}{}^{i}=u_{H}{}^{i}(x_{n},y_{m})\), \(S_{k}^{i}=S^{i}(x_{n},y_{m},i\Delta t)\), \(S_{1k}^{i}=S_{1}(x_{n},y_{m})\), \(S_{2k}^{i}=S_{2}(x_{n},y_{m})\). Here, \(1\leq k=n+1+m(H+1)\leq G=(H+1)^{2}\), \(0\leq n,m\leq H\), \(0\leq i\leq I=T/\Delta t\).