A new numerical method for heat equation subject to integral speci…cations

(1)

A new numerical method for heat equation subject to integral speci…cations

H.M. Jaradat

^a

, M.M.M. Jaradat

^b;

, Zead Mustafa

^b

O. Alsayyed

^c

aDepartment of Mathematics, Al al-Bayt University, Jordan

Department of Mathematics and Applied Sciences, Dhofar University, Salalah, Oman

bDepartment of Mathematics, Statistics and Physics, Qatar University, Doha, Qatar

cDepartment of Mathematics, Hashemite University, Jordan

Abstract

In this research a numerical technique is developed for the one-dimensional heat equation that combine classical and integral boundary conditions. The combined Lablace transform, high-precision Quadrature schemes and Stehfest inversion algorithm have been propsed for the numerical solution of the problem. A Lablace transform method is introduced for the solution of considered equation, de…nite integrals approximated by high-precison Quadrature schemes. To invert the equation numerically back into the time domain we apply Stehfest inversion algorithm. The accuracy and computational e¢ ciency of the proposed method are veri…ed with the help of the numerical examples.

Key words: Heat equation; Nonlocal Boundary Value Problems; Laplace Inversion; High-Precision Quadrature Schemes; Stehfest inversion algorithm.

Introduction

In 1963, nonlocal boundary equation have been presented by Cannon [5], and Batten [3], independently. Then, parabolic initial-boundary problems with nonlocal integral conditions for parabolic equations were investigated by Kamynin

Corresponding author: Tel:+97455793995.

Email address: [email protected] (M.M.M. Jaradat).

(2)

[24] and Ionkin [23]. Over the previous years nonclassical problems for partial di¤erential equations have been widely used for a description a number of phenomena in modern physics and technology. Nonclassical problems with nonlocal conditions include relations between boundary values of an unknown solution and its derivatives and their values at internal points of a domain.

Nonlocal problems with integral conditions which are naturally generalization of discrete nonlocal conditions can be considered as mathematical models the processes with inaccessible boundary.

This paper is focused on the numerical solution of the following di¤usion equation

@u

@t

@²u

@x² +u=q(x; t); x2(0;1); 0< t < T; (1) with initial condition

u(x;0) =f(x); x2(0;1);0< t T; (2) and integral conditions

Z 1

0

u(x; t)dx=g₁(t); 0< t < T; (3)

Z 1

0

b(x)u(x; t)dx=g2(t); 0< t < T; (4) whereq(x; t); f(x); g₁(t); g₂(t); b(x)are known functions andT is a given con- stant. The mathematical modeling of this type of problems are encountered in heat transmission theory, in thermoelasticity, and in plasma physics [32,33]

and can be reduced to the nonlocal problems. Therefor, partial di¤erential equations with nonlocal boundary conditions have received much attention in last 20 years. However, most of the articles were directed to the second order parabolic equations, particularly to heat conduction equations. Recently parabolic equations with nonlocal boundary conditions have been treated ex- tensively by …nite di¤erence methods, …nite element procedures, boundary element techniques, spectral schemes, Adomian decomposition method, and the semidiscretization procedures [12–15]. The numerical techniques developed in [16] are based on three-level explicit …nite di¤erence procedures. Ang [1] developed a numerical technique for solution of the studied model, in [30]

a di¤erent approach is used by using combined …nite di¤erence and spectral methods for solving the hyperbolic equation with integral condition. The proof of the existence, uniqueness and continuous dependence of the strong solution upon the data for an initial-boundary value problem and integral conditions for this problem is studied by Bouziani [4]. The theoretical discussion of these case of equations can be found in [6,7]. The famous work of Lin [26] was one of the …rst to the solution of similar parabolic inverse problems. Dehghan applied some …nite-di¤erence schemes [17,18] and a shifted Tau method [19] for

(3)

solving similar problem. Similar form of this case of parabolic equations have been considered by various authors [8,9,21].

The purpose of present article is to give a method of solution to problem (1)- (4) using Laplace transform technique. Often the analytical inverse transform is too di¢ cult to …nd or evaluate in closed form. Numerical inversion methods are then used to overcome this di¢ culty. There are many approximate Laplace inversion algorithms. In this paper we will use the Stehfest inversion algorithm [34] in order to invert the Laplace transform. We use High-Precision Quadrature Schemes for numerical approximation of the integrals appear in the method.

The rest of this paper is organized as follow, in the next section, we introduce the High-Precision Quadrature Schemes which can be used to approximate de…nite integrals. Laplace transform technique to solve problem (1)- (4) is uses in Section 2. Moreover, we will get analytical solution representation using Stehfest inversion algorithm. Numerical results are given in Section 3.

