A Spectral Method for Convection-Diffusion Equations

Peng Guo^1*, Qin Wang², Zhengang Zhao³

¹The Business School, Shanghai Dianji University, Shanghai, China.
²College of Transport & Communications, Shanghai Maritime University, Shanghai, China.
³Department of Fundamental Courses, Shanghai Customs College, Shanghai, China.

DOI: 10.4236/am.2022.1312061                

Abstract

In the practical problems such as nuclear waste pollution and seawater intrusion etc., many problems are reduced to solving the convection-diffusion equation, so the research of convection-diffusion equation is of great value. In this work, a spectral method is presented for solving one and two dimensional convection-diffusion equation with source term. The finite difference method is also used to solve the convection diffusion equation. The numerical experiments show that the spectral method is more efficient than other methods for solving the convection-diffusion equation.

Keywords

Convection-Diffusion Equation, Central Finite Difference Method, Upwind Difference Method, Chebyshev, Spectral Method

1. Introduction

In recent decades, with the increasing attention of human beings to environmental problems, environmental pollution problems have become an increasing object of concern. Such as groundwater pollutants transfer, diffusion problem of pollutants in the ocean; absorption of chemical substances in riverbeds, distribution of pollutants in nuclear pollution, long-range propagation of pollutants in the atmosphere, etc., all of these phenomenon can be modelling by the convection-diffusion equation. The study of convection-diffusion equation can provide theoretical support for pollution prediction and the development trend of pollutants, therefore the numerical analysis of convection-diffusion equation has important theoretical and practical value.

In general, nearly all the partial differential equations are difficult or cannot give the analytical solutions. So many scholars have proposed various numerical algorithms to cope with such problems, the traditional numerical algorithms including finite difference method, finite element method, finite volume method, especially for finite difference method, many scholars have proposed various compact schemes, i.e., using equal number of nodes with different coefficients to construct some new schemes with high precision. These numerical methods have greatly enriched the numerical solution of the convection-diffusion.

Next, we will focus on the following convection-diffusion equation,

$-\epsilon \Delta u+\left(\alpha ,\beta \right)\cdot \nabla u+\gamma u=f,x\in \Omega .$ (1)

The above Equation (1) is the convection-diffusion equation, where the second derivative i.e. $\Delta u$ is the diffusion term, the first derivative term $\nabla u$ is the convective term, f is the source term, $\epsilon$ is the coefficient of diffusion term, $\alpha$ and $\beta$ are constants.

In general, it is difficult to give an analytical solution for the convection-diffusion equation, so the numerical solution becomes a good way to cope with it. Many scholars proposed many different methods such as finite difference method, finite element method, finite volume method, etc. for solving the convective diffusion equation. Dennis and Hundson [1] proposed a 4th-order compact finite difference method for the Navier-stokes type convection-diffusion equation. Lele [2] proposed a compact finite difference format with pseudo-spectral resolution. Fu and Ma [3] proposed an upwind compact difference method. Cockburn and Shu [4] constructed the nonlinear compact format with fourth-order accuracy, etc. For some other methods and discussions of the convection-diffusion equation we can refer to [5] - [16].

Compared with the methods mentioned earlier, the spectral method has a higher accuracy in solving partial differential equations [17] [18]. In particular, for convection-dominated convective diffusion equations, the spectral method is effective in overcoming the numerical oscillation phenomenon in addition to giving numerical solutions with higher accuracy. In this paper we will consider the Chebyshev spectral method to give numerical solutions of two convective diffusion equations. For more applications of the Chebyshev spectral method and the related theory, we can refer to [17] - [23].

2. Preliminaries

Some basic contents including Chebyshev polynomials and Chebyshev points will be introduced in this part. The Chebyshev points will be used to construct the differentiation matrices, which is the key point to obtain high-precision solutions.

Definition 1. [17] When the weight function $\rho \left(x\right)=\frac{1}{\sqrt{1-x^2}}$, the Chebyshev polynomials can be expressed as follows

$T_n\left(x\right)=\cos \left(n\arccos \left(x\right)\right),$ (2)

where $\left|x\right|\le 1$.

Figure 1 shows the Chebyshev polynomials when $n=1,2,3,4$. Using the trigonometric constancy relation, we also have the following Chebyshev polynomial recurrence relation,

