# Paper 15-Identification Problem of Source Term of A Reaction Diffusion Equation by editorijacsa

VIEWS: 2 PAGES: 6

• pg 1
```									                                                                                      (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  yik1 )  ri  a j yikj11
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 yi11  (1  pi  ri a1  ci ) yi1  ( pi  ri a2 ) yi11  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 yik1  (1  pi  ri a1  ci ) yik 1  ( pi  ri a2 ) yik1  ri  a j yikj11
1
%
1
(13)                                                                                                                                                                                                                                                  j 3

k
where
 yik    j ( yik 1i  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 i11  (1  pi  ra1  ci ) i1  ( pi  ra2 ) i11
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 ik1  (1  pi  ra1  ci ) ik 1  ( pi  ra2 ) ik1
1               i                          i        1

we obtain  | a |  | b |  | c |  | d || a  b  c  d | .                                                    i 1                                k 1
ri  a j  ikj11  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 ,  m1 )' ,
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

1i  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 |  |  l11 |)  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 ) |  l11 | 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  l11  (1  pl  rl a1  ci ) l1  ( pl  rl a2 ) l11
we easily know that                    d         j
 1n.                                                     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
1i  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 |  ls11 | (1  pl  rl a1  ci ) |  ls 1 |
l 1
( pl  rl a2 ) |  ls11 | rl  a j |  lsj11 |
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  ls11  (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 ) ls11  rl  a j  ls11 |
j
( pl  rl a2 ) | el11 | rl  a j | el1 j 1 |
j 3                                                                                                  j 3

s 1                                                                                                         %
| re  (1  pl  rl a1  ci )el1  ( pl  rl a2 )el11
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)   01 ( 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 s1  d2 s2  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
%
rei11  (1  pi  ra1  ci )ei1  ( pi  ra2 )ei11
i                  i                      i
  s1 ( 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]
%
reik1  (1  pi  ra1  ci )eik 1  ( pi  ra2 )eik1                                                  1
i   1               i                         i       1
, thus there is a constant   0 such that
i 1                                   k 1
1
ri  a j eikj11  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 ||   k1 ( 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

1i  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

```
To top