A Comparative Numerical Study of Parabolic Partial Integro-Differential Equation Arising from Convection-Diffusion

Kamil Khan,Arshed Ali,＊,Fazal-i-Haq,Iltaf Hussain and Nudrat Amir

1Department of Mathematics,Islamia College,Peshawar,25000,Pakistan

2Department of Mathematics,Statistics and Computer Science,The University of Agriculture,Peshawar,25000,Pakistan

3Department of Basic Sciences&Islamiat,University of Engineering&Technology,Peshawar,25000,Pakistan

4Department of Basic Sciences&Humanities,CECOS University of Information Technology and Emerging Sciences,Peshawar,25000,Pakistan

ABSTRACT This article studies the development of two numerical techniques for solving convection-diffusion type partial integro-differential equation(PIDE)with a weakly singular kernel.Cubic trigonometric B-spline(CTBS)functions are used for interpolation in both methods.The first method is CTBS based collocation method which reduces the PIDE toan algebraic tridiagonal system of linear equations.The other method is CTBSbased differential quadrature method which converts the PIDE to a system of ODEs by computing spatial derivatives as weighted sum of function values.An efficient tridiagonal solver is used for the solution of the linear system obtained in the first method as well as for determination of weighting coefficients in the second method.An explicit scheme is employed as time integrator to solve the system of ODEs obtained in the second method.The methods are tested with three nonhomogeneous problems for their validation.Stability,computational efficiency and numerical convergence of the methods are analyzed.Comparison of errors in approximations produced by the present methods versus different values of discretization parameters and convection-diffusion coefficients are made.Convection and diffusion dominant cases are discussed in terms of Peclet number.The results are also compared with cubic B-spline collocation method.

KEYWORDS Partial integro-differential equation;convection-diffusion;collocation method;differential quadrature;cubic trigonometric B-spline functions;weakly singular kernel

1 Introduction

Partial integro-differential equations have been applied to model different physical phenomena in science and engineering such as heat conduction [1],reactor dynamics [2],flow in fractured biomaterials [3],electricity swaptions [4],financial option pricing [5],viscoelasticity [6],population dynamics [7],and convection-diffusion [8-11].

Convection-diffusion integro-differential equation has been used to model several important physical phenomena such as heat and mass transfer,flows in porous media,current density in fluids,and pollutant transport in atmosphere,streams,rivers and oceans (see [8,10,11] and the references therein).

PIDE models have been rarely solved analytically and their general solution is only obtained under restrictive conditions [12].Due to such restrictive conditions the models become over simplified and also limit their physical relevance.Therefore,alternative approaches have been obtained for the solution of PIDEs including finite difference methods [13,14],finite element method [15],radial basis function collocation methods [16],spectral method [17],quasi-wavelet method [18],and spline methods [19-21].

This paper is devoted to comparative study of two efficient numerical techniques for solution of convection-diffusion parabolic PIDE [8,10] defined by:

with initial condition

boundary conditions

where the integral term is known as memory term,g(x),g1(t),g2(t) are known functions,andK(q,t)=(t−q)−μ,0<μ<1,is weakly singular kernel.The parameterscanddare positive constants which represent quantities of the convection and diffusion processes respectively.Different numerical techniques are developed to obtain the solution of the PIDE (1).Siddiqi et al.[8]employed cubic b-spline functions to spatial derivatives and Euler backward formula to time derivative to solve the PIDE (1).Ali et al.[9] constructed quartic B-spline collocation technique to obtain the solution of the PIDE (1).Fahim et al.[10] have solved (1) using product trapezoidal integration rule and sinc collocation method.Recently,Al-Humedi et al.[11] have solved (1) using cubic B-spline Galerkin method with quadratic weight function.

