Numerical Method Using Cubic Trigonometric B-Spline Technique for Nonclassical Diffusion Problems

Abstract

A new two-time level implicit technique based on cubic trigonometric B-spline is proposed for the approximate solution of a nonclassical diffusion problem with nonlocal boundary constraints. The standard finite difference approach is applied to discretize the time derivative while cubic trigonometric B-spline is utilized as an interpolating function in the space dimension. The technique is shown to be unconditionally stable using the von Neumann method. Several numerical examples are discussed to exhibit the feasibility and capability of the technique. The 𝐿₂ and 𝐿_∞ error norms are also computed at different times for different space size steps to assess the performance of the proposed technique. The technique requires smaller computational time than several other methods and the numerical results are found to be in good agreement with known solutions and with existing schemes in the literature.

1. Introduction

This study deals with the numerical solution of a nonclassical diffusion problem with two nonlocal boundary constraints using cubic trigonometric B-splines. This problem arises in several branches of science. In particular, electrochemistry [1], heat conduction process [2], thermoelasticity [3], plasma physics [4], semiconductor modeling [5], biotechnology [6], control theory, and inverse problems [7]. The analysis, development, and implementation of numerical methods for the solution of such diffusion problems have received wide attention in the literature.

Consider an insulated rod of length 𝐿 located on the 𝑥-axis of the interval [0,𝐿]. Let the rod have a source of heat. Let 𝑢(𝑥,𝑡) denote the temperature in the insulated rod with ends held at constant temperatures 𝑇₁ and 𝑇₂, and the initial temperature distribution along the rod is 𝑔₁(𝑥). The problem is to study the flow of heat in the rod and in this paper the partial differential equation governing the flow of heat in the rod is given by the diffusion equation with specification of energy: 𝜕𝑢()𝜕𝑡𝑥,𝑡=𝛼²𝜕²𝑢𝜕𝑥²(𝑥,𝑡)+𝑞(𝑥,𝑡)0≤𝑥≤𝐿,0≤𝑡≤𝑇 (1) with the initial constraint 𝑢(𝑥,𝑡=0)=𝑔₁(𝑥)0≤𝑥≤𝐿 (2) and the nonlocal boundary constraints 𝜉₁𝑢(𝑥=0,𝑡)+𝜉₂𝑢_𝑥()=𝑥=0,𝑡_𝐿0𝑔₂(𝑥)𝑢(𝑥,𝑡)𝑑𝑥+ℎ₁(𝑡)=𝑇₁,𝜉₃𝑢(𝑥=𝐿,𝑡)+𝜉₄𝑢_𝑥(𝑥=𝐿,𝑡)=^𝐿₀𝑔₃(𝑥)𝑢(𝑥,𝑡)𝑑𝑥+ℎ₂(𝑡)=𝑇₂0<𝑡≤𝑇, (3) where 𝜉_𝑖, 𝑖=1,2,3,4 are known constants, 𝑞,𝑔_𝑖(𝑖=1,2,3) are known continuous functions, 𝛼 is the thermal diffusivity of the rod, and ℎ_𝑖(𝑖=1,2) are prescribed functions. This problem has been studied by Dehghan [8], Martín-Vaquero and Vigo-Aguiar [9], Li and Wu [10], and Golbabai and Javidi [11]. Several physical circumstances might be modeled by equation and constraints (1)–(3) and several examples of application in physics with comprehensive derivations of the above mentioned problem can be found in [3, 12–14].

There are several numerical methods in the literature that have been developed for solving the proposed problem (1) subject to initial and nonlocal boundary constraints (2)-(3). The methods were based, for instance, on the forward Euler method, the backward Euler approach or the Crank-Nicolson scheme [15, 16], Laplace transformation [17], and so forth. Dehghan [8] presented four finite difference approaches, namely, the BTCS (backward time centred space) scheme, the implicit (3,3) Crandall’s formula, the 3-point FTCS (forward time centred space) two-level scheme, and the Dufort-Frankel three-level approach for the numerical solution of parabolic equation with nonlocal specification. Martín-Vaquero and Vigo-Aguiar [9] improved the order of convergence of the implicit (3,3) Crandall’s formula proposed by Dehghan [8] and also improved the accuracy of the method. Martín-Vaquero and Vigo-Aguiar [18] developed an algorithm for the solution of the heat conduction equations with nonlocal constraints which reduced the CPU time and enhanced the accuracy of (3, 3) Crandall’s formula proposed in [8]. Li and Wu [10] proposed an algorithm which was based on the transverse method of lines (TMOL) which can reduce a nonclassical diffusion equation to a series of ordinary differential equations (ODEs). Subsequently, the authors in [10] used an analytic reproducing kernel technique to solve ODEs with integral boundary constraints. Dehghan [19, 20], jointly with Tatari and Dehghan [21], has proposed several efficient techniques for the numerical solutions of partial differential equations subject to nonlocal boundary constraints. Golbabai and Javidi [11] introduced a Chebyshev spectral collocation method (CSCM) based on Chebyshev polynomials for solving a parabolic problem subject to nonlocal boundary constraints. For more details on other numerical methods for the solution of a one-dimensional heat equation subject to nonlocal boundary constraints in the literature, see [22–31].

