A Numerical Method of High Accuracy for Linear Parabolic Partial Differential Equations

Abstract

We report a new numerical algorithm for solving one-dimensional linear parabolic partial differential equations (PDEs). The algorithm employs optimal quadratic spline collocation (QSC) for the space discretization and two-stage Gauss method for the time discretization. The new algorithm results in errors of fourth order at the gridpoints of both the space partition and the time partition, and large time steps are allowed to save computational cost. The stability of the new algorithm is analyzed for a model problem. Numerical experiments are carried out to confirm the theoretical order of accuracy and demonstrate the effectiveness of the new algorithm.

1. Introduction

Quadratic spline collocation (QSC) is a kind of numerical methods for solving systems of differential equations, which gives rise to an approximate solution in the quadratic spline space. An interesting property of the QSC method is that the optimal order of convergence can be obtained by adding appropriate perturbations to the spatial differential operators. Such an idea is used in the smooth cubic spline collocation for a special linear initial value problem [1]. Then a modified version of the cubic spline collocation is formed in [2] for more general problems. On this basis, the cubic B-spline scaling function has been studied and applied widely for a variety of problems [3–8]. The optimal QSC method is applied for the solution of boundary value problems (BVPs) in [9–11] and extended to solve elliptic partial differential equations (PDEs) in [12], as well as parabolic PDEs in [13, 14].

In [13, 14], the optimal QSC is combined with the Crank-Nicolson (CN) method to formulate several revised one-step QSC-CN algorithms for linear parabolic PDEs. These algorithms give fourth-order accuracy at the midpoints and gridpoints of the space partition and second-order accuracy at the gridpoints of the time partition, while they require solving a tridiagonal linear system at each time step. However, the time step of the QSC-CN algorithms should be small enough to achieve the overall high performance. In addition, the oscillations are likely to be introduced by the application of the CN method. In fact, we can employ the optimal QSC method directly for the parabolic PDE and use high-order numerical methods for the resulting collocation equations, such as the two-stage Gauss method [15, 16].

Based on QSC and the two-stage Gauss method, we can propose a QSC-TG algorithm in this paper for linear parabolic PDEs. The new algorithm gives high-order accuracy at the gridpoints of both the space partition and the time partition. Because large time steps are allowed by the two-stage Gauss method, we do not need too much computations to reach a comparable order of accuracy. Moreover, the QSC-TG algorithm has competitive stability properties superior to the QSC-CN algorithms presented in [13, 14] and immunes to oscillations.

The remainder of this paper is organized as follows. In Section 2, we formulate the QSC-TG algorithm for one-dimensional linear parabolic PDEs and present the corresponding order of accuracy. In Section 3, we analyze numerically the stability properties of the new algorithm. Numerical experiments are given in Section 4 to illustrate the effectiveness of the new algorithm.

2. The QSC-TG Algorithm for Linear Parabolic PDEs

We consider the numerical solution of the one-dimensional linear parabolic PDE 𝜕𝑢()𝜕𝜕𝑡=𝑝𝑥,𝑡²𝑢𝜕𝑥²+𝑞(𝑥,𝑡)𝜕𝑢()()𝜕𝑥+𝑟𝑥,𝑡𝑢+𝑓𝑥,𝑡,𝑎<𝑥<𝑏,0<𝑡<𝑇, (2.1) subjecting to the initial condition 𝑢(𝑥,0)=𝛾(𝑥),𝑎≤𝑥≤𝑏 (2.2) and the boundary conditions 𝑢(𝑎,𝑡)=𝛼(𝑡),𝑢(𝑏,𝑡)=𝛽(𝑡),0≤𝑡≤𝑇. (2.3) For simplicity, we consider in this paper the homogeneous Dirichlet boundary conditions as follows: 𝑢(𝑎,𝑡)=0,𝑢(𝑏,𝑡)=0,0≤𝑡≤𝑇, (2.4) where 𝑝,𝑞,𝑟,𝑓,and𝛾 are given functions, [𝑎,𝑏] is the space domain, [0,𝑇] is the time interval, and the function 𝑢(𝑥,𝑡) is to be computed. Here, we assume that 𝑝(𝑥,𝑡)>0.

