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

Abstract

Numerical solutions of one-dimensional heat and advection-diffusion equations are obtained by collocation method based on cubic B-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.

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, 𝑢(𝑥,𝑡), at a point 𝑥(𝑎≤𝑥≤𝑏) in a one-dimensional moving fluid with a constant speed 𝛽 and diffusion coefficient 𝛼 in 𝑥 direction at time 𝑡(𝑡≥0) is given by the one-dimensional advection-diffusion equation, which is in the form𝑢_𝑡+𝛽𝑢_𝑥=𝛼𝑢_𝑥𝑥,𝑎≤𝑥≤𝑏,𝑡≥0,(1.1) subject to the initial condition𝑢(𝑥,0)=𝜙(𝑥),𝑥∈[𝑎,𝑏],(1.2) and the boundary conditions𝑢(𝑏,𝑡)=𝑢(𝑔𝑎₁,𝑡)=𝑔₀(𝑡),(1.3a)(𝑡),𝑡∈[0,T],(1.3b)where 𝑔₀ and 𝑔₁ 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-degree 𝐵-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 𝐶¹-spline collocation method for the solution with fourth-order accuracy in both space and time variables, 𝑂(ℎ⁴,𝑘⁴). 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 cubic 𝐵-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 𝐵-spline would be applied to interpolate the solutions at time 𝑡. 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, cubic 𝐵-splines are used to construct the numerical solutions to solve the problems. Consider a partition of [𝑎,𝑏] that is equally divided by knots 𝑥_𝑖 into 𝑛 subinterval [𝑥_𝑖,𝑥_𝑖+1], where 𝑖=0,1,…,𝑛−1 such that 𝑎=𝑥₀<𝑥₁<⋯<𝑥_𝑛=𝑏. Hence, an approximation 𝑈(𝑥,𝑡) to the exact solution 𝑢(𝑥,𝑡) based on collocation approach can be expressed as [4]𝑈(𝑥,𝑡)=^𝑛−1_𝑖₌₋₃𝐶_𝑖(𝑡)𝐵_3,𝑖(𝑥),(2.1) where 𝐶_𝑖(𝑡) are time-dependent quantities to be determined and 𝐵_3,𝑖(𝑥) are third-degree 𝐵-spline functions which are defined by the relationship [5]𝐵_3,𝑖(𝑥)=16ℎ³⎧⎪⎪⎩⎪⎨⎪ℎ𝑥₃−𝑥_𝑖³,𝑥∈𝑥_𝑖,𝑥_𝑖+1,ℎ₃+3ℎ²𝑥−𝑥_𝑖+1+3ℎ𝑥−𝑥_𝑖+1²−3𝑥−𝑥_𝑖+1³,𝑥∈𝑥_𝑖+1,𝑥_𝑖+2,𝑥+_𝑖+43ℎ²𝑥_𝑖+3−𝑥+3ℎ𝑥_𝑖+3−𝑥²−3𝑥_𝑖+3−𝑥³,𝑥∈𝑥_𝑖+2,𝑥_𝑖+3,−𝑥³,𝑥∈𝑥_𝑖+3,𝑥_𝑖+4,(2.2) where ℎ=(𝑏−𝑎)/𝑛. The approximation 𝑈^𝑘𝑖 at the point (𝑥_𝑖,𝑡_𝑘) over the subinterval [𝑥_𝑖,𝑥_𝑖+1] can be simplified into𝑈^𝑘𝑖=^𝑖−1_𝑗_=𝑖−3𝐶^𝑘𝑗𝐵_3,𝑗(𝑥),(2.3) where 𝑖=0,1,…,𝑛. To obtain the approximations of the solutions, the values of 𝐵_3,𝑖(𝑥) 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 of 𝐵_𝑖, 𝐵^𝑖, and 𝐵^𝑖.

The approximations of the solutions of (1.1) at 𝑡_𝑗+1th time level can be considered by [6]:𝑈_𝑡^𝑘𝑖+(1−𝜃)𝑓^𝑘𝑖+𝜃𝑓^𝑖𝑘+1=0,(2.4)

where 𝑓^𝑘𝑖=𝛽(𝑈_𝑥)^𝑘𝑖−𝛼(𝑈_𝑥𝑥)^𝑘𝑖 and the superscripts 𝑘 and 𝑘+1 are successive time levels, 𝑘=0,1,2,…. Now, discretizing the time derivative by a first-order accurate forward difference scheme and rearranging the equation, we obtain𝑈^𝑖𝑘+1+𝜃Δ𝑡𝑓^𝑖𝑘+1=𝑈^𝑘𝑖−(1−𝜃)Δ𝑡𝑓^𝑘𝑖,(2.5) where Δ𝑡 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𝑈^𝑖𝑘+1+0.5Δ𝑡𝑓^𝑖𝑘+1=𝑈^𝑘𝑖−0.5Δ𝑡𝑓^𝑘𝑖(2.6) for 𝑖=0,1,…,𝑛 at each level of time. Therefore, a linear system of order (𝑛+1) is obtained with (𝑛+3) unknowns 𝐂^𝑘+1=(𝐶^𝑘+1−3,𝐶^𝑘+1−2,…,𝐶^{𝑘+1𝑛−1}) at the level time 𝑡=𝑡_𝑘+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𝑈𝑈^0𝑘+1=𝑔₀𝑡_𝑘+1,(2.7a)_𝑛𝑘+1=𝑔₁𝑡_𝑘+1.(2.7b)Equations (2.6), (2.7a)-(2.7b) lead to a (𝑛+3)×(𝑛+3) tridiagonal matrix system, which can be solved by the Thomas algorithm. Once the initial vector 𝐂⁰ has been calculated from the initial conditions [7], the approximation solution 𝑈^𝑖𝑘+1 at each level of time 𝑡_𝑘+1 can be determined by the vector 𝐂^𝑘+1 which is found by solving the recurrence relation repeatedly.

