EFFICIENT ALGORITHMS FOR SOLVING NONLINEAR VOLTERRA INTEGRO-DIFFERENTIAL EQUATIONS WITH TWO POINT BOUNDARY CONDITIONS

. In this paper a collection of efficient algorithms are described for solving an algebraic system with a symmetric Toeplitz coecient matrix. Systems of this form arise when approximating the solution of boundary value Volterra integro-differential equations with finite difference methods. In the nonlinear case, an iterative procedure is required and is incorporated into the algorithms presented. Numerical examples illustrate the results.


INTRODUCTION
The form of the system to be solved in this paper is AY F(Y) (I I) where F(Y) is a nonlinear function in an unknown vector Y and A is a symmetric Toeplitz matrix.In [1,4] circulant matrices and synunetric band matrices are discussed in connection with solving a linear system.In [7] there is an application of a fast algorithm to a system of the form (1.1).In [5] this algorithm involves the solution of a second order Volterra integro-differential equation.The background for that work and the work considered here is based on the p.113].Consider two square matrices A and B and two column vectors u and v related by A B nv T If the inverse ofB exists and vTB'*u, 1, then A "t B " + B'muvTB " I/(I-vTB'u) In this paper an application involving the numerical solution of a second order boundary value problem of the Volterra type will be discussed.In particular, the problem model will be nonlinear in the unknown and the computer implementation will employ the Sherman-Morrison formula.Two numerical examples will be given to compare this method with an efficient form ofthe LU method and a variant ofthe method described in [5] and [7].

I)
fand K are uniformly continuous in each variable ii) for the function f and for all (x,y,z), (x,y,z) and (x,y,z) in Rz, ]ffx,y,z)-f(x,y,z)l ly-YI lffx,y,z) ffx,y,z)l lzzl iii) for the function K and for all (x,t,y) and (x,t,) in R, IK(x,t,y) K(x,t,y)l s l=Jly-YI and iv) the functions fr' fz and K r are continuous and satisfy I(x,t,y)< 0 for all (x,t,y)6R and (x,y,z)6R2.
f(x,y,z)>0, fz(x,y,z)a0 and A proof ofthe uniqueness of a solution for problems satisfying the above can be found in Shaw [6].
To develop a discrete solution for (2.1) let IN={X.x--ih, I=0(1)N, Nh=a} be a partition of I=[0,a].A general k-step method of solution is given by where gi Ya*i h2 -. h %x(..,j,y?, > ,, "-o 0 j-O (2.2) The coefficients {wj} denote the weights for the choice ofa quadrature rule.I,n addition, there are special rules required in connection with starting values.These rules have weights w, msmax{k,s} and s is related to the order ofthe method.Note that when k=-2 the two boundary conditions provide the required auxiliary conditions and the (N+1)x(N+1) system can then be solved.In the linear case only one solve is required but in the nonlinear case an iterative scheme must be employed.
For the remainder ofthis section we consider a model ofthe form ya)(x) y(x) + g(x,y(x),z(x)) with boundary values y(0)=b and y(a)fb.Applying a two-step method to (2.3) gives rise to an (IN+l)x(N+1) system.Furthermore, if the method for solving the differential equation is denoted by (p,o) with then a particular method, known as Numerov's method, has (ao, a l,a 2)=(1, -2,1) and [3o, [3 ,13z)=(1/12,10/12,1/12).The quadrature rule Q selected as part of this particular method is the fourth order Newton-Gregory rule.Conditions for the convergence oftwo part methods ((p, o),Q) can be found in [6].To present the method in a way which is more appropriate for the next section, remove the terms involving Yo and Ys to the fight hand side to give an (N-1)x(N-1) system.This system can be expressed as Of particular interest is that C is a diagonally dominant symmetric Toeplitz matrix.