Finally, Section 4 will be our conclusions of this paper.

1 The Quadrature Scheme

In order to introduce the method used, some preliminary explanations are needed. Let a …nite interval (a; b) be given, as well as the positive integer numbersmandn. We seth= (b a)=nandh_j =a+jh (andB_2i for Bernoulli numbers). From Euler-Maclaurin formula, assuming that the function has at least 2m+ 2 continuous derivatives, it follows

I =

Z b a

f(x)dx=h

Xn j=0

f(x_j) h

2(f(a) +f(b))

Xm i=1

h²ⁱB_2i

(2i)! f⁽²ⁱ ¹⁾(b) f⁽²ⁱ ¹⁾(a) E;

where

E = h^2m+2(b a)B_2m+2f^(2m+2)( ) (2m+ 2)!

for some 2(a; b). For more details one can refer to [2,35–37].

Transforming the integral of f(x) on the interval [ 1;1] to an integral on ( 1;1) can be done using a change of variablex=g(t). In this case we can write, forh >0,

(4)

I =

Z 1

1

f(x)dx=

Z ₁

1

f(g(t))g⁰(t)dt=h

X1 j= 1

w_jf(x_j) +E =I_h+E; (5)

where x_j = g(h_j) and w_j = g⁰(h_j). We truncate the in…nite summation (5) into a …nite one as

I_h^(N) =h

N⁺

X

j= N

w_jf(x_j):

HereN =N +N⁺+ 1 is the number of the sampling points actually used.

For the integral

I =

Z 1

1

f(x)dx;

the transformation

x=g(t) = tanh

2sinht gives a double exponential formula

I_h = 2h

X1 j= 1

f tanh

2sinhhj coshhj cosh² ₂ sinhhj :

Note that the abscissasx_j and the weightsw_j can be computed only one time for a givenh, and then used for other problems. Typically one selectsh = 2 ^m, for some m. It is found that m = 12is more than su¢ cient to evaluate most integrals to 500-digit accuracy . One typically proceeds one “level”at a time, where level k usesh= 2 ^k, starting with level one and continuing until either a fully accurate result has been obtained or the …nal (m-th) level has been completed. For more details refer to [2,36].

2 Analysis of the method

2.1 Method based on Laplace transform

First we take a Laplace transform on both sides of (1)-(4) with respect to t, we get

d²

dx²U(x; s) (s+ 1)U(x; s) = f(x) Q(x; s); (6)

Z 1

0

U(x; s)dx=G₂(s); (7)

Z 1

0

b(x)U(x; s)dx=G₂(s); (8)

(5)

whereU(x; s) =L[u(x; t)],Q(x; s) = L[q(x; t)],G₁(s) = L[g₁(t)]andG₂(s) = L[g₂(t)]. We have now a boundary value problem governed by a second order inhomogeneous ordinary di¤erential equation. The general solution of (6) can be given by

U(x; s) =C₁(s)e^p^s+1x+C₂(s)e ^p^s+1x 1 ps+ 1

Z x 0

[f( ) +Q( ; s)] (9) sinh(p

s+ 1 (x ))d :

Substituting (9) into the integral conditions (7)-(8), we have C1(s)^he^p^s+1 1ⁱ C2(s)^he ^p^s+1 1ⁱ=p

s+ 1G1(s) (10)

+ 1

ps+ 1

Z 1

0

[f( ) +Q( ; s)][cosh(p

s+ 1 (1 )) 1]d :

and

C1(s)

Z 1

0

b(x)e^p^s+1xdx +C2(s)

Z 1

0

b(x)e ^p^s+1xdx =G2(s) (11)

+ 1

ps+ 1

Z 1

0

[f( ) +Q( ; s)]

Z 1

b(x) sinh(p

s+ 1 (x ))dx d :

Solving (10)-(11) for C₁(s)and C₂(s) we have

C₁(s) C₂(s)

!

=

0

B@a₁₁(s) a₁₂(s) a₂₁(s) a₂₂(s)

1 CA

1

b₁(s) b₂(s)

!

; (12)

where

a₁₁(s) = e^p^s+1 1; a₁₂(s) = 1 e ^p^s+1; (13) b₁(s) = p

s+ 1G₁(s) + 1 ps+ 1

Z 1

0

[f( ) +Q( ; s)] (14) [cosh(p

s+ 1 (1 )) 1]d ; a₂₁(s) =

Z 1

0

b(x)e^p^s+1xdx; a₂₂(s) =