For convenience, we denote the spatial differential operator ℒ by ()𝜕ℒ𝑢=𝑝𝑥,𝑡²𝑢𝜕𝑥²+𝑞(𝑥,𝑡)𝜕𝑢()𝜕𝑥+𝑟𝑥,𝑡𝑢. (2.5) Then (2.1) is equivalent to the equation 𝜕𝑢()𝜕𝑡=ℒ𝑢+𝑓𝑥,𝑡,𝑎<𝑥<𝑏,0<𝑡<𝑇. (2.6)

2.1. QSC for Parabolic PDEs

Following [13, 14], we consider a uniform partition of the space domain Δ=𝑎=𝑥₀<𝑥₁<⋯<𝑥_𝐽=𝑏 (2.7) with mesh size Δ𝑥=(𝑏−𝑎)/𝐽. Denote 𝐏²([𝑥_𝑗,𝑥_𝑗+1]) as the set of the quadratic polynomials on [𝑥_𝑗,𝑥_𝑗+1], and let 𝐏²_Δ,1=𝑣∈𝐂¹([𝑎,𝑏])∣𝑣∈𝐏²𝑥_𝑗,𝑥_𝑗+1,𝑗=0,1,…,𝐽−1, (2.8) where 𝐂¹([𝑎,𝑏]) is the set of the functions with continuous first-order derivatives on [𝑎,𝑏].

Because of the three coefficients to be determined in each quadratic polynomial and the continuity of the first-order derivative at the gridpoints, the dimension of the space 𝐏²_Δ,1 is 𝐽+2. Furthermore, let {𝜙₀,𝜙₁,…,𝜙_𝐽+1} be a set of piecewise polynomial basis functions of the space 𝐏²_Δ,1.

The approximate solution in space 𝐏²_Δ,1 to system (2.1), (2.2), and (2.4) can be written as 𝑢_Δ(𝑥,𝑡)=^𝐽+1_𝑗=0𝑐_𝑗(𝑡)𝜙_𝑗(𝑥), (2.9) where 𝑐_𝑗(𝑡),𝑗=0,1,…,𝐽+1, are degrees of freedom. For any fixed value of 𝑡, 𝐽+2 relations are needed to specify the approximate solution 𝑢_Δ(𝑥,𝑡). Obviously, two relations of them can be obtained from the boundary conditions. Therefore, we have to choose other 𝐽 points on (𝑎,𝑏), from which the 𝐽 relations can be found. The 𝐽 points together with the two boundary points are called collocation points, which can be denoted by 𝜏₀=𝑎;𝜏_𝑗=𝑥_𝑗−1+𝑥_𝑗2,1≤𝑗≤𝐽;𝜏_𝐽+1=𝑏. (2.10) With the collocation points, we can obtain the following relations: 𝜕𝑢_Δ𝜏_𝑗,𝑡𝜕𝑡=ℒ𝑢_Δ𝜏_𝑗,𝑡+𝑓𝜏_𝑗,𝑢,𝑡_Δ(𝑎,𝑡)=𝑢_Δ(𝑏,𝑡)=0, (2.11) where 𝑢_Δ(𝜏_𝑗,0) is the value at the point 𝜏_𝑗 on the interpolation of 𝛾 in 𝐏²_Δ,1. For some specially chosen basis functions for 𝐏²_Δ,1, the relations (2.11) have simple forms.