$T_{n+1}\left(x\right)=2x{T}_n\left(x\right)-T_{n-1}\left(x\right),n\ge 1.$ (3)

This Chebyshev polynomial recurrence relation (3) is equivalent to the explicit expression (2). The first several terms of the Chebyshev recurrence relation are as follows

$\{\begin{array}{l}{T}_0\left(x\right)=1,\hfill \\ T_1\left(x\right)=x,\hfill \\ T_2\left(x\right)=2x^2-1,\hfill \\ T_3\left(x\right)=4x^3-3x.\hfill \\ ⋮\hfill \end{array}$ (4)

Clearly, when $x_k=\cos \frac{k\pi }{n},k=0,1,\cdots ,n$, the Chebyshev polynomials equal to zero. All these points are called Chebyshev points. For some finite difference schemes, the equidistant nodes are used widely because of their simplicity and convenience. But for the equidistant nodes, with the number of nodes increases, a major disadvantage is the Runge phenomenon. Especially near the end of the interval, the Runge phenomenon is very obvious. Therefore, to avoid the Runge phenomenon, equidistant nodes are not a suitable choice, and non-equidistant nodes are a way to cope with this problem. Chebyshev points are dense near the ends of the interval and relatively sparse in the middle, so interpolation with Chebyshev nodes can better avoid the Runge phenomenon.

Figure 1. Chebyshev polynomials, n = 1, 2, 3, 4.

Property 1.

${\displaystyle {\int }_{-1}^1}\rho \left(x\right)T_m\left(x\right)T_n\left(x\right)\text{d}x=\{\begin{array}{ll}0,\hfill & m\ne n,\hfill \\ \frac{\pi }{2},\hfill & m=n=0,\hfill \\ \pi ,\hfill & m=n\ne 0.\hfill \end{array}$ (5)

where the weight function is $\rho \left(x\right)=\frac{1}{\sqrt{1-x^2}}$.

Since the studied convective-diffusion equation contains derivative terms, it is necessary to consider the Chebyshev derivative matrix to approximate the derivative terms, and next we give the specific representation of the Chebyshev derivative matrix.

$D=\left(\begin{array}{ccccc}-\frac{2N^2+1}{6}& \cdots & 2\frac{{\left(-1\right)}^i}{-1-x_j}& \cdots & \frac{1}{2}{\left(-1\right)}^N\\ ⋮& \ddots & ⋮& \frac{{\left(-1\right)}^{i+j}}{x_i-x_j}& ⋮\\ -\frac{1}{2}\frac{{\left(-1\right)}^i}{1+x_i}& \cdots & \frac{-x_j}{2\left(1-x_j^2\right)}& \cdots & -\frac{1}{2}\frac{{\left(-1\right)}^i}{1-x_i}\\ ⋮& \frac{{\left(-1\right)}^{i+j}}{x_i-x_j}& ⋮& \ddots & ⋮\\ \frac{1}{2}{\left(-1\right)}^N& \cdots & 2\frac{{\left(-1\right)}^{N+j}}{1-x_j}& \cdots & \frac{2N^2+1}{6}\end{array}\right)$ (6)

where D is the Chebyshev matrix of first order derivative. For each element of the Chebyshev first order derivative matrix, the following is shown

$\{\begin{array}{l}{D}_{0,0}=-\frac{2N^2+1}{6},\\ D_{i,i}=\frac{-x_i}{2\left(1-x_i^2\right)},i=1,2,\cdots ,N-1,\\ D_{i,j}=\frac{c_i{\left(-1\right)}^{i+j}}{c_j\left(x_i-x_j\right)},i\ne j,i,j=0,1,\cdots ,N,\\ D_{N,N}=\frac{2N^2+1}{6},\end{array}$ (7)

where

$c_i=\{\begin{array}{l}2,i=0,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}N,\\ 1,\text{else}.\end{array}$ (8)

For the second order and higher order derivatives, we approximate them by $D^k\mathrm{,}k\ge 2$ i.e. we have the following approximate representation,

$\left(\begin{array}{c}{u}^{\left(k\right)}\left(x_1\right)\\ u^{\left(k\right)}\left(x_2\right)\\ ⋮\\ u^{\left(k\right)}\left(x_n\right)\end{array}\right)=D^k\left(\begin{array}{c}u\left(x_1\right)\\ u\left(x_2\right)\\ ⋮\\ u\left(x_n\right)\end{array}\right).$ (9)

There are some other polynomials such as Legendre polynomials, Laguerre polynomials, Jacobi polynomials etc. For a specific discussion of special polynomials and their related properties, one can refer to [17] [18].

Theorem 1. [24] Suppose u is analytical, $\rho$ is sum of semimajor and semiminor axis lengths, P is the interpolation function on the Chebyshev points, Then,

$\left|P_N\left(x\right)-u\left(x\right)\right|=O\left({\rho }^{-N}\right).$ (10)

One of the advantages of spectral methods is that they have high accuracy. Theoretically speaking, Chebyshev spectral method can achieve any order of accuracy. In fact, the accuracy of the Chebyshev spectrum method depends on the smoothness of the solution of the original problem. If the solution u of the original problem is an analytic function in the complex plane, then the convergence efficiency completely depends on the analytic ellipse of u.

3. Nunmerical Algorithms for Two Types of Convection-Diffusion Equations

In the following part, we will discuss the numerical solutions of several different types of convection-diffusion equations. For each convection-diffusion equation, we will use the finite difference method and the Chebyshev spectral method to get the numerical solution, respectively.

3.1. One Dimensional Convection-Diffusion Equation

In this part we will consider the following one dimensional convection-diffusion equation,

$\{\begin{array}{l}-\epsilon u_{xx}+\alpha u_x=f\left(x\right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}x\in \Omega ,\\ {u|}_{\partial \Omega }=g.\end{array}$ (11)

where the diffusion coefficient $\epsilon$ and convection coefficient $\alpha$ are constants. $f\left(x\right)$ is the cource term.

Using the central finite difference method, the second order derivative term $u_{xx}$ has the following approximate format

${u_{xx}|}_{x_i}=\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}+O\left(h^2\right).$ (12)

For the first order derivative term $u_x$, there are many different ways to approximate it, such as the central finite difference method, the upwind finite difference method, and the modified finite difference method, etc. Finally, we will discuss several different numerical algorithms and compare the advantages and disadvantages of these algorithms.

Example 1.

$\{\begin{array}{l}-\epsilon u_{xx}+u_x=1,\\ u\left(0\right)=0,\\ u\left(1\right)=0.\end{array}$ (13)

In this case, we will give the corresponding numerical solution by finite difference method. In the interval $\left[\mathrm{0,1}\right]$ we will use $n+1$ nodes and give a equidistant discretization, i.e. each subinterval the length $h=\frac{1}{n}$. Then each point of the interval can be expressed as $x_i=\frac{i}{n},i=1,2,\cdots ,n$.

For the first order derivative term $u_x$, we use

${u_x|}_{x_i}=\frac{u_{i+1}-u_{i-1}}{2h}+O\left(h^2\right).$ (14)

Then for this problem, the following central finite difference scheme is obtained.

$-\epsilon \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}+\frac{u_{i+1}-u_{i-1}}{2h}=1,$ (15)

The trunction error is $O\left(h^2\right)$.

Let’s denote $\frac{2\epsilon }{h^2}=\alpha$, $-\frac{\epsilon }{h^2}+\frac{1}{2h}=\beta$, $-\frac{\epsilon }{h^2}-\frac{1}{2h}=\gamma$, then we can write the discrete form of convection-diffusion equation as the following expression,

$\left(\begin{array}{ccccccc}\alpha & \beta & 0& \cdots & 0& 0& 0\\ \gamma & \alpha & \beta & \cdots & 0& 0& 0\\ 0& \gamma & \alpha & \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & \gamma & \alpha & \beta \\ 0& 0& 0& \cdots & 0& \gamma & \alpha \end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ ⋮\\ u_{n-2}\\ u_{n-1}\end{array}\right)=\left(\begin{array}{c}1-\gamma u_0\\ 1\\ 1\\ ⋮\\ 1\\ 1-\beta u_n\end{array}\right)\mathrm{.}$ (16)

Figures 2-5 show the numerical and exact solutions of the convection diffusion equation for different $\epsilon$ respectively. Figure 2 and Figure 3 show that when $\epsilon >1$, the numerical solution agrees very well with the exact solution, and

Figure 2. The numerical solution and exact solution for $\epsilon =100$.

Figure 3. The numerical solution and exact solution for $\epsilon =2$.

Figure 4. The numerical solution and exact solution for $\epsilon =0.05$.

Figure 5. The numerical solution and exact solution for $\epsilon =0.01$.

there are no significant errors from the plots. Figure 4 and Figure 5 show that the numerical and exact solutions agree well with relatively small errors in the left half of the interval, but the errors become larger in the right side of the interval, especially near the right endpoint of the interval $\left[\mathrm{0,1}\right]$.

As shown in Table 1, the absolute error between the exact and numerical solutions when $\epsilon$ takes different values, i.e. $\epsilon =100,2,\frac{1}{20},\frac{1}{100}$. When the coefficient $\epsilon$ is less than 1, the numerical solution of the convective diffusion equation deviates significantly from the exact solution at the discrete points on the right side of the interval, and the numerical format has a significant oscillation at these points. When $\epsilon$ is less than 1, the smaller $\epsilon$ is, the larger the error is.

The traditional differential format produces numerical oscillations when the coefficient is less than one, especially when the coefficient is much less than one. This is mainly because, when the coefficient is much less than 1, the equation at this point becomes a convective dominance problem.

Table 1. Absolute error with different $\epsilon$.

In order to overcome numerical oscillations, scholars have proposed upwind finite difference method (16), corrected finite difference method (17), characteristic finite difference method and other methods based on the traditional finite difference method.

Upwind finite difference method,

$-\epsilon \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}+\frac{u_i-u_{i-1}}{h}=1.$ (17)

Corrected finite difference method,

$-\epsilon \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}+\left[\lambda \frac{u_i-u_{i-1}}{h}+\left(1-\lambda \right)\frac{u_{i+1}-u_{i-1}}{2h}\right]=1.$ (18)

For the corrected finite difference method, is a combination of central difference method and upwind difference method with weight $\lambda$ and $1-\lambda$ separately, where $0\le \lambda \le 1$.

Similarly, the traditional finite element method also has numerical oscillations for the convection dominant equation. Many scholars have proposed a series of correction techniques to combine the finite element method with other methods, and proposed the characteristic finite element method, characteristic hybrid element method, etc.

Figure 6 and Figure 7 show the numerical solutions of various different algorithms when $\epsilon$ takes different parameters. For the central differential format, the numerical oscillation is obvious when $\epsilon$ is smaller. The upwind differential format and the modified upwind differential format can overcome the numerical oscillation phenomenon, but the accuracy of these methods is not good enough. Therefore, an algorithm that can avoid numerical oscillations and at the same time obtain high numerical accuracy is necessary for solving the convection-diffusion equation.

Table 2 shows the absolute errors between the numerical solution and the exact solution obtained with different numerical algorithms. Error₁, Error₂, and

Figure 6. The numerical solution and exact solution for $\epsilon =\frac{1}{100}$.

Figure 7. The numerical solution and exact solution for $\epsilon =\frac{1}{200}$.

Table 2. Absolute errors for different methods.

Error₃ in the table are the absolute errors between the analytical solution and the numerical solution obtained in the central difference format, upwind difference format, and the corrected upwind difference format, respectively. For the same points as shown in Table 2, especially at the points near the right end of the interval, Errors₂ and Error₃ are significantly smaller than Error₁, which also means that there is no numerical oscillation for the corresponding numerical format.

Algorithm 1: Chebyshev spectral method for one dimentional convection-diffusion equation

For the convection-diffusion equation, we can obtain the following numerical format, using the Chebyshev spectral method,

$-\epsilon D^{2}u+Du=f,$ (19)

where D is the Chebyshev matrix, $f={\left(\mathrm{1,1,}\cdots \mathrm{,1}\right)}^{\text{T}}$. If we denote L as the coefficient matrix $-\epsilon D^2+D$, then we have

$Lu=f.$ (20)

Next, we consider the boundary conditions $u\left(0\right)=1,u\left(1\right)=0$. We need to make a small modification to the coefficient L. After replacing the first and last row of the matrix L with $\left(\mathrm{1,0,}\cdots \mathrm{,0,0}\right)\mathrm{,}\left(\mathrm{0,0,}\cdots \mathrm{,0,1}\right)$ respectively, we denote the new coefficient matrix as $\hat{L}$. For the right-hand side term f of the Equation (18), we also need to modify the first and last terms to 0, we denote the modified f as $\hat{f}$. Finally, we get the numerical format for the convection-diffusion equation with Chebyshev spectral method,

$\hat{L}u=\hat{f}.$ (21)

Figure 8 shows the numerical solutions given by the Chebyshev spectral method and several of the previously discussed finite difference methods. It is clear that the Chebyshev spectral method does not have numerical oscillations while being able to obtain higher accuracy. This is mainly because the Chebyshev points are not uniformly distributed and are more concentrated at the two ends of the interval, so there is no numerical oscillation near the interval endpoints.

In Table 3, Error₁, Error₂, Error₃ and Error_ch are the errors of cerntral difference method, upwind difference method, corrected difference method and Chebyshev spectral method respectively and 21 nodes are used for each format. From Table 3, it is clear that the numerical solution calculated with the central difference format has the largest error and the numerical solution calculated by the Chebyshev spectral method has the smallest error.

Table 4 shows the absolute errors between the numerical results calculated by the Chebyshev spectral method and exact solutions when different parameters are taken. It is obvious that when the coefficients are fixed, the larger the number of nodes, the smaller the error, meanwhile, when the coefficients are very small there is no numerical oscillation. When the node N is large enough, the absolute error values in the table are already very small and exceed the minimum value of the system, so some of the data values in the table are 0.

Figure 8. Numerical results with different methods for $\epsilon =\frac{1}{100}$.

Table 3. Maximum error of different methods. ( $\epsilon =\frac{1}{100}$ ).

Table 4. Maximum errors of Chebyshev spectral method with different nodes N and $\epsilon$.

3.2. Two Dimensional Convection-Diffusion Equation

In this part we will focus on the two dimensional convection-diffusion equation. The specific equation format is as follows,

$-\epsilon \left(u_{xx}+u_{yy}\right)+\alpha \left(x\mathrm{,}y\right)u_x+\beta \left(x\mathrm{,}y\right)u_y+u=f\left(x\mathrm{,}y\right)\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>}\left(x\mathrm{,}y\right)\in \Omega \mathrm{,}$ (22)