The study of B-spline functions is a key element in computer-aided geometric design [32–35]. It has also attracted attention in the literature [36–51] to the numerical solution of various differential equations [38–40]. This is because they have important geometric properties and features that make them amenable to more detailed analysis. Numerical methods based on B-spline functions of various degrees have been utilized for solving initial and boundary value problems. As examples, a cubic B-spline collocation method was used to solve a nonlinear diffusion equation subject to certain initial and Dirichlet boundary constraints [41], a finite element method based on bivariate splines has been used for solving parabolic partial differential equation [42], and the combination of finite difference approach and cubic B-spline method was applied for the solution of a one-dimensional heat equation subject to local boundary constraints [43, 44]. Goh et al. [45] presented a comparison of cubic B-spline and extended cubic uniform B-spline based collocation methods for solving a one-dimensional heat equation with a nonlocal initial constraint and concluded that extended cubic uniform B-spline with an appropriate value of parameters gives better results than the cubic B-spline. A finite difference scheme based on cubic B-spline was also used for solving the one-dimensional wave equation [43], advection-diffusion equation [44], one-dimensional coupled viscous Burgers’ equation [47], system of strongly coupled reaction-diffusion equations [48], and one-dimensional hyperbolic problems [49].

In our present paper, a new two-time level implicit technique is developed to approximate the solution of the nonclassical diffusion problem (1) subject to initial constraints in (2) and nonlocal boundary constraints in (1)–(3). The technique is based on the cubic trigonometric B-spline functions. A finite difference approach and 𝜃-weighted scheme are applied for the time and space discretization, respectively. Some researchers have considered the ordinary B-spline collocation method for solving the heat equation subject to local and nonlocal boundary constraints but, so far as we are aware, not with the cubic trigonometric B-spline collocation method. Cubic trigonometric B-spline is used as an interpolating function in the space dimension. The unconditional stability property of the method is proved by von Neumann method. The feasibility of the method is shown by test problems with 𝑘=𝑠ℎ², 𝑠=1,2,4,5 instead of smaller time step size 𝑘=0.4ℎ² and the approximated solutions are found to be in good agreement with the known exact solutions.

The outline of this study is as follows. A numerical solution of nonclassical diffusion problem is presented in Section 2. In Section 3, the cubic trigonometric B-spline is utilized as an interpolating function in the space dimension. The von Neumann approach is used to prove the stability of the method in Section 4. Numerical examples are considered in Section 5 to show the achievability of the proposed method. Finally, the concluding remarks of this study are given in Section 6.

2. Solution of Nonclassical Diffusion Problem

Consider a uniform mesh Ω with grid points (𝑥_𝑖,𝑡_𝑛) to discretize the grid region Δ=[𝑎,𝑏]×[0,𝑇] with 𝑥_𝑖=𝑎+𝑖ℎ, 𝑖=0,1,2,…,𝑁 and 𝑡_𝑛=𝑛𝑘, 𝑛=0,1,2,3,…,𝑀, 𝑀𝑘=𝑇. Here the quantities ℎ and 𝑘 are mesh space size and time step size, respectively. The time derivative can be approximated by using the standard finite difference formula: 𝜕𝑢^𝑛=𝑢𝜕𝑡^𝑛+1−𝑢^𝑛𝑘. (4) Using the approximation of (4), (1) becomes 𝑢^𝑛+1−𝑢^𝑛𝑘=𝛼²𝜕²𝑢^𝑛𝜕𝑥²+𝑞𝑥_𝑖,𝑡_𝑛+1. (5) Using 𝜃-weighted technique, the space derivatives of (5) can be written as 𝑢^𝑛+1−𝑢^𝑛𝑘=𝜃𝛼²𝜕²𝑢^𝑛+1𝜕𝑥²+(1−𝜃)𝛼²𝜕²𝑢^𝑛𝜕𝑥²+𝑞𝑥_𝑖,𝑡_𝑛+1, (6) where 0≤𝜃≤1 and the subscripts 𝑛 and 𝑛+1 are successive time levels. It is noted that the system becomes an explicit scheme when 𝜃=0, a fully implicit scheme when 𝜃=1, and a Crank-Nicolson scheme when 𝜃=1/2 [43, 49]. In this paper, we use the Crank-Nicolson approach. Hence, (6) becomes 𝑢^𝑛+1−𝑢^𝑛𝑘=12𝛼²𝜕²𝑢^𝑛+1𝜕𝑥²+12𝛼²𝜕²𝑢^𝑛𝜕𝑥²+𝑞𝑥_𝑖,𝑡_𝑛+1. (7) After simplification, (7) leads to 2𝑢^𝑛+1−𝑘𝛼²𝑢^𝑛+1_𝑥𝑥=2𝑢^𝑛+𝑘𝛼²𝑢^𝑛_𝑥𝑥+2𝑘𝑞𝑥_𝑖,𝑡_𝑛+1. (8) The space derivatives are approximated by using cubic trigonometric B-spline and are discussed in the next section.