ALGORITHM DEVELOPMENT
To solve system (2.4) an iteration scheme must be employed.Using {i} to denote the = iterate gives There are several ways to solve system (3.1).For example, by efficiently using LU decomposition on the Toeplitz tridiagonal symmetric system the operation count is about 9n, where n is the size of the system.For repeated applications with the same coefficient matrix C and different fight hand sides the count is 5n.The form ofthe LU method used can be described as follows: 2) 3) 4) This method of efficiem LU decomposition will be referred to as Method 1.
Working with C, note that {c,I>{CO{ +lcxl, %=5 and the division by -c gives a matrix which is tridiagonal symmetric and diagonally dominant.The nonzero elements in each row are given by (-1, ., -1) with .=2+(10h2/12).Denoting this matrix by D gives Dyi.
_I 012 Let u=(-g)e] and v--e, where e=(1,O O)T.Perturb D such that 13 is given by Then it is easily seen that Dffibuv'.Method 2 uses an LU decomposition of f) from [4] as follows: For ease of notation, let K =h2(BG(yi)+s(Yi))+R represent the s iterate ofthe right hand side of equation (3.2).Let aj=('J-gJZ'2)d where d=g'(1-'a'+)) ", m=N-1.The steps in the method are given by 1) factor b into its LU decomposition and determine the aj, j= 1(1)N- 2) solve the system Z'= K where Z t= () 3) evaluate Y' (j) where ,+' -9 , J= (I)N- 4) if IY /' -lil > tolerance, repeat from step 2 There are just a few calculations for step and approximately 4n calculations for step 2. The remaining steps require about 5n operations for a total of9n.For repeated applications with the same 13 the operation count is about 6n.
As a third approach, Method 3 utilizes the method given in [2,p.113].To obtain an approximate solution the steps are as follows: factor15 imo its LU decomposition solve 13 W u solve 13 Z Y calculate the constant 13--vZ/(1-vTW) and set Y Z+I3W if IY /' -1 > tolerance, repeat from step 3 As before, step essentially requires just a few operations and step 2 requires 4n operations.Since 1,0 0)3, step 4 requires 2n operations.The remaining steps give a total of 9n operations that reduces to 6n with repeated iterations.
To close this section, we will look at the convergence of the three methods.For the nonlinear problems, each method is iterative in nature.For the direct LU method, the application of the method is equivalent to solving Y+'= A" (hX(BG(Y) + S(Y)) + R).
Subtracting the same equation with the exact solution ofthe discrete system being used gives us where L is a Lipschitz constant and convergence takes place provided h2LI,IIB|<I.Note that A is a continuous function of h which tends to J as h-,0.The matrix -A is monotone [3,p.360] and hence is invertible.The inverse of a monotone matrix has nonnegative elements.Since -J and -A are monotone and (-A)-(-J)0 we see that (-J")-(-Aa)0 [3,Theorem 7.5, p.362].It follows, therefore, that 0s(-A'l)s(-J).
From [3] it is also known that IFI N2/8 and therefore Il is bounded.
The other two methods are somewhat similar so only the last one will be analyzed.There are two parts to the iteration process.In the first part, an approximation is obtained.By adding an appropriate correction, the next iterate Y is given by Y/ =b "' K + = 15 "vb -' K%. (3.4) Again, exact values for the solution of the discrete problem are substituted and the result is subtracted from (3.4) to give In particular, vTb "* h2B(G(Y)-G(Y)+S(Y)-S(Y)) is just a constant because ofthe definition of v. Hence, equation (3.5) can be written as E i+' (I -* + x I -1 nv 1 -1 )h2B(G(yi).G(Y)+S(yi).S(y))

NUMERICAL RESULTS
Two Volterra integro-differential equations of the form (2.1) are used to illustrate the methods described.Method is the efficient form ofthe LU method, Method 2 is a variant of the fast algorithm given in [4] and [6] and Method 3 is the implementation ofthe Sherman-Morrison formula [see, 2].The iteration count for all methods is given along with the average CPU time in seconds.AlL programs were nan on a SUN SPARC 20 in double precision arithmetic.Another approach to solving system (3.1) is given in the paper by Yah and Chung [7].Their modification reduces the operation count for step 3 from 5n to 4n+2t where represents the number of operations in the correction term approximation (tsn).The actual value of depends on the degree of the diagonal dominance ofthe coefficient matrix C. To describe their algorithm, one can replace the approximation to the correction term in step 3 ofMethod 2 by y i.1 Z Let b=t "].The vector p is given by pfCu,b b',0 0) T where the number of components in p is dependent on [.1 and a tolerance %.In particular, % is chosen such that tog(ll 2) t() > togClbl)

Mathematical Problems in Engineering
Special Issue on Time-Dependent Billiards

Call for Papers
This subject has been extensively studied in the past years for one-, two-, and three-dimensional space.Additionally, such dynamical systems can exhibit a very important and still unexplained phenomenon, called as the Fermi acceleration phenomenon.Basically, the phenomenon of Fermi acceleration (FA) is a process in which a classical particle can acquire unbounded energy from collisions with a heavy moving wall.This phenomenon was originally proposed by Enrico Fermi in 1949 as a possible explanation of the origin of the large energies of the cosmic particles.His original model was then modified and considered under different approaches and using many versions.Moreover, applications of FA have been of a large broad interest in many different fields of science including plasma physics, astrophysics, atomic physics, optics, and time-dependent billiard problems and they are useful for controlling chaos in Engineering and dynamical systems exhibiting chaos (both conservative and dissipative chaos).We intend to publish in this special issue papers reporting research on time-dependent billiards.The topic includes both conservative and dissipative dynamics.Papers discussing dynamical properties, statistical and mathematical results, stability investigation of the phase space structure, the phenomenon of Fermi acceleration, conditions for having suppression of Fermi acceleration, and computational and numerical methods for exploring these structures and applications are welcome.
To be acceptable for publication in the special issue of Mathematical Problems in Engineering, papers must make significant, original, and correct contributions to one or more of the topics above mentioned.Mathematical papers regarding the topics above are also welcome.