The cubic trigonometric B-spline methods are numerically accurate,computationally fast,consistent,and has ability to get the approximate solution at any point in the domain with more accuracy as compared to the conventional finite difference method which approximates solution at only selected points in the domain.In recent years solution of several partial differential equations have been found using these methods [22-28] due to special geometric features of trigonometric splines [22,29-31].In 1964,Schoenberg [29] introduced piecewise trigonometric spline functions and founded the theory of locally supported trigonometric splines,named as trigonometric B-splines.It was derived that trigonometric spline can be written as linear combination of trigonometric B-splines.Lycche et al.[32] derived a recurrence relation using trigonometric B-splines divided differences for stable computation of trigonometric splines.They also formulated derivatives of trigonometric B-splines and provided an integral representation of trigonometric divided differences with trigonometric B-spline as kernel.Later on several other researchers have enhanced the theory of trigonometric B-splines by further developing many important properties such asC2continuity,nonnegativity,partition of unity,smoothness,curve control and its analysis,and banded interpolation matrix (see [22,30,31] and the references therein).The first application of cubic trigonometric B-spline functions appeared in the literature in 2010 [33] when these functions were used for the solution of Two-point boundary value problem.In 2014,CTBS functions were extended for the solution of initial and boundary-value hyperbolic problems [22].Recently,CBTS functions have been employed for the solution of PIDE problem [19,34].The CTBS functions are also employed for the solution different problems which include Hunter Saxton equation [25],time fractional diffusion-wave equation [26],Fisher’s reaction-diffusion equations [28],non-conservative linear transport problems [27],Nonclassical diffusion problems [35],Hyperbolic telegraph equation [36],Fisher’s equations [37],and coupled Burgers’equations [38].

In 1937,Frazer et al.[39] initiated the use of collocation method.This method is a special variant of the method of weighted residuals which uses dirac delta functions as weighting functions to minimize the residual.The main advantage of the collocation method is that as the number of collocation points increases,more points from the domain satisfy the functional equation.As a result approximate solution obtained by this method approaches to exact solution.Another advantage of this method is that using the dirac delta function makes computation of residual easy.The method is coupled with other methods such as Galerkin method,least squares method,Newton’s method,and finite difference method for the solution of various problems(see [22,40,41] and the references therein).In this paper,the first technique is developed by coupling CTBS functions with collocation procedure (CTBS-CT) for the solution of problem (1).The main advantage of this technique is that it reduces (1) to a tridiagonal system of linear equations which can be solved by an efficient tridiagonal solver.

The differential quadrature (DQ) method was projected as strong alternative of conventional methods i.e.,finite difference and finite element methods.This method produces accurate numerical results with little computational effort as it uses smaller set of grid points while the conventional methods need relatively larger sets [42,43].The DQ method was first introduced by Bellman et al.[44] in 1971,for approximating derivative of a sufficiently smooth function.The basic structure of DQ method is that it uses a weighted linear sum of functional values to approximate the function derivatives.The weighting coefficients play a key role since accuracy of the approximation depends on the accuracy of these coefficients.Initially Bellman and his co-researchers have developed two methods for determination of weighting coefficients [45].However,the methods lead to ill-conditioned algebraic system for larger set of grid points.Shu [46]presented a systematic theoretical analysis and applications of earlier work about this method,and introduced two explicit approaches (i) a recurrence formula,and (ii) matrix multiplication,for computing the weighting coefficients of higher order derivatives to avoid the ill-conditioned system.The selection of grid points also effects accuracy of the DQ method and better accuracy was achieved using non-uniform and scattered grid points,and ghost points [34,43,47,48].The DQ method is extended for the solution of various problems including Fisher’s reaction-diffusion equations [28],coupled viscous Burgers’equations [38],space-dependent inverse heat problems [48],Kawahara equation [49],and shock wave simulations [50].

Recently,Korkmaz et al.[27] pioneered DQ method along with CTBS functions for the solution of second order advection-diffusion equation.Tasmir et al.[28] presented modified CTBS functions based DQ method for second order nonlinear Fisher’s reaction-diffusion equations.Kumar Singh et al.[38] extended the method for the solution of one and two dimensional coupled viscous Burgers’equations.Arora et al.[51] developed a hybrid CTBS based DQ method for the solution of nonlinear Burgers’equations.In this paper,the other technique is developed by coupling CTBS functions with differential quadrature (CTBS-DQ) approach.The key feature of this approach is that it reduces (1) to a system of ODEs which can be solved by any ODE solver.