3. Cubic Trigonometric B-Spline Technique

In this section, we discuss the cubic trigonometric B-spline collocation method (CuTBSM) for the numerical solution of the nonclassical diffusion equation (1). Consider a mesh 𝑎≤𝑥≤𝑏 which is equally divided by knots 𝑥_𝑖 into 𝑁 subintervals [𝑥_𝑖,𝑥_𝑖+1], 𝑖=0,1,2,…,𝑁−1 where 𝑎=𝑥₀<𝑥₁<⋯<𝑥_𝑁=𝑏. Our approach for the nonclassical diffusion equation using collocation method with cubic trigonometric B-spline is to seek an approximate solution as [37] 𝑈(𝑥,𝑡)=^𝑁−1_𝑖=−3𝐶_𝑖(𝑡)𝑇𝐵_𝑖(𝑥), (9) where 𝐶_𝑖(𝑡) are to be determined for the approximated solutions 𝑈(𝑥,𝑡) to the exact solutions 𝑢(𝑥,𝑡), at the point (𝑥_𝑖,𝑡_𝑛). 𝑇𝐵_𝑖(𝑥) are twice continuously differentiable piecewise cubic trigonometric B-spline basis functions over the mesh defined by [49–51] 𝑇𝐵_𝑖(𝑥)=1𝜔𝑝³𝑥_𝑖,𝑥∈𝑥_𝑖,𝑥_𝑖+1𝑝𝑥_𝑖𝑝𝑥_𝑖𝑞𝑥_𝑖+2+𝑞𝑥_𝑖+3𝑝𝑥_𝑖+1+𝑞𝑥_𝑖+4𝑝²𝑥_𝑖+1,𝑥∈𝑥_𝑖+1,𝑥_𝑖+2𝑞𝑥_𝑖+4𝑝𝑥_𝑖+1𝑞𝑥_𝑖+3+𝑞𝑥_𝑖+4𝑝𝑥_𝑖+2+𝑝𝑥_𝑖𝑞²𝑥_𝑖+3,𝑥∈𝑥_𝑖+2,𝑥_𝑖+3𝑞³𝑥_𝑖+4,𝑥∈𝑥_𝑖+3,𝑥_𝑖+4, (10) where 𝑝𝑥_𝑖=sin𝑥−𝑥_𝑖2,𝑞𝑥_𝑖𝑥=sin^𝑖−𝑥2,ℎ𝜔=sin2sin(ℎ)sin3ℎ2 (11) and where ℎ=(𝑏−𝑎)/𝑛. The approximations 𝑈^𝑛_𝑖 at the point (𝑥_𝑖,𝑡_𝑛) over subinterval [𝑥_𝑖,𝑥_𝑖+1] can be defined as 𝑈^𝑛_𝑖=^𝑖−1_{𝑗=𝑖−3}𝐶^𝑛_𝑘𝑇𝐵_𝑗(𝑥). (12) In order to obtain the approximations to the solutions, the values of 𝑇𝐵_𝑖(𝑥) and its derivatives at nodal points are required and these derivatives are tabulated in Table 1, where 𝑎₁=sin²(ℎ/2)(ℎ)(),𝑎sinsin3ℎ/2₂=2(ℎ),𝑎1+2cos₃3=−(),𝑎4sin3ℎ/2₄=3(),𝑎4sin3ℎ/2₅=3(1+3cos(ℎ))16sin²(()(,𝑎ℎ/2)(2cosℎ/2+cos3ℎ/2))₆=−3cos²(ℎ/2)sin²((ℎ.ℎ/2)(2+4cos)) (13) Using approximate functions (10) and (12) and following Mittal and Arora [46], the values at the knots of 𝑈^𝑛_𝑖 and their derivatives up to second order are determined in terms of time parameters 𝐶^𝑛_𝑗 as 𝑢^𝑛_𝑖=𝑎₁𝐶^𝑛_𝑖−3+𝑎₂𝐶^𝑛_𝑖−2+𝑎₁𝐶^𝑛_𝑖−1,𝑢_𝑥^𝑛_𝑖=𝑎₃𝐶^𝑛_𝑖−3+𝑎₄𝐶^𝑛_𝑖−1,𝑢_𝑥𝑥^𝑛_𝑖=𝑎₅𝐶^𝑛_𝑖−3+𝑎₆𝐶^𝑛_𝑖−2+𝑎₅𝐶^𝑛_𝑖−1. (14) Substituting (12) into (8) gives the following equation: 2^𝑖−1_{𝑗=𝑖−3}𝐶_𝑗^𝑛+1𝑇𝐵_𝑗𝑥_𝑖−𝑘𝛼^2𝑖−1_{𝑗=𝑖−3}𝐶_𝑗^𝑛+1𝑇𝐵_𝑗^′′𝑥_𝑖=2^𝑖−1_{𝑗=𝑖−3}𝐶^𝑛_𝑗𝑇𝐵_𝑗𝑥_𝑖+𝑘𝛼^2𝑖−1_{𝑗=𝑖−3}𝐶^𝑛_𝑗𝑇𝐵_𝑗^′′𝑥_𝑖+2𝑘𝑞𝑥_𝑖,𝑡_𝑛+1. (15) The system thus obtained on simplifying (15) consists of 𝑁+1 linear equations in 𝑁+3 unknowns 𝐶^𝑛+1=(𝐶^𝑛+1₋₃,𝐶^𝑛+1₋₂,𝐶^𝑛+1₋₁,…,𝐶^𝑛+1_𝑁−1) at the time level 𝑡=𝑡_𝑛+1. Equation (9) is applied to the boundary constraints (2) and (3) for two additional linear equations to obtain a unique solution of the resulting system: 𝜉₁𝑈0,𝑡_𝑛+1+𝜉₂𝑈_𝑥0,𝑡_𝑛+1=^𝐿₀𝑔₂(𝑥)𝑢𝑥,𝑡_𝑛+1𝑑𝑥+ℎ₁𝑡_𝑛+1,𝜉₃𝑈𝐿,𝑡_𝑛+1+𝜉₄𝑈_𝑥𝐿,𝑡_𝑛+1=^𝐿₀𝑔₃(𝑥)𝑢𝑥,𝑡_𝑛+1𝑑𝑥+ℎ₂𝑡_𝑛+1. (16)(17) From (15), (16), and (17), the system can be written in the matrix vector form as follows: 𝑀𝐶^𝑛+1=𝑁𝐶^𝑛+𝑏, (18) where 𝐶^𝑛+1=𝐶^𝑛+1₋₃,𝐶^𝑛+1₋₂,𝐶^𝑛+1₋₁,…,𝐶^𝑛+1_𝑁−1^𝑇,𝐶^𝑛=𝐶^𝑛₋₃,𝐶^𝑛₋₂,𝐶^𝑛₋₁,…,𝐶^𝑛_𝑁−1^𝑇,𝑏=𝛼₁𝑡_𝑛+1,2𝑘𝑞𝑥₀,𝑡_𝑛+1,2𝑘𝑞𝑥₁,𝑡_𝑛+1,…,2𝑘𝑞𝑥_𝑁,𝑡_𝑛+1,𝛽₁𝑡_𝑛+1^𝑇,𝑛=0,1,2,…,𝑀 (19) and 𝑀 and 𝑁 are 𝑁+3-dimensional matrix given by 𝑝𝑀=⁰𝑞₀𝑝₀₀𝑝0⋯00₁𝑞₁𝑝₁0⋯000𝑝₁𝑞₁𝑝₁0⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮00⋯⋯𝑝₁𝑞₁𝑝₁00⋯⋯𝑝₀𝑞₀𝑝₀₀,𝑝𝑁=0000⋯00₂𝑞₂𝑝₂0⋯000𝑝₂𝑞₂𝑝₂0⋯0⋮⋮⋮⋮⋮⋮⋮00⋯⋯𝑝₂𝑞₂𝑝₂,00⋯⋯000 (20) where 𝑝₀=𝜉₁𝑎₁+𝜉₂𝑎₃,𝑞₀=𝜉₁𝑎₁,𝑝₀₀=𝜉₁𝑎₁+𝜉₂𝑎₄,𝑝₁=2𝑎₁−𝑘𝛼²𝑎₅,𝑞₁=2𝑎₂−𝑘𝛼²𝑎₆,𝑝₂=2𝑎₁+𝑘𝛼²𝑎₅,𝑞₂=2𝑎₂+𝑘𝛼²𝑎₆,𝛼₁𝑡_𝑛+1=^𝐿₀𝑔₂(𝑥)𝑢𝑥,𝑡_𝑛+1𝑑𝑥+ℎ₁𝑡_𝑛+1,𝛽₁𝑡_𝑛+1=^𝐿₀𝑔₃(𝑥)𝑢𝑥,𝑡_𝑛+1𝑑𝑥+ℎ₂𝑡_𝑛+1. (21) Thus, the system (18) becomes a matrix system of dimension (𝑁+3)×(𝑁+3) which is a tridiagonal system that can be solved by the Thomas Algorithm [16].

