UNIVERSITY OF 
 
 ILLINOIS LIBRARY 
 
 AI UR3AINA- CHAMPA 'GM 
 
 ENGINEERING 
 
NOTICE: Return or renew all Library Materials! The Minimum Fee tor 
 each LOM BOOK i, $50.00. J^O? ]m 
 
 The person charging this material is'twpbnsible for 
 its return to the library from which it was withdrawn 
 on or before the Latest Date stamped below. 
 
 Theft, mutilation, and underlining of books are reasons for discipli- 
 nary aclloifcand may result in dismissal front fta University. 
 To renfew c&JI Telephone Center, 333-8400 . 
 
 i..NI',LkSlTv OF ILLINOIS MBRAFfy;'/f'4U^ANA-CHAMPAIGN 
 
Digitized by the Internet Archive 
 
 in 2012 with funding from 
 
 University of Illinois Urbana-Champaign 
 
 http://archive.org/details/investigationofo75belf 
 
ENGINEERING LIBRARY 
 UNIVERSITY OF ILLINOIS 
 
 UR8ANA, ILLINOIS 
 
 enter for Advanced Computation 
 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 
 
 URBANA, ILLINOIS 61801 
 
 CAC Document No. 75 
 
 Investigation of One-Step Methods for 
 Integro-Differential Equations 
 
 "by 
 
 Geneva G. Belford 
 and 
 Surender Kumar Kenue 
 
 May 1973 
 
CAC DOCUMENT NO. 75 
 
 INVESTIGATION OF ONE-STEP METHODS 
 
 FOR 
 
 INTEGRO-DIFFERENTIAL EQUATIONS 
 
 by 
 
 Geneva G. Belford 
 and 
 Surender Kumar Kenue 
 
 Center for Advanced Computation 
 University of Illinois at Urbana-Champaign 
 Urbana, Illinois 6l801 
 
 May 1973 
 
 This work was supported in part "by the Advanced 
 Research Projects Agency of the Department of 
 Defense and was monitored by the U.S. Army Research 
 Office-Durham under Contract No. DAHC0U-72-C-0001. 
 
voiNEERIMfi U8RAW 
 
 ABSTRACT 
 
 One-step methods for solving integro-differential equations 
 are studied from the point of view of desiring that the method give 
 good accuracy when the true solution asymptotically goes to zero. A 
 test equation is proposed and absolute stability is investigated. 
 
TABLE OF CONTENTS 
 
 Page 
 
 1. INTRODUCTION 1 
 
 2. STABILITY OF EULER'S METHOD 2 
 
 3. GENERALIZED EULER METHODS 8 
 
 k. REFERENCES 11 
 
• 
 
1. INTRODUCTION 
 
 Behavior of a number of physical and biological systems may 
 he described by integro-differential equations of the form 
 
 (1) y (x) = F (x, y(x), J K (x, t, y(t)) dt ) . 
 
 o 
 
 A familiar example is the population equation [2] 
 
 2 , X 
 
 (2) y (x) = ay(x) - by (x) - y(x) J K(x-t) y (t) dt , 
 
 in which y represents a population and x represents time. The integral 
 
 factor takes into account "heredity", or in general the past history of 
 
 the system. Thus in applications to mechanics such integral factors 
 
 may provide for component fatigue. 
 
 In spite of their importance, only very recently has work 
 
 been done on the development of efficient and accurate numerical methods 
 
 for solving these equations. The simplest methods proposed [l, 3] are 
 
 the one-step methods, analogous to Euler's method and its generalizations, 
 
 in which the left hand side of (l) is replaced by a first-order 
 
 difference y ( x+h ) - y ( x ) and the right side is evaluated at x (or, 
 h 
 
 more generally, is some convex linear combination of evaluations at 
 x and x+h). Beginning with a given initial value y(o), one then carries 
 out the familiar step-by-step process of determining an approximate 
 solution at points nh, n = 1, 2, ... . 
 
 Multistep methods have also been proposed [U, 6], but in the 
 interests of simplicity we have restricted ourselves to one-step methods 
 in this preliminary study. In Section 2 we look at some stability 
 questions for Euler's method and in Section 3 we consider possible 
 advantages of using a modified Euler's method. 
 
