VIEWS: 2 PAGES: 6 CATEGORY: Research POSTED ON: 9/17/2011
(IJACSA) International Journal of Advanced Computer Science and Applications, Vol. 2, No. 8, 2011 Identiﬁcation Problem of Source Term of A Reaction Diffusion Equation Bo Zhang Linyi University at Feixian Feixian, Shandong, P.R.China y y ( x, t ) y ( x, t ) Abstract—This paper will give the numerical difference scheme u ( x) v( x) c ( x ) y ( x, t ) with Dirichlet boundary condition， and prove stability and t x x convergence of the difference scheme, final numerical experiment results also confirm effectiveness of the algorithm. b( x, t ), x , t 0, (2) Keywords- Fractional derivative; Numerical difference scheme; The gradient regularization method. Dy( x, t ) g ( x, t ), x , t 0, (3) I. INTRODUCTION Source term inversion in groundwater pollution is a class Ey( x, 0) f ( x), x , of inverse problems[7,8], and is also important ﬁeld of (4) inverse problem research. Scholars have done a large amount of work and obtained many results. Reference [2] presented a where 0 1,1 2, D represents matrix operator new gradient regularization algorithm to solve an inverse with boundary condition, E represents matrix operator with problem of determining source terms in one-dimension initial condition, y ( x, t ) and b( x) represent undetermined solute transport with ﬁnal observations, and reference [3] source term function and undetermined vector function proposed a implicit method to solve a class of space-time respectively. fractional order diffusion equations with variable coefficient. Inverse problem of this problem is to determine unknown However, when fractional derivative replaces second vector function b( x) by Eqs. (2)-(4) and the additional derivative in diffusion equations, there is anomalous condition below diffusion phenomenon. In this paper, we give the numerical difference scheme in source term identiﬁcation with Dirichlet y( x, t ) |x T ( x). boundary condition, and we prove the stability and convergence of the difference scheme, also verify the (5) practicality and effectiveness of the algorithm through II. NUMERICAL DIFFERENCE SCHEME numerical experiment. Considering the following space-time reaction diffusion In this paper, we use difference scheme to solve the equations forward problem, and when solving the inverse problem, we use the gradient regularization method based on Tikhonov y ( x, t ) y ( x, t ) y ( x, t ) regularization strategy. Here, the additional information for u ( x) v( x) source term identiﬁcation is set as the ﬁnal observations, and t x x suppose that the source term function are only concerned c( x) y( x, t ) b( x, t ), with the space variable and has nothing to do with the time variable. (6) In fact, the solute transport model can be described by the y( x, 0) f ( x), following equation[1] (7) y 2 y ( x, t ) y ( x, t ) y(0, t ) g1 (t ), y( L, t ) g 2 (t ). u ( x) v( x) c ( x ) y ( x, t ) (8) t x 2 x b( x, t ), x [0, L], t [0, T ], Where u( x) 0, v( x) 0, c( x) 0 are continuous (1) functions on [0, L], g1 (t ) 0, g 2 (t ) 0 are continuous by introducing fractional derivative and adding initial functions on [0, T ], b( x, t ) is continuous on [0, L] boundary conditions[5], Eqs. (1) will be modified as the y ( x, t ) y ( x, t ) following problem [0, T ], (0,1), (1, 2), , are t x 76 | P a g e www.ijacsa.thesai.org (IJACSA) International Journal of Advanced Computer Science and Applications, Vol. 2, No. 8, 2011 Caputo fractional derivative and Riemann-Liouville 1 k y ( xi , tk 1 j ) y ( xi , tk j ) fractional derivative respectively[5] (2 ) [( j 1)1 j1 ] j 0 y ( x, t ) 1 y ( x, ) (1 ) 0 t (t ) d , u ( xi ) i 1 t h a y( x ( j 1)h, t j i k 1 ) - c( xi ) y ( xi , tk 1 ) j 0 (9) y ( xi , tk 1 ) - y ( xi -1 , tk 1 ) y ( x, t ) 1 2 y ( , t ) b( xi , tk 1 ) O( h). t -v( xi ) d . h x (2 ) x 2 0 (x ) 1 (10) % Let ui u ( xi ), vi v( xi ), ci c( xi ), ci (2 ) 0, Suppose that Eqs. (6)-(8) has a unique and smooth bik b( xi , tk ), j ( j 1) 1 enough solution. T / n is time step, x h L / m is j1 , j 0,1, 2,L , n; space step, tk k (k 0,1, 2,L , n), ri xi ih(i 0,1, 2,L , m). For the time fractional ui vi derivative, we usually adopt following finite difference (2 ) 0, pi (2 ) 0, g1k g1 (tk ), g 2 k approximation h h g 2 (tk ), yi represents numerical solution of Eqs. (6)-(8), k y ( xi , tk 1 ) then we have t k i 1 1 k y ( xi , t j 1 ) y ( xi , t j ) ( j 1) d j ( yik 1 j yik j ) pi ( yik 1 yik1 ) ri a j yikj11 1 (1 ) j (tk 1 ) j 0 j 0 j 0 1 k y ( xi , tk 1 j ) y ( xi , tk j ) % % ci yik 1 ci bik 1 . (2 ) [( j 1)1 (14) j 0 j1 ] O( ), Since local truncation error is O( h), thus the (11) difference scheme is consistent[10]. Eqs. (14) will be repalced by y ( xi , tk 1 ) y ( xi , tk 1 ) y ( xi 1 , tk 1 ) O(h). When k 0, x h (12) i 1 ri yi11 (1 pi ri a1 ci ) yi1 ( pi ri a2 ) yi11 ri a j yi1 j 1 % y ( x, t ) j 3 For , by using Grünwald’s improved formula x % y ci bi1 , 0 i [4], we have that (15) y ( xi , tk 1 ) 1 i 1 when k 0, x h a y( x ( j 1)h, t j i k 1 ) O( h), j 0 i 1 ri yik1 (1 pi ri a1 ci ) yik 1 ( pi ri a2 ) yik1 ri a j yikj11 1 % 1 (13) j 3 k where yik j ( yik 1i yik j ) (2 21 ) yik k yi0 ci bik 1 % a0 1, a1 , a j (1) j , j j 0 k 1 ( 1)L ( j 1) / ( j !), j j 1, 2,3,L . y k j i [2( j 1)1 ( j 2)1 j1 ] j 1 Substituting Eqs. (11)-(13) into Eqs. (6)-(8), we obtain k 1 d y yik j d j 1 k yi0 ci bik 1 , k 1 i % j 1 (16) 77 | P a g e www.ijacsa.thesai.org (IJACSA) International Journal of Advanced Computer Science and Applications, Vol. 2, No. 8, 2011 where d j j 1 j , i 1, 2,L , m 1; k 1, 2,L , n 1. When k 0, III. STABILITY AND CONVERGENCE OF THE DIFFERENCE % ri i11 (1 pi ra1 ci ) i1 ( pi ra2 ) i11 i i SCHEME i 1 Lemma 2.1 For arbitrary real number a, b, c, d , we have ri a j i1 j 1 i0 , j 3 | a | | b | | c | | d || a b c d | . Proof. From | b || a b c d a c d | when k 0, | a b c d | | a | | c | | d |, % ri ik1 (1 pi ra1 ci ) ik 1 ( pi ra2 ) ik1 1 i i 1 we obtain | a | | b | | c | | d || a b c d | . i 1 k 1 ri a j ikj11 d1 ik d j 1 ik j k i0 . Lemma 2.2 (1) a j 0( j 2). (2) For any positive j 3 j 1 N a 0. E k (1k , 2 ,L , m1 )' , k k integer N ，we have j Let we prove j 0 || E || || E || with mathematical induction in the k 0 Proof. (1) Note that a j (1) j and 1 2, we j following. know that a j 0( j 2). When k 1, let | l | max | i | . Note that rl 0, 1 1 1i m 1 pl 0, a0 1, a1 0, we have from Lemma 2.2 that (2) Since (1 x) x , x [1,1], let j j x 1, l 1 j 0 || E1 || | l1 || l1 | pl (| l1 | | l11 |) rl a j | l1 | then a j 0. From Lemma 2.2(1), we have that % j 0 j 0 rl | | (1 pl rl a1 ci ) | l1 | 1 l 1 N a j a j 0. l 1 ( pl rl a2 ) | l11 | rl a j | l1 j 1 |. j 0 j N 1 j 3 n d l 1 1 n ; (2)d j 0, j 1 j . Lemma2.3 (1) j 1 j Note that pl rl a2 0, rl a j % 0,1 pl rl a1 ci 0, j 3 and from Lemma 2.1, we further obtain that Proof. (1)From d j j 1 j , j ( j 1) 1 j1 , n % | l1 || rl l11 (1 pl rl a1 ci ) l1 ( pl rl a2 ) l11 we easily know that d j 1n. l 1 j 1 rl a j l1 j 1 | | l0 ||| E 0 || . 1 j 3 (2) Let h( x) ( x 1) x1 ( x 1), then h '( x) Thus || E || | l ||| E || . 1 1 0 (1 )[( x 1) x ] 0, so h( x) is decreasing function, d j j 1 j h( j 1) h( j ) 0. Therefore, Assume that we always have || E || || E || when k 0 we have that d j 0, j 1 j . k s, let | ls 1 | max | is 1 |, then when k s 1, we 1i m 1 A. Stability of the difference scheme have Theorem 2.1 Implicit difference schemes defined by Eqs. % (15)-(16) are unconditionally steady for initial value[10]. | ls 1 | rl | ls11 | (1 pl rl a1 ci ) | ls 1 | l 1 ( pl rl a2 ) | ls11 | rl a j | lsj11 | k Proof. Suppose that %, yi represent solutions of Eqs. yi k j 3 k (15)-(16) for initial value f1 ( x), f 2 ( x) respectively, and bi % | rl ls11 (1 pl rl a1 ci ) ls 1 k is accurate vale, then error % yi satisfies k k yi i 78 | P a g e www.ijacsa.thesai.org (IJACSA) International Journal of Advanced Computer Science and Applications, Vol. 2, No. 8, 2011 l 1 l 1 ( pl rl a2 ) ls11 rl a j ls11 | j ( pl rl a2 ) | el11 | rl a j | el1 j 1 | j 3 j 3 s 1 % | re (1 pl rl a1 ci )el1 ( pl rl a2 )el11 1 | d1 ls d j 1 ls j s l0 | l l 1 j 1 l 1 s 1 rl a j el1 j 1 | | Rl1 | d1 || E s || d j 1 || E s j || s || E 0 || j 3 j 1 ( 1 h) 01 ( 1 h). s 1 (d1 d j 1 s ) || E 0 || Assume that || e || s 1 ( s 1 1 h) when k s, s 1 j 1 let | el s 1 | max | ei |, then when k s 1, we have that 1 i m 1 (1 s s ) || E 0 || || E 0 || . Consequently, the desired result follows. || es 1 || | els 1 | B. Convergence of the difference scheme s 1 Suppose that y ( xi , tk ) is exact solution of the d1 || e s || d j 1 || e s j || ( 1 h) j 1 differential equation at ( xi , tk ). Let ei y ( xi , tk ) yi , k k ek (e1k , e2 ,L , k em 1 )' , k then yik y( xi , tk ) eik , (1 d1 s1 d2 s2 L d s 0 1 ) ( 1 h), 1 1 as in Lemma 2.3, we have found that j s ( j s), so 1 1 substituting it into Eqs. (15)-(16), and note e 0, we have 0 that we futher obtain s 1 When k 0, || e s 1 || s ( di 1 s ) ( h) 1 1 i 0 % rei11 (1 pi ra1 ci )ei1 ( pi ra2 )ei11 i i i s1 ( 1 h). i 1 ri a j ei1 j 1 Ri1 , The desired result follows. k j 3 1 1 1 Since lim lim lim when k 0, k k k k k k k [(1 1 / k ) 1 1] % reik1 (1 pi ra1 ci )eik 1 ( pi ra2 )eik1 1 i 1 i i 1 , thus there is a constant 0 such that i 1 k 1 1 ri a j eikj11 d1eik d j 1eik j Rik 1 , || E k || k ( 1 h) (k ) ( h). j 3 j 1 where | Ri | ( h), is a positive constant and it k 1 Consequently, we can obtain the following result when k T . has nothing to do with h, . Theorem 2.3 Suppose that y ( xi , tk ) denotes exact Theorem 2.2 There is a constant 0 such that k solution at ( xi , tk ), yi is numerical solution of implicit % difference scheme, then there exists a constant T 0 || ek || k1 ( 1 h), k 1, 2,L , n, 1 such that where has nothing to do with h, , and || e || k % | y( xi , tk ) yik | ( h), 1 i m,1 k n. max | eik | . 1 i m 1 IV. SOURCE TERM INVERSION Proof. When k 1, let | el | max | ei |, then || e || 11 1 1i m 1 The inverse problem, which is composed of Eqs. (2)-(5), 1 is to solve nonlinear operator equation | e |, we have from Lemma 2.1 that A[b( x)] ( x). l % | e | rl | e | (1 pl rl a1 ci ) | el1 | 1 1 l l 1 79 | P a g e www.ijacsa.thesai.org (IJACSA) International Journal of Advanced Computer Science and Applications, Vol. 2, No. 8, 2011 Suppose that y(b( x); x, t ) denotes solution of Eqs. (2)- K0 ( AT A I )1 AT (Q P), * (4) for b( x), b0 ( x) denotes a function near b ( x), where (20) n where b0 ( x) k ( x) K i 0 i T 0 ( x), and b* ( x) denotes exact i 1 y (b0 ( x); x1 , T ) yT (x1 ) solution of Eqs. (2)-(4), b( x) K , 1 ( x), 2 ( x), L , are a y (b ( x); x , T ) y (x ) P , Q , 0 2 T 2 group of primary functions on K , then a tiny disturbing M M quantity of b0 ( x) is y (b0 ( x); xM , T ) yT (xM ) n b0 ( x) ki0 i ( x) K 0T ( x)， A (aij ) M N , aij y (b0 ( x); xi , T ). i 1 k j (17) Substituting Eqs. (20) into Eqs. (17), we can obtain where ( x) ( 1 ( x), 2 ( x),L , n ( x))T , b0 ( x) and a new approximation of the exact solution, K T (k1 , k2 ,L , kn ) R n . b1 ( x) b0 ( x) b0 ( x). Assume that y (b1 ( x); x, t ) denotes the solution of initial Repeating the above, until satisfies the precision requirement. boundary value problem for b1 ( x), where b1 ( x) b0 ( x) b0 ( x), using Taylor formula, then we have[6] V. NUMERICAL E XPERIMENTS In order to verify the effectiveness of the algorithm in the y(b0 ( x) b0 ( x); x, t ) source term identification, we do the following numerical experiment[7]. For simplicity, we set part variables as y(b0 ( x); x, t ) T y(b0 ( x); x, t ) K 0 o(|| b0 ( x) ||), K 0 follows by using the Tikhonov regularization method, solving b( x) c( x) 0.05, u( x) 292, v( x) 365, L 4000, T 11, is converted into K 0 , and K 0 can be determined by local m 20, n 11, [0.01, 0.01, 0.01]. minimum of the following function[9] Where is the increment of K when computing matrix A from Eqs. (20). Moreover, we always take polynomial G[ K0 ]=||y(b0 ( x) b0 ( x); x, t ) y( x, T ) ||2 ( [0,T ]) L 2 ' function space as primary function space in the following computation, and setting initial boundary condition as S[ K0 ] =||y(b0 ( x); x, t ) y( x, T ) follows T y(b0 ( x); x, t ) K0 ||2 ( [0,T ]) S[ K0 ], K 0 L 2 ' yT ( x) 0.057 x 45.6, 0 x 4000, (18) g1 (t ) 0.724t 2 45.6, 0 t 11, where x , denotes regularization coefficient, ' S[ K 0 ] denotes steady functional of K 0 . g 2 (t ) 2.2t 2 273.4, 0 t 11. Let b( x) 1 x in the Eqs. (6)-(8), and substituting Assume that there are discrete points xm (m 1, 2,L , initial boundary condition, we can obtain y( X , T ) by M ) on ' , S[ K0 ]= K0 K0 , then T sloving forward problem. And as the additional data, we can do inversion calculation by using the above algorithm. Let G[ K0 ]= K0 AT A K0 2 K0 AT ( P Q) T T initial iteration vector K0 [1,1,1], and iterative termination condition b( x) 1e 4, then we obtain inversion results ( P Q)T ( P Q) K0 K0 . T under different regularization coefficient(see TABLEⅠ let ), (19) theta 1e 3, we get inversion results under different initial It is easy to prove that solving the local minimum values value(see TABLE Ⅱ). of Eqs. (19) is equivalent to slove ( A A I ) K0 A T T (Q P), so we have TABLE Ⅰ. THE INVERSION RESULTS UNDER DIFFERENT REGULARIZATION COEFFICIENT Iteration times Results 80 | P a g e www.ijacsa.thesai.org (IJACSA) International Journal of Advanced Computer Science and Applications, Vol. 2, No. 8, 2011 1e-1 4 1.0000 -1.0000 0 To better simulate the errors generated by actual data, 1e-2 3 1.0000 -1.0000 0 and verify the effectiveness of the algorithm, we choose the 1e-3 3 1.0000 -1.0000 0 disturbance error V V (1 ), where [1,1], and 1e-4 3 1.0000 -1.0000 0 0 is error level. TABLE Ⅱ. THE INVERSION RESULTS UNDER DIFFERENT INITIAL VALUE According to the above algorithm, we do 8 times Iteration numerical experiments, and obtain the inversion results under K0 Results times different error level (see TABLE Ⅲ ). Besides, the -100 -100 -100 4 1.0000 -1.0000 0 -10 -10 -10 3 1.0000 -1.0000 0 comparison of inversion results and exact solution can see 1 1 1 3 1.0000 -1.0000 0 Figure 1 when =0.01. 10 10 10 3 1.0000 -1.0000 0 100 100 100 4 1.0000 -1.0000 0 Figure 1. The comparison of inversion results and exact solutions TABLE Ⅲ. THE INVERSION RESULTS UNDER D IFFERENT ERROR LEVEL Times =0.01 =0.05 =0.1 1 0.9941 -0.9997 -0.0000 1.4988 -1.0243 0.0000 0.3838 -0.9700 -0.0000 2 1.1639 -1.0080 0.0000 0.3221 -0.9670 -0.0000 2.3115 -1.0639 0.0000 3 1.1098 -1.0053 0.0000 0.8026 -0.9903 -0.0000 -1.0527 -0.9000 -0.0000 4 0.9818 -0.9991 -0.0000 1.9119 -1.0444 0.0000 -0.5123 -0.9264 -0.0000 5 0.7984 -0.9902 -0.0000 1.8730 -1.0425 0.0000 -0.2448 -0.9394 -0.0000 6 1.1346 -1.0066 0.0000 0.8121 -0.9909 -0.0000 -0.2617 -0.9386 -0.0000 7 0.9768 -0.9989 -0.0000 1.8243 -1.0401 0.0000 1.4347 -1.0212 0.0000 8 1.0483 -1.0024 0.0000 0.0742 -0.9549 -0.0000 0.0459 -0.9535 -0.0000 mean value 0.9941 -1.0013 0.0000 1.1399 -1.0067 0.0000 0.2631 -0.9641 -0.0000 Through the above numerical experiment, we find that the Vol.172,pp. 65-77, 2004. inversion results and exact solution are almost the same, and [5] Podlubny I, Fractional differential equations, Academic Press, San this shows that the above algorithm is feasible and very Diego,1999. effective. [6] Jinqing Liu, Gongsheng Li and Yu Ma, Gradient—regulation method for determining a pollution source term in groundwater, Journal of REFERENCES Shandong University of Technology(Natural Science Edition), Vol. [1] Andreas Kirsch, An introduction to the mathematical theory of inverse 21,No.2,pp.17-21,2007. problems, Springer, Karlsruhe, 1996. [7] Gongsheng Li, Yongji Tan and Xiaoqin Wang, Inverse problem method [2] Gongsheng Li, Jinqing Liu, et al., A new gradient regularization on determining magnitude of groundwater pollution sources, algorithm for source term inversion in 1D solute transportation with ﬁnal Mathematica Applicata, Vol.18, No.1, pp.92-98, 2005. observations, Appl. Math. Comput, Vol.196,pp. 646-660, 2008. [8] Nazheng Sun, Inverse problem in groundwater modeling. Kluwer, [3] Yang Zhang, A ﬁnite method for fractional partial differential Dordrecht, 1994. equation,Applied Mathematics and Computation, Vol.215, pp. 524-529, [9] Tingyan Xiao, Shengen Yu and Yanfei Wang, Numerial solution of 2009. inverse problem, Beijing: Science Press, 2003. [4] Meerschaert M M and Tadjeran C, Finite difference approximations for [10] Wensheng Zhang, The finite difference method of partial differential fractional advection-dispersion ﬂow equations, J. Comput. Appl. Math., equation in science calculation, Beijing: Higher Education Press, 2006. 81 | P a g e www.ijacsa.thesai.org