Table 1 

Values TB_𝑖(𝑥) and its derivatives.

3.1. Initial State Vector 𝐶⁰

After the initial vectors 𝐶⁰ have been computed from the initial constraints, the approximate solutions 𝑈_𝑖^𝑛+1 at a particular time level can be calculated repeatedly by solving the recurrence relation (15) [40].

The initial vectors 𝐶⁰ can be obtained from the initial condition and boundary values of the derivatives of the initial condition as follows [40, 49]: 𝑈⁰_𝑖𝑥=𝑔^′₁𝑥_𝑖𝑈,𝑖=0,_0𝑖=𝑔₁𝑥_𝑖𝑈,𝑖=0,1,2,…,𝑁,_0𝑖𝑥=𝑔^′₁𝑥_𝑖,𝑖=𝑁. (22) Thus (22) yields a (𝑁+3)×(𝑁+3) matrix system, of the form 𝐴𝐶^𝑛=𝑑, (23) where 𝐶^𝑛=𝐶^𝑛₋₃,𝐶^𝑛₋₂,𝐶^𝑛₋₁,…,𝐶^𝑛_𝑁−1^𝑇,𝑑=𝑔^′₁𝑥₀,𝑔₁𝑥₀,…,𝑔₁𝑥_𝑁,𝑔^′₁𝑥_𝑁^𝑇𝑎,𝑛=0,𝐴=³0𝑎₄𝑎0⋯00₁𝑎₂𝑎₁0⋯000𝑎₁𝑎₂𝑎₁0⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮00⋯⋯𝑎₁𝑎₂𝑎₁00⋯⋯𝑎₃0𝑎₄. (24) The solution of this system can be found by the use of the Thomas Algorithm.