2. STABILITY OF EULER'S METHOD 
 
 In a recent paper, [l], Cryer and Tavernini discussed an 
 Euler' s -method approach to the computation of solutions to Volterra 
 functional differential equations. Under very general hypotheses they 
 show that the method is stable in the sense that the errors are bounded. 
 Just as for ordinary differential equations, however, one would expect 
 that a stronger stability condition is required in order to compute 
 valid solutions to a functional differential equation with a decreasing 
 exponential solution. We here investigate this question in a particularly 
 simple context. * 
 
 Since it is impossible to investigate the stability properties 
 of a method as it would apply to any equation, it is common in studying 
 methods for ordinary differential equations to introduce the notion of 
 a "test equation" — a simple equation for which' the method is easily 
 analyzed to yield information on how the method will work for a large 
 class of "similar" equations. The standardly used test equation for 
 ordinary differential equations is 
 (3) y (x) + ay = (a>o) 
 
 with initial condition y(o) = 1« Nothing this simple will provide a 
 valid test for an integro-differential equation method, since an equation 
 involving both an integral and a derivative is essentially analogous to 
 a second-order differential equation. With this in mind, we decided on 
 the test equation 
 
 y + ay + (a 2 /U) / y(t) dt = 
 (h) 
 
 y (o) = 1. 
 
For a>0, this has the solution 
 
 (5) y = exp (- ax/2) [l - ax/2]. 
 
 Although this is not the simple decreasing exponential obtained for 
 
 (3), it does have the desired testing behavior of approaching zero 
 
 asymptotically as x-* + °°. 
 
 To solve (h) numerically, we define x = nh (n = 0, 1, ...) 
 
 n 
 
 and y(x ) = y . Using Euler's method together with the trapezoidal 
 rule for the quadrature, we then obtain the following numerical 
 approximation to (h) . 
 
 y = (l - ah)y 
 
 (6) y 2 = 6y 1 - (y/2)y 
 
 
 
 n-1 
 
 y n+l = 6y n " Y I y j " (Y/2) y (n>l) 
 
 where 6=l-ah-ah/8 
 
 2 2 
 and y = a- h /h. 
 
 In order to study error propagation, we assume that y n is exact, and 
 that an error E is introduced into y . This error then propagates 
 according to 
 
 (7) E 2 = 6E , 
 
 n-1 
 
 E ^ = 6E - y I E . 
 n+1 n L J 
 
