Algorithm with guaranteed accuracy for computing a solution to an initial value problem for linear difference equations

Applied Mathematics

Algorithm with guaranteed accuracy for

computing a solution to linear dierence

equations

Vladimir Vaskevich

1?

, Haydar Bulgak

2

, Cengiz Cinar

2 1 Sobolev Institute of Mathematics, SB RAS, Novosibirsk, Russia

e-mail:vask@math.nsc.ru

2 Research Centre of Applied Mathematics, Selcuk University, Konya, Turkey

e-mail:bulgk@karatay1.cc.selcuk.edu.tr

Received: August 31, 2000

Summary.

Consider an initial value problem for simultaneous lin-ear dierence equations x(n+ 1) = Ax(n) +f(n), x(0) = a, with

A the N N rational matrix, ff(n)g a sequence of N-dimensional

rational vectors, andaa rationalN-dimensional vector. The problem has a unique solution but to compute fx(n)g in the interval 0M]

with M a nonnegative integer, we approximate the reals and carry out the elementary arithmetic operations in special way. By means of this algorithm, we solve the initial value problem for a discrete asymptotically stable matrixAwith guaranteed accuracy.

Key words:

initial value problems, simultaneous linear dierence equations, algorithms with guaranteed accuracy, discrete asymptoti-cally stable matrices

Mathematics Subject Classication (1991): 65F30, 65G10

1. Introduction

Let , P;, P+, and k be integers  > 1. In Forsythe, Malcolm,

and Moler 7], the special nite subsets

F

(P;P+k) of rational ? This work was supported in part by TUBITAK within the framework of

numbers were introduced as follows

F

(P;P+k) = f0gfz2Rjz= p k X j=1 mj ;j wherepmj 2Zp2 P ;P+] mj 2 0;1]m 1 >0 g:

The sets

F

(P;P+k) are called formats. In this paper, we use

the following properties of formats (see, e.g., 10], 2], 3], and 9]). The number 1(

F

) = 

P+(1 ;

;k) is the greatest member of

F



0(

F

) = 

P;;1 is the least positive element of

F

. If 

1(

F

) =  1;k,

then there are no members of

F

in (11+1(

F

)) but the sum 1+1(

F

)

belongs to

F

. If

F

(P;P+k) 

G

(P  ;P  +k ), then

G

is called

athin formatwith respect to

F

. In this case we also write that

F

is a

rough formatwith respect to

G

. Before the arithmetic operations can be performed on R, we have to choose a format for approximations of the reals and in the sequel we refer to the chosen format as the standard format. Given a real number z and a format

F

, we dene z]F as the member of

F

which is the nearest measured by

jjto z.

To each format

F

, we assign also

F

-approximations a]F ( a

i]F) of

a real vector aand A]F ( A

ij]F) of a real matrixA.

Given an integer M and a set

x

= (x(0)x(1):::x(M)) of real vectorsx(n)2R

N, we use the norm jj

x

jj M  max 0k M jjx(k)jj 2

with a double bar notationsjjjj

2 to designate the Euclidean length

of a vector inRN and the spectral norm of a matrix. LetAbe a real

matrix. We dene the quality !(A) of discrete time stability by the formula !(A) = sup T>0jjx(0)jj 2 =1 T X n=0 jjx(n)jj 2 2

with x(n) a solution to the following simultaneous linear dierence equations

x(n+ 1) =Ax(n) n= 01::::

IfAis a discrete asymptotically stable matrix, then the value of!(A) is nite (see, e.g., 1]). We also use the Euclidean norm kAk

E of A dened as followskAk E = ( N P ik =1 ja ik j 2)1=2.

Let A = (Aij) be the N N rational matrix and let a = (ai)

be a rational N-dimensional vector. We also consider a set

f

= (f(0)f(1):::f(M;1)), wheref(n) is an N-dimensional rational

vector. We assume that the entries ofA,f(n), and aare members of the standard format

F

 i.e.,A= A]F,f(n) = f(n)]F, and a= a]F.

Given a positive integerM, we consider the initial value problem (1:1) x(n+ 1) =Ax(n) +f(n) x(0) =a n= 01:::M;1:

Problem (1.1) has a unique solution (see 6, p. 124] and 11, p. 767]) (1:2) x(n) =Ana+ n;1 X k =0 An;1;kf(k) n= 12:::M: If !(A)1=2 jjajj 2+ 2!(A) 3=2 jj

f

jj M 

1(

F

), then the entries of x(n)