4. Stability

In this section, the von Neumann stability method is applied for investigating the stability of the proposed scheme. This approach has been used by many researchers [18, 40, 44, 45, 47–49]. Substituting the approximate solution 𝑈 and its derivatives at knots with 𝑞(𝑥,𝑡)=0, into (6) yields a difference equation with variables 𝐶_𝑚 given by 𝑎₁−𝑘𝜃𝛼²𝑎₅𝐶^𝑛+1_𝑚−3+𝑎₂−𝑘𝜃𝛼²𝑎₆𝐶^𝑛+1_𝑚−2+𝑎₁−𝑘𝜃𝛼²𝑎₅𝐶^𝑛+1_𝑚−1=𝑎₁+𝑘(1−𝜃)𝛼²𝑎₅𝐶^𝑛_𝑚−3+𝑎₂+𝑘(1−𝜃)𝛼²𝑎₆𝐶^𝑛_𝑚−2+𝑎₁+𝑘(1−𝜃)𝛼²𝑎₅𝐶^𝑛_𝑚−1. (25) Substituting the values of 𝑎_𝑖,𝑖=1,2,5,6 into (25) we obtain 𝑝₁−𝑘²𝜃𝛼²𝑝₂𝐶^𝑛+1_𝑚−3+𝑝₃+𝑘²𝜃𝛼²𝑝₄𝐶^𝑛+1_𝑚−2+𝑝₁−𝑘²𝜃𝛼²𝑝₂𝐶^𝑛+1_𝑚−1=𝑝₁+𝑘²(1−𝜃)𝛼²𝑝₂𝐶^𝑛_𝑚−3+𝑝₃−𝑘²(1−𝜃)𝛼²𝑝₄𝐶^𝑛_𝑚−2+𝑝₁+𝑘²(1−𝜃)𝛼²𝑝₂𝐶^𝑛_𝑚−1, (26) where 𝑝₁=16sin²ℎ2ℎ2cos2+cos3ℎ2×ℎ1+2cos2,𝑝₂=3sin(ℎ)sin3ℎ2ℎ1+3cos2cosec²ℎ2×ℎ1+2cos2,𝑝₃=32sin(ℎ)sin3ℎ2ℎ2cos2+cos3ℎ2,𝑝₄=24cot²ℎ2sinℎsin3ℎ2ℎ2cos2+cos3ℎ2. (27) Simplifying it leads to 𝑤₁𝐶^𝑛+1_𝑚−3+𝑤₂𝐶^𝑛+1_𝑚−2+𝑤₁𝐶^𝑛+1_𝑚−1=𝑤₃𝐶^𝑛_𝑚−3+𝑤₄𝐶^𝑛_𝑚−2+𝑤₃𝐶^𝑛_𝑚−1, (28) where 𝑤₁=𝑝₁−𝑘𝜃𝛼²𝑝₂,𝑤₂=𝑝₃+𝑘𝜃𝛼²𝑝₄,𝑤₃=𝑝₁+𝑘(1−𝜃)𝛼²𝑝₂,𝑤₄=𝑝₃−𝑘(1−𝜃)𝛼²𝑝₄. (29) Now on inserting the trial solutions (one Fourier mode out of the full solution) at a given point 𝑥_𝑚𝐶^𝑛_𝑚=𝛿^𝑛exp(𝑖𝑚𝜂ℎ) into (28) and rearranging the equations, 𝜂 is the mode number, ℎ is the element size, and 𝑖=√−1, we get 𝑤₁𝛿^𝑛+1𝑒^{𝑖𝜂(𝑚−3)ℎ}+𝑤₂𝛿^𝑛+1𝑒^{𝑖𝜂(𝑚−2)ℎ}+𝑤₁𝛿^𝑛+1𝑒^{𝑖𝜂(𝑚−1)ℎ}=𝑤₃𝛿^𝑛𝑒^{𝑖𝜂(𝑚−3)ℎ}+𝑤₄𝛿^𝑛𝑒^{𝑖𝜂(𝑚−2)ℎ}+𝑤₃𝛿^𝑛𝑒^{𝑖𝜂(𝑚−1)ℎ}. (30) Dividing (30) by 𝛿^𝑛𝑒^{𝑖𝜂(𝑚−2)ℎ} and rearranging the equation, we get 𝛿𝑤₂+𝑤₁cos𝜂ℎ=𝑤₄+𝑤₃cos𝜂ℎ. (31) Let 𝐴=𝑝₃+𝑝₁cos𝜂ℎ,𝐵=𝑘𝛼²𝑝₄−𝑝₂cos𝜂ℎ. (32) Therefore, (31) can be written as 𝛿(𝐴+𝜃𝐵)−(𝐴−(1−𝜃)𝐵)=0. (33) Equation (33) can be rewritten as 𝛿𝑋₂+𝑖𝑌−𝑋₁+𝑖𝑌=0, (34) where𝑋₁=(()𝐵),𝑋𝐴−1−𝜃₂=(𝐴+𝜃𝐵),𝑌=0. (35) For stability, the maximum modulus of the eigenvalues of the matrix has to be less than or equal to one [47]. Since 𝐴>0, 𝐵>0, and 0≤𝜃≤1, we always have 𝜉²=𝑋²₁+𝑌²𝑋²₂+𝑌²≤1. (36) Thus, from (36), the proposed scheme for nonclassical diffusion equation (with term, 𝑞(𝑥,𝑡)=0) is unconditionally stable and it is also unconditionally stable with a general term 𝑞(𝑥,𝑡), since the modulus of the eigenvalues must be less than one [47]. We recall Duhamel’s principle ([15], chapter 9); a scheme is stable for equation 𝑃_𝑘,ℎ𝑣=𝑓 if it is stable for the equation 𝑃_𝑘,ℎ𝑣=0. This means that there are no constraints on grid size ℎ and step size in time level 𝑘, but we should prefer those values of ℎ and 𝑘 for which we obtain the best accuracy of the scheme.

