- Email: [email protected]
Numerical solution of stochastic linear heat conduction problem by using new algorithms
Applied Mathematics and Computation 163 (2005) 97–106 www.elsevier.com/locate/amc
Numerical solution of stochastic linear heat conduction problem by using new algorithms K. Maleknejad *, F. Mirzaee Department of Mathematics, Iran University of Science and Technology, Narmak, Tehran 16844, Iran
Abstract We consider a stochastic linear heat conduction problem so, it is reduced to a special integral equation of the second kind. Our aim, in this paper is to give stable numerical algorithms for approximating solution of integral equation and heat conduction. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Stochastic heat conduction problem; Stochastic integral equations; Numerical stochastic
1. Introduction This paper deals with equations of the form 2 r u K 2 u ¼ W_ ðx; tÞ; ou ¼ 0; x 2 Ge Rd ; on S
ð1Þ
where W_ is a space–time noise and S is the boundary of Ge (region exterior to S and Gi is region interior) while n stands for the unit normally outwards to S, we suppose S be a smooth or regular boundary and K be a given number. Let l be be the measure dl dt on Rdþ1 , moreover a r-finite measure on Rd . Suppose l W be a Gaussian additive set function the Borel sets of Rdþ1 such that W ðAÞ is a
*
Corresponding author. E-mail addresses: [email protected] (K. Maleknejad), [email protected] (F. Mirzaee). 0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.10.054
98
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
ðAÞ such that Gaussian random variable of mean zero and variance l A \ B ¼ U, then W ðAÞ and W ðBÞ are independant. We call W a white noise based on v or we say v is underlying measure of W . For x ¼ ðx1 ; x2 ; . . . ; xd Þ 2 Rd define: Wxt ¼ W fðx1 ; 1Þ ðxd ; þ1Þ ð0; tÞg: Now, W is not differentiable, but its derivative exists as a distribution. Let odþ1 W W_ ¼ ox1
oxx;td ot. System (1) is very interested when d ¼ 3, since the problem embodies the solution of the heat conduction problem of an infinite expands of material containing cavity S on which the temperature gradient is zero. Transient heat conduction with white noise occurs in numerous engineering applications [1,2,10], and it is important to appreciate the different methods for dealing with it. We know, if geometrical complexities and/or the form of boundary conditions preclude their use, recourse must be made to finite-difference and finite-element methods, with the digital computer, such methods may be used to solve any conduction problems, regardless of complexity. Now, we propose new methods for doing and then many advantages will be obtained. Outline in this paper is as follows: in Section 2, we consider heat conduction problem for case d ¼ 3 and so that it reduces to an integral equation of second kind. In Section 3, we want to solve integral equation by circulant integral operators then those will be constructed. Finally, we survey the stochastic numerical methods by using preconditioned conjugated gradient (PCG) then, various algorithms are given.
2. Integral equation formulation and smoothness of solution We assume that u can be represented as the sum of a volume potential and a single-layer potential as follows: 0
4puðx Þ ¼
Z
W_ ðn0 ÞEðx0 ; n0 ÞdV þ Ge
Z
rðn0 ÞEðx0 ; n0 ÞdS;
ð2Þ
S
where the first integral is called the volume potential or Newtonian potential and the second integral is the simple or single-layer potential [10,14,17]. r is0 an n0 jÞ unknown source density and so, x0 ¼ ðx; tÞ, n0 ¼ ðn; tÞ and Eðx0 ; n0 Þ ¼ expðKjx jx0 n0 j is a fundamental solution of the Helniholtz equation, such that x, n 2 S. The next step is to take the normal derivative of both sides of Eq. (2), let x approach S, and use the single-layer potential properties and Newtonian potential.
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
99
Remark. By using the properties of single-layer we consider the derivative of u taken in direction of a line normal to the surfaces S in the outward direction form S. Then, we have Z ou cosðn x; nÞ ¼ 2prðpÞ þ rðQÞ dS ð3Þ 2 on pþ jn xj S and Z ou cosðn x; nÞ ¼ 2prðpÞ þ rðQÞ dS; 2 on p jn xj S
ð4Þ
where pþ and p signify that approach S from Gi , and Ge respectively such that, both x and n are on S, from these, we obtain the jump, the normal derivative, of 1 ou u across S; r ¼ 4p ðon jp ou j Þwhen on pþ Z r ds : ð5Þ u¼ jx nj S The result is satisfied by (3) and (4) we have Z Z 0 0 oEðx0 ; n0 Þ 0 oEðx ; n Þ 0 _ dV 2prðx Þ þ rðn0 Þ dS 0¼ W ðn Þ on on Ge S
ð6Þ
which is Fredholm integral equation of the second kind. Let S be chosen so that, (6) leads to an important integral equation as follows: Z rðxÞ ¼ aðx; yÞkðx yÞrðyÞdy þ w; ð7Þ G0
where G0 ¼ fx 2 Rdþ1 ; 0 6 xj < 1; j ¼ 1; 2; . . . ; d; d þ 1g, we assume integral operator be self-adjoint and also the operator is positive definite then a 2 C m ðG0 Þ; k 2 C ðm1Þ ð0; 1Þ;
mP1
ð8Þ
and the estimate jk ðm1Þ ðlÞj 6 Cjljb
ð0 < b < m;
C ¼ constÞ
ð9Þ
is valid for ð1; 1Þ n f0g. Then, we can suppose w be a wiener process to define on a probability space ðX; F; P Þ. Furthermore, we can easily write for noninteger b, this equality yields [4,5,9] jk ðiÞ ðtÞj 6 Ci ð1 þ jtj
bþm1k
Þ;
i ¼ 0; 1; . . . ; m 1:
ð10Þ
In particular, for b < m 1 the factor kðtÞ of the kernel can have only a jump discontinuity at t ¼ 0. For b > m 1 it can have an integrable power
100
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
singularity at t ¼ 0. For integer, b estimates (10) are valid for all derivatives except for the derivative of the order i ¼ m 1 b, for which (9) yields jk ðiÞ ðtÞj 6 Ci ðln jtj þ 1Þ;
i ¼ m 1 b:
ð11Þ
Let us introduce the Banach function space (where it isðmÞas a solution space for ðtÞj equipped r). Es ¼ Esb;m ¼ fr 2 C½0; 1Þ \ C m ð0; 1Þj sup0
ðtÞj with the norm krkEr ¼ max0
jrðkÞ ðtÞj 6 Ck ½tbþmk þ ðs tÞbþmk k ¼ j þ 1; . . . ; m: For integer bð1 6 b 6 m 1Þ, the inclusion r 2 Es implies that r 2 C j1 ½0; 1Þ \ C m ð0; 1Þ j ¼ m b P 1 and we have, jrðjÞ ðtÞj 6 Cj ðj ln tj þ ln j bþmk ðs tÞj þ 1Þ, jrðkÞ ðtÞj 6 Ck ðtbþmk þ ðs tÞ Þ k ¼ j þ 1; . . . ; m. Clearly, the smoothness of the solution to (7) is described by the following theorem (and [5,7,9]). Theorem. Let conditions (8) and (9) be satisfied and let w 2 Es . If Eq. (7) has an integrable solution r, then r 2 Es . Now, we study approximate solutions to Eq. (7) by using circulant integral operators in Es .
3. Circulant integral operator 3.1. Preliminary We consider Eq. (7) by the projection method [4,9] where the solution rðxÞ is approximated by rs ðxÞ of above equation Z rs ðxÞ ¼ aðx; yÞkðx; yÞrs ðyÞ þ w; ð12Þ 00
G
where G00 ¼ fx 2 G0 j 0 6 xj < s; j ¼ 1; 2; . . . ; d; d þ 1g. lims!1 krs rkLP ðG00 Þ ¼ 0 for 1 6 p 6 1. If we define R aðx; yÞkðx; yÞ½ dy; x 2 G00 ; G00 As ½ ¼ 0 x 62 G00 ;
It
is
clear
that
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
101
then As is a compact operator and As : LP ðG00 Þ ! Esb;m . Therefore the spectrum of the operator I As is clustered around 1 and moreover solving (12) by iterative methods such as the Conjugate Gradient (CG) method will be less expensive than direct methods [9,12,13]. For speeding up the convergence rate of the CG method is to apply a preconditioner to (12). Thus instead ðI þ Cs Þ1 ðI As Þrs ¼ ðI þ Cs Þ1 w:
ð13Þ
We will call the operator Cs a preconditioner for the operator As (or a circulant integral operator). A good preconditioner Cs is an operator that is close to As in some norm and the operator equation ðI þ Cs Þrs ¼ w
ð14Þ R
is easier to solve than (13). Let G ¼ fCs jCs ½ ¼ G00 hs ðx yÞ½ dy; x 2 G00 g be a class of preconditioned (where hs is periodic function in G00 ). Our scale is for choosing Cs as follows: (i) minimize kAs Cs k2 (here, we obtain ‘‘Optimal Circlant Integral’’); 1 2 (ii) minimize kI ðI þ Cs Þ ðI As Þk (here, we obtain ‘‘Super-Optimal Integral’’) on G with Hilbert–Schmidt norm. In [6,18] were constructed the circulant optimal and other circulant integral preconditioners when in Eq. (12) aðx; yÞ ¼ 1 and k was nonsingular, but here, we suppose aðx; yÞ 6¼ 1 and kernel function has singularity. 3.2. Construct optimal circulant integral If we call first minimizer the optimal preconditioner and denote it by pðAs Þ then, we construct the optimal circulant integral preconditioners pðAs Þ for integral operators As . The preconditioner pðAs Þ is defined to be the circulant integral operator that minimizes the Hubert–Schmidt norm Z sZ 2 kAs Hs k ¼ ðaðx; yÞkðx yÞ hs ðx yÞÞðaðx; yÞkðx yÞ 0
00
G
hs ðx yÞÞdy dx
ð15Þ
over all circulant integral operators Cs or over G. We first give the expression of the kernel function of pðAs Þ. Since a, k 2 Lp ðG00 Þ, we can write by using Fourier expansions [6] aðx; yÞKðx yÞ ¼
1 X m;n¼1
vm;n um ðxÞ un ðyÞ;
ð16Þ
102
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
where um is one of the eigenfunctions the operator Cs with eigenvalues as follows: ( m 2 Z; um ðxÞ ¼ p1ffis e2pimx=s ; pffiffiffi pffiffiffi R ð17Þ um ðxÞdx; m 2 Z: km ¼ sðhs ; um Þ ¼ s G00 hs ðxÞ (Remark. Then (14) can be solved efficiently by using the Fourier transforms.) Moreover, Z Z vm;n ¼ aðx; yÞkðx yÞun ðyÞ um ðxÞdx dy ¼ ðAs un ; um Þs ; ð18Þ G00
G00
where m, n 2 Z. By means of Fourier expansion, we can write hs ðx yÞ ¼
1 X
km um ðxÞ um ðyÞ;
x; y 2 G00 :
m¼1
Combining this with (18) and using the orthogonality of un , we can rephrase the distance (15) as 2
kAs Cs k ¼
þ1 X
1 X
2
jvm;m km j þ
m¼1
2
jvm;n j :
m;n¼1 m6¼n
Clearly, the expression becomes minimal if and only if km ¼ vm;m ¼ ðAs um ; um Þs for all integers m. 3.3. Construct super-optimal circlant integral This section will be an extended idea in [18]. We consider another type of circulant integral operators which are obtained by minimizing the Hilbert– Schmidt norm kI ðI þ Cs Þ1 ðI As Þk2
ð19Þ 1
over all circulant integral operators of G such that ðI þ Cs Þ exists. Let Cs be a circulant integral operator with kernel function and eigenvalues km given in (17). If I þ Cs P is invertible, then 1 þ kn 6¼ 0 for all integers n. Moreover, since hs 1 2 is in L2 ðG00 Þ, n¼1 ðkn Þ < 1 and therefore j1 þ kn j P 1=2 for all jnj sufficiently large. On the other hard, I þ Cs has inverse ðI þ Cs Þ1 ¼ I Ks , where Ks is also a circulant integral operator with kernel function ks ðx yÞ ¼
þ1 X n¼1
kn un ðxÞ un ðyÞ: 1 þ kn
ð20Þ
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
103
Because, the kernel function of Cs Kr at the point ðt sÞ is given by Z þ1 þ1 X X kn km um ðxÞ um ðy 0 Þ un ðy 0 Þun ðyÞdy 0 00 1 þ kn G m¼1 n¼1 ¼
1 X n¼1
k2n un ðnÞ un ðyÞ 1 þ kn
so, the function ks is a s-periodic function and it is straight forward to check that the kernel function of Cs Ks Cs Ks ¼ ðI þ Cs ÞðI Ks Þ I. Then, problem of minimizing the norm (19) becomes the problem of minimizing kI ðI Ks ÞðI þ As Þk over integral operator Ks . We asP all circulant 2 sume vn;m ¼ ðAs un ; um Þ and nm ¼ þ1 n¼1 jvn;m j . Then, kI ðI Ks ÞðI þ As Þk ¼ kKs þ As Ks As k and by using (20) the kernel function of Ks Ks As þ As at point ðx; yÞ is given by 1 X dm;n km km ; vm;n vm;n um ðxÞuðyÞ þ 1 þ km 1 þ km m;n¼1
1 X dm;n km vm;n ¼ un ðyÞ; um ðxÞ 1 þ km m;n¼1 where dm;n denotes the Kronecker symbol. By the definition of Hilbert–Schmidt norm 2 þ1 X dm;n km vm;n ; kI ðI Cs Þ1 ðI þ As Þk2 ¼ 1þk m m;n¼1 it is clear, that the above expression is minimized if and only if the term is minimized for all integers m, however, by differentiating this quotient with respect to the real and imaginary parts of km , we see that nm þvm;m minimum will be obtained if we set km ¼ 1þ for all m, hence we have obvm;m tained kernel function of super-optimal circulant Cs as follows: ! 1 X nm þ vm;m Cs ðx yÞ ¼ um ðyÞ: um ðxÞ 1 þ vm;m m¼1 Pþ1 jvm;n j2 þj1þvm;m j2 Moreover, since 1 þ km ¼ n¼1 1þvm;m 6¼ 0, we see that I þ Cs is invertible. jkm j2 km vm;n km vm;n þnm j1þkm j2
4. A survey of stochastic numerical methods by using PCG algorithms In (6), we saw that one of the integrals is a stochastic integral, so we can write [3,15,16]:
104
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
wN ;k ¼
N X
bkj : xj :n00 :
k ¼ 1; . . . ; N ;
j¼1
where wN ;k is a simulation of approximate vector or random vector for w and here collocations of random variables fx1 ; . . . ; xN g correspond to the process rE rn R 0 0 from W ¼ Ge W_ ðn0 Þ oEðxon;n Þ dV so, independent standard Gaussian random variables fn001 ; . . . ; n00N g are normalized increments of the martingale g at jump N moments and ðbkj Þk;j¼1 are determined matrix, moreover the symbol: xj :n00 denotes the wick product [16] of the random variables xj . Then, we can generate wN ;k , by using [3,8,15,16,19]. Now, we consider these algorithms, which we call AlglR, Alg2R and Alg3R. In all those algorithms, the integral, are discretized by the rectangular quadrature rule, that the discretization is done at integer grid points 1; . . . ; N where N ¼ s > 1. Thus, integral equation (12) is replaced by a system of linear equations. The algorithm AlglR is the basic CG algorithm [9,11–13] without a preconditioner. In algorithm Alg2R and Alg3R the preconditioners are obtained respectively by optimal circulant integral in Section 3.2 and superoptimal circulant integral in Section 3.3. For stopping criterion is krk k=kr0 k < 107 , where rk is the residual vector of PCG, CG after k iteration. We show AlglS, Alg2S and Alg3S when algorithms the integral are discretized by the Simpson’s 1/3 quadrature rule. On below, we first test the performance of integral equations with normally distributed random variables as white noise, that is means and variance 2 EðDwN ;i Þ ¼ 0, EððDwN ;i Þ Þ ¼ ni for the each iteration. Tables 1 and 2 show, the numbers of iteration required for convergence for different choice of algorithms for Examples 1 and 2. We see from the tables that the numbers of iterations of the preconditioned systems are significantly less than that of the nonpreconditioned system. Moreover, we see that the iterations performance of optimal circulant operators be smaller than super-optimal circulant operators. All computations were computed using programs written in the symbolic language mathematical version 2.1.
Table 1 Iteration counts N
Alg1R
Alg3R
Alg2R
Alg1S
Alg3S
Alg2S
20 40 60 80 100
18 27 42 56 83
9 12 14 21 23
8 11 12 14 19
18 25 40 51 79
10 13 14 20 22
9 11 13 15 18
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
105
Table 2 Iteration counts N
Alg1R
Alg3R
Alg2R
Alg1S
Alg3S
Alg2S
20 40 60 80 100
15 23 37 43 78
11 12 19 24 32
10 10 13 17 27
17 22 31 37 75
10 11 17 22 31
9 9 11 14 23
Example 1. In the singular integral equation of the second kind Z 1 Z 1 rffiffiffiffiffiffiffiffiffiffiffirffiffiffiffiffiffiffiffiffiffi 1þx 1þt rðx; tÞ ¼ logðjx tjÞrðx; tÞdx dt þ w; 1x 1t 1 1
ð21Þ
¼ jx þ tj so, in the two-dimensional case, the random where we Rassume in w, oE on 1 R1 _ variable 1 1 W ðxÞ oE dx dt is Gaussian centered with variance 105 and on covariance 106 , so it can easily be simulated (by the classical Monte Carlo method). The numerical solutions is given at the point t ¼ 0:25 in Table 1. Remark. We can approximate gðx; tÞ ¼ ð1 þ xÞð1 tÞrðx; tÞ by a finite Chebyshev series gðx; tÞ ffi
p X p X k¼0
akl Tl ðxÞTk ðtÞ:
l¼0
Then rðx; tÞ ffi
p X p X k¼0
l¼0
Z
1
Z
1
akl 1
1
Tl ðxÞTk ðlÞ logðjx tjÞdx dt pffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 1 t 2
and using 1 p
Z
1 1
logðjx tjÞTk ðtÞ dt pffiffiffiffiffiffiffiffiffiffiffiffi ¼ 1 t2
log 2 Tk ðxÞ=k
k¼0 kP1
and the orthogonality of the Chebyshev polynomials with respect to the weight 1=2 function ð1 x2 Þ , the integral (21) can be determined in closed form. On the theoretical grounds we would expect good results. ¼ ðx þ tÞ2 which variance 103 Example 2. In Example 1, we consider oE on and convariance is 104 , so at the point t ¼ 0:25, we have given resultants in Table 2.
106
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
References [1] S.K. Biswas, N.U. Ahmed, Stabilization of systems governed by the wave equation in the presence of distributed white noise, IEEE Trans. Autom. Contr. AC. 30 (10) (1985). [2] S.H. Crandall, W.Q. Zhu, Random vibration. A survey of recent development, J. Appl. Mech. 50 (1980) 953–962. [3] A.A. Dorogavtsev, An extended stochastic integral for smooth functionals of white noise, Ukrain. Math. J. 41 (1989) 1252–1258. [4] I. Goberg, I. Feldman, Convolution equation and projection methods for their solution, Trans. Math. Monogr. (1974) 41. [5] I. Goberg, N. Krupnik, in: One-Dimensional Linear Singular Integral Equations, Volumes I, II, Birkh€ auser Verlag, Basel, 1992. [6] I. Goberg et al., Fast preconditioned conjugate gradient algorithms, SIAM J. Matrix Anal. Appl. 31 (1994) 429–443. [7] M.A. Goberg, Numerical Solution of Integral Equations, Plenum Press, New York, 1990. [8] P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, SpringerVerlag, Berlin, 1995, Second corrected printing. [9] R. Kress, Linear Integral Equations, Springer-Verlag, Berlin, 1989. [10] L.S. Langston, Heat transfer from multidimensional objects using one-dimensional solutions for heat loss, Int. J. Heat Mass Transfer 25 (1982) 149–150. [11] K. Maleknejad, Numerical experimentation with time-series methods for convolution integral equation, J. Computat. Appl. Math. 31 (1990) 211–225. [12] K. Maleknejad, D. Rostami, Numerical solution of linear operator integral equation of the first kind by using methods of PCG, in: Lectures Note on Numerical Analysis, Polvdiv, Bulgarria, 4, 1995, pp. 111–122. [13] K. Maleknejad, D. Rostami, The PNCGC methods for solving some nonlinear boundary integral equation, Bull. of the Iranian Mathematics Society 26 (1) (2000) 7–29. [14] G.E. Myers, Analytical Methods in Conduction Heat Transfer, McGraw-Hill, New York, 1971. [15] T. Sekiguchi, V. Shiota, L2 -Theory of noncausal stochastic integrals, Math Rep. Toyama Univ. 8 (1985) 119–195. [16] A.V. Skorokhod, On a generalization of the stochastic integral, Theory Random Process. 7 (1979) 118–127. [17] P.J. Schneider, Conduction Heat Transfer, Addison-Wesley, Reading, MA, 1955. [18] E.E. Tyrtyshnikov, V.V. Strela, Which circulant preconditioner is better, Math. Comp. 65 (1996) 135–150. [19] R.J. Williams, K.L. Chung, Introduction to Stochastic Integration, Birkh€auser, Boston, 1990.
Numerical solution of stochastic linear heat conduction problem by using new algorithms K. Maleknejad *, F. Mirzaee Department of Mathematics, Iran University of Science and Technology, Narmak, Tehran 16844, Iran
Abstract We consider a stochastic linear heat conduction problem so, it is reduced to a special integral equation of the second kind. Our aim, in this paper is to give stable numerical algorithms for approximating solution of integral equation and heat conduction. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Stochastic heat conduction problem; Stochastic integral equations; Numerical stochastic
1. Introduction This paper deals with equations of the form 2 r u K 2 u ¼ W_ ðx; tÞ; ou ¼ 0; x 2 Ge Rd ; on S
ð1Þ
where W_ is a space–time noise and S is the boundary of Ge (region exterior to S and Gi is region interior) while n stands for the unit normally outwards to S, we suppose S be a smooth or regular boundary and K be a given number. Let l be be the measure dl dt on Rdþ1 , moreover a r-finite measure on Rd . Suppose l W be a Gaussian additive set function the Borel sets of Rdþ1 such that W ðAÞ is a
*
Corresponding author. E-mail addresses: [email protected] (K. Maleknejad), [email protected] (F. Mirzaee). 0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.10.054
98
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
ðAÞ such that Gaussian random variable of mean zero and variance l A \ B ¼ U, then W ðAÞ and W ðBÞ are independant. We call W a white noise based on v or we say v is underlying measure of W . For x ¼ ðx1 ; x2 ; . . . ; xd Þ 2 Rd define: Wxt ¼ W fðx1 ; 1Þ ðxd ; þ1Þ ð0; tÞg: Now, W is not differentiable, but its derivative exists as a distribution. Let odþ1 W W_ ¼ ox1
oxx;td ot. System (1) is very interested when d ¼ 3, since the problem embodies the solution of the heat conduction problem of an infinite expands of material containing cavity S on which the temperature gradient is zero. Transient heat conduction with white noise occurs in numerous engineering applications [1,2,10], and it is important to appreciate the different methods for dealing with it. We know, if geometrical complexities and/or the form of boundary conditions preclude their use, recourse must be made to finite-difference and finite-element methods, with the digital computer, such methods may be used to solve any conduction problems, regardless of complexity. Now, we propose new methods for doing and then many advantages will be obtained. Outline in this paper is as follows: in Section 2, we consider heat conduction problem for case d ¼ 3 and so that it reduces to an integral equation of second kind. In Section 3, we want to solve integral equation by circulant integral operators then those will be constructed. Finally, we survey the stochastic numerical methods by using preconditioned conjugated gradient (PCG) then, various algorithms are given.
2. Integral equation formulation and smoothness of solution We assume that u can be represented as the sum of a volume potential and a single-layer potential as follows: 0
4puðx Þ ¼
Z
W_ ðn0 ÞEðx0 ; n0 ÞdV þ Ge
Z
rðn0 ÞEðx0 ; n0 ÞdS;
ð2Þ
S
where the first integral is called the volume potential or Newtonian potential and the second integral is the simple or single-layer potential [10,14,17]. r is0 an n0 jÞ unknown source density and so, x0 ¼ ðx; tÞ, n0 ¼ ðn; tÞ and Eðx0 ; n0 Þ ¼ expðKjx jx0 n0 j is a fundamental solution of the Helniholtz equation, such that x, n 2 S. The next step is to take the normal derivative of both sides of Eq. (2), let x approach S, and use the single-layer potential properties and Newtonian potential.
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
99
Remark. By using the properties of single-layer we consider the derivative of u taken in direction of a line normal to the surfaces S in the outward direction form S. Then, we have Z ou cosðn x; nÞ ¼ 2prðpÞ þ rðQÞ dS ð3Þ 2 on pþ jn xj S and Z ou cosðn x; nÞ ¼ 2prðpÞ þ rðQÞ dS; 2 on p jn xj S
ð4Þ
where pþ and p signify that approach S from Gi , and Ge respectively such that, both x and n are on S, from these, we obtain the jump, the normal derivative, of 1 ou u across S; r ¼ 4p ðon jp ou j Þwhen on pþ Z r ds : ð5Þ u¼ jx nj S The result is satisfied by (3) and (4) we have Z Z 0 0 oEðx0 ; n0 Þ 0 oEðx ; n Þ 0 _ dV 2prðx Þ þ rðn0 Þ dS 0¼ W ðn Þ on on Ge S
ð6Þ
which is Fredholm integral equation of the second kind. Let S be chosen so that, (6) leads to an important integral equation as follows: Z rðxÞ ¼ aðx; yÞkðx yÞrðyÞdy þ w; ð7Þ G0
where G0 ¼ fx 2 Rdþ1 ; 0 6 xj < 1; j ¼ 1; 2; . . . ; d; d þ 1g, we assume integral operator be self-adjoint and also the operator is positive definite then a 2 C m ðG0 Þ; k 2 C ðm1Þ ð0; 1Þ;
mP1
ð8Þ
and the estimate jk ðm1Þ ðlÞj 6 Cjljb
ð0 < b < m;
C ¼ constÞ
ð9Þ
is valid for ð1; 1Þ n f0g. Then, we can suppose w be a wiener process to define on a probability space ðX; F; P Þ. Furthermore, we can easily write for noninteger b, this equality yields [4,5,9] jk ðiÞ ðtÞj 6 Ci ð1 þ jtj
bþm1k
Þ;
i ¼ 0; 1; . . . ; m 1:
ð10Þ
In particular, for b < m 1 the factor kðtÞ of the kernel can have only a jump discontinuity at t ¼ 0. For b > m 1 it can have an integrable power
100
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
singularity at t ¼ 0. For integer, b estimates (10) are valid for all derivatives except for the derivative of the order i ¼ m 1 b, for which (9) yields jk ðiÞ ðtÞj 6 Ci ðln jtj þ 1Þ;
i ¼ m 1 b:
ð11Þ
Let us introduce the Banach function space (where it isðmÞas a solution space for ðtÞj equipped r). Es ¼ Esb;m ¼ fr 2 C½0; 1Þ \ C m ð0; 1Þj sup0
ðtÞj with the norm krkEr ¼ max0
jrðkÞ ðtÞj 6 Ck ½tbþmk þ ðs tÞbþmk k ¼ j þ 1; . . . ; m: For integer bð1 6 b 6 m 1Þ, the inclusion r 2 Es implies that r 2 C j1 ½0; 1Þ \ C m ð0; 1Þ j ¼ m b P 1 and we have, jrðjÞ ðtÞj 6 Cj ðj ln tj þ ln j bþmk ðs tÞj þ 1Þ, jrðkÞ ðtÞj 6 Ck ðtbþmk þ ðs tÞ Þ k ¼ j þ 1; . . . ; m. Clearly, the smoothness of the solution to (7) is described by the following theorem (and [5,7,9]). Theorem. Let conditions (8) and (9) be satisfied and let w 2 Es . If Eq. (7) has an integrable solution r, then r 2 Es . Now, we study approximate solutions to Eq. (7) by using circulant integral operators in Es .
3. Circulant integral operator 3.1. Preliminary We consider Eq. (7) by the projection method [4,9] where the solution rðxÞ is approximated by rs ðxÞ of above equation Z rs ðxÞ ¼ aðx; yÞkðx; yÞrs ðyÞ þ w; ð12Þ 00
G
where G00 ¼ fx 2 G0 j 0 6 xj < s; j ¼ 1; 2; . . . ; d; d þ 1g. lims!1 krs rkLP ðG00 Þ ¼ 0 for 1 6 p 6 1. If we define R aðx; yÞkðx; yÞ½ dy; x 2 G00 ; G00 As ½ ¼ 0 x 62 G00 ;
It
is
clear
that
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
101
then As is a compact operator and As : LP ðG00 Þ ! Esb;m . Therefore the spectrum of the operator I As is clustered around 1 and moreover solving (12) by iterative methods such as the Conjugate Gradient (CG) method will be less expensive than direct methods [9,12,13]. For speeding up the convergence rate of the CG method is to apply a preconditioner to (12). Thus instead ðI þ Cs Þ1 ðI As Þrs ¼ ðI þ Cs Þ1 w:
ð13Þ
We will call the operator Cs a preconditioner for the operator As (or a circulant integral operator). A good preconditioner Cs is an operator that is close to As in some norm and the operator equation ðI þ Cs Þrs ¼ w
ð14Þ R
is easier to solve than (13). Let G ¼ fCs jCs ½ ¼ G00 hs ðx yÞ½ dy; x 2 G00 g be a class of preconditioned (where hs is periodic function in G00 ). Our scale is for choosing Cs as follows: (i) minimize kAs Cs k2 (here, we obtain ‘‘Optimal Circlant Integral’’); 1 2 (ii) minimize kI ðI þ Cs Þ ðI As Þk (here, we obtain ‘‘Super-Optimal Integral’’) on G with Hilbert–Schmidt norm. In [6,18] were constructed the circulant optimal and other circulant integral preconditioners when in Eq. (12) aðx; yÞ ¼ 1 and k was nonsingular, but here, we suppose aðx; yÞ 6¼ 1 and kernel function has singularity. 3.2. Construct optimal circulant integral If we call first minimizer the optimal preconditioner and denote it by pðAs Þ then, we construct the optimal circulant integral preconditioners pðAs Þ for integral operators As . The preconditioner pðAs Þ is defined to be the circulant integral operator that minimizes the Hubert–Schmidt norm Z sZ 2 kAs Hs k ¼ ðaðx; yÞkðx yÞ hs ðx yÞÞðaðx; yÞkðx yÞ 0
00
G
hs ðx yÞÞdy dx
ð15Þ
over all circulant integral operators Cs or over G. We first give the expression of the kernel function of pðAs Þ. Since a, k 2 Lp ðG00 Þ, we can write by using Fourier expansions [6] aðx; yÞKðx yÞ ¼
1 X m;n¼1
vm;n um ðxÞ un ðyÞ;
ð16Þ
102
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
where um is one of the eigenfunctions the operator Cs with eigenvalues as follows: ( m 2 Z; um ðxÞ ¼ p1ffis e2pimx=s ; pffiffiffi pffiffiffi R ð17Þ um ðxÞdx; m 2 Z: km ¼ sðhs ; um Þ ¼ s G00 hs ðxÞ (Remark. Then (14) can be solved efficiently by using the Fourier transforms.) Moreover, Z Z vm;n ¼ aðx; yÞkðx yÞun ðyÞ um ðxÞdx dy ¼ ðAs un ; um Þs ; ð18Þ G00
G00
where m, n 2 Z. By means of Fourier expansion, we can write hs ðx yÞ ¼
1 X
km um ðxÞ um ðyÞ;
x; y 2 G00 :
m¼1
Combining this with (18) and using the orthogonality of un , we can rephrase the distance (15) as 2
kAs Cs k ¼
þ1 X
1 X
2
jvm;m km j þ
m¼1
2
jvm;n j :
m;n¼1 m6¼n
Clearly, the expression becomes minimal if and only if km ¼ vm;m ¼ ðAs um ; um Þs for all integers m. 3.3. Construct super-optimal circlant integral This section will be an extended idea in [18]. We consider another type of circulant integral operators which are obtained by minimizing the Hilbert– Schmidt norm kI ðI þ Cs Þ1 ðI As Þk2
ð19Þ 1
over all circulant integral operators of G such that ðI þ Cs Þ exists. Let Cs be a circulant integral operator with kernel function and eigenvalues km given in (17). If I þ Cs P is invertible, then 1 þ kn 6¼ 0 for all integers n. Moreover, since hs 1 2 is in L2 ðG00 Þ, n¼1 ðkn Þ < 1 and therefore j1 þ kn j P 1=2 for all jnj sufficiently large. On the other hard, I þ Cs has inverse ðI þ Cs Þ1 ¼ I Ks , where Ks is also a circulant integral operator with kernel function ks ðx yÞ ¼
þ1 X n¼1
kn un ðxÞ un ðyÞ: 1 þ kn
ð20Þ
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
103
Because, the kernel function of Cs Kr at the point ðt sÞ is given by Z þ1 þ1 X X kn km um ðxÞ um ðy 0 Þ un ðy 0 Þun ðyÞdy 0 00 1 þ kn G m¼1 n¼1 ¼
1 X n¼1
k2n un ðnÞ un ðyÞ 1 þ kn
so, the function ks is a s-periodic function and it is straight forward to check that the kernel function of Cs Ks Cs Ks ¼ ðI þ Cs ÞðI Ks Þ I. Then, problem of minimizing the norm (19) becomes the problem of minimizing kI ðI Ks ÞðI þ As Þk over integral operator Ks . We asP all circulant 2 sume vn;m ¼ ðAs un ; um Þ and nm ¼ þ1 n¼1 jvn;m j . Then, kI ðI Ks ÞðI þ As Þk ¼ kKs þ As Ks As k and by using (20) the kernel function of Ks Ks As þ As at point ðx; yÞ is given by 1 X dm;n km km ; vm;n vm;n um ðxÞuðyÞ þ 1 þ km 1 þ km m;n¼1
1 X dm;n km vm;n ¼ un ðyÞ; um ðxÞ 1 þ km m;n¼1 where dm;n denotes the Kronecker symbol. By the definition of Hilbert–Schmidt norm 2 þ1 X dm;n km vm;n ; kI ðI Cs Þ1 ðI þ As Þk2 ¼ 1þk m m;n¼1 it is clear, that the above expression is minimized if and only if the term is minimized for all integers m, however, by differentiating this quotient with respect to the real and imaginary parts of km , we see that nm þvm;m minimum will be obtained if we set km ¼ 1þ for all m, hence we have obvm;m tained kernel function of super-optimal circulant Cs as follows: ! 1 X nm þ vm;m Cs ðx yÞ ¼ um ðyÞ: um ðxÞ 1 þ vm;m m¼1 Pþ1 jvm;n j2 þj1þvm;m j2 Moreover, since 1 þ km ¼ n¼1 1þvm;m 6¼ 0, we see that I þ Cs is invertible. jkm j2 km vm;n km vm;n þnm j1þkm j2
4. A survey of stochastic numerical methods by using PCG algorithms In (6), we saw that one of the integrals is a stochastic integral, so we can write [3,15,16]:
104
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
wN ;k ¼
N X
bkj : xj :n00 :
k ¼ 1; . . . ; N ;
j¼1
where wN ;k is a simulation of approximate vector or random vector for w and here collocations of random variables fx1 ; . . . ; xN g correspond to the process rE rn R 0 0 from W ¼ Ge W_ ðn0 Þ oEðxon;n Þ dV so, independent standard Gaussian random variables fn001 ; . . . ; n00N g are normalized increments of the martingale g at jump N moments and ðbkj Þk;j¼1 are determined matrix, moreover the symbol: xj :n00 denotes the wick product [16] of the random variables xj . Then, we can generate wN ;k , by using [3,8,15,16,19]. Now, we consider these algorithms, which we call AlglR, Alg2R and Alg3R. In all those algorithms, the integral, are discretized by the rectangular quadrature rule, that the discretization is done at integer grid points 1; . . . ; N where N ¼ s > 1. Thus, integral equation (12) is replaced by a system of linear equations. The algorithm AlglR is the basic CG algorithm [9,11–13] without a preconditioner. In algorithm Alg2R and Alg3R the preconditioners are obtained respectively by optimal circulant integral in Section 3.2 and superoptimal circulant integral in Section 3.3. For stopping criterion is krk k=kr0 k < 107 , where rk is the residual vector of PCG, CG after k iteration. We show AlglS, Alg2S and Alg3S when algorithms the integral are discretized by the Simpson’s 1/3 quadrature rule. On below, we first test the performance of integral equations with normally distributed random variables as white noise, that is means and variance 2 EðDwN ;i Þ ¼ 0, EððDwN ;i Þ Þ ¼ ni for the each iteration. Tables 1 and 2 show, the numbers of iteration required for convergence for different choice of algorithms for Examples 1 and 2. We see from the tables that the numbers of iterations of the preconditioned systems are significantly less than that of the nonpreconditioned system. Moreover, we see that the iterations performance of optimal circulant operators be smaller than super-optimal circulant operators. All computations were computed using programs written in the symbolic language mathematical version 2.1.
Table 1 Iteration counts N
Alg1R
Alg3R
Alg2R
Alg1S
Alg3S
Alg2S
20 40 60 80 100
18 27 42 56 83
9 12 14 21 23
8 11 12 14 19
18 25 40 51 79
10 13 14 20 22
9 11 13 15 18
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
105
Table 2 Iteration counts N
Alg1R
Alg3R
Alg2R
Alg1S
Alg3S
Alg2S
20 40 60 80 100
15 23 37 43 78
11 12 19 24 32
10 10 13 17 27
17 22 31 37 75
10 11 17 22 31
9 9 11 14 23
Example 1. In the singular integral equation of the second kind Z 1 Z 1 rffiffiffiffiffiffiffiffiffiffiffirffiffiffiffiffiffiffiffiffiffi 1þx 1þt rðx; tÞ ¼ logðjx tjÞrðx; tÞdx dt þ w; 1x 1t 1 1
ð21Þ
¼ jx þ tj so, in the two-dimensional case, the random where we Rassume in w, oE on 1 R1 _ variable 1 1 W ðxÞ oE dx dt is Gaussian centered with variance 105 and on covariance 106 , so it can easily be simulated (by the classical Monte Carlo method). The numerical solutions is given at the point t ¼ 0:25 in Table 1. Remark. We can approximate gðx; tÞ ¼ ð1 þ xÞð1 tÞrðx; tÞ by a finite Chebyshev series gðx; tÞ ffi
p X p X k¼0
akl Tl ðxÞTk ðtÞ:
l¼0
Then rðx; tÞ ffi
p X p X k¼0
l¼0
Z
1
Z
1
akl 1
1
Tl ðxÞTk ðlÞ logðjx tjÞdx dt pffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 1 t 2
and using 1 p
Z
1 1
logðjx tjÞTk ðtÞ dt pffiffiffiffiffiffiffiffiffiffiffiffi ¼ 1 t2
log 2 Tk ðxÞ=k
k¼0 kP1
and the orthogonality of the Chebyshev polynomials with respect to the weight 1=2 function ð1 x2 Þ , the integral (21) can be determined in closed form. On the theoretical grounds we would expect good results. ¼ ðx þ tÞ2 which variance 103 Example 2. In Example 1, we consider oE on and convariance is 104 , so at the point t ¼ 0:25, we have given resultants in Table 2.
106
K. Maleknejad, F. Mirzaee / Appl. Math. Comput. 163 (2005) 97–106
References [1] S.K. Biswas, N.U. Ahmed, Stabilization of systems governed by the wave equation in the presence of distributed white noise, IEEE Trans. Autom. Contr. AC. 30 (10) (1985). [2] S.H. Crandall, W.Q. Zhu, Random vibration. A survey of recent development, J. Appl. Mech. 50 (1980) 953–962. [3] A.A. Dorogavtsev, An extended stochastic integral for smooth functionals of white noise, Ukrain. Math. J. 41 (1989) 1252–1258. [4] I. Goberg, I. Feldman, Convolution equation and projection methods for their solution, Trans. Math. Monogr. (1974) 41. [5] I. Goberg, N. Krupnik, in: One-Dimensional Linear Singular Integral Equations, Volumes I, II, Birkh€ auser Verlag, Basel, 1992. [6] I. Goberg et al., Fast preconditioned conjugate gradient algorithms, SIAM J. Matrix Anal. Appl. 31 (1994) 429–443. [7] M.A. Goberg, Numerical Solution of Integral Equations, Plenum Press, New York, 1990. [8] P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, SpringerVerlag, Berlin, 1995, Second corrected printing. [9] R. Kress, Linear Integral Equations, Springer-Verlag, Berlin, 1989. [10] L.S. Langston, Heat transfer from multidimensional objects using one-dimensional solutions for heat loss, Int. J. Heat Mass Transfer 25 (1982) 149–150. [11] K. Maleknejad, Numerical experimentation with time-series methods for convolution integral equation, J. Computat. Appl. Math. 31 (1990) 211–225. [12] K. Maleknejad, D. Rostami, Numerical solution of linear operator integral equation of the first kind by using methods of PCG, in: Lectures Note on Numerical Analysis, Polvdiv, Bulgarria, 4, 1995, pp. 111–122. [13] K. Maleknejad, D. Rostami, The PNCGC methods for solving some nonlinear boundary integral equation, Bull. of the Iranian Mathematics Society 26 (1) (2000) 7–29. [14] G.E. Myers, Analytical Methods in Conduction Heat Transfer, McGraw-Hill, New York, 1971. [15] T. Sekiguchi, V. Shiota, L2 -Theory of noncausal stochastic integrals, Math Rep. Toyama Univ. 8 (1985) 119–195. [16] A.V. Skorokhod, On a generalization of the stochastic integral, Theory Random Process. 7 (1979) 118–127. [17] P.J. Schneider, Conduction Heat Transfer, Addison-Wesley, Reading, MA, 1955. [18] E.E. Tyrtyshnikov, V.V. Strela, Which circulant preconditioner is better, Math. Comp. 65 (1996) 135–150. [19] R.J. Williams, K.L. Chung, Introduction to Stochastic Integration, Birkh€auser, Boston, 1990.