Upto the best of our knowledge a comparative study of PIDE (1) using the proposed techniques has not been done.Moreover,solving PIDEs of the form (1) numerically have two major challenges in addition to discretization of temporal and spatial derivatives:the singular kernel and approximation of the memory term which also affect stability and accuracy of a numerical method.

Remaining work of the paper is arranged as follows:Section 2 shows development of the proposed methods through combination of CTBS functions with collocation and differential quadrature procedures.Section 3 analyzes stability of the present schemes.Section 4 shows numerical results of both methods including error norms,computational efficiency,conditioning,spectral radii,convergence and comparison with an existing method in order to establish the current approaches.Section 5 concludes the findings and outcomes of the paper.

2 The Proposed Methods

In order to develop,the proposed methods,consider the problem (1)-(3).The spatial domainΓis divided intoNsub-intervalsΓn=[xn−1,xn],n=1,2,...,N,of equal lengthby the nodes xn,n∈I={0,1,...,N}with a=x0and b=xN.Let tl=lδt,l=0,1,2,...,M,whereδt is temporal step.

2.1 CTBS-CT

Taking t=tl+1in Eq.(1) and the temporal derivative is approximated by Euler backward formula which reduces Eq.(1) to the following form:

The integral part of Eq.(4) is approximated as below ([9]):

whereuldenotesu(x,tl+1) andandbs→0 ass→∞.

Putting Eq.(5) in Eq.(4),we get

wherefl+1=f(x,tl+1).

Simplifying Eq.(6),we get

Now collocating Eq.(7) at x=xr,r∈I,we obtain

where=u(xr,tl+1).Next using CTBS functions Sr(x),given in [22,23] to approximateu(xr,t).

Let

where

Lemma[22,33]:The values of Sr(x),S′r(x) and S′′r(x) at the node xsare obtained as:

and

Then using Eq.(9),the functionuand its derivatives at therth node (x=xr,r∈I)can be expressed as

Thus by using (10) in (8) and combining like terms,we get

whereα1=[(1−b0δt)δ1+cδtδ3−dδtδ5],α2=[(1−b0δt)δ2−dδtδ6],α3=[(1−b0δt)δ1+cδtδ4−dδtδ5],and

Forr=0 andN,Eq.(11) reduces torespectively.After eliminating the parametersandusing (10) and the boundary conditions (3),we get the linear system ofN+1 equations inN+1 unknowns whose matrix form is as follows:

where

Cl+1=and

Eq.(12) can be re-written as

where D=A−1B and K=A−1F.The system Eq.(12) is solved using the efficient tridiagonal solver“Thomas Algorithm.”

2.2 CTBS-DQ

The DQ approach provideskth order derivative of the functionu(x,t) from the values ofu(x,t) at xr,r∈I1={1,2,...,N},as

where,k=1,2,...arekth order weighting coefficients (WCs) which are determined by test functions.Different basis functions have been employed as test functions for determination of the WCs such as CTBS functions [27,38],polynomials [42],radial basis functions [48],cubic B-spline functions [49],and sinc functions [50].In CTBS-DQ method we take the following test functions which are defined as modified CTBS functions [51]:

which form a basis in regionΓ.

Thus for each basis functionϖl(xr),Eq.(14) gives

which leads to the following matrix form:

The weighting coefficients,r,s∈I1,are obtained by solving the system (15) through the well known efficient tridiagonal solver “Thomas algorithm.” Now using (14) in (1) and taking x=xrto obtain the differential quadrature based technique,

The explicit time marching process for Eq.(16) is obtained as follows:

where ul=[ulr:r∈I1]T.Using Eq.(5) in Eq.(17) and combining the like terms,we get

Eq.(18) leads to the following matrix form:

where

3 Stability Analysis

In this section we assess stability of the present schemes (13) and (19).Leten,unandrepresent error,exact solution and approximate solution atnth time level respectively,thenen=un−.The error equation for the methods (13) and (19) is given by