The initial vector 𝐂⁰ can be obtained from the initial condition and boundary values of the derivatives of the initial condition as the following expressions [6]:(1)(𝑈^0𝑖)_𝑥=𝜙^(𝑥_𝑖), 𝑖=0,(2)𝑈^0𝑖=𝜙(𝑥_𝑖), 𝑖=0,1,…,𝑛,(3)(𝑈^0𝑖)_𝑥=𝜙^(𝑥_𝑖), 𝑖=𝑛.

This yields a (𝑛+3)×(𝑛+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 𝑥_𝑚𝐶^𝑘𝑚=𝛿^𝑘exp(𝑖𝜂𝑚ℎ),(3.1) where 𝑖=√−1and 𝜂 is the mode number. By substituting (2.3) into (2.5) and rearranging the equation, it leads to𝑝₁𝐶^{𝑘+1𝑚−3}+𝑝₂𝐶^{𝑘+1𝑚−2}+𝑝₃𝐶^{𝑘+1𝑚−1}=𝑝₄𝐶^𝑘𝑚−3+𝑝₅𝐶^𝑘𝑚−2+𝑝₆𝐶^𝑘𝑚−1,(3.2) where 𝑝₁=16+𝜃2Δℎ𝑡𝛽−𝜃Δℎ²𝑡𝛼𝑝₂,=46+2𝜃ℎΔ²𝑡𝛼𝑝₃,=16−𝜃2Δℎ𝑡𝛽−𝜃Δℎ²𝑡𝛼𝑝₄,=16−(1−2𝜃ℎ)Δ𝑡𝛽+(1−ℎ𝜃²)Δ𝑡𝛼𝑝₅,=46−2(1−ℎ𝜃²)Δ𝑡𝛼𝑝₆,=16+(1−2𝜃ℎ)Δ𝑡𝛽+(1−ℎ𝜃²)Δ𝑡𝛼.(3.3) Inserting the trial solution (3.1) into (3.2) and simplifying the equation give𝛿=𝐴+𝑖𝐵𝐶+𝑖𝐷,(3.4) where𝐴=13(2+cos𝜂ℎ)−2(1−ℎ𝜃²)Δ𝑡𝛼𝐵=(1−𝜃)Δ𝑡𝛽(1−cos𝜂ℎ),𝐶=1ℎsin𝜂ℎ,3(2+cos𝜂ℎ)+2𝜃ℎΔ²𝑡𝛼𝐷=−𝜃Δ𝑡𝛽(1−cos𝜂ℎ),ℎsin𝜂ℎ.(3.5) 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𝑎²+𝑏²≤𝑐²+𝑑²or||𝛿||²=𝑎²+𝑏²𝑐²+𝑑²≤1.(3.6) 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]:𝑢_𝑡=𝑢_𝑥𝑥,0<𝑥<1,𝑡>0,(4.1)

with initial and boundary conditions𝑢(𝑥,0)=sin(𝜋𝑥),𝑢(0,𝑡)=𝑢(1,𝑡)=0.(4.2) The exact solution is known to be 𝑢(𝑥,𝑡)=exp(−𝜋²𝑡)sin(𝜋𝑥). This problem is tested by different values of ℎ and Δ𝑡 to show the capability of the presented method for solving one-dimensional heat equation. The final time is chosen as 𝑇=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.

4.2. Problem 2

Consider the advection-diffusion equation in (1.1) with 𝛽=1, 𝛼=0.1, as follows [2]:𝑢_𝑡+𝑢_𝑥=0.1𝑢_𝑥𝑥,0<𝑥<1,𝑡>0,(4.3)

where the initial condition is given by𝑢(𝑥,0)=exp(5𝑥)cos𝜋2𝑥+0.25sin𝜋2𝑥(4.4)

and the exact solution𝑢(𝑥,𝑡)=exp5𝑥−𝑡2exp−4𝜋0²𝑡cos𝜋2𝑥+0.25sin𝜋2𝑥.(4.5) The boundary conditions at 𝑥=0 and 𝑥=1 can be obtained from the exact solution. Table 3 shows the absolute errors of the approximations at the grid points when 𝑇=2. It can be noticed that the present method is comparable with cubic 𝐶¹-spline collocation method. The approximations of the solutions over a time period 𝑡∈[0,2] along 𝑥 is depicted in Figure 1.

Table 3 

Absolute error obtained with Δ𝑡=2ℎ at 𝑇=2 for problem 2.

Figure 1 

Spatial-time approximations for problem 2 with ℎ=0.05 and Δ𝑡=0.5ℎ.

5. Conclusions

A numerical method based on collocation of cubic 𝐵-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 𝐵-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.

Acknowledgment

The authors would like to acknowledge with thanks the financial support from Malaysian Government in the form of Fundamental Research Grant Scheme (FRGS) of number 203/PMATHS/6711150.