Let 𝜙(𝑥)=12𝑥²(),0≤𝑥≤1,−2𝑥−1²()()+2𝑥−1+1,1≤𝑥≤2,3−𝑥²,2≤𝑥≤3,0,elsewhere. (2.12) We choose a set of quadratic B-splines 𝜙_𝑗(𝑥)=𝜙𝑥−𝑎Δ𝑥−𝑗+2,𝑗=0,1,…,𝐽+1 (2.13) as the basis functions of the space 𝐏²_Δ,1. Following [13], we denote by 𝐏²_Δ,1 the space of quadratic splines satisfying homogeneous Dirichlet boundary conditions. The dimension of the space 𝐏²_Δ,1 is 𝐽, and a set of its basis functions is 𝜙₁=𝜙₁−𝜙₀;𝜙_𝑗=𝜙_𝑗,𝑗=2,…,𝐽−1;𝜙_𝐽=𝜙_𝐽−𝜙_𝐽+1. (2.14) Then the quadratic spline approximate solution of system (2.1), (2.2), and (2.4) can be written as 𝑢_Δ(𝑥,𝑡)=^𝐽_𝑗=1𝑐_𝑗(𝑡)𝜙_𝑗(𝑥). (2.15) The values of the quadratic spline basis functions and their derivatives at the collocation points are 𝜙_𝑗(𝑥)=18,𝑥=𝜏_𝑗−1,34,𝑥=𝜏_𝑗,18,𝑥=𝜏_𝑗+1,𝜙^′_𝑗(𝑥)=12Δ𝑥,𝑥=𝜏_𝑗−1,0,𝑥=𝜏_𝑗,−12Δ𝑥,𝑥=𝜏_𝑗+1,𝜙_𝑗^′′(𝑥)=1Δ𝑥²,𝑥=𝜏_𝑗−1,−2Δ𝑥²,𝑥=𝜏_𝑗,1Δ𝑥²,𝑥=𝜏_𝑗+1, (2.16) respectively.

Denote diag{} as a diagonal matrix with the diagonal elements listed in the brackets. Let 𝐷_𝑝(𝑡)=diag𝑝𝜏₁,𝑡,…,𝑝𝜏_𝐽,𝐷,𝑡_𝑞(𝑡)=diag𝑞𝜏₁,𝑡,…,𝑞𝜏_𝐽,𝐷,𝑡_𝑟(𝑡)=diag𝑟𝜏₁,𝑡,…,𝑟𝜏_𝐽,𝑡_→,𝑓(𝑡)=𝑓𝜏₁,𝑡,…,𝑓𝜏_𝐽,𝑡^𝖳. (2.17) The relations (2.11) lead to the following system of ODEs: 𝑄₀𝑑𝑐(𝑡)=1𝑑𝑡Δ𝑥²𝐷_𝑝(𝑡)𝑄₂+Δ𝑥2𝐷_𝑞(𝑡)𝑄₁+Δ𝑥²𝐷_𝑟(𝑡)𝑄₀𝑐(𝑡)+^→𝑓(𝑡)[],𝑐(0),𝑡∈0,𝑇=𝑐₀, (2.18) where 𝑐(𝑡)=[𝑐₁(𝑡),…,𝑐_𝐽(𝑡)]^𝖳, 𝑄₀=1851𝟎161⋱⋱⋱161𝟎15_𝐽×𝐽,𝑄₁=11𝟎−101⋱⋱⋱−101𝟎−1−1_𝐽×𝐽,𝑄₂=−31𝟎1−21⋱⋱⋱1−21𝟎1−3_𝐽×𝐽. (2.19) The vector 𝑐₀∈𝐑^𝐽 satisfies 𝑄₀𝑐₀=^→𝛾, where ^→𝛾 is the interpolation of 𝛾(𝑥) at collocation points.