where the diffusion coefficient $\epsilon$ is a constant, $\alpha \left(x\mathrm{,}y\right)\mathrm{,}\beta \left(x\mathrm{,}y\right)$ are the convection coefficients, and $f\left(x\mathrm{,}y\right)$ is the source term. The convection-diffusion satisfies Dirichlet boundary condition, i.e. ${u|}_{\partial \Omega }=g\left(x\mathrm{,}y\right)$.

Both in x and y directions we use the central difference method respectively, we have

${u_{xx}|}_{x_i,y_j}=\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}+O\left(h^2\right),$ (23)

${u_{yy}|}_{x_i,y_j}=\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2}+O\left(h^2\right),$ (24)

${u_x|}_{x_i,y_j}=\frac{u_{i+1,j}-u_{i-1,j}}{2h}+O\left(h^2\right),$ (25)

${u_y|}_{x_i,y_j}=\frac{u_{i,j+1}-u_{i,j-1}}{2h}+O\left(h^2\right).$ (26)

Then we have the following numerical scheme.

$\begin{array}{l}-\epsilon \left(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}+\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2}\right)\\ +{\alpha }_{i,j}\frac{u_{i+1,j}-u_{i-1,j}}{2h}+{\beta }_{i,j}\frac{u_{i,j+1}-u_{i,j-1}}{2h}+u_{i,j}=f_{i,j}.\end{array}$ (27)