Letting A = ah, one defines the region of absolute stability of the 
 
 method as the set of complex values of A for which E •*■ as n ■*■ **». 
 
 n 
 
 Since the order of the difference equation (T) changes with 
 
 n, we use an indirect approach to determine the behavior of E . Define 
 
 n 
 
 the accumulated error 
 n 
 
 T = Ye,. 
 
 n L j 
 
 Then, from (T), T satisfies the simple 3-term recursion 
 
 (8) T = (1+6) T - (y+6) T (n>l), 
 
 n+I n n-1 — 
 
 with T ■ 0, T = E . 
 
 The solution to equation (5) is 
 
 (9) T = E l 
 
 n 
 
 [(l +5 ) 2 - My+6)] 1/2 
 
 {f, - 4} 
 
 where 
 
 (10) 
 
 ((1+5) ± [(1+6) 2 - My+5)] 1/2 J 
 = 1 - (A/2) - (A 2 /l6) ± (A/l+)[A+A 2 /l6] l/2 . 
 
 q ± 2 = 1/2 \ (1+6) ± [(1+6)" - U(y+6)] 
 
 It is clear that if both q_ and |q_ are less than one, then T ■*» as 
 
 1 n l ' ' ^2 ' n 
 
 n-»- °° and E must approach zero also. If we now write E = T - T , in 
 n n n n-1 
 
 terms of q_ and q and consider the various possibilities, it becomes 
 
 clear that the conditions I q_ |<1, I q_ I <1 are also necessary for E ■+ 0. 
 
 1 TL ' ' 2 ' n 
 
 Thus the region of absolute stability is the domain in which both | q. | 
 and |q-| are less than one. Looking at real values of A, we see that 
 the condition for this is 0<A<2. Computer trials for A in the vicinity 
 of 2 confirm that A = 2 is the dividing point between error growth and 
 
error decay. To determine the situation for complex A, we have computed 
 \n |, |q I on a grid of values in the complex plane. The region on 
 which they are both less than one is only slightly distorted from being 
 the disk | X - 1 1 < 1 (See the Figure). 
 
 For comparison, recall that the region of absolute stability 
 for Euler's method applied to the ordinary differential equation (3) is 
 precisely the disk |A-1|<1, where again A = ah. Intuitively, one might 
 feel that addition of an integral operator term and the accompanying 
 summing of errors (see (?)) would have a serious effect on the region of 
 absolute stability. Our analysis shows that this is not the case. 
 
 A very closely related question is whether the numerical 
 solution y approaches zero as n -> °° as it ought. This question may be 
 investigated by a very similar procedure. Define 
 n 
 
 (11) \ = Uy 
 
 5=0 
 
 Then, using (6) with y = 1, we obtain 
 
 (12) Y n+1 = (1+5 )Y n - (Y+<5)Y n _ 1 + (y/2) (n>l) 
 
 with initial conditions Y = 1, Y. = 2 - A. This is an inhomogeneous 
 
 o 1 
 
 difference equation with solution 
 
 (13) Y n " 2l^T «3-2X-q 2 V - (3-2X- 9l )q 2 n } + 1/2. 
 
 where q and q are given by (10). If both |q | and | q ? | are less than 
 
 one, Y -+■ 1/2 as n -> °°, and the convergence of the series 
 n 
 
 oo 
 
 I y. implies that y . -> o as j ■> °°. Thus we see that, as commonly 
 3=o J J 
 
Figure . The broken line denotes the "boundary of the region of 
 absolute stability in the complex X plane. For comparison, the solid 
 line is the circle \\ - ll =1. 
 
occurs in the case of ordinary differential equations, the numerical 
 solution has the proper asymptotic behavior within the region of 
 absolute stability. 
 
 In summary, it appears that the simple Euler's method is 
 perhaps not as useful as a universal method as Cryer and Tavernini 
 seem to suggest. The lack of absolute stability in some situations 
 is a serious problem. Furthermore, there seems to be no straightforward 
 way to determine a parameter A which, for a general integro- differential 
 equation, plays roughly the same role as the ah in the test equation. 
 That is, there is much more work to be done on prediction of whether 
 Euler's method is likely to be absolutely stable for a particular 
 equation and a particular step size h. 
 
3. GENERALIZED EULER METHODS 
 
 In this section we briefly consider the class of 
 one-step methods which may be written 
 
 (Ik) y n+l 
 
 - = (l-y) F .. + yF , 
 n+1 n 
 
 where 0<y< 1 and F represents the right side of (l) evaluated at x . 
 n n 
 
 (When F involves an integral, the "evaluation" involves an approximate 
 quadrature, of course.) For y = 1 we have the simple Euler's method 
 discussed in the last section. When y<l the method becomes implicit, 
 but the increase in computational complexity may be offset by improved 
 stability properties (i.e. stability for larger values of h). In the 
 case of ordinary differential equations this is certainly true. 
 Applying (lM to the test equation (3) it is a trivial exercise to show 
 that the method is "A-stable" (absolutely stable for all X such that 
 ReX>0) whenever 0<y<— . 
 
 It seems very reasonable that this result will also 
 extend to integro-differential equations. In an attempt to do this, the 
 approach of Section 2, with the same test equation, has been applied to 
 the generalized methods. One readily obtains a stability condition of 
 the form |q |, |q_|<l; the difficulty arises in trying to determine for 
 what values of X these conditions hold. Recall that in the last section 
 the boundary of the domain of absolute stability was found by computing 
 values of |q |, |q p |. Such a computational approach has no hope of 
 demonstrating absolute stability for half of the complex plane. Some 
 other method must be devised. When X is real, however, it is possible 
 
to show analytically that absolute stability holds for all A>0 when 
 0<y<-. 
 
 An advantage of having a class of methods depending upon a 
 parameter (y) is that this parameter may sometimes "be adjusted so as to 
 improve the numerical solution in some sense. Thus Liniger and 
 Willoughby [5] have used the free parameter y in order to fit the 
 numerical solution to exponentially-decreasing behavior and thus obtain 
 improved solutions to "stiff" ordinary differential equations. 
 
 We have done some experimentation on this sort of idea for 
 the case of integro-differential equations. A difficulty is that the 
 solution (5) to the test equation (h) is initially almost a straight 
 line, the exponential dominance appearing later. Hence it is hard to 
 see how a y might be chosen to mimic the step-by-step behavior of the 
 true solution throughout. Believing it of overriding importance to 
 begin with a good numerical solution, we have opted to use y to attempt 
 to get an exact solution at the point x = h. 
 
 Evaluating the exact solution (5) at h we get 
 
 (15) y(h) = (l-A/2) exp(-A/2), 
 
 while the numerical method (again with the trapezoidal rule used for 
 the quadrature) yields 
 
 f^\ - (l-Ay-A 2 (l-y)/8) 
 
 1 + A(l-y) + A (l-y)/8 
 
 Equating the expressions in (15) and (l6), we find that 
 
 (IT) (l-A/2) exp(-A/2) - 1 + A 
 
 A{(l-A/8) - (l-A/2) (l+A/8) exp(-A/2)} 
 
 = 1-y 
 
10 
 
 As A increases, y rapidly decreases to zero (when X is about k) and 
 turns negative. Hence only a very small range of A values produce 
 meaningful values of the parameter y from equation (IT). There seems, 
 therefore, to be very little value in pursuing this approach. It does 
 appear to indicate, however, that if some member of the class of 
 generalized Euler formulas (lU) is to be chosen for routine use, the 
 best choice might well be y = (backward Euler), both for initial 
 accuracy and long-term stability. It should be noted that all members 
 of the class should have 0(h ) local truncation error and thus should 
 be first order. (See [l] for a detailed error analysis for Euler's 
 method (y = l).) It might well be worthwhile to develop higher order 
 single-step methods for integro-differential equations. 
 
11 
 
 REFERENCES 
 
 1. Cryer, C. W. and Tavernini, L. s The numerical solution of Volterra 
 functional differential equations by Euler's method, SIAM 
 J. Numer. Anal. 9 (1972), pp. 105-129. 
 
 2. Davis, H. T. , Introduction to Nonlinear Differential and Integral 
 Equations, United States Atomic Energy Commission, I960, 
 Chapter 13. 
 
 Day, J. T. , Note on the Numerical Solution of Integro-differential 
 Equations, The Computer Journal, 9 (1967) pp. 39^-5. 
 
 Kemper, G. A. , Linear Multistep Methods for a Class of Functional 
 Differential Equations, Numer. Math. 19 (1972), pp. 361-372. 
 
 5. Liniger, W. and Willoughby, R. A., Efficient Integration Methods 
 
 for Stiff Systems of Ordinary Differential Equations, 
 SIAM J. Numer. Anal. 7 (1970), pp. 1+7-66. 
 
 6. Linz, P. , Linear Multistep Methods for Volterra Integro-differential 
 Equations, J. Assoc. Comp. Mach., l6 (1969), pp. 295-301. 
 
UNCLASSIFIED 
 
 Security Classification 
 
 DOCUMENT CONTROL DATA -R&D 
 
 f Security elmmalllemllon ol till: toay of abaraacf and inatea* 
 
 itiWw a»K ba mmtmrmd arbaw Ina ovurmll tmport la elm filled 
 
 1. originating activity fCarporala author) 
 
 Center for Advanced Computation 
 University of Illinois at Urb ana- Champaign 
 Urbana, Illinois 6l801 
 
 la. REPORT SECURITY C L A * SI F I c * T I Of 
 
 UNCLASSIFIED 
 
 2b. GROUP 
 
 J REPORT TITLE 
 
 Investigation of One-Step Methods for Integro-Differential Equations 
 
 4. DESCRIPTIVE NOTES (Type o/ report mtd rncluelre dmtea) 
 
 Research Report 
 
 ». AUTHORISI (Firm I mm, mlddlm Initial. Imml name) 
 
 Geneva G. Belford and Surender Kumar Kenue 
 
 • REPORT DATE 
 
 May 1973 
 
 7*. TOTAL MO. OF PACES 
 
 15 
 
 7b. NO. OF REFS 
 
 6 
 
 •a. CONTRACT OR CRANT NO. 
 
 DAHC04-72-C-0001 
 
 b. PROJECT NO. 
 
 AREA Order No. 1899 
 
 KS» 
 
 CAC Document No. 75 
 
 OTHER REPORT NOIJI (Any emSer number* tfiaf may be amelfned 
 thle report) 
 
 10 DISTRIBUTION STATEMENT 
 
 Copies may "be requested from the address given in (l) above, 
 
 II. SUPPLEMENTARY NOTES 
 
 None 
 
 12. SPONSORING MILITARY ACTIVITY 
 
 U.S. Army Research Office - Durham 
 Duke Station, Durham, North Carolina 
 
 13 ABSTRAC T 
 
 One-step methods for solving integro-differential equations 
 are studied from the point of view of desiring that the method give 
 good accuracy when the true solution asymptotically goes to zero. A 
 test equation is proposed and absolute stability is investigated. 
 
 DD ,'°?..1473 
 
 UNCLASSIFIED 
 
 Security Classification 
 
UNCLASSIFIED 
 
 Security Classification 
 
 KEY WOROS 
 
 no LI WT 
 
 Integro-differential equations 
 Euler's method 
 Absolute stability 
 
 UNCLASSIFIED 
 
 Security Classification