AN INTERVAL DIFFERENCE METHOD FOR SOLVING THE WAVE EQUATION

(1)

__________________________________________

* Poznan University of Technology.

Barbara SZYSZKA*

AN INTERVAL DIFFERENCE METHOD FOR SOLVING THE WAVE EQUATION

In the paper a difference interval method for solving the wave equation together the initial- boundary value problems is presented. Using an interval method together floating-point interval arithmetic guarantee, that obtained interval solutions contain all numerical errors. Additionally, each exact solution is included into interval solution. In numerical experiments it is guarantee contain all numerical errors in obtained interval solutions. Taken into consideration is the central discretization method with regard to space and time. An initial condition is approximated by the third-order Taylor polynomial with local truncation error of order O(h⁴). In the paper new formula, which described discretization of the initial condition, is proposed. Therefore more exact solutions are obtained then in the previous considerations.

1. INTRODUCTION

The studies on Partial Differential Equations (PDE) in context to interval methods on floating point interval arithmetic are continued by author. The difference interval methods together initial-boundary conditions for solving the wave equation were presented at GAMM [5], MACMESE [7], ICNAAM [8] and PPAM [9] Conferences. The main point of construction interval methods is to contain error of method into obtained solutions. Additionally, floating-point interval arithmetic is used for computer calculations. Thus solutions, in form of intervals, contain all possible numerical errors.

The methods presented in [5] and [8] concern on the initial condition error of order O(h) and the estimations of errors are depended on a parameter of the wave equation (v), but such estimations were correct only for small value v1. Another way estimation error of method was presented in [9] for central- and in [7] for central-backward difference interval method. In both papers the estimation of error does not depend on parameter v. Additionally, an initial condition was considered with local truncation error of order O(h⁴).

In this paper the central difference interval method for solving wave equation together initial-boundary conditions is considered. The estimation of method error is similar to method presented in [9], but an initial condition (with local truncation error of order O(h⁴)) is approximated by another formula. In numerical experiments solutions obtained while using conventional and interval methods are compared.

(2)

2. THE WAVE EQUATION

Taken into consideration is the string, which is an example of the one dimensional wave equation [2], [4] and [10].

Assuming that the string is well flexible and homogeneous (mass of string per unit length  is a constant). The tension T of the string is constant and larger than the force of gravity (no other external forces act on the string). Supposing damping effects are neglected and the amplitude is not too large and each point of the string can move only in the vertical direction.

If the string is stretched between two points (see Fig. 1) x0 and x L and function uu

 

x,t denotes the amplitude of the string's displacement in the time t, then the wave equation is in the form

   

, , 0 ,

2 2 2 2

2 







t t x u x

t x

v u (1)

for 0x L and t 0, where v ² T/ = const.

Since the ends of the string are secured to the x -axis, function u satisfies (for 0



t ) the following Dirichlet boundary conditions . 0 ) , ( , 0 ) , 0