where D is the amplification matrix.The schemes (13) and (19) will be stable if ‖D‖≤1 (see [52]and Lax-Richtmyer condition for stability).This condition is equivalent toρ(D)≤1,whereρ(D)represents spectral radius of the matrix D [53] and is defined asρ(D)=max1≤i≤N|λi|andλiis an eigenvalue of D.The values ofλidepend onN,δt,c,dandμ.Computational values ofρ(D) for different values of the parametersN,δt,c,dandμare provided in the next section which show that the schemes (13) and (19) satisfy the condition of stability.

4 Test Problems

In this section three examples are taken from the literature [8] to validate the present approaches (13) and (19) for the solution of the problem (1)-(3) and comparison is also made with results obtained by cubic B-spline collocation method (CBSC) [8].The methods are examined viaE∞,E2error norms [8].The schemes are also tested for convection and diffusion dominated problems via the Peclet numberHigh values ofPlead to convection dominated problem while lower values ofPlead to diffusion dominated problem [8].Computations are carried out with Corei3,2.4 GHz processor,4 GB Ram,Matlab7.5 and computer run time (RT) is provided in seconds.

4.1 Problem 1

We take Eqs.(1)-(3) with convection-diffusion coefficientsc=0.1,d=0.1 andelsewhere mentioned,and choosef(x,t) such that the exact solution is given by [8]

Table 1:Error norms for different values of convection and diffusion coefficients using δt=10−4,M=100,N=100

Table 2:Error norms vs.N using δt=10−4,M=100

Table 3:Error norms at various time levels M using δt=10−4 and N=100

Table 4:Convergence in space using δt=10−5 and M=100

Table 5:Convergence in time using N=100 and t=0.1

Table 6:Error norms vs. the parameter μ using δt=10−4,N=100 and M=100

Table 7:Comparison of error norms with CBSC method [8],δt=10−5

Figure 1:Numerical solutions obtained by (a) CTBS-DQ (b) CTBS-CT,using N=50,δt=10−4,M=1000,c=0.1,d=0.1,μ=1/4

Figure 2:Error in numerical solutions obtained by (a) CTBS-DQ (b) CTBS-CT,using N=50,δt=10−4,M=1000,c=0.1,d=0.1,μ=1/4

Figure 3:Numerical solutions obtained by (a) CTBS-DQ (b) CTBS-CT,over the time interval[0,1] using N=50,δt=10−4,c=0.1,d=0.1,μ=1/4

The conditionsg,g1,g2are found from Eq.(21).The error normsE∞,E2are obtained using the present methods (13) and (19) along with different values of the parameters number of nodesN,time stepδt,time levelM,convection and diffusion coefficientsc,d,andμwhich are reported in Tabs.1-6.The effect of the convection and diffusion coefficient parametersc,don the accuracy of both the methods are shown in Tab.1.Computational efficiency RT of the methods are given in Tabs.2 and 3 which shows that both the proposed methods are highly efficient.Numerical rate of convergence with respect toNandδt are provided in Tabs.4 and 5,and it can be seen that the CTBS-CT has higher rate of convergence than the CTBS-DQ.However,the rate of convergence of CTBS-DQ becomes higher than that of CTBS-CT for larger values of spatial nodesN.From Tabs.4-6,it can be noted that numerical values ofρ(D) for both methods satisfy the stability condition given in Section 3,whereasκ(A) regarding CTBS-CT remains reasonably small.From Tab.6,it can be noted that as the parameterμdecreases accuracy of the CTBS-DQ slightly increases as compared to the CTBS-CT.Furthermore,better accuracy of CTBS-CT than CBSC method [8] is evident from Tab.7.The solution of (1) for cases of convection dominant and diffusion dominant problems are also provided in this table.It can also be seen from these tables that the CTBS-CT gives better results than the CTBS-DQ.Fig.1 represents solutions obtained using CTBS-DQ and CTBS-CT atM=1000.Fig.2 displays errors in approximation obtained by the present methods.Fig.3 shows the CTBS based solutions at different time levels.Fig.4 represents condition numbers of system matrices of the methods.Fig.5 compares convergence in space for the proposed methods.

