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

Authors: Jun Liu; Yan Wang
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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 Section2, 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(2.1)∂u∂t=p(x,t)∂2u∂x2+q(x,t)∂u∂x+r(x,t)u+f(x,t),a<x<b,0<t<T,
subjecting to the initial condition
(2.2)u(x,0)=γ(x),a≤x≤b
and the boundary conditions
(2.3)u(a,t)=α(t),u(b,t)=β(t),0≤t≤T.
For simplicity, we consider in this paper the homogeneous Dirichlet boundary conditions as follows:
(2.4)u(a,t)=0,u(b,t)=0,0≤t≤T,
where p,q,r,f,andγ are given functions, [a,b] is the space domain, [0,T] is the time interval, and the function u(x,t) is to be computed. Here, we assume that p(x,t)>0.For convenience, we denote the spatial differential operatorℒ by
(2.5)Lu=p(x,t)∂2u∂x2+q(x,t)∂u∂x+r(x,t)u.
Then (2.1) is equivalent to the equation
(2.6)∂u∂t=Lu+f(x,t),a<x<b,0<t<T.


### 2.1. QSC for Parabolic PDEs

Following [13, 14], we consider a uniform partition of the space domain
(2.7)Δ={a=x0<x1<⋯<xJ=b}
with mesh size Δx=(b-a)/J. Denote P2([xj,xj+1]) as the set of the quadratic polynomials on [xj,xj+1], and let
(2.8)PΔ,12={v∈C1([a,b])∣v∈P2([xj,xj+1]),j=0,1,…,J-1},
where C1([a,b]) is the set of the functions with continuous first-order derivatives on [a,b].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 spacePΔ,12 is J+2. Furthermore, let {ϕ0,ϕ1,…,ϕJ+1} be a set of piecewise polynomial basis functions of the space PΔ,12.The approximate solution in spacePΔ,12 to system (2.1), (2.2), and (2.4) can be written as
(2.9)uΔ(x,t)=∑j=0J+1‍cj(t)ϕj(x),
where cj(t),j=0,1,…,J+1, are degrees of freedom. For any fixed value of t, J+2 relations are needed to specify the approximate solution uΔ(x,t). Obviously, two relations of them can be obtained from the boundary conditions. Therefore, we have to choose other J points on (a,b), from which the J relations can be found. The J points together with the two boundary points are called collocation points, which can be denoted by
(2.10){τ0=a;τj=xj-1+xj2,1≤j≤J;τJ+1=b}.
With the collocation points, we can obtain the following relations:
(2.11)∂uΔ(τj,t)∂t=LuΔ(τj,t)+f(τj,t),uΔ(a,t)=uΔ(b,t)=0,
where uΔ(τj,0) is the value at the point τj on the interpolation of γ in PΔ,12. For some specially chosen basis functions for PΔ,12, the relations (2.11) have simple forms.Let(2.12)ϕ(x)=12{x2,0≤x≤1,-2(x-1)2+2(x-1)+1,1≤x≤2,(3-x)2,2≤x≤3,0,elsewhere.
We choose a set of quadratic B-splines
(2.13)ϕj(x)=ϕ(x-aΔx-j+2),j=0,1,…,J+1
as the basis functions of the space PΔ,12. Following [13], we denote by P~Δ,12 the space of quadratic splines satisfying homogeneous Dirichlet boundary conditions. The dimension of the space P~Δ,12 is J, and a set of its basis functions is
(2.14){ϕ~1=ϕ1-ϕ0;ϕ~j=ϕj,j=2,…,J-1;ϕ~J=ϕJ-ϕJ+1}.
Then the quadratic spline approximate solution of system (2.1), (2.2), and (2.4) can be written as
(2.15)uΔ(x,t)=∑j=1J‍cj(t)ϕ~j(x).
The values of the quadratic spline basis functions and their derivatives at the collocation points are
(2.16)ϕ~j(x)={18,x=τj-1,34,x=τj,18,x=τj+1,ϕ~j′(x)={12Δx,x=τj-1,0,x=τj,-12Δx,x=τj+1,ϕ~j′′(x)={1Δx2,x=τj-1,-2Δx2,x=τj,1Δx2,x=τj+1,
respectively.Denotediag{} as a diagonal matrix with the diagonal elements listed in the brackets. Let
(2.17)Dp(t)=diag⁡{p(τ1,t),…,p(τJ,t)},Dq(t)=diag⁡{q(τ1,t),…,q(τJ,t)},Dr(t)=diag⁡{r(τ1,t),…,r(τJ,t)},f→(t)=(f(τ1,t),…,f(τJ,t))T.
The relations (2.11) lead to the following system of ODEs:
(2.18)Q0dc(t)dt=1Δx2[Dp(t)Q2+Δx2Dq(t)Q1+Δx2Dr(t)Q0]c(t)+f→(t),t∈[0,T],c(0)=c0,
where c(t)=[c1(t),…,cJ(t)]T,
(2.19)Q0=18(510161⋱⋱⋱161015)J×J,Q1=(110-101⋱⋱⋱-1010-1-1)J×J,Q2=(-3101-21⋱⋱⋱1-2101-3)J×J.
The vector c0∈RJ satisfies Q0c0=γ→, where γ→ is the interpolation of γ(x) 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:
(2.20)Q0dc(t)dt=1Δx2[Dp(t)(Q2+124Qxx)+Δx2Dq(t)(Q1-124Qx)+Δx2Dr(t)Q0]c(t)+f→(t),c(0)=c0,
where
(2.21)Qxx=(-1116-146-10⋯0-56-4100⋯01-46-410⋯0⋱⋱⋱⋱⋱0⋯01-46-410⋯001-46-50⋯0-16-1416-11)J×J,Qx=(7-2-44-10⋯030-2100⋯0-120-210⋯0⋱⋱⋱⋱⋱0⋯0-120-210⋯00-120-30⋯01-442-7)J×J.By the system (2.20) from the optimal QSC, the discretization error 𝒪(Δx4) 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.22)L(t)=Dp(t)Q2+Δx2Dq(t)Q1+Δx2Dr(t)Q0,P(t)=124Dp(t)Qxx-124Δx2Dq(t)Qx,
and system (2.20) can be rewritten as
(2.23)Q0dc(t)dt=1Δx2(L(t)+P(t))c(t)+f→(t),t∈[0,T],c(0)=c0.
We employ a kind of two-stage implicit Runge-Kutta methods for system (2.23), with the following scheme:
(2.24)Cn+1=Cn+12(K1+K2),Q0K1=ΔtΔx2(L(tn+ω1Δt)+P(tn+ω1Δt))(Cn+14K1+ω3K2)+Δtf→(tn+ω1Δt),Q0K2=ΔtΔx2(L(tn+ω2Δt)+P(tn+ω2Δt))(Cn+ω4K1+14K2)+Δtf→(tn+ω2Δt),
where ω1=1/2+3/6,ω2=1/2-3/6,ω3=1/4+3/6,ω4=1/4-3/6,Δt=T/N, and Cn is an approximation to c(nΔt). 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{Cn}n=1N from scheme (2.24), we can obtain the approximate solution of system (2.1), (2.2), and (2.4) by {Un|Un=Q0Cn}n=1N. 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 𝒪(Δt4+Δx4).


## 2.1. QSC for Parabolic PDEs

Following [13, 14], we consider a uniform partition of the space domain
(2.7)Δ={a=x0<x1<⋯<xJ=b}
with mesh size Δx=(b-a)/J. Denote P2([xj,xj+1]) as the set of the quadratic polynomials on [xj,xj+1], and let
(2.8)PΔ,12={v∈C1([a,b])∣v∈P2([xj,xj+1]),j=0,1,…,J-1},
where C1([a,b]) is the set of the functions with continuous first-order derivatives on [a,b].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 spacePΔ,12 is J+2. Furthermore, let {ϕ0,ϕ1,…,ϕJ+1} be a set of piecewise polynomial basis functions of the space PΔ,12.The approximate solution in spacePΔ,12 to system (2.1), (2.2), and (2.4) can be written as
(2.9)uΔ(x,t)=∑j=0J+1‍cj(t)ϕj(x),
where cj(t),j=0,1,…,J+1, are degrees of freedom. For any fixed value of t, J+2 relations are needed to specify the approximate solution uΔ(x,t). Obviously, two relations of them can be obtained from the boundary conditions. Therefore, we have to choose other J points on (a,b), from which the J relations can be found. The J points together with the two boundary points are called collocation points, which can be denoted by
(2.10){τ0=a;τj=xj-1+xj2,1≤j≤J;τJ+1=b}.
With the collocation points, we can obtain the following relations:
(2.11)∂uΔ(τj,t)∂t=LuΔ(τj,t)+f(τj,t),uΔ(a,t)=uΔ(b,t)=0,
where uΔ(τj,0) is the value at the point τj on the interpolation of γ in PΔ,12. For some specially chosen basis functions for PΔ,12, the relations (2.11) have simple forms.Let(2.12)ϕ(x)=12{x2,0≤x≤1,-2(x-1)2+2(x-1)+1,1≤x≤2,(3-x)2,2≤x≤3,0,elsewhere.
We choose a set of quadratic B-splines
(2.13)ϕj(x)=ϕ(x-aΔx-j+2),j=0,1,…,J+1
as the basis functions of the space PΔ,12. Following [13], we denote by P~Δ,12 the space of quadratic splines satisfying homogeneous Dirichlet boundary conditions. The dimension of the space P~Δ,12 is J, and a set of its basis functions is
(2.14){ϕ~1=ϕ1-ϕ0;ϕ~j=ϕj,j=2,…,J-1;ϕ~J=ϕJ-ϕJ+1}.
Then the quadratic spline approximate solution of system (2.1), (2.2), and (2.4) can be written as
(2.15)uΔ(x,t)=∑j=1J‍cj(t)ϕ~j(x).
The values of the quadratic spline basis functions and their derivatives at the collocation points are
(2.16)ϕ~j(x)={18,x=τj-1,34,x=τj,18,x=τj+1,ϕ~j′(x)={12Δx,x=τj-1,0,x=τj,-12Δx,x=τj+1,ϕ~j′′(x)={1Δx2,x=τj-1,-2Δx2,x=τj,1Δx2,x=τj+1,
respectively.Denotediag{} as a diagonal matrix with the diagonal elements listed in the brackets. Let
(2.17)Dp(t)=diag⁡{p(τ1,t),…,p(τJ,t)},Dq(t)=diag⁡{q(τ1,t),…,q(τJ,t)},Dr(t)=diag⁡{r(τ1,t),…,r(τJ,t)},f→(t)=(f(τ1,t),…,f(τJ,t))T.
The relations (2.11) lead to the following system of ODEs:
(2.18)Q0dc(t)dt=1Δx2[Dp(t)Q2+Δx2Dq(t)Q1+Δx2Dr(t)Q0]c(t)+f→(t),t∈[0,T],c(0)=c0,
where c(t)=[c1(t),…,cJ(t)]T,
(2.19)Q0=18(510161⋱⋱⋱161015)J×J,Q1=(110-101⋱⋱⋱-1010-1-1)J×J,Q2=(-3101-21⋱⋱⋱1-2101-3)J×J.
The vector c0∈RJ satisfies Q0c0=γ→, where γ→ is the interpolation of γ(x) 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:
(2.20)Q0dc(t)dt=1Δx2[Dp(t)(Q2+124Qxx)+Δx2Dq(t)(Q1-124Qx)+Δx2Dr(t)Q0]c(t)+f→(t),c(0)=c0,
where
(2.21)Qxx=(-1116-146-10⋯0-56-4100⋯01-46-410⋯0⋱⋱⋱⋱⋱0⋯01-46-410⋯001-46-50⋯0-16-1416-11)J×J,Qx=(7-2-44-10⋯030-2100⋯0-120-210⋯0⋱⋱⋱⋱⋱0⋯0-120-210⋯00-120-30⋯01-442-7)J×J.By the system (2.20) from the optimal QSC, the discretization error 𝒪(Δx4) 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.22)L(t)=Dp(t)Q2+Δx2Dq(t)Q1+Δx2Dr(t)Q0,P(t)=124Dp(t)Qxx-124Δx2Dq(t)Qx,
and system (2.20) can be rewritten as
(2.23)Q0dc(t)dt=1Δx2(L(t)+P(t))c(t)+f→(t),t∈[0,T],c(0)=c0.
We employ a kind of two-stage implicit Runge-Kutta methods for system (2.23), with the following scheme:
(2.24)Cn+1=Cn+12(K1+K2),Q0K1=ΔtΔx2(L(tn+ω1Δt)+P(tn+ω1Δt))(Cn+14K1+ω3K2)+Δtf→(tn+ω1Δt),Q0K2=ΔtΔx2(L(tn+ω2Δt)+P(tn+ω2Δt))(Cn+ω4K1+14K2)+Δtf→(tn+ω2Δt),
where ω1=1/2+3/6,ω2=1/2-3/6,ω3=1/4+3/6,ω4=1/4-3/6,Δt=T/N, and Cn is an approximation to c(nΔt). 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{Cn}n=1N from scheme (2.24), we can obtain the approximate solution of system (2.1), (2.2), and (2.4) by {Un|Un=Q0Cn}n=1N. 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 𝒪(Δt4+Δx4).


## 3. Stability of QSC-TG

In this section, we consider the stability of QSC-TG for a model problem as follows:(3.1)∂u∂t=p∂2u∂x2,a<x<b,0<t<T,u(x,0)=γ(x),a≤x≤b,u(a,t)=0,u(b,t)=0,0≤t≤T,
where p is a positive constant, γ is a given function, [a,b] is the space domain, [0,T] is the time interval, and u(x,t) is the unknown function.For system (3.1), the optimal QSC gives rise to the following collocation equation:
(3.2)Q0dc(t)dt=p1Δx2(Q2+124Qxx)c(t),t∈[0,T],Q0c(0)=γ→,
where γ→ is the interpolation of γ(x) at the collocation points. The two-stage Gauss method for system (3.2) is
(3.3)Cn+1=Cn+12(K1+K2),Q0K1=pΔtΔx2(Q2+124Qxx)(Cn+14K1+ω3K2),Q0K2=pΔtΔx2(Q2+124Qxx)(Cn+ω4K1+14K2),
where Q0, Q2, and Qxx are defined in (2.20). Denote σ=p(Δt/Δx2) and Q¯=σ(Q2+(1/24)Qxx). Substitute K1 and K2 of (3.3) into the first equation of (3.3), then we have
(3.4)Cn+1=Cn+12(II)(Q0-14Q¯-ω3Q¯-ω4Q¯Q0-14Q¯)-1(Q¯00Q¯)(CnCn)≜QCn,
where the matrix Q is of size J×J, which can be regarded as the iteration matrix for the scheme described by (3.4). We denote by Qi the ith power of Q and follow the stability analysis presented in [14]. The stability of (3.4) is guaranteed if limΔx→0(maxi=0,1,…,T/Δt∥Qi∥∞) is bounded independently of Δx.For (3.4), we consider the quantities ∥Qi∥∞, with p=1, Δx=(b-a)/J, and Δt=σΔx2, for several values of σ and J. Figure 1 shows how ∥Qi∥∞ behaves as i increases, with σ=20, for several values of J. It can be observed that the quantities ∥Qi∥∞ are bounded by a constant which is independent of J, and thus limΔx→0(maxi=0,1,…,T/Δt∥Qi∥∞) is bounded independently of Δx. Therefore, the scheme of QSC-TG for the model problem (3.1) is stable without any restriction on the time step size.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.
(a)
(b)
(c)
(d)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
(3.5)A¯c¯i+1=B¯c¯i,
with the matrices
(3.6)A¯=Q0-p2ΔtΔx2Q2,B¯=Q0+p2ΔtΔx2Q2+p24ΔtΔx2Q¯xx,
where
(3.7)Q¯xx=(0000-56-411-46-41⋱⋱⋱⋱⋱⋱⋱⋱⋱⋱1-46-411-46-50000).
Let Q¯=A¯-1B¯, 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 ∥Qi∥∞ and ∥Q¯i∥∞, with the same values of σ and J, are also shown in Figure 1. It can be seen clearly that, for fixed value of J, the upper bound of the quantities ∥Qi∥∞ is smaller than the upper bound of the quantities ∥Q¯i∥∞.In Figure2, we compare the spectral radii of the iteration matrix Q with the spectral radii of the iteration matrix Q¯, versus σ, for different values of J. We can observe that the spectral radii of the matrix Q are smaller than those of Q¯ for all values of σ. It seems that the QSC-TG algorithm is better than QSC-CN0 as far as stability is concerned.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).
(a)
(b)
(c)
(d)


## 4. Numerical Experiments

We compute a linear parabolic PDE as follows:(4.1)∂∂tu(x,t)=p∂2∂x2u(x,t)+(2pπ2-1)e-t/2sin⁡(πx),u(0,t)=u(1,t)=0,0≤t≤5,u(x,0)=2sin⁡(πx),0≤x≤1,
where p is a constant. The exact solution of system (4.1) is
(4.2)u(x,t)=2e-t/2sin⁡(πx),0≤x≤1,0≤t≤5.To perform the QSC-TG algorithms, we employ a uniform partition(4.3){0=x0<x1<⋯<xJ=1}
for the space domain [0,1] with mesh size Δx=1/J and choose the collocation points
(4.4){τj=(j-12)Δx,j=1,2,…,J}.
The resulting collocation equation can be written as
(4.5)Q0dc(t)dt=pΔx2(Q2+124Qxx)c(t)+(2pπ2-1)e-t/2g→,t∈[0,5],c(0)=2Q0g→,
where the vector g→=(sin(πτ1),sin(πτ2),…,sin(πτJ))T and the matrices Q0, Q2, and Qxx, of size J×J, are the same as the matrices in system (2.20).By the QSC-TG algorithm, we can obtain an approximate solutionuΔ to system (4.1). The resulting error is measured by
(4.6)max⁡1≤j≤J,0≤i≤N|uΔ(j,i)-u(τj,ti)|,
where uΔ(j,i) denotes the (j,i)th entry of uΔ, and u(τj,ti) is the true solution at (τj,ti).Case 1 (whenp=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Δt=1/1024, which is small enough to ensure the errors are introduced by QSC, the observed errors with different values of J and the corresponding orders of accuracy are shown in Table 1. Similarly, if we choose J=1024, which is big enough to ensure the errors are introduced by the time integration, the observed errors with different values of Δt and the corresponding orders of accuracy are shown in Table 2. The results in Tables 1 and 2 confirm the discretization errors of 𝒪(Δt4+Δx4) 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Δt/Δx2=100 behaves worse than that with Δt/Δx2=98. In fact, the approximate solutions by QSC-CN0 contain spurious oscillations if the value of Δt/Δx2 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 p=1, for several values of J, where Δt=1/1024.
J
Error
Order
J
Error
Order
8
8.62
e
-
005
—
64
6.91
e
-
009
4.05
16
2.25
e
-
006
5.26
128
4.30
e
-
010
4.01
32
1.14
e
-
007
4.30
256
2.76
e
-
011
3.96Table 2
Observed errors and orders of accuracy by QSC-TG for system (4.1) with p=1, for several values of time steps, where Δx=1/1024.
Δ
t
Error
Order
Δ
t
Error
Order
1/4
1.84
e
-
005
—
1/32
4.58
e
-
009
3.99
1/8
1.16
e
-
006
3.99
1/64
3.29
e
-
010
3.80
1/16
7.26
e
-
008
4.00
1/128
6.49
e
-
011
2.34Table 3
Observed errors and time cost by QSC-CN0 and QSC-TG for system (4.1) with p=1, for several values of time steps and mesh sizes.
Δ
x
QSC-CN0
QSC-TG
Δ
t
=
98
Δ
x
2
Δ
t
=
100
Δ
x
2
Δ
t
=
Δ
x
Error
Time (s)
Error
Time (s)
Error
Time (s)
1
/
64
1.04
e
-
008
0.025
2.44
e
-
008
0.025
7.19
e
-
009
0.61
1
/
128
3.39
e
-
010
0.34
1.38
e
-
009
0.33
4.48
e
-
010
0.80
1
/
256
1.09
e
-
011
10.25
8.25
e
-
011
9.87
2.87
e
-
011
4.07
1
/
512
3.66
e
-
012
316.13
8.49
e
-
012
312.41
4.74
e
-
012
21.59Case 2 (whenp=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 p=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 p=0.1, for several values of time steps and mesh sizes.
Δ
x
QSC-CN0
QSC-TG
Δ
t
=
8
Δ
x
2
Δ
t
=
20
Δ
x
2
Δ
t
=
2.5
Δ
x
Error
Time (s)
Error
Time (s)
Error
Time (s)
1
/
64
1.02
e
-
008
1.55
1.27
e
-
007
0.27
4.54
e
-
009
0.075
1
/
128
7.17
e
-
010
45.92
7.83
e
-
009
7.54
2.83
e
-
010
0.26
1
/
256
4.75
e
-
011
1494.55
4.88
e
-
010
236.46
1.81
e
-
011
1.26
1
/
512
—
—
—
—
2.91
e
-
012
6.02


## 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.



---

*Source: 497365-2012-07-31.xml*