APPLICATION OF THE BOUNDARY ELEMENT METHOD USING DISCRETIZATION IN TIME
FOR NUMERICAL SOLUTION OF HYPERBOLIC EQUATION
Maria Lupa, Ewa Ładyga
Institute of Mathematics, Czestochowa University of Technology, Poland lupa@imi.pcz.pl
Abstract. The hyperbolic equation (1D problem) supplemented by adequate boundary and initial conditions is considered. This equation is solved using the combined variant of the boundary element method. The problem is also solved in analytical way. The comparison of the results obtained by means of these two methods confirms the effectiveness and accuracy of the BEM.
1. Formulation of the problem The following equation is considered
2 2 2
2 ( , ) ( , )
) , (
x t x U t
t x U t
t x U
∂
=∂
∂ +∂
∂
∂ (1)
where U(x, t) is an unknown function, x is the spatial co-ordinate and t is the time.
The equation (1) is supplemented by the boundary conditions
0 ) , 1 ( : 1 , 0
0 ) , 0 ( : 0 , 0
=
=
>
=
=
>
t U x
t
t U x
t (2)
and the initial ones
) 0 , : (
0 , 1 0
0 )
0 , ( : 0 , 1 0
0 0
∂ =
= ∂
<
<
>
=
=
<
<
=
t t
t x t U
x
U x U t
x
(3)
This type of boundary and initial conditions allows to solve the problem ana- lytically and in this way the results obtained by means of the boundary element method using discretization in time can be compared with the analytical solution.
2. Boundary element method
To solve the equation (1), the BEM using discretization in time is applied [1, 2]. So, the time grid
∞
<
<
<
<
<
<
<
=t0 t1 K tf−2 tf−1 tfK tF
0 (4)
with constant step ∆t=tf −tf−1 is introduced.
For the time interval
[
tf−2,tf]
the following approximations of time derivative can be taken into accountt t x U t x U t
t x
U f f
t
t f ∆
) , ( ) , ( )
,
( −1
=
= −
∂
∂ (5)
and
t t x U t x U t
t x
U f f
t
t f 2∆
) , ( ) , ( )
,
( −2
=
= −
∂
∂ (6)
or
t
t x U t
x U t
x U t
t x
U f f f
t
t f 2∆
) , ( ) , ( 4 ) , ( 3 )
,
( −1 −2
=
+
= −
∂
∂ (7)
The second time derivative is approximated in following way
2
2 1
2 2
) (
) , ( ) , ( 2 ) , ( ) , (
t
t x U t
x U t
x U t
t x
U f f t
∆
−
− +
= −
∂
∂ (8)
Let β =1∆t and Uf =U(x, f∆t). At the f-th time step t= f∆t (f ≥2) the equa- tion (1) can be approximately rewritten as
2 0
1 2
2 − + − =
∂
∂ f f− f−
f
U C U B U x A
U (9)
where for the first variant (equation (5))
2 2
2 β, 2β β, β
β + = + =
= B C
A (10)
for the second variant (equation (6))
β β β
β
β 2
, 1 2 2 ,
1 2 2
2+ = = −
= B C
A (11)
and for the third variant (equation (7))
β β β
β β
β 2
), 1 (
2 2 ,
3 2 2
2+ = + = +
= B C
A (12)
For equation (9) the weighted residual criterion is applied [1]
∫
=
− + −
∂
∂ − − ∗
1
0
2 1
2 2
0 ) , ( x dx U
U C U B U x A
U f f f
f ξ (13)
where ξ∈(0, 1) is the observation point, U * (ξ, x) is the fundamental solution and this function should fulfil the equation
) , ( ) , ) (
,
( *
2 2
x x
U x A
x
U ξ − ξ =−δ ξ
∂
∂ ∗
(14)
where δ(ξ, x) is the Dirac function. It can be check that the following function
( )
exp( )2 , 1
* x A
A x
U ξ = − −ξ (15)
fulfills the equation (14).
Additionally, the function q* (ξ, x) resulting from fundamental solution is defined
x x x U
q ∂
−∂
= ( , )
) , (
*
* ξ ξ (16)
and it can be calculated analytically
) (
2 exp ) ( ) sgn ,
*(
A x x
x
q ξ = −ξ − −ξ (17)
where sgn(⋅) is the sign function.
Integrating twice by parts the first component of equation (13) and taking into account the property (14) of fundamental solution one obtains
∫
− =+
+
−
∂
− ∂
∂
∂
−
−
=
=
=
= 1
0
* 2 1
1
0 1 *
0
*
0 ) , ( ) (
) , ) (
, ) (
, (
x d x U U C U B
t x U
x U U
x x U U
f f
f x
x f
x
x f
ξ ξ ξ ξ
(18)
or
) , ( ) , 0 ( ) 0 , ( ) , 1 ( ) 1 , (
) , 0 ( ) 0 , ( ) , 1 ( ) 1 , ( ) , (
*
*
*
*
f f
f
f f
f
t P t U q
t U q
t q U
t q U
t U
ξ ξ
ξ
ξ ξ
ξ
+
−
=
=
−
+ (19)
where q(x,tf)=−∂U(x,tf)/∂x and
x d x U U C U
B t
P f =
∫
1 f− − f−0
* 2
1 ) ( , )
( ) ,
(ξ ξ (20)
For ξ → 0 + one obtains
) , 0 ( ) , 0 ( ) 0 , 0 ( ) , 1 ( ) 1 , 0 (
) , 0 ( ) 0 , 0 ( ) , 1 ( ) 1 , 0 ( ) , 0 (
*
*
*
*
f f
f
f f
f
t P t U q
t U q
t q U
t q U
t U
+
−
=
=
− +
+
+ (21)
and for ξ → 1- one has
) , 1 ( ) , 0 ( ) 0 , 1 ( ) , 1 ( ) 1 , 1 (
) , 0 ( ) 0 , 1 ( ) , 1 ( ) 1 , 1 ( ) , 1 (
*
*
*
*
f f
f
f f
f
t P t U q
t U q
t q U
t q U t U
+
−
=
=
− +
−
− (22)
The system of equations (21), (22) can be written in the matrix form
+
−
−
−
−
=
−
−
−
−
+ +
) , 1 (
) , 0 ( ) , 1 (
) , 0 ( 1 ) 1 , 1 ( ) 0 , 1 (
) 1 , 0 ( 1 ) 0 , 0 (
) , 1 (
) , 0 ( ) 1 , 1 ( ) 0 , 1 (
) 1 , 0 ( ) 0 , 0 (
*
*
*
*
*
*
*
*
f
t P
t P t
U t U q
q
q q
t q
t q U
U
U U
f f
f f
f
(23)
or
+
=
) , 1 (
) , 0 ( ) , 1 (
) , 0 ( )
, 1 (
) , 0 (
22 21
12 11 22
21 12 11
f f f
f f
f
t P
t P t
U t U H H
H H t
q t q G G
G G
(24) where G11=−G22=−1/(2 A), G12=−G21=−1/(2 A)exp(− A),
2 /
22 1
11 =H =−
H , H12=H21=1/2exp(− A)
This system of equations allows to find the boundary values U(0, t f ), U(1, t f ).
Next, the values of U at the internal nodes ξ∈(0, 1) are calculated using the for- mula (c.f. equation (19))
) , ( ) , 0 ( ) 0 , ( ) , 1 ( ) 1 , (
) , 0 ( ) 0 , ( ) , 1 ( ) 1 , ( ) , (
*
*
*
*
f f
f
f f
f
t P t q U
t q U
t U q
t U q
t U
ξ ξ
ξ
ξ ξ
ξ
+ +
−
+
−
= (25)
this means (c.f. formulas (15) and (17))
) , ( ) , 1 ( ] ) 1 ( 2 exp[
) 1 , 0 ( ) 2 exp(
1
) , 1 ( ] ) 1 ( 2exp[
) 1 , 0 ( ) 2exp(
) 1 , (
f f
f
f f
f
t P t q A A
t q A A
t U A t
U A t
U
ξ ξ
ξ
ξ ξ
ξ
+
−
−
−
− +
+
−
− +
−
=
(26)
3. Analytical solution
To solve the problem (1), (2), (3) analytically, the Fourier method [3, 4] is ap- plied. So, one assumes that
∑
∞=
=
0
) , ( )
, (
n
n x t U t
x
U (27)
where
) ( ) ( ) ,
(x t X x T t
Un = n n (28)
and then
) ( ) ) (
, ), (
( ) ) (
, (
) ( ) ) (
, ), (
( ) ) (
, (
'' 2
' 2
'' 2
' 2
t T x X x
t x t U
T x x X
t x U
t T x X t
t x t U
T x t X
t x U
n n n
n n n
n n n
n n n
∂ =
= ∂
∂
∂
∂ =
= ∂
∂
∂
(29)
Putting (29) into (1) one obtains
) ( ) ( ) ( ) ( ) ( )
(x T' t X x T '' t X'' x T t
Xn n + n n = n n (30)
this means
) (
) ( ) (
) ( ) (
)
( '' ''
'
x X
x X t T
t T t T
t T
n n n
n n
n + = (31)
It is assumed that
'' 2 ''
'
) (
) ( ) (
) ( ) (
) (
n n
n n
n n
n
x X
x X t T
t T t T
t
T + = =−λ (32)
where
λ
n ≠0 are the constants.If the functions Xn (x), Tn (t) will be the solutions of equations
'' 2 '
) (
) ( ) (
) (
n n
n n
n
t T
t T t T
t
T + =−λ , 2
''
) (
) (
n n
n
x X
x
X =−λ (33)
then these functions will fulfill the equation (31). The equations (33) can be writ- ten in the form
0 ) ( )
( )
( ' 2
'' t +T t + T t =
Tn n λn n (34)
and
0 ) ( )
( 2
'' x + X x =
Xn λn n (35)
The solution of equation (35) is following
x B
x A
x
Xn( )= ncosλn + nsinλn (36)
Taking into account the boundary conditions (2) one has 0
) 0
( = n=
n A
X , Xn(1)=Bnsinλn=0→λn=nπ (37) this means
x n B x
Xn( )= nsin π (38)
The equation (34) can be written in the form
0 ) ( )
( )
( ' 2 2
'' t +T t +n T t =
Tn n π n (39)
Because
0 4
1− n2π2< (40)
so the solution of equation (39) is the following
) sin cos
( 2) ( exp )
( t C t D t
t
T n = − n αn + n αn (41)
where
2 1 4 2 2−
= π
αn n (42)
Finally, the function (27) has the form
∑
∞=
+
−
=
1
) sin cos
( sin 2) ( exp ) , (
n
n n n
n t F t
E x t n
t x
U π α α (43)
where: En=BnCn, Fn=Bn Dn are the constants.
Now, the initial conditions (3) should be taken into account, this means
0 1
sin )
0 ,
(x E n x U
U
n
n =
=
∑
∞=
π (44)
and
0 sin
2 ) ( 1 )
, (
0 1
= +
−
∂ =
∂
∑
∞= n=
n n n t
x n F t E
t x
U α π (45)
For the arguments x∈ [-1, 1] function U(x, 0) can be extended on the uneven function
=
∈
=
−
∈
−
−
=
=
1 ,
0
) 1 , 0 ( ,
0 ,
0
) 0 , 1 ( ,
1 ,
0
) 0 , (
0 0
x x U
x x U
x
x
U (46)
Taking into account the expansion of this function into a Fourier series one obtains
[ ]
∫
= − −=
1
0
0
0 2 1 ( 1)
sin
2 n
n n
dx U x n U
E π π (47)
From condition (45) results that the zero function is expanded into the Fourier series. In such case
2 0
1 + =
− En αn Fn (48)
this means
[
n]
n
n n
F = U0 1−(−1) α
π (49)
Finally, one obtains
∑
∞=
− +
− −
=
1
0 1 sin )
cos 2 ( ) sin
1 ( ) 1 ( 2 exp )
, (
n
n n n n
t t
x n n
t t U
x
U α
α α
π π (50)
4. Results of computations
Application of the boundary element method using discretization in time re- quires a proper assumption of time step ∆t. Additionally, the integral (20) should be determined with sufficient accuracy. To calculate this integral the domain [0, 1]
is divided into equal M sub-domains and six-point Gauss integral formula is used.
All computations have been done under the assumption that U0 =1.
In Figure 1 the curves at the point x=0.5 for different approximations of time derivative (c.f. formulas (5), (6), (7)) are shown. It is visible that the results are practically the same. The calculations have been done for ∆t=0.02, M = 100 and for these values very good agreement with analytical solution has been observed.
Figure 2 illustrates the distribution of function U for times 1, 2, 3, 4 and 5 s.
Fig. 1. Course of function U at x= 0.5
It should be pointed out that the influence of time step ∆t on the results of com- putations is big. In Figure 3 the curves at the point x=0.25 and for
05 . 0 , 02 . 0 , 01 .
=0
∆t (M = 100) are shown.
Summing up, the BEM using discretization in time constitutes the effective numerical method of hyperbolic equation solution but it requires a proper choice of time step ∆t and number of internal cells M.
Fig. 2. Distribution of function U
Fig. 3. Curves for different ∆t at the point x= 0.25
References
[1] Majchrzak E., Boundary element method in heat transfer, Publ. of the Czest. Univ. of Technol- ogy, Czestochowa 2001 (in Polish).
[2] Majchrzak E., Ładyga E., Mendakiewicz J., Piasecka Belkhayat A., Different variants of the boundary element method for parabolic equations, Scientific Research of the Institute of Mathe- matics and Computer Science 2004, 1(3), 127-132.
[3] Kącki E., Partial differential equations, WNT, Warsaw 1995 (in Polish).
[4] Lupa M., Analytical solution of Cattaneo equation, Scientific Research of the Institute of Ma- thematics and Computer Science 2007, 1(6), 127-132.