Figure 4:Condition numbers of system matrix A and interpolation matrix G corresponding to the methods (12) and (18) respectively,vs. the number of grid points N,for δt=10−4,c=0.1,d=0.1,μ=1/4

4.2 Problem 2

We take the Eqs.(1)-(3) with the convection-diffusion coefficientsc=0.5,d=0.005,andμ=1/3.f(x,t) is chosen so that the exact solution of (1)-(3) is given by [8]

Figure 5:Convergence of the approximate solution in space using δt=10−4,c=0.1,d=0.1,μ=1/4

Table 8:Error norms using M=100,μ=1/3

Table 9:Comparison of error norms with CBSC method [8] for N=100,P=1

The conditionsg,g1,g2are found from Eq.(22).The error norms for different values ofN,δt andMare computed which are reported in Tabs.8,9 alongwith the results of CBSC method [8].It can be noted from Tab.8 that errors produced by both the proposed methods decrease as number of nodesNincreases.Also both the methods provide considerably better accuracy than the method [8].Fig.6 represents solutions obtained using CTBS-DQ and CTBS-CT atM=1000.Fig.7 depicts the solution at different time levels.

Figure 6:Numerical solutions obtained by (a) CTBS-DQ (b) CTBS-CT,using N=50,δt=10−4,M=1000,c=0.5,d=0.005,μ=1/3

Figure 7:Numerical solutions obtained by (a) CTBS-DQ (b) CTBS-CT,over time interval [0,0.1]using N=50,δt=10−4,c=0.5,d=0.005,μ=1/3

4.3 Problem 3

For this example we take the convection-diffusion coefficientsc=0.5,d=0.001 andμ=1/4 for the problem (1)-(3) andf(x,t) is chosen so that the exact solution is given by [8]

Table 10:Error norms using M=100

Table 11:Comparison of error norms with CBSC method [8] using N=100,δt=10−5,P=5

Figure 8:Numerical solutions obtained by CTBS-CT,using N=50,δt=10−4,c=0.5,d=0.001,μ=1/4 at t=0,0.1

The conditionsg,g1,g2are found from Eq.(23).The error norms for different values ofN,δt andMare computed which are provided in Tabs.10 and 11.It can be noted from Tabs.10 and 11 that CTBS-CT provided better results than CTBS-DQ.The results in Tab.11 correspond to convection dominated problem forP=5.Moreover,the results of the present methods are also compared with the method CBSC [8].Fig.8 shows solutions obtained using CTBS-CT at t=0,0.1.Fig.9 shows the solution obtained by CTBS-DQ at various time levels.

Figure 9:Numerical solutions obtained by CTBS-DQ at different times in the interval [0,0.1]using N=50,δt=10−4,c=0.5,d=0.001,μ=1/4

5 Conclusion

Two cubic trigonometric B-spline functions based methods are constructed for comparative study of approximate solution of a parabolic type partial integro-differential equation with a weakly singular kernel.The methods are validated using three test problems and the results are compared via error norms.The methods are computationally efficient and improved accuracy is obtained for relatively larger number of spatial nodes.Stability of the methods is shown via spectral radius.It is also observed that condition numbers of the system matrix obtained from cubic trigonometric B-spline collocation method does not increase with increase in the number of spatial nodes and hence the computation remained stable.Both methods do not require much smaller time step to attain high accuracy.In most cases,it is found that cubic trigonometric B-spline collocation method provides better accuracy as compared to cubic trigonometric B-spline differential quadrature.Also,in two examples the cubic trigonometric B-spline collocation technique provided better accuracy than cubic B-spline collocation method.Both techniques have the ability to solve convection dominated and diffusion dominated problems.Due to excellent agreement of these methods with the exact solution,the proposed techniques are efficient in employing to get approximate solution of PIDEs.

Acknowledgement:The authors are thankful to the reviewers for their helpful comments which improved the work of this paper.

Funding Statement:The author(s) received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.