5. Results and Discussions

In this section, the cubic trigonometric B-spline collocation method is employed to obtain the numerical solutions for one-dimensional nonclassical diffusion problem with nonlocal boundary constraints given in (1)–(3). Two numerical examples are discussed in this section to exhibit the capability and efficiency of the proposed trigonometric spline method. Numerical results are compared with existing methods in the literature and with the exact solution at the different nodal points 𝑥_𝑖 for some time levels 𝑡_𝑛 using some particular space step size ℎ and time step 𝑘. In order to calculate the maximum errors and relative 𝐿₂ error norms of the proposed method numerically, we use the following formulas: 𝐿_∞=max_𝑖𝑈_num𝑥_𝑖,𝑇−𝑢_exact𝑥_𝑖,𝐿,𝑇₂=𝑈_num(𝑥,𝑡)−𝑢_exact(𝑥,𝑡)𝑢_exact(.𝑥,𝑡) (37)

Example 1. Consider the nonclassical diffusion problem ((1)–(3)), with 𝑞(𝑥,𝑡)=−𝑒^{−(𝑥+sin𝑡)}()𝑔1+cos𝑡,0≤𝑥≤1,0<𝑡≤1,₁(𝑥)=𝑒^−𝑥𝜉,0≤𝑥≤1,_𝑖(𝑖=1,3)=1,𝜉_𝑖()𝑔𝑖=2,4=0,₂(𝑥)=𝑒𝑒−2𝑥,𝑔₃(𝑥)=2ℎ𝑒+sin1−cos1cos𝑥,𝐿=1,_𝑖(𝑡)(𝑖=1,2)=0,0<𝑡≤1. (38) This test problem is from Dehghan [8] and Li and Wu [10] and the known solution is 𝑢(𝑥,𝑡)=𝑒^{−(𝑥+sin𝑡)}. We compare the maximum errors with TMOL [10] when they are considered with space size ℎ=0.05 and errors are recorded with several values of time step 𝑘, given in Table 2 and also shown in Figure 1. It is worth noting that the results obtained using CuTBSM are more accurate as compared to TMOL [10]. We also compare the relative errors of numerical value of 𝑢(0.6,0.1) with different space step ℎ=0.05,0.025,0.01 and with time step size 𝑘=𝑠ℎ², 𝑠=1,2,4,5, instead of 𝑘=0.4ℎ² which was used in the BTCS [8], the implicit (3,3) Crandall’s formula [8], the 3-point FTCS [8], and the Dufort-Frankel three-level approach [8] and TMOL [10], and these results are tabulated in Table 3. It is clearly shown from this table that the obtained results by using CuTBSM are more precise as compared to methods in [8, 10]. Figure 2 shows the approximate solution and exact solution for this example at different time levels with ℎ=0.05 and 𝑘=0.01.