( t  u Lt 

u (2)

The initial position and velocity of the string are adequately given by the Cauchy initial conditions

 

, ( ), ), ( ) 0 , (

0

t x t x u

x x

u

t





 





(3) for 0x L, where  and  are given functions.

3. THE CENTRAL DIFFERENCE METHOD

To set up the difference method the space x and time t are divided into n and m equal parts adequately, where  x0

n

h L , k t 0. Thus, the space-time grid with the mesh points



x ,i tj



, i0 ,1, ,n, j0 ,1, ,m, is received, where

,

, t j k

h i

x_i   _j   (4)

The

 

2

2 ,

x t x u _i _j



 and

 

2

2 ,

t t x u _i _j



 into the Taylor series are expanded to obtain the central difference formulas. Then putting them into the wave equation (1) and substituting approximation u_ij for function u



xi,tj



, the central difference method is obtained as follows

(3)



i j i j



ij ij M j

i u u u u e

u     

 _, ² _₁_, _₁_, _, _₁ _, _₁

2 1) (

2  , (5)

where

h

 k

  . (6)

The truncation error of method –e is given as _M

   

, , 12

,

12 ⁴

2 4 2 2 4

4 4

x t h u

v k t

x k u

e_M ⁱ ^j ⁱ ^j



 



   

(7) for _i



x_i_₁,x_i_₁



, _j



t_j_1,t_j_1



, and i1,2,,n1 and j1,2,,m1.

The Dirichlet conditions (2), for j1,2,,m are in the form ,

, 0

,

0 _j u_n _j 

u (8)

To obtain the second initial condition in (3) function u in a third-degree Taylor polynomial about a point t₀ is used and becomes

         

, , 6

, 2

, ,

, ₃ ⁰

3 3 2

0 2 2 0 0

1 C

i i

i e

t t x u k t

t x u k t

t x k u t x u t x

u 





 



 







for i1,2,,n1, and if uC⁴



0,L



.

eC is a truncation error of the second initial condition (3) and is defined as

 

4 4

4 ,~

24 t

x u

e_C k ⁱ ⁱ





 

, (9)

where t₀ ~_i t₁ for i1,2,,n1.

In this case, the third order approximation of the derivative

 

0

,

 

 t t

t x

u is used, because the order of error e_C (9) is the same as error e (7), and estimation error _M of e will be true for estimation error of _M e . Then only one estimation is needed _C for both truncation errors.

If  and  exists and the initial conditions (3) are held by the wave equation (1) (see [1]), then for i1,2,,n1 and t₀ 0 we can write

  

ⁱ

  

ⁱ

 

_i

t

i v x

dx x v d x

t x v u t

t x

u  

 

 

 





2 2 2 2 2

0 2 2 0 2

0

2 , ,

, (10)

and

     

.

, ,

,

2 2 2 2 2 2 2

0 0 2

2 2 0 0 2

2

0 3

i i

i

t i

x dx v

x v d

x v x

t t x u v x

t t x u t t

t x u

 

   



 

 



 







 







 



 







 





 (11)

(4)

Summarizing, for i1,2,,n1 and j0, the Cauchy conditions (3) become

6 , 2

,

2 3 2 2 1

, 0 ,

C i i

i i i

i i

e k v

k v k u

u



















(12)

where the truncation error e_C is defined by formula (9).

4. THE CENTRAL DIFFERENCE INTERVAL METHOD

In conventional methods approximation errors e and _M e are usually ignored, _C but in the case of interval methods these errors must be estimated and involved into interval solution [3].

First, the estimations of errors e (7) and _M e (9) are needed. In the paper [9] _C the way of these estimations is presented in detail and leads to the formulas

  

^M ^M



H

E_M K 1 ,

12

2 2 2









 , (13)



M M



V

E_C K ,

24

2 2





 , (14)

for truncation errors e and _M e analogously. The value M is described by the _C formula

 

^,

2

4 5 max

. 1

1 , 1 , , 1 , 1

, 1 , 1 1 , 1 1 , 1 1 , 1 1 , , 2 , 1,2, , 1 2 1

2



















 











j i j i j i j i

j i j i j i j i j i m

j n

i

u u u u

u u

k u M h



(15)

where the terms u_i_,_j are calculated by the conventional central difference method.

Finally, for i1,2,...,n1 and j1,2,,m1, the central difference interval method we can write as follows



Γ²1



Ui,j Γ²



Ui_₁,j Ui_₁,j



Ui,j_₁Ui,j_₁EM

2 (16)

where E is given by formula (13). _M

The Dirichlet conditions (8), for j1,2,,m, are interval numbers given as ],

0 , 0