The spline collocation approximate solution from system (2.18) has second-order accuracy. Similar to the way to get optimal spline collocation approximation for BVPs in [9, 13], the optimal-order approximation to the system (2.1), (2.2), and (2.4) can be obtained by the following system of ODEs with extra perturbations: 𝑄₀𝑑𝑐(𝑡)=1𝑑𝑡Δ𝑥²𝐷_𝑝(𝑡)𝑄₂+1𝑄24^𝑥𝑥+Δ𝑥2𝐷_𝑞(𝑡)𝑄₁−1𝑄24^𝑥+Δ𝑥²𝐷_𝑟(𝑡)𝑄₀𝑐(𝑡)+^→𝑓(𝑡),𝑐(0)=𝑐₀, (2.20) where 𝑄_𝑥𝑥=−1116−146−10⋯0−56−4100⋯01−46−410⋯0⋱⋱⋱⋱⋱0⋯01−46−410⋯001−46−50⋯0−16−1416−11_𝐽×𝐽,𝑄_𝑥=7−2−44−10⋯030−2100⋯0−120−210⋯0⋱⋱⋱⋱⋱0⋯0−120−210⋯00−120−30⋯01−442−7_𝐽×𝐽. (2.21)

By the system (2.20) from the optimal QSC, the discretization error 𝒪(Δ𝑥⁴) at the midpoints and gridpoints of the uniform space partition can be given.

2.2. The Two-Stage Gauss Method for the Collocation Equation

Denote the matrices 𝐿(𝑡)=𝐷_𝑝(𝑡)𝑄₂+Δ𝑥2𝐷_𝑞(𝑡)𝑄₁+Δ𝑥²𝐷_𝑟(𝑡)𝑄₀,𝑃(𝑡)=1𝐷24^𝑝(𝑡)𝑄_𝑥𝑥−124Δ𝑥2𝐷_𝑞(𝑡)𝑄_𝑥, (2.22) and system (2.20) can be rewritten as 𝑄₀𝑑𝑐(𝑡)=1𝑑𝑡Δ𝑥²(𝐿(𝑡)+𝑃(𝑡))𝑐(𝑡)+^→𝑓(𝑡)[],𝑐(0),𝑡∈0,𝑇=𝑐₀. (2.23) We employ a kind of two-stage implicit Runge-Kutta methods for system (2.23), with the following scheme: 𝐶_𝑛+1=𝐶_𝑛+12𝐾₁+𝐾₂,𝑄₀𝐾₁=Δ𝑡Δ𝑥²𝐿𝑡_𝑛+𝜔₁Δ𝑡+𝑃𝑡_𝑛+𝜔₁Δ𝑡𝐶_𝑛+14𝐾₁+𝜔₃𝐾₂+Δ𝑡^→𝑓𝑡_𝑛+𝜔₁,𝑄Δ𝑡₀𝐾₂=Δ𝑡Δ𝑥²𝐿𝑡_𝑛+𝜔₂Δ𝑡+𝑃𝑡_𝑛+𝜔₂Δ𝑡𝐶_𝑛+𝜔₄𝐾₁+14𝐾₂+Δ𝑡^→𝑓𝑡_𝑛+𝜔₂Δ𝑡, (2.24) where 𝜔₁=1/2+√3/6,𝜔₂=1/2−√3/6,𝜔₃=1/4+√3/6,𝜔₄=1/4−√3/6,Δ𝑡=𝑇/𝑁, and 𝐶_𝑛 is an approximation to 𝑐(𝑛Δ𝑡). The Runge-Kutta method (2.24), which is based on Gauss-Legendre quadrature, is also called the two-stage Gauss method, with the fourth order of accuracy.

Based on the results {𝐶_𝑛}^𝑁_𝑛=1 from scheme (2.24), we can obtain the approximate solution of system (2.1), (2.2), and (2.4) by {𝑈_𝑛|𝑈_𝑛=𝑄₀𝐶_𝑛}^𝑁_𝑛=1. The hybrid algorithm by QSC and the two-stage Gauss method (2.24) is called QSC-TG algorithm in this paper, which gives rise to discretization errors of 𝒪(Δ𝑡⁴+Δ𝑥⁴).

3. Stability of QSC-TG