Table 2 

Absolute errors for Example 1 with ℎ=0.05 and several values of time step by present method (CuTBSM) and compare with TMOL [10].

Table 3 

Relative errors of numerical method at 𝑢 (0.6, 0.1) for Example 1 with 𝑘=0.4ℎ² and various values of space step.

Figure 1 

Absolute errors of numerical values at different time levels for Example 1 by using ℎ=0.05 and 𝑘=0.001.

Figure 2 

A comparison of numerical solution and known solution at different time levels for Example 1 with ℎ=0.05 and 𝑘=0.01.

Example 2. We consider another numerical test problem, with 𝑞()=𝑥𝑥,𝑡−2²+𝑡+1(𝑡+1)³𝑔0≤𝑥≤1,0<𝑡≤1,₁(𝑥)=𝑥²𝜉,0≤𝑥≤1,_𝑖(𝑖=1,3)=1,𝜉_𝑖(𝑖=2,4)=0,𝑔₂(𝑥)𝑔=𝑥,₃(𝑥)ℎ=𝑥,𝐿=1,0≤𝑥≤1,₁(𝑡)=−14(𝑡+1)²,ℎ₂(𝑡)=34(𝑡+1)²0<𝑡≤1. (39) The exact solution of this equation is 𝑢(𝑥,𝑡)=𝑥²/(𝑡+1)² and this test problem has been taken from [8, 10, 11]. This problem is tested using different values of ℎ and 𝑘 to show the capability of the present method for solving nonclassical diffusion equation ((1)–(3)). The final time is taken 𝑇=1.0. The maximum errors of the numerical method are calculated at different time levels with different time step size and it is observed that they are more accurate as compared to TMOL [10] and Chebyshev spectral collocation method (CSCM) based on Chebyshev polynomials [11]. The numerical errors are tabulated in Tables 4 and 5 and are also depicted graphically in Figure 3. The relative errors of numerical value 𝑢(0.6,1.0) with different space step ℎ=0.05,0.025,0.01 and with time step size 𝑘=𝑠ℎ², 𝑠=1,2,4,5, instead of 𝑘=0.4ℎ² which was used in the BTCS [8], the implicit (3,3) Crandall’s formula [8], the 3-point FTCS [8], and the Dufort-Frankel three-level approach [8] and TMOL [10] and they are recorded in Table 6. It is worth noting that numerical results are much better than the methods in [8, 10, 11]. A comparison of numerical solutions at different time levels with known solution is presented graphically in the Figure 4.