Z 1

0

b(x)e ^p^s+1xdx; (15) b₂(s) = G₂(s) + 1

ps+ 1

Z 1

0

[f( ) +Q( ; s)] (16)

Z 1

b(x) sinh(p

s+ 1 (x ))dx]d :

Thus, to …nd out the solution in Laplace domain one has to evaluate all the integrals appear in (13)-(16). Using High-Precision Quadrature Schemes we

(6)

have the following approximations of the above integrals

Z 1

0

K(x)e ^p^s+1xdx= 1 2

Z 1

1

K 1

2(x+ 1) e ^p^s+1(¹₂^(x+1))dx ' 2h

X1 j= 1

w_jK 1

2(x_j+ 1) e ¹²^p^s+1(x^j⁺¹⁾;

Z 1

0

F( ; s)[cosh(p

s+ 1 (1 )) 1]d = 1 2

Z 1

1

F(1

2( + 1); s) [cosh(p

s+ 1 1 1

2( + 1) ) 1]d

= 1 2

Z 1

1

F(1

2( + 1); s) [cosh(p

s+ 1 1

2(1 ) ) 1]d ' 2h

X1 j= 1

w_jF 1

2(x_j + 1); s x (cosh(p

s+ 1 1

2(1 x_j) ) 1) and

Z 1

0

F( ; s)

Z 1

K(x) sinh(p

s+ 1 (x ))dx d ' 2h

X1 j= 1

w_jF 1

2(x_j + 1); s 1 ¹₂(x_j+ 1)

2

!

h

X1 i= 1

w_iK 1 ¹₂(x_j + 1)

2 x_i+1 + ¹₂(x_j + 1) 2

!

sinh(p

s+ 1 1 ¹₂(x_j+ 1)

2 x_i+ 1 + ¹₂(x_j+ 1) 2

1

2(x_j + 1)

!

;

where the nodesx_j and the weightsw_j are given by x_j = tanh

2sinhhj , w_j = coshhj

cosh² ₂ sinhhj : (17)

2.2 Laplace Inversion Schemes

The Gaver-Stehfest algorithm for numerical inversion of Laplace transform was developed in the late 1960s. Due to its simplicity and good performance it is becoming increasingly more popular in such diverse areas as geophysics, operations research, economics, …nancial and actuarial mathematics, computational physics, engeneering and chemistry [20,22].

(7)

Assume that f : (0;1) ! R is a locally integrable function, such that its Laplace transform

L[f(t)] =F(s) =

Z ₁

0

f(t)e ^stdt;

is …nite for all s >0. The problem consists in recovering the original function f(t) given that we knowF(s). This problem has numerous applications, and it has attracted a lot of attention from researchers over the last fty years (see [10] for an up-to-date exposition of this area). The exact inversion is normally di¢ cult to carry out, so approximate inversion techniques are used. There are many approximate Laplace inversion algorithms. (For more details see [28,29,31,38]).

The numerical inversion of Laplace transform arises in many areas of science and engineering. Stehfest [34] derived the Gaver-Stehfest algorithm for the numerical inversion of Laplace transforms. For most of the more interesting problems, however, numerical inverting often has numerical accuracy problems [11,25,27]. As such, small rounding errors in computation may signi…cantly o¤set the results, rendering these algorithms impractical to apply.

The Gaver-Stehfest method uses the summation:

f(t) = ln 2 t

XN n=1

nF nln 2 t

!

: (18)

The n coe¢ cients only depend on the number of expansion terms, N (which must be even), they are:

n = ( 1)^n+N=2

minfXn;N=2g

k=b(n+1)=2c

k^N=2(2k)!

(N=2 k)!k! (k 1)! (n k)! (2k n)!: (19) The n coe¢ cients become very large and alternate in sign when increasing n. The precision of the Stehfest inversion method depends on the Stehfest numberN. Indeed, one can see in equation 18 that the inversion is based on a summation ofN weighted values. The default Stehfest number is often chosen in the range6 N 18.

Taking this into account, we can obtain the solution to problem (1)-(4) as:

u(x; t) = ln 2 t

XN n=1

n[C₁ep_n_{ln 2}

t +1x+C₂e p_{nln 2}

t +1x

pt pnln 2 +t

Z x 0

[f( ) +Q( ;nln 2

t )] sinh(

snln 2

t + 1 (x ))d ];

(8)

where n is given by (19).

3 Applications

In this section, we illustrate e¢ ciency and accuracy of the presented method by the following numerical examples.

Example 1 Consider the heat equationn

@u

@t

@²u

@x² +u= 2t+t²+x; x2(0;1); 0< t < T;

with initial condition

u(x;0) =x; x 2(0;1);0< t T;

and integral conditions

Z 1

0

u(x; t)dx= 1

2 +t²; 0< t < T;

Z 1

0

xu(x; t)dx= 1 3+ 1

2t²; 0< t < T:

We can verify that the exact solution of this problem is u(x; t) = x+t². The absolute errors in the approximation are shown in Fig. 1.

@u

@t

@²u

@x²+u= 11+6t+11t² x 2xt xt² 4x² 8tx² 4x²t²; x2(0;1); 0< t < T;

u(x;0) = 3 x 4x²; x 2(0;1);0< t T;

Z 1

0

u(x; t)dx= 7

6(1 +t²); 0< t < T;

Z 1

0

(1 + 2x)u(x; t)dx= 3 2 +3

2t²; 0< t < T:

(9)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10^-8

10^-7 10^-6

x

error

Fig. 1. Absolute errors between numerical and exact solution in Example 1 for N = 16; t= 1; x2[0;1]:

The exact solution of this problem is u(x; t) = (1 +t²)(3 x 4x²). The absolute errors in the approximation are shown in Fig. 2.

@u

@t

@²u

@x² +u= (10 2x)e^t; x2(0;1); 0< t < T;

u(x;0) = 5 x; x2(0;1);0< t T;

Z 1

0

u(x; t)dx= 9

2e^t; 0< t < T;

Z 1

0

xu(x; t)dx= 13

6 e^t; 0< t < T:

The exact solution of this problem isu(x; t) = 5 x)e^t. The absolute errors in the approximation are shown in Fig. ??.

(10)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10^-7

10^-6 10^-5

x

error

Fig. 2. Absolute errors between numerical and exact solution in Example 2 for N = 16; t= 2; x2[0;1]:

4 Conclusions

In this article we presented a computational method for solving the parabolic heat equation with an integral condition. A Laplace transform method is introduced for the solution of considered equation. Then high-precision Quadrature schemes are used to approximate the resulting de…nite integrals. Stehfest inversion algorithm is applied to invert the equation numerically back into the time domain. The numerical results show that our new technique described in this paper is an accurate and reliable analytical technique worked very well for the studied problem. The new technique can be extended to high dimensional parabolic equations with integral conditions.

References

[1] W.T. Ang, A numerical method for the wave equation subject to a non-local conservation condition, Appl. Numer. Math., 56 (2006) 1054–1060.

[2] D.H. Bailey, X.S. Li, K. Jeyabalan, A Comparison of Three High-Precision Quadrature Programs, Experimental Mathematics, 14 (2005) 317-329.

(11)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10^-7

10^-6 10^-5 10^-4

x

error

Fig. 3. Absolute errors between numerical and exact solution in Example 3 for N = 16; t= 0:4; x2[0;1]:

[3] J.G. Batten, Second-order correct boundary conditions for the numerical solution of the mixed boundary problem for parabolic equations. Math Comput., 17(1963) 405-13.

[4] A. Bouziani and N.E. Benouar, Problème mixte avec conditions intégrales pour uneclasse d’équations paraboliques, Comptes Rendus de l’Académie des Sciences, Paris, t.321, Série I, (1995) 1177-1182.

[5] J. R. Cannon, The solution of the heat equation subject to the speci…cation of energy. Quart Appl Math., 21(1963) 155-60.

[6] J.R. Cannon, and Y. Lin, Determination of parameter p(t) in some quasilinear parabolic di¤erential equations, Inverse Problems., 4(1988) 35-45.

[7] J.R. Cannon, and Y. Lin, Determination of parameter p(t) in Holder classes for some semilinear parabolic equations, Inverse Problems 4(1988) 595-606.

[8] J.R. Cannon, Y. Lin, A Galerkin procedure for di¤usion equations with boundary integral conditions, Int. J. Eng. Sci., 28(1990) 579-587.

[9] J.R. Cannon, S. Prez-Esteva, J.A. van der Hoek, A Galerkin, procedure for the di¤usion equation subject to the speci…cation of mass, SIAM J. Numer. Anal., 24(1987) 499-515.

[10] A.M. Cohen. Numerical methods for Laplace transform inversion, volume 5 of Numerical Methods and Algorithms. Springer, 2007.

(12)

[11] B. Davies, B. Martin, Numerical Inversion of the Laplace Transform, A Survey and Comparison of Methods, Journal of Computational Physics 33(1) (1979) 132.

[12] M. Dehghan, Numerical solution of a parabolic equation with non-local boundary speci…cations, Appl. Math. Computat., 145 (2003) 185–194.

[13] M. Dehghan, The use of Adomian decomposition method for solving the one dimentional parabolic equation with non-local boundary spesi…cation, Int. J.

Comput. Math.,81 (2004) 25–34.

[14] M. Dehghan, E¢ cient techniques for the second-order parabolic equation subject to nonlocal speci…cations, Appl. Numer. Math., 52 (2005) 39–62.

[15] M. Dehghan, Solution of partial integro-di¤erential equation arising from viscoelasticity, Int. J. Comput. Math., 83 (2006) 123–129.

[16] M. Dehghan, On the solution of an initial-boundary value problem that combines Neumann and integral condition for the wave equation, Numer.

Methods Partial Di¤erential Equations, 21 (2005) 24–40.

[17] M. Dehghan, Numerical solution of one-dimential parabolic inverse problem, Appl. Math. Comp., 136(2003) 333-344.

[18] M. Dehghan, Parameter Determination in a Partial Di¤erential Equation from the Overspeci…ed Data, Math. Computer Modelling., 41(2005) 197-213.

[19] M. Dehghan, A. Saadatmandi, A tau method for the one-dimensional parabolic inverse problem subject to temperature overspeci…cation, Comput Math Appl., 52(2006) 9339.

[20] D.P. Gaver, Observing stochastic processes and approximate transform inversion, Oper. Res, 14 (1966) 444-459.

[21] A.B. Gumel, On the numerical solution of the di¤usion equation subject to the speci…cation of mass. J Aust Math Soc Ser B., 40(1999) 47583.

[22] H. Hassanzadeh, M. Pooladi-Darvish, Comparison of di¤erent numerical Laplace inversion methods for engineering applications, Appl. Math. Comp., 189 (2007) 1966-1981.

[23] N.I. Ionkin, Solutions of boundary value problem in heat conductions theory with nonlocal boundary conditions, Di¤erents. Uravn., 13(2)(1977) 294–304.

[24] L.I. Kamynin, A boundary value problem in the theory of the heat conduction with nonclassical boundary condition, Z. Vychisl. Mat. Fiz., 4(6)(1964) 1006- 1024.

[25] K.L. Kuhlman, Review of Inverse Laplace Transform Algorithms for Laplace- Space Numerical Approaches, Numerical algorithms 63(2) (2013) 339-355.

[26] Y. Lin, An inverse problem for a class of quasilinear parabolic equation, SIAM.

J. Math. Anal., 22(1991) 146-156.

(13)

[27] A. Murli, M. Rizzardi, Algorithm 682 Talbots Method for the Laplace Inversion Problem, ACM Transactions on Mathematical Software 16(2) (1990) 158-168.

[28] R. Piessens, New quadrature formulas for the numerical inversion of the Laplace transform, BIT, 9 (1969) 351-361.

[29] R. Piessens, Gaussian quadrature formulas for the numerical inversion of the Laplace transform. J. Engr. Math., 5 (1971) 1-9.

[30] M. Ramezani, M. Dehghan, M. Razzaghi, Combined Finite Di¤erence and Spectral Methods for the Numerical Solution of Hyperbolic Equation with an Integral Condition, Numerical Methods for Partial Di¤erential Equations, 24 (2007) 1-8.

[31] T. Sakurai, Numerical inversion of the Laplace transform of functions with discontinuities. Adv. Appl. Prob., 36 (2004) 616-642.

[32] A.A. Samarskii, Some problems in di¤erential equations theory, Di¤erents.

Uravn.,16 (1980), 1925-1935.

[33] P. Shi, M. Shilor, Design of contact patterns in one dimensional thermoelasticity, in Theorical Aspects of industrial design, Society of Industrial and Appl Math, Philadelphia 1992.

[34] H. Stehfest, Numerical inversion of the Laplace transform, Comm. ACM, 13(1) (1970) 47-49.

[35] H. Takahasi, M. Mori, Double exponential formulas for numerical integration, Pubi. RIMS Kyoto Univ., 9 (1974) 721-741.

[36] H. Takahasi, M. Mori, Error estimation in the numerical integration of analytic functions, Rep. Comput. Centre Univ. Tokyo, 3 (1970) 41-108.

[37] H. Takahasi, M. Mori,

Quadrature formulas obtained by variable transformation, Numer. Math., 21 (1973) 206-219.

[38] A. Talbot, The accurate inversion of Laplace transforms. J. Inst. Maths.

Applics., 23 (1979) 97-120.