In this section, we consider the stability of QSC-TG for a model problem as follows: 𝜕𝑢𝜕𝜕𝑡=𝑝²𝑢𝜕𝑥²𝑢()(𝑥)𝑢()(),𝑎<𝑥<𝑏,0<𝑡<𝑇,𝑥,0=𝛾,𝑎≤𝑥≤𝑏,𝑎,𝑡=0,𝑢𝑏,𝑡=0,0≤𝑡≤𝑇, (3.1) where 𝑝 is a positive constant, 𝛾 is a given function, [𝑎,𝑏] is the space domain, [0,𝑇] is the time interval, and 𝑢(𝑥,𝑡) is the unknown function.

For system (3.1), the optimal QSC gives rise to the following collocation equation: 𝑄₀𝑑𝑐(𝑡)1𝑑𝑡=𝑝Δ𝑥²𝑄₂+1𝑄24^𝑥𝑥𝑐(𝑡)[],𝑄,𝑡∈0,𝑇₀𝑐(0)=^→𝛾, (3.2) where ^→𝛾 is the interpolation of 𝛾(𝑥) at the collocation points. The two-stage Gauss method for system (3.2) is 𝐶_𝑛+1=𝐶_𝑛+12𝐾₁+𝐾₂,𝑄₀𝐾₁=𝑝Δ𝑡Δ𝑥²𝑄₂+1𝑄24^𝑥𝑥𝐶_𝑛+14𝐾₁+𝜔₃𝐾₂,𝑄₀𝐾₂=𝑝Δ𝑡Δ𝑥²𝑄₂+1𝑄24^𝑥𝑥𝐶_𝑛+𝜔₄𝐾₁+14𝐾₂, (3.3) where 𝑄₀, 𝑄₂, and 𝑄_𝑥𝑥 are defined in (2.20). Denote 𝜎=𝑝(Δ𝑡/Δ𝑥²) and 𝑄=𝜎(𝑄₂+(1/24)𝑄_𝑥𝑥). Substitute 𝐾₁ and 𝐾₂ of (3.3) into the first equation of (3.3), then we have 𝐶_𝑛+1=𝐶_𝑛+12𝑄𝐼𝐼⁰−14𝑄−𝜔₃𝑄−𝜔₄𝑄𝑄₀−14Q⁻¹0𝑄0𝑄𝐶_𝑛𝐶_𝑛≜𝑄𝐶_𝑛, (3.4) where the matrix 𝑄 is of size 𝐽×𝐽, which can be regarded as the iteration matrix for the scheme described by (3.4). We denote by 𝑄^𝑖 the 𝑖th power of 𝑄 and follow the stability analysis presented in [14]. The stability of (3.4) is guaranteed if lim_Δ𝑥→0(max_{𝑖=0,1,…,𝑇/Δ𝑡}‖𝑄^𝑖‖_∞) is bounded independently of Δ𝑥.

For (3.4), we consider the quantities ‖𝑄^𝑖‖_∞, with 𝑝=1, Δ𝑥=(𝑏−𝑎)/𝐽, and Δ𝑡=𝜎Δ𝑥², for several values of 𝜎 and 𝐽. Figure 1 shows how ‖𝑄^𝑖‖_∞ behaves as 𝑖 increases, with 𝜎=20, for several values of 𝐽. It can be observed that the quantities ‖𝑄^𝑖‖_∞ are bounded by a constant which is independent of 𝐽, and thus lim_Δ𝑥→0(max_{𝑖=0,1,…,𝑇/Δ𝑡}‖𝑄^𝑖‖_∞) is bounded independently of Δ𝑥. Therefore, the scheme of QSC-TG for the model problem (3.1) is stable without any restriction on the time step size.

(a)

(b)

(c)

(d)

Figure 1 

The infinity norm of the powers of the iteration matrix for the scheme of QSC-TG and that for QSC-CN0 applied to the model problem (3.1); 𝜎=20.