, [

,

0 _j U_n _j 

U (17)

and Cauchy conditions (12) for i1,2,...,n1 yield

6 , 2

,

2 3 2

2 1

, 0 ,

C i i

i i

i i i

E K V

K V K

U U



 



 

















(18) where E is defined by expression (14). _C

(5)

Above relations are satisfied for the mesh points (x_i,t_j)(X_i,T_j). Functions U,

,  and values H, K ,  are adequately interval extensions [6] of u, , , and h, k,



respectively.

To find the interval solution at the mesh points of the grid we have to solve the



ⁿ^¹

 

^ ^m^¹



system of interval linear equations for the



ⁿ^¹

 

^ ^m^¹



unknowns Uij. The formula (16) together Dirichlet (17) and Cauchy (18) conditions leads to obtaining the following system of interval linear equations









































m

m R

R R

U U U

I A I I A I

I A I





3 2 3

2

(19)

where I is identity matrix and A is the tridiagonal matrix in the form

 1

) 1 (   

































n n

A

) 1 - 2(

2 2



 (20)

Vectors U_j and R_j, for j2,3,,m, are adequately vectors of unknowns and constants and can be presented as follows



j j n j



^T

j U U U

U  ₁_, , ₂_, ,, _₁_, , Rj 



R₁_,j,R₂_,j,,Rn_₁_,j



^T (21) Elements of vector R_j, which contain the truncation errors E_M and E_C, are obtained by the following formulas

 

   



(( 1) ) (( 1) )



( ),

) ( ) ( ) 1 ( 2 )

) 1 ((

) ) 1 ((

: 2 , , 3 , 2

), ( )

2 ( ) 2 ( )

( ) ( ) 1 ( 2

2

2 2

2 ,

2 2

2 , 1

iH E

H i K H i

E iH K iH E

H i K H i E

R

n i

H E H K H E

H K H E

R

C

C C

M i

C C

M

















































































 for

 

. :

, , 5 , 4 , 1 , , 2 , 1

, ) ( ) ( :

1 , , 2 , 1

), ) 1 ((

) ) 1 ((

) 1 ( 2

) ) 2 ((

, 3 , 2

2 2

, 1

M j i

C M

i

C C M

n

E R m

j n i

E iH K iH E

R n

i

H n E

H n K H n

E H n K H n E

R































































 



 for

for

To solve the interval linear systems of equations (19) interval algorithm based on Gaussian elimination with complete pivoting is used.

All calculations are made on the computer on floating-point interval arithmetic in Delphi Pascal.

(6)

5. FLOATING-POINT INTERVAL ARITHMETIC

Computer implementations in floating-point interval arithmetic give solutions, in form of intervals, which contain initial-data errors, data representation errors and rounding errors. Presented central difference interval method contain error of method. So, interval methods together floating-point interval arithmetic give solutions, which contain all possible numerical errors. All calculations have been made using IntervalArithmetic unit, written by prof. Andrzej Marciniak, in the Delphi Pascal language. This unit allowed representing an input data and making all calculations on floating-point interval arithmetic.

For example real number 10 is represented as following interval number .



0.0999999999999999,0.1000000000000001



1 . 0 

where the diameter of this interval number is equals d=6.7710^-21.

6. NUMERICAL EXPERIMENTS

In an electric transmission line (the exercise is put in [1]) of length L that carries alternating current of high frequency, the voltage V and current i are described by

, 0 , 0

,

, 0 , 0

,

2 2 2

2

2 2 2

2







 







 

 



t L x t

LC i x

i

t L t x

LC V x

V

x x

where L=0.3 [henries/ft] (the inductance per unit length) and C=0.1 [farads/ft] (the capacitance per unit length). The line Lx=200 is feet long

The voltage V and current i satisfy following Dirichlet and Cauchy conditions

       

   



,0



0,



,0



0, 0 200.

, 200 0

200, cos 5 . 5 0 , 200,

sin 110 0 ,

, 0 , 0 , 200 ,

0 ,

0 , 200 ,

0



 

 











x t x

x i t V

x x x

x i x

V

t t

i t i t

V t V



The voltage V and current i at time t=0.5 are found.

To obtain the interval solutions the value M for estimation errors is needed. For voltage V the value M 0,00033, and for current i the value M 0,000017. Solutions of voltage V and current i, calculated for point (x,t)=(100.00, 0.25), (which have the same values in this example), are presented in Table 1.

The maximal widths of interval solutions for voltage V and current i at all grid points are presented on the Fig.1.

(7)

Table 1. Solutions for point (100.00, 0.25) and estimations of errors for voltage V and current i

U(100.00,0.25)

n=m real INT.left INT.right INT.width EC EM

10 109,9719606 109,9710809 109,9727476 1,67E-03 9E-11 3E-05 20 109,9717870 109,9715338 109,9720286 4,95E-04 5E-12 2E-06 30 109,9717547 109,9716298 109,9718762 2,46E-04 1E-12 3E-07 40 109,9717434 109,9716635 109,9718219 1,58E-04 3E-13 1E-07 50 109,9717382 109,9716773 109,9717984 1,21E-04 1E-13 5E-08 60 109,9717354 109,9716812 109,9717891 1,08E-04 ^7E-14 2E-08 70 109,9717337 109,9716765 109,9717905 1,14E-04 ^4E-14 1E-08 80 109,9717325 109,9652352 109,9782297 1,30E-02 ^2E-14 7E-09 90 109,9717318 -10,2657582 230,2092216 2,40E+02 ^1E-14 4E-09

Fig. 1. The maximal widths of interval solutions for voltage V and current i

7. CONCLUSIONS

Interval methods for solving wave equation on floating-point interval arithmetic give solutions, in form of intervals, which contain all possible numerical errors.

Presented interval method on floating point interval arithmetic guarantee correct digits in solutions. For example the in the interval solution:

X = [ 1.5707016101111265, 1.5707621013768904],

correct digits 1.5707 are guaranteed, while adequate real solution obtained by conventional method is x = 1.57073185574401, and we do not know, how many digits are corrected.

The main point of an interval method is to make good quality error estimation of discretization method. If the estimation is poor, the widths of interval-solutions cannot be satisfied. Truncation errors of method and initial condition are the same order in the presented method and both errors are depended on the same value M (it

(8)

is easier to make only one estimation). Additionally, proposed initial condition approximation affords to consider the larger value v in the wave equation until now.

Computer memory and calculation time (for large system of interval linear equations) are usually disadvantages of interval methods together floating point interval arithmetic.

REFERENCES

[1] Burden R. L., Faires J. D., Numerical Analysis third Edition, PWS, Boston, 1985.

[2] Davis J. L., Finite difference methods in dynamics of continuous media, Macmillan Publishing Company, New York, 1986.

[3] Kalmykov S. A., Šokin Yu. I., Yuldašew Z. H., Methods of Interval Analysis, Nauka, Novosibirsk 1986 (in Russian).

[4] Kącki E., Partial Differential Equations in physics and techniques problems, WNT Warszawa 1995 (in Polish).

[5] Marciniak A., Szyszka B., An interval method for solving some partial differential equations of second order, Proceedings on CD: 80th Annual Meeting of the International Association of Applied Mathematics and Mechanics, GAMM 2009.

[6] Moore R. E., Methods and Applications of Interval Analysis, SIAM, 1979.

[7] Szyszka B., An Interval Difference Method for Solving Hyperbolic Partial Differential Equations, Mathematical Methods and Techniques in Engineering and Environmental Science, MACMESE Conference, p. 330-334, WSEAS Press, 2011.

[8] Szyszka, B., Central difference interval method for solving the wave equation.

ICNAAM (International Conference on Numerical Analysis and Applied Mathematics) 2010. AIP Conference Proceedings, vol. 1281, pp. 2173-2176, 2010.

[9] Szyszka B., The Central Difference Interval Method for Solving the Wave Equation, Lecture Notes in Computer Science, Springer-Verlag, Berlin, (in print).

[10] Zill D. G., Differential Equations with Boundary-Value Problems, PWS Publishers, Boston, 1986.