In the above numerical format ${\alpha }_{i,j}=\alpha \left(x_i,y_j\right)$, ${\beta }_{i,j}=\beta \left(x_i,y_j\right)$, $f_{i,j}=f\left(x_i,y_j\right)$. Both in x and y directions, we use the same step size h.

Example 2.

$\{\begin{array}{l}-\epsilon u_{xx}+u_{yy}+\frac{1}{1+y}{u}_x=f\left(x,y\right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\left(x,y\right)\in \left[0,1\right]\times \left[0,1\right],\\ u\left(x,0\right)={\text{e}}^{-x}+2^{-\frac{1}{\epsilon }},\\ u\left(x,1\right)={\text{e}}^{1-x}+2,\\ u\left(0,y\right)={\text{e}}^y+2^{-\frac{1}{\epsilon }}{\left(1+y\right)}^{1+\frac{1}{\epsilon }},\\ u\left(1,y\right)={\text{e}}^{y-1}+2^{-\frac{1}{\epsilon }}{\left(1+y\right)}^{1+\frac{1}{\epsilon }}.\end{array}$ (28)

On the right hand side of the above equation $f\left(x,y\right)=\left(-2\epsilon -\frac{1}{1+y}\right){\text{e}}^{y-x}-2^{-\frac{1}{\epsilon }}\left(1+\frac{1}{\epsilon }\right){\left(1+y\right)}^{\frac{1}{\epsilon }-1}$.

The numerical solution of Example 2 using the central finite difference algorithm (26) is shown in Figure 9, and the analytic solution of the corresponding problem is shown in Figure 10. It is clearly that when $\epsilon$ is small, there is a significant difference between the numerical solution and the exact solution at the points near the end of the definition domain. Figure 12 shows the absolute errors of the numerical and exact solutions at the nodes, which are smaller when x and y are near the left endpoint 0 and have become quite large when x and y are near the right endpoint 1.

Figure 11 and Figure 12 show the errors between the numerical solution given by the central difference method and the exact solution when the diffusion coefficient $\epsilon$ takes different values. When the diffusion coefficient $\epsilon$ is large, the numerical solution given by the central difference algorithm agrees better with the exact solution, but when the diffusion coefficient is small, the error between the numerical solution and the exact solution is larger.

Figure 9. Numerical solution $\hat{u}\left(x\mathrm{,}y\right)$ with $\epsilon =\frac{1}{100},N_x=N_y=21$.

Figure 10. Exact solution $u\left(x\mathrm{,}y\right)$ with $\epsilon =\frac{1}{100}$.

Figure 11. The absolute error of central difference method for $\epsilon =10,N_x=N_y=21$.

Figure 12. The absolute error of central difference method for $\epsilon =\frac{1}{100},N_x=N_y=21$.

Table 5 shows the maximum error in the definition domain when the nodes and diffusion coefficients take different values. From Table 5 we can clearly see that when the diffusion coefficient $\epsilon \ge 1$, the maximum error is relatively small regardless of the number of nodes, but when the diffusion coefficient $\epsilon$ is small, as shown in the table when $\epsilon =\frac{1}{100}$ or $\epsilon =\frac{1}{200}$, the maximum error has reached an unacceptable level even if the number of nodes is large.

Algorithm 2: Chebyshev spectral method for two dimenditional convection- diffusion equation.

For the Chebyshev derivative matrix D, as discussed in the previous section, we will not expand the discussion in detail here. If we arrange the solution u into a column vector along the y direction. Then $u_{xx}+u_{yy}$ can be expressed as $\left(D^2\otimes I+I\otimes D^2\right)u$, where $\otimes$ is the Kronecker product. The first order term $\frac{1}{1+y}{u}_x$ can be writen as $diag\left(\frac{1}{1+y_0},\cdots ,\frac{1}{1+y_n}\right)\left(D\otimes I\right)U$. Similarly, for the right end term f of the equation we rearrange it and denote it as F. Then for two dimenditional convection-diffusion Equation (27) we have the following Chebyshev spectral numerical format.

$\left(-\epsilon \left(D^2\otimes I+I\otimes D^2\right)+diag\left(\frac{1}{1+y_0}\mathrm{,}\cdots \mathrm{,}\frac{1}{1+y_n}\right)\left(D\otimes I\right)\right)U=F\mathrm{.}$ (29)

If we denote $-\epsilon \left(D^2\otimes I+I\otimes D^2\right)+diag\left(\frac{1}{1+y_0}\mathrm{,}\cdots \mathrm{,}\frac{1}{1+y_n}\right)\left(D\otimes I\right)$ as L, then the above format can be simply written as

$LU=F.$ (30)

However, this Format (28) or (29) only considers the internal nodes, and the boundary conditions are not yet considered. If we take the boundary conditions into account, then we need to make a small modification to (28) or (29), i.e., replace the coefficient matrix and the corresponding row of the right terminal term with the corresponding relationship at the boundary. If we denote the modified coefficient matrix and the right terminal term as $\hat{L}$ and $\hat{F}$, respectively, then we have,

$\hat{L}U=\hat{F}.$ (31)

We use the Chebyshev spectral method given above to calculate the numerical solution of Example 2. Figure 13 and Figure 14 show the errors when $\epsilon$ is taken as 10 and $\frac{1}{100}$ respectively, and 21 nodes are taken in the x and y directions.

Figure 13 shows that when $\epsilon =10$, the maximum absolute error is less than 1.5 × 10⁻⁹. When the diffusion coefficient is $\epsilon =\frac{1}{100}$, the absolute error of the numerical solution given by the spectral method becomes larger but is still less than 4 × 10⁴. This result is much better than the numerical solution given by the central difference format, and there is no numerical oscillation.

Table 6 gives the absolute errors of the numerical solutions obtained by the Chebyshev spectral method for example 2 when different parameters are taken.

Figure 13. The absolute error of Chebyshev spectral method for $\epsilon =10,N_x=N_y=21$.

Table 5. Absolute error of central difference method with different parameters.

Figure 14. The absolute error of Chebyshev spectral method for $\epsilon =\frac{1}{100},N_x=N_y=21$.

Table 6. Absolute error of Chebyshev spectral method with different parameters.

The errors are all better than the finite difference method when the diffusion coefficient is smaller, which also shows that the spectral method is better adapted than the finite difference method for solving the convective dominance equation.

4. Conclusions

In this work, we mainly focus on the one and two dimensional convection-diffusion equation. For the one dimensional convection-diffusion equation, we give the Chebyshev spectral format, in addition the central difference format, the upwind format, and the correction format are also given. Through specific numerical case, it is obvious that the finite difference format has numerical oscillations in the convection dominant case and the accuracy of the differential format is not very satisfactory. In contrast, the numerical solution given by the Chebyshev spectral method is in good agreement with the exact solution and has no numerical oscillations. For the two-dimensional convection-diffusion equation, we mainly discuss the numerical format corresponding to the two dimensional Chebyshev spectral method, and also give the central difference format. Through specific numerical results, it is obvious that the Chebyshev spectral method not only obtains higher accuracy, but also has significant advantages over the difference method in suppressing numerical oscillations.

The Chebyshev spectral method mentioned above is generally only adapted to rectangular regions, and it is difficult to extend it to general regions. Therefore, the spectral method is combined with other methods, such as the finite element method, to extend its adaptation to the general area, which is very useful when dealing with practical problem models, not only to give more accurate solutions, but also to save the computational workload significantly.

Acknowledgements

This research was funded by the Humanity and Social Science Youth foundation of Ministry of Education (No. 18YJC630120), the Applied Mathematics of Shanghai Dianji University (No. 16JCXK02) and the Key Projects of Shanghai Soft Science Research Program (No. 18692106600), the Natural Science Foundation of Shanghai (No. 19ZR1422000).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Dennis, S.C.R. and Hundson, J.D. (1989) Compact h4 Finite-Difference Approximations to Operators of Navier-Stokes Type. Journal of Computational Physics, 85, 390-416. https://doi.org/10.1016/0021-9991(89)90156-3
[2]	Lele, S.K. (1992) Compact Finite Difference Schemes with Spectral-Like Resolution. Journal of Computational Physics, 103, 16-42. https://doi.org/10.1016/0021-9991(92)90324-R
[3]	Ma, Y.W. and Fu, D.X. (1995) Super Compact Finite Difference Methods with Uniform and Nonuniform Grid System. Proceedings of 6th International Symposium on Computational Fluid Dynamics, Lake Tahoe, Nevada, 1435-1440.
[4]	Cockburn, B. and Shu, C.-W. (1994) Nonlinearly Stable Compact Schemes for Shock Calculations. SIAM Journal on Numerical Analysis, 31, 607-630. https://doi.org/10.1137/0731033
[5]	Bassi, F. and Rebay, S. (2002) Numerical Evaluation of Two Discontinuous Galerkin Methods for the Compressible Navier-Stokes Equations. International Journal for Numerical Methods in Fluids, 40, 197-207. https://doi.org/10.1002/fld.338
[6]	Cockburn, B. and Shu, C.W. (1998) The Local Discontinuous Galerkin Method for Time-Dependent Convection Diffusion System. SIAM Journal on Numerical Analysis, 35, 2440-2463. https://doi.org/10.1137/S0036142997316712
[7]	Li, L.Y. and Yin, Z. (2017) Numerical Simulation of Groundwater Pollution Problems Based on Convection Diffusion Equation. American Journal of Computational Mathematics, 7, 350-370. https://doi.org/10.4236/ajcm.2017.73025
[8]	Ghorbanzadeh, M. and Kerayechian, A. (2011) Analysis of an Il’in Scheme for a System of Singularly Perturbed Convection-Diffusion Equations. Applied Mathematics, 2, 866-873. https://doi.org/10.4236/am.2011.27116
[9]	Spotz, W.F. and Carey, G.F. (2001) Extension of High-Order Compact Schemes to Time-Dependent Problems. Numerical Methods for Partial Differential Equations, 17, 657-672. https://doi.org/10.1002/num.1032
[10]	Patidar, K.C. (2007) High Order Parameter Uniform Numerical Method for Singular Perturbation Problems. Applied Mathematics and Computation, 188, 720-733. https://doi.org/10.1016/j.amc.2006.10.040
[11]	Liu, Q. and Shi, D. (2020) New Error Analysis of a Second Order BDF Scheme for Unsteady Natural Convection Problem. Applied Numerical Mathematics, 154, 243-259. https://doi.org/10.1016/j.apnum.2020.04.004
[12]	Andargie, A. and Reddy, Y.N. (2007) Fitted Fourth-Order Tridiagonal Finite Difference Method for Singular Perturbation Problems. Applied Mathematics and Computation, 192, 90-100. https://doi.org/10.1016/j.amc.2007.02.123
[13]	Li, R., Wu, Q. and Zhu, S. (2019) Proper Orthogonal Decomposition with Supg-Stabilized Isogeometric Analysis for Reduced Order Modelling of Unsteady Convection-Dominated Convection-Diffusion-Reaction Problems. Journal of Computational Physics, 387, 280-302. https://doi.org/10.1016/j.jcp.2019.02.051
[14]	Yang, J.Q. and Zhao, B.X. (2021) Numerical Investigation of Double-Diffusive Convection in Rectangular Cavities with Different Aspect Ratio I: High-Accuracy Numerical Method. Computers & Mathematics with Applications, 94, 155-169. https://doi.org/10.1016/j.camwa.2021.05.002
[15]	Tong, F.H., Feng, X.L. and Li, Z.L. (2021) Fourth Order Compact FD Methods for Convection Diffusion Equations with Variable Coefficients. Applied Mathematics Letters, 121, Article ID: 107413. https://doi.org/10.1016/j.aml.2021.107413
[16]	Sweilam, N.H., Moharram, H.M. and Ahmed, S. (2012) On the Parallel Iterative Finite Difference Algorithm for 2-D Poisson’s Equation with MPI Cluster. 2012 8th International Conference on Informatics and Systems (INFOS), Giza, 14-16 May 2012, MM-78-MM-85.
[17]	Jie, S. and Tao, T. (2006) Spectral and High-Order Methods with Applications. Science Press, Beijing.
[18]	Jie, S., Tao, T. and Wang, L.L. (2011) Spectral Methods: Algorithms, Analysis and Applications, Springer-Verlag, Berlin.
[19]	Gu, Z. and Chen, Y. (2014) Chebyshev Spectral-Collocation Method for a Class of Weakly Singular Volterra Integral Equations with Proportional Delay. Journal of Numerical Mathematics, 22, 311-342. https://doi.org/10.1515/jnma-2014-0014
[20]	Kumar, S. (2022) Numerical Solution of Fuzzy Fractional Diffusion Equation by Chebyshev Spectral Method. Numerical Methods for Partial Differential Equations, 38, 490-508.
[21]	Rauch, R., Trummer, M.R. and Williams, J.F. (2021) A Spectral Collocation Method for Mixed Functional Differential Equations. Applied Numerical Mathematics, 161, 101-110. https://doi.org/10.1016/j.apnum.2020.10.011
[22]	Trefethen, L.N. (2000) Spectral Methods in MATLAB. SIAM, Philadelphia. https://doi.org/10.1137/1.9780898719598
[23]	Kosloff, D. (2012) Solution of the Equations of Dynamic Elasticity by a Chebychev Spectral Method. Geophysics, 55, Article 734. https://doi.org/10.1190/1.1442885
[24]	Berrut, J.P. and Klein, G. (2014) Linear Barycentric Rational Quadrature. Journal of Computational and Applied Mathematics, 259, 95-107. https://doi.org/10.1016/j.cam.2013.03.044