To illustrate the advantages on the stability of QSC-TG, we recall the QSC-CN0 algorithm presented in [13, 14] for the model problem (3.1). Applying the QSC-CN0 algorithm to the model problem (3.1) gives rise to the linear equations of the form 𝐴𝑐_𝑖+1=𝐵𝑐_𝑖, (3.5) with the matrices 𝐴=𝑄₀−𝑝2Δ𝑡Δ𝑥²𝑄₂,𝐵=𝑄₀+𝑝2Δ𝑡Δ𝑥²𝑄₂+𝑝24Δ𝑡Δ𝑥²𝑄_𝑥𝑥, (3.6) where 𝑄_𝑥𝑥=000𝟎−56−411−46−41⋱⋱⋱⋱⋱⋱⋱⋱⋱⋱1−46−411−46−5𝟎000. (3.7) Let 𝑄=𝐴⁻¹𝐵, which is the iteration matrix of the scheme described by (3.5). It has been analyzed in [13, 14] that the QSC-CN0 algorithm for the model problem (3.1) is stable without any restriction on the time step size. The results by the comparison between the behaviors of ‖𝑄^𝑖‖_∞ and ‖𝑄^𝑖‖_∞, with the same values of 𝜎 and 𝐽, are also shown in Figure 1. It can be seen clearly that, for fixed value of 𝐽, the upper bound of the quantities ‖𝑄^𝑖‖_∞ is smaller than the upper bound of the quantities ‖𝑄^𝑖‖_∞.

In Figure 2, we compare the spectral radii of the iteration matrix 𝑄 with the spectral radii of the iteration matrix 𝑄, versus 𝜎, for different values of 𝐽. We can observe that the spectral radii of the matrix 𝑄 are smaller than those of 𝑄 for all values of 𝜎. It seems that the QSC-TG algorithm is better than QSC-CN0 as far as stability is concerned.

(a)

(b)

(c)

(d)

Figure 2 

The spectral radii of the iteration matrices for the scheme of QSC-TG and these for QSC-CN0 applied to the model problem (3.1).

4. Numerical Experiments

We compute a linear parabolic PDE as follows: 𝜕𝑢()𝜕𝜕𝑡𝑥,𝑡=𝑝²𝜕𝑥²𝑢(𝑥,𝑡)+2𝑝𝜋²−1𝑒^−𝑡/2(),𝑢()()𝑢()()sin𝜋𝑥0,𝑡=𝑢1,𝑡=0,0≤𝑡≤5,𝑥,0=2sin𝜋𝑥,0≤𝑥≤1, (4.1) where 𝑝 is a constant. The exact solution of system (4.1) is 𝑢(𝑥,𝑡)=2𝑒^−𝑡/2sin(𝜋𝑥),0≤𝑥≤1,0≤𝑡≤5. (4.2)

To perform the QSC-TG algorithms, we employ a uniform partition 0=𝑥₀<𝑥₁<⋯<𝑥_𝐽=1 (4.3) for the space domain [0,1] with mesh size Δ𝑥=1/𝐽 and choose the collocation points 𝜏_𝑗=1𝑗−2Δ𝑥,𝑗=1,2,…,𝐽. (4.4) The resulting collocation equation can be written as 𝑄₀𝑑𝑐(𝑡)=𝑝𝑑𝑡Δ𝑥²𝑄₂+1𝑄24^𝑥𝑥𝑐(𝑡)+2𝑝𝜋²−1𝑒^→−𝑡/2𝑔[],𝑐(0),𝑡∈0,5=2𝑄₀^→𝑔, (4.5) where the vector ^→𝑔=(sin(𝜋𝜏₁),sin(𝜋𝜏₂),…,sin(𝜋𝜏_𝐽))^T and the matrices 𝑄₀, 𝑄₂, and 𝑄_𝑥𝑥, of size 𝐽×𝐽, are the same as the matrices in system (2.20).