Table 4 

Absolute errors for Example 2 with ℎ=0.05 and several values of time step by present method and TMOL [10].

Table 5 

Absolute errors of numerical value of 𝑢 for Example 2 with different space step size, 𝑘=0.0001, at selected time levels and compare with CSCM [11].

Table 6 

Relative errors of numerical method at 𝑢 (0.6, 1.0) for Example 2 with various values of space step.

Figure 3 

Absolute errors of numerical values at different time levels for Example 2 by using ℎ=0.05 and 𝑘=0.001.

Figure 4 

A comparison of numerical solution and known solution at different time levels for Example 2 with ℎ=0.05 and 𝑘=0.01.

Example 3. Finally, we consider nonclassical diffusion problem ((1)-(3)), with 𝑞(𝑥,𝑡)=−𝑒^𝑥(1+𝑡)²1+𝑡²²𝑔,0≤𝑥≤1,0<𝑡≤1,₁(𝑥)=𝑒^𝑥𝜉,0≤𝑥≤1,_𝑖(𝑖=1,3)=1,𝜉_𝑖(𝑖=2,4)=0,𝑔₂(𝑥)=(𝑥+1)𝑒,𝑔₃(𝑥)ℎ=0,𝐿=1,0≤𝑥≤1,₁(𝑡)=0,ℎ₂(𝑡)=𝑒1+𝑡²0<𝑡≤1. (40) This test problem has been taken from [15] and its exact solution is 𝑢(𝑥,𝑡)=𝑒^𝑥/(1+𝑡²). The final time is taken as 𝑇=1. The maximum errors of the proposed scheme are considered at 𝑇=1 with different time step sizes that are depicted graphically in Figure 5. The relative errors of numerical value 𝑢(0.6,1.0) are calculated with different space size step with 𝑘=𝑠ℎ², 𝑠=1,2,4,5, and they are given in Table 7. It is worth noting that the numerical results are found to be in good agreement with exact solutions. A comparison of numerical solutions at different time levels with known solution is presented graphically in Figure 6.

Table 7 

Relative errors of numerical method at 𝑢 (0.6, 1.0) for Example 3 with different time step size and various values of space step.

Figure 5 

Absolute errors of numerical values at 𝑇=1 for Example 3 with ℎ=0.01 and different time step size.

Figure 6 

A comparison of numerical solution and known solution at different time levels for Example 3 with ℎ=0.05 and 𝑘=0.01.

6. Concluding Remarks

In this paper, a new two-time level implicit scheme based on cubic trigonometric B-spline has been used to solve the nonclassical diffusion problem with known initial and with nonlocal boundary constraints instead of the usual boundary constraints. A usual finite difference discretization is used for time derivatives and cubic trigonometric B-spline is applied for space derivatives. It is noted that the accuracy of solution may reduce as time increases due to the time truncation errors of time derivative term [47]. The cubic trigonometric B-spline method used in this paper is simple and straightforward to apply. An advantage of using the cubic trigonometric B-spline method outlined in this paper is that it produces a spline function on each new time line which can be used to obtain the solutions at any intermediate point in the space direction whereas the finite difference approach yields the solution only at the selected points. The CuTBSM has approximated the solution with more accurate results for time step size 𝑘=𝑠ℎ², 𝑠=1,2,4,5, as compared to some finite difference schemes with smaller time step size 𝑘=0.4ℎ² such as BTCS, Crandall’s formula, FTCS, the Dufort-Frankel scheme, and TMOL based on reproducing kernel. The proposed method is shown to be unconditionally stable. It is also evident from the examples that the approximate solution is very close to the exact solution.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors are grateful to the anonymous reviewers for their helpful, valuable comments, and suggestions to improve the paper. This study was fully supported by FRGS Grant no. 203/PMATHS/6711324 from the School of Mathematical Sciences, Universiti Sains Malaysia, Penang, Malaysia. The first author was supported by Post Doctorate Fellowship from School of Mathematical Sciences, Universiti Sains Malaysia, Penang, Malaysia, during a part of the time in which the research was carried out.