# Cubic B-Spline Collocation Method for One-Dimensional Heat and Advection-Diffusion Equations

Authors: Joan Goh; Ahmad Abd. Majid; Ahmad Izani Md. Ismail
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Numerical solutions of one-dimensional heat and advection-diffusion equations are obtained by collocation method based on cubicB-spline. Usual finite difference scheme is used for time and space integrations. Cubic B-spline is applied as interpolation function. The stability analysis of the scheme is examined by the Von Neumann approach. The efficiency of the method is illustrated by some test problems. The numerical results are found to be in good agreement with the exact solution.


---

## Body

## 1. Introduction

The combination of advection and diffusion is important for mass transport in fluids. It is well known that the volumetric concentration of a pollutant,u(x,t), at a point x(a≤x≤b) in a one-dimensional moving fluid with a constant speed β and diffusion coefficient α in x direction at time t(t≥0) is given by the one-dimensional advection-diffusion equation, which is in the form(1.1)ut+βux=αuxx,a≤x≤b,t≥0,
subject to the initial condition(1.2)u(x,0)=ϕ(x),x∈[a,b],
and the boundary conditions(1.3a)u(a,t)=g0(t),(1.3b)u(b,t)=g1(t),t∈[0,T],where g0 and g1 are assumed to be smooth functions. It should be noted that, when β=0, the advection-diffusion equation will be reduced to the one-dimensional heat equation in the case of thermal diffusion.Advection-diffusion equation arises very frequently in transferring mass, heat, energy, and vorticity in chemistry and engineering. Thus, it has been of interest to many authors. A third-degreeB-spline function has been used by Caglar et al. for solving one-dimensional heat equation with a nonlocal initial condition [1]. Mohebbi and Dehghan [2] have presented a fourth-order compact finite difference approximation and cubic C1-spline collocation method for the solution with fourth-order accuracy in both space and time variables, O(h4,k4). In [3], Dag and Saka concluded that collocation scheme is easy to implement compared to other numerical methods with giving a better result.In this paper, a combination of finite difference approach and cubicB-spline method would be considered to solve the one-dimensional heat and advection-diffusion equation. Forward finite difference approach would be used for discretizing the derivative of time, while cubic B-spline would be applied to interpolate the solutions at time t. Von Neumann approach would be used to prove the unconditionally stable property of the method. Finally, the approximated solutions and the numerical errors would be presented to demonstrate the efficiency of the method.


## 2. Collocation Method