By the QSC-TG algorithm, we can obtain an approximate solution 𝑢_Δ to system (4.1). The resulting error is measured by max_{1≤𝑗≤𝐽,0≤𝑖≤𝑁}𝑢_Δ𝑗,𝑖−𝑢𝜏_𝑗,𝑡_𝑖, (4.6) where 𝑢_Δ(𝑗,𝑖) denotes the (𝑗,𝑖)th entry of 𝑢_Δ, and 𝑢(𝜏_𝑗,𝑡_𝑖) is the true solution at (𝜏_𝑗,𝑡_𝑖).

Case 1 (when 𝑝=1). We first investigate the contributions of the optimal quadratic spline collocation and the two-stage Gauss method, respectively. If we choose the time step Δ𝑡=1/1024, which is small enough to ensure the errors are introduced by QSC, the observed errors with different values of 𝐽 and the corresponding orders of accuracy are shown in Table 1. Similarly, if we choose 𝐽=1024, which is big enough to ensure the errors are introduced by the time integration, the observed errors with different values of Δ𝑡 and the corresponding orders of accuracy are shown in Table 2. The results in Tables 1 and 2 confirm the discretization errors of 𝒪(Δ𝑡⁴+Δ𝑥⁴) for the QSC-TG algorithm.
To show the advantages of the QSC-TG algorithm, we employ the QSC-CN0 algorithm presented in [13, 14] for comparison, which has been proved to be much efficient and unconditionally stable. In Table 3, we present the errors and the time cost (measured in seconds) by the QSC-CN0 algorithm and the QSC-TG algorithm for system (4.1). We can see that the QSC-TG algorithm needs much less running time when achieving a desired high accuracy.
Furthermore, we notice that the QSC-CN0 algorithm with Δ𝑡/Δ𝑥²=100 behaves worse than that with Δ𝑡/Δ𝑥²=98. In fact, the approximate solutions by QSC-CN0 contain spurious oscillations if the value of Δ𝑡/Δ𝑥² is large [17]. The QSC-TG algorithm has no such restriction and immunes to oscillations.

Table 1 

Observed errors and orders of accuracy by QSC-TG for system (4.1) with 𝑝=1, for several values of 𝐽, where Δ𝑡=1/1024.

Table 2 

Observed errors and orders of accuracy by QSC-TG for system (4.1) with 𝑝=1, for several values of time steps, where Δ𝑥=1/1024.

Table 3 

Observed errors and time cost by QSC-CN0 and QSC-TG for system (4.1) with 𝑝=1, for several values of time steps and mesh sizes.

Case 2 (when 𝑝=0.1). To show the advantages of QSC-TG more clearly, we implement the QSC-TG and QSC-CN0 again for system (4.1) with 𝑝=0.1. The observed errors and the running time are compared in Table 4. We can see that the QSC-TG algorithm costs much less running time than QSC-CN0, while it reaches much better accuracy.

Table 4 

Observed errors and time cost by QSC-CN0 and QSC-TG for system (4.1) with 𝑝=0.1, for several values of time steps and mesh sizes.

5. Conclusions

We have proposed a QSC-TG algorithm for solving linear one-dimensional parabolic PDEs. The space discretization is dealt with by the optimal quadratic spline collocation, and the time discretization is treated by the two-stage Gauss method. High order of accuracy both in space and time discretizations can be achieved. The QSC-TG algorithm has been confirmed numerically to be unconditional stable. The results of the numerical experiments show that the QSC-TG algorithm costs much less running time than the very efficient QSC-CN0 algorithm, which is presented in [13, 14], when solving the same system and achieving the same high accuracy.

Acknowledgments

This work was supported by the Natural Science Foundation of China (NSFC) under Grant 11071192 and the International Science and Technology Cooperation Program of China under Grant 2010DFA14700. This work was partly supported by NSFC under Grant 11101127 (Y. Wang).