are in the interval (;

1(

F

)1(

F

)).

Since the entries ofA and f(n) are members of

F

, it follows that they also are the members of

G

, a thin format with respect to

F

. Given

F

, we choose

G

such that

(1:3) N1(

G

) 1 4 1(

G

)! 3=2(A)  1(

F

) 41(

G

)! 2(A) 1 (N+1)1(

G

)! 3=2(A) kAk E 1 16 120(

G

)  !(A)(N+1) 3=2  0(

F

):

Then we use Algorithm 1 to compute approximations y(n) of the solutionx(n) to (1.1).

Algorithm 1

Initial Value Problem for Linear Dierence Equation

M andNfgiven positive integersg

G=G(P;P+k) andG=G(2P;;12P+2k)finitial formatsg A= A] G= ( A ij) fgiven a matrixA2G NN g f= f] G= ( fj) anda= a] G= ( aj)fgiven vectorsf2G N and a2G N g y(0) =a for n= 1 toM do computey(n) =Ay(n;1) +f(n;1) as follows. z=y(n;1) for j= 1 toN do S 0j= f j ff

j is the rational entry of

f(n;1)g for k= 1 toNdo

compute the product A j k z k fA j k and z

k are rational numbers g

replaceAj kzkby the memberk j = Aj kzk] G of G computeh k j = S k ;1j+  k j fS k ;1j and 

k j are rational numbers g replaceh k j by the member S k j =  h k j]G of G  endfor

replace the rational numberSN j by the memberHj= SN j] Gof G. endfor y(n) =HfwithH= (H 1 :::H N) g

We emphasize that arithmetic operations in Algorithm 1 are per-formed on the members of

G

, a thin format with respect to

F

.

It is important to know how well y(n) agrees with x(n). To this end, we establish the following

Theorem 1.

If (1.3) holds, then

(1:4) jj

x

;

y

jj M 2  1(

F

) + 8(N+ 1)1(

G

)! 3=2(A) kAk E  jj

x

jj M + 161(

G

)! 3=2(A)(N+ 1) jj

f

jj M;1+0(

F

)

where

x

= (x(n))is the solution to (1.1) and

y

= (y(n))is the vector computed by Algorithm 1.

If the entries ofA are not rational numbers, then we refer to 5].

2. The estimates of the error

Proof of Theorem 1. Consider the set

g

= (g(0)g(1):::g(M;1))

of the vectorsg(n) =y(n+1);Ay(n);f(n). It is not hard to show

(2:1) jj

x

;

y

jj M 2! 3=2(A) jj

g

jj M;1: To estimatejj

g

jj

M;1from above, we turn to Algorithm 2.

Algorithm 2

The inner product (xy)

Nfa positive integerg G=G(P ; P + k) andG = G ( 2P ; ;12P + 2k)finitial formatsg x= x] G= ( xk) andy= y] G= ( yk)fgiven vectorsx2G N and y2G N g

compute the inner product (xy) as follows. S0= 0

for k= 1 toN do

compute the productx k y k fx kand y

kare rational numbers g

replacexkyk by the memberk= xkyk] G of G computeh k= S k ;1+  k fS k ;1and 

k are rational numbers g

replacehk by the memberSk= hk] G of

G endfor

replace the rational numberSN by the memberH= SN] Gof G. Put xy] GG = H. Given xy2 G N, we use xy] GG

 to designate the real number H

If jjxjj 2  1(

G

)=2 and jjyjj 2 

1(

G

)=2, then the following

inequality (2:2)   (xy); xy] GG      1(

G

) j(xy)j+ 0(

G

) + N1(

G

) 1;N 1(

G

)=2 N X j=1 jx j jjy j j+ N0(

G

) 1;N 1(

G

)=2 holds 8, p. 445]. By (1.3), we infer 1 1;N 1(

G

)=2 1 1;N 1(

G

)=2 2:

This inequality, together with (2.2), yields

(2:3)   (xy); xy] GG      1(

G

) j(xy)j+ 0(

G

) + 2N1(

G

) N X j=1 jx j jjy j j+ 2N 0(

G

):

Theith entry ofAy(n)+f(n) equals the inner product of the vector

af(i) = (a

i1ai2:::aiNfi) 2

G

N+1

and the vector

y(N+1)(n) = (y(n)1) 2

G

N+1:

By the denition ofy(n), we also have

yi(n+ 1) = af

(i)y(N+1)(n)]

GG i= 12:::N:

Hence, the equality (2:4) jjg(n)jj 2 2= N X i=1   (af (i)y(N+1)(n)) ; af (i)y(N+1)(n)] GG     2

holds. From (2.3) and (2.4) it follows that (2:5) jjg(n)jj 2  1(

G

) kAy(n) +f(n)k 2+ (n) where (2:6) (n) = 2(N + 1)1(

G

)  kAk E ky(n)k 2+ kf(n)k 2  + 2(N + 1)3=2 0(

G

) + p N0(

G

) n= 12:::M ;1:

Apply rst (2.5) and next (1.1). Then, by easy calculations, we infer jjg(n)jj 2  1(

G

) kx(n+ 1)k 2+1(

G

) kAk 2 ky(n);x(n)k 2+ (n):

Whence and from (2.1) we have (2:7) jj

x

;

y

jj M 2 1(

G

)! 3=2(A) jj

x

jj M + 21(

G

)! 3=2(A) kAk 2 jj

x

;

y

jj M+ 2! 3=2(A) max 0nM;1 (n):

Bearing in mind thatkAk 2 ! 1=2(A) 5] and inserting (1.3) in (2.7), we get (2:8) jj

x

;

y

jj M  1(

F

) jj

x

jj M+ 4! 3=2(A) max 0nM;1 (n):

By denition (2.6) of (n), we conclude that max 0nM;1 (n) 2(N+1) 1(

G

)  kAk E k

x

k M+ k

f

k M;1  +p N0(

G

) + 2(N+ 1)1(

G

) kAk E k

x

;

y

k M + 2(N+ 1) 3=2 0(

G

):

Inserting this estimate in (2.8), we obtain (2:9)  1;8(N + 1) 1(

G

)! 3=2(A) kAk E  k

x

;

y

k M  1(

F

) + 8(N+ 1)1(

G

)! 3=2(A) kAk E  k

x

k M + 8(N+ 1)1(

G

)! 3=2(A) k

f

k M;1+ 4 p N!3=2(A) 0(

G

) + 8(N+ 1)3=2!3=2(A) 0(

G

):

Since (1.3) holds, it follows that 1;8(N + 1) 1(

G

)! 3=2(A) kAk E 1=2 4!3=2(A)  p N0(

G

) + 2(N+ 1) 3=2 0(

G

)   0(

F

):

Substituting these inequalities in (2.9), we have (2:_{10) 12}jj

x

;

y

jj M  1(

F

)+81(

G

)! 3=2(A) kAk E(N+1)  jj

x

jj M + 81(

G

)! 3=2(A)(N + 1) k

f

k M;1+0(

F

):

References

1. Akin, O. and Bulgak, H. (1998): Linear Dierence Equations and Stability Theory in Turkish], Selcuk University, Research Centre of Applied Mathe-matics, Konya.

2. Bulgak, A. (1998): Checking of a well-conditioning of an interval matrix, Siberian J. of Dierential Equations 3, No. 4, 75{79.

3. Bulgak, H. (1999): Pseudoeigenvalues, spectral portrait of a matrix and their connections with dierent criteria of stability, in Error Control and Adaptivity in Scientic Computing, (Bulgak, H. and Zenger, C., Eds.), Kluwer Academic Publishers, 95{124.

4. Bulgakov, A. Ya., and Godunov, S. K. (1988): Circle dichotomy of the matrix spectrum, Siberian Math. J.29, 734{744.

5. Bulgak, A., Bulgak, H., and Vaskevich, V. L. (2000): Computing an initial value problem for systems of linear dierence equations with error estimates, in Cubature Formulae and their Applications, (Noskov M. V., Ed.), Krasno-yarsk, 238{258.

6. Elaydi, S. N. (1996):An Introduction to Dierence Equations, Springer, New York.

7. Forsythe, G. E., Malcolm, M. A., and Moler, C. B. (1977):Computer Methods for Mathematical Computations, Englewood Clis, N. J., Prentice-Hall. 8. Godunov, S. K., Antonov, A. G., Kiriluk, O. P., and Kostin V. I. (1993):

Gua-ranteed Accuracy in Numerical Linear Algebra, Kluwer Academic Publisher, Dordrecht.

9. Golub, G. H., and Van Loan, C. F. (1989):Matrix Computations, The John Hopkins University Press.

10. Kulisch, U. W., and Miranker, W. L. (1986): The arithmetic of the digital computer: a new approach, SIAM Review28, 1{40.

11. Myskis, A. D. (1979):Advanced Mathematics for Engineers. Special Courses, Mir Publishers, Moscow.