In this paper, cubicB-splines are used to construct the numerical solutions to solve the problems. Consider a partition of [a,b] that is equally divided by knots xi into n subinterval [xi,xi+1], where i=0,1,…,n-1 such that a=x0<x1<⋯<xn=b. Hence, an approximation U(x,t) to the exact solution u(x,t) based on collocation approach can be expressed as [4](2.1)U(x,t)=∑i=-3n-1Ci(t)B3,i(x),
where Ci(t) are time-dependent quantities to be determined and B3,i(x) are third-degree B-spline functions which are defined by the relationship [5](2.2)B3,i(x)=16h3{(x-xi)3,x∈[xi,xi+1],h3+3h2(x-xi+1)+3h(x-xi+1)2-3(x-xi+1)3,x∈[xi+1,xi+2],h3+3h2(xi+3-x)+3h(xi+3-x)2-3(xi+3-x)3,x∈[xi+2,xi+3],(xi+4-x)3,x∈[xi+3,xi+4],
where h=(b-a)/n. The approximation Uik at the point (xi,tk) over the subinterval [xi,xi+1] can be simplified into(2.3)Uik=∑j=i-3i-1CjkB3,j(x),
where i=0,1,…,n. To obtain the approximations of the solutions, the values of B3,i(x) and its derivatives at the knots are needed. Since the values vanish at all other knots, they are omitted from Table 1.Table 1
Values ofBi, Bi′, and Bi′′.
xixi+1xi+2xi+3xi+4Bi01646160Bi′012h0-12h0Bi′′01h2-2h21h20The approximations of the solutions of (1.1) at tj+1th time level can be considered by [6]:(2.4)(Ut)ik+(1-θ)fik+θfik+1=0,wherefik=β(Ux)ik-α(Uxx)ik and the superscripts k and k+1 are successive time levels, k=0,1,2,…. Now, discretizing the time derivative by a first-order accurate forward difference scheme and rearranging the equation, we obtain(2.5)Uik+1+θΔtfik+1=Uik-(1-θ)Δtfik,
where Δt is the time step. Note that the system becomes an explicit scheme when θ=0, a fully implicit scheme when θ=1, and a mixed scheme of Crank-Nicolson when θ=0.5 [6]. Here, Crank-Nicolson approach is used. Hence, (2.5) takes the form(2.6)Uik+1+0.5Δtfik+1=Uik-0.5Δtfik
for i=0,1,…,n at each level of time. Therefore, a linear system of order (n+1) is obtained with (n+3) unknowns Ck+1=(C-3k+1,C-2k+1,…,Cn-1k+1) at the level time t=tk+1. To solve the system, two additional linear equations are needed. Thus, (2.3) is applied to the boundary conditions (1.3a)-(1.3b) to obtain(2.7a)U0k+1=g0(tk+1),(2.7b)Unk+1=g1(tk+1).Equations (2.6), (2.7a)-(2.7b) lead to a (n+3)×(n+3) tridiagonal matrix system, which can be solved by the Thomas algorithm. Once the initial vector C0 has been calculated from the initial conditions [7], the approximation solution Uik+1 at each level of time tk+1 can be determined by the vector Ck+1 which is found by solving the recurrence relation repeatedly.The initial vectorC0 can be obtained from the initial condition and boundary values of the derivatives of the initial condition as the following expressions [6]:(1)
(Ui0)x=ϕ′(xi),   i=0,(2)
Ui0=ϕ(xi),   i=0,1,…,n,(3)
(Ui0)x=ϕ′(xi),  i=n.This yields a(n+3)×(n+3) matrix system where the solution can be found by Thomas algorithms.


## 3. Stability Analysis

Von Neumann stability method is applied for analyzing the stability of the proposed scheme. This type of stability analysis had been used by many researchers [3, 8–10]. Consider the trial solution (one Fourier mode out of the full solution) at a given point xm(3.1)Cmk=δkexp⁡(iηmh),
where i=-1and η is the mode number. By substituting (2.3) into (2.5) and rearranging the equation, it leads to(3.2)p1Cm-3k+1+p2Cm-2k+1+p3Cm-1k+1=p4Cm-3k+p5Cm-2k+p6Cm-1k,
where
(3.3)p1=16+θΔtβ2h-θΔtαh2,p2=46+2θΔtαh2,p3=16-θΔtβ2h-θΔtαh2,p4=16-(1-θ)Δtβ2h+(1-θ)Δtαh2,p5=46-2(1-θ)Δtαh2,p6=16+(1-θ)Δtβ2h+(1-θ)Δtαh2.
Inserting the trial solution (3.1) into (3.2) and simplifying the equation give(3.4)δ=A+iBC+iD,
where(3.5)A=13(2+cos⁡ηh)-2(1-θ)Δtαh2(1-cos⁡ηh),B=(1-θ)Δtβhsin⁡ηh,C=13(2+cos⁡ηh)+2θΔtαh2(1-cos⁡ηh),D=-θΔtβhsin⁡ηh.
If the amplification factor |δ|≤1, then the proposed scheme is stable, or else the approximations grow in amplitude and become unstable. As θ=0.5 is used in the proposed scheme, thus substitute the θ value into (3.4) and after some algebraic manipulation, it can be noticed that(3.6)a2+b2≤c2+d2or|δ|2=a2+b2c2+d2≤1.
Thus, this had been proved that the presented numerical scheme for the advection-diffusion equation is unconditionally stable.


## 4. Numerical Results

### 4.1. Problem 1

Suppose the heat equation is as follows [11]:(4.1)ut=uxx,0<x<1,t>0,with initial and boundary conditions(4.2)u(x,0)=sin⁡(πx),u(0,t)=u(1,t)=0.
The exact solution is known to be u(x,t)=exp(-π2t)sin(πx). This problem is tested by different values of h and Δt to show the capability of the presented method for solving one-dimensional heat equation. The final time is chosen as T=1. The maximum absolute errors of the method are compared with those obtained by Crank-Nicolson (CN) scheme and compact boundary value method (CBVM) in [11]. The numerical errors are presented in Table 2. Although the fourth-order compact boundary value method gives a much more better solution, the present method is still well compared with the Crank-Nicolson scheme.Table 2
Maximum absolute error obtained for problem 1.
h=ΔtCN [11]CBVM [11]Present method1/51.1×10-12.8×10-21.4145×10-11/103.0×10-23.8×10-33.7195×10-21/206.9×10-32.7×10-48.4588×10-31/401.7×10-31.3×10-52.0698×10-31/804.2×10-45.1×10-75.1473×10-4


### 4.2. Problem 2

Consider the advection-diffusion equation in (1.1) with β=1, α=0.1, as follows [2]:(4.3)ut+ux=0.1uxx,0<x<1,t>0,where the initial condition is given by(4.4)u(x,0)=exp⁡(5x)[cos⁡(π2x)+0.25sin⁡(π2x)]and the exact solution(4.5)u(x,t)=exp⁡(5(x-t2))exp⁡(-π240t)[cos⁡(π2x)+0.25sin⁡(π2x)].
The boundary conditions at x=0 and x=1 can be obtained from the exact solution. Table 3 shows the absolute errors of the approximations at the grid points when T=2. It can be noticed that the present method is comparable with cubic C1-spline collocation method. The approximations of the solutions over a time period t∈[0,2] along x is depicted in Figure 1.Table 3
Absolute error obtained withΔt=2h at T=2 for problem 2.
Grid pointh=0.02h=0.01C1-spline [2]Present methodC1-spline [2]Present method0.17.1744×10-68.2212×10-61.8035×10-62.0556×10-60.21.1019×10-52.2566×10-52.7685×10-65.6432×10-60.31.6596×10-54.5188×10-54.1679×10-61.1298×10-50.42.4579×10-57.7748×10-56.1705×10-61.9435×10-50.53.5871×10-51.2011×10-49.0026×10-63.0020×10-50.65.1637×10-51.6809×10-41.2955×10-54.2001×10-50.77.3208×10-52.1002×10-41.8360×10-55.2464×10-50.81.0163×10-42.2264×10-42.5476×10-55.5602×10-50.91.3624×10-41.6833×10-43.4134×10-54.2039×10-5Figure 1
Spatial-time approximations for problem 2 withh=0.05 and Δt=0.5h.


## 4.1. Problem 1

Suppose the heat equation is as follows [11]:(4.1)ut=uxx,0<x<1,t>0,with initial and boundary conditions(4.2)u(x,0)=sin⁡(πx),u(0,t)=u(1,t)=0.
The exact solution is known to be u(x,t)=exp(-π2t)sin(πx). This problem is tested by different values of h and Δt to show the capability of the presented method for solving one-dimensional heat equation. The final time is chosen as T=1. The maximum absolute errors of the method are compared with those obtained by Crank-Nicolson (CN) scheme and compact boundary value method (CBVM) in [11]. The numerical errors are presented in Table 2. Although the fourth-order compact boundary value method gives a much more better solution, the present method is still well compared with the Crank-Nicolson scheme.Table 2
Maximum absolute error obtained for problem 1.
h=ΔtCN [11]CBVM [11]Present method1/51.1×10-12.8×10-21.4145×10-11/103.0×10-23.8×10-33.7195×10-21/206.9×10-32.7×10-48.4588×10-31/401.7×10-31.3×10-52.0698×10-31/804.2×10-45.1×10-75.1473×10-4


## 4.2. Problem 2

Consider the advection-diffusion equation in (1.1) with β=1, α=0.1, as follows [2]:(4.3)ut+ux=0.1uxx,0<x<1,t>0,where the initial condition is given by(4.4)u(x,0)=exp⁡(5x)[cos⁡(π2x)+0.25sin⁡(π2x)]and the exact solution(4.5)u(x,t)=exp⁡(5(x-t2))exp⁡(-π240t)[cos⁡(π2x)+0.25sin⁡(π2x)].
The boundary conditions at x=0 and x=1 can be obtained from the exact solution. Table 3 shows the absolute errors of the approximations at the grid points when T=2. It can be noticed that the present method is comparable with cubic C1-spline collocation method. The approximations of the solutions over a time period t∈[0,2] along x is depicted in Figure 1.Table 3
Absolute error obtained withΔt=2h at T=2 for problem 2.
Grid pointh=0.02h=0.01C1-spline [2]Present methodC1-spline [2]Present method0.17.1744×10-68.2212×10-61.8035×10-62.0556×10-60.21.1019×10-52.2566×10-52.7685×10-65.6432×10-60.31.6596×10-54.5188×10-54.1679×10-61.1298×10-50.42.4579×10-57.7748×10-56.1705×10-61.9435×10-50.53.5871×10-51.2011×10-49.0026×10-63.0020×10-50.65.1637×10-51.6809×10-41.2955×10-54.2001×10-50.77.3208×10-52.1002×10-41.8360×10-55.2464×10-50.81.0163×10-42.2264×10-42.5476×10-55.5602×10-50.91.3624×10-41.6833×10-43.4134×10-54.2039×10-5Figure 1
Spatial-time approximations for problem 2 withh=0.05 and Δt=0.5h.


## 5. Conclusions

A numerical method based on collocation of cubicB-spline had been described in the previous section for solving one-dimensional heat and advection-diffusion equations. A finite difference scheme had been used for discretizing time derivatives and cubic B-spline for interpolating the solutions at each time level. From the test problems, the obtained results show that the presented method is capable for solving one-dimensional heat and advection-diffusion equations accurately with a promised stability.



---

*Source: 458701-2012-05-13.xml*