Recurrence Relations for Ordinary Differential Equations 

Newton's identities can be used to compute the coefficients of multistep recurrence formulae for simulating the solutions of ordinary linear differential equations of high order. This approach gives the exact rootmatched recurrence formula without having to determine the complex roots of the characteristic polynomials. 

1.0 Introduction 

Let x(t) and y(t) denote continuous functions of time related according to the Nth order differential equation 

_{} 

where the coefficients a_{k} and b_{k} are constants. Simple linear relationships of this form are widely used in numerical simulations of heat transfer, chemical reactions, fluid flows, mechanical dynamics, Markov processes, etc. [1] (Relationships of this form are also used in signal filtering and controlsystem compensation, but in those applications the frequency response is usually emphasized over the exact timedomain response.) A variety of methods (cf [2,3]) have been used to numerically generate the solutions y(t) of such equations for arbitrary forcing functions x(t). In the digital environment equation (1) is often simulated by means of a multistep linear recurrence relation of the form 

_{} 

where X_{n} and Y_{n} are the column vectors 

_{} 

for some constant Dt (called the stepsize of the simulation), and g and h are constant row vectors 

_{} 

Given the coefficients a_{k} and b_{k} of (1), our objective is to compute, as quickly and efficiently as possible, the components of h and g. 

It is well known [4] that if the recurrence is required to match the exact homogeneous response for y(t) when x(t) is constantly zero, then the components of h must be proportional to the elementary symmetric functions [8] (with alternating signs) of the quantities _{}, i = 1, 2, .., N, where a_{i} are the roots of the characteristic polynomial of the left side of (1). The usual procedure for determining this "rootmatched" recurrence is to solve the characteristic polynomial for the roots a_{i}, and then compute the elementary symmetric functions of these roots [5]. For high order relationships the most troublesome aspect of this approach is the determination of all the complex roots of the characteristic polynomials, especially when there are multiple roots [6]. The most common way of simplifying the calculations is to use a substitution method, e.g., bilinear substitution [7], but such methods do not reproduce the exact homogeneous response. The purpose of this note is to describe how, using Newton's Identities, we can proceed directly from the coefficients of (1) to the exact rootmatched recurrence coefficients without solving the characteristic polynomials. 

2.0 Determining the Root Matched Recurrence 

For any integer k, let w(k) denote the sum of the kth powers of the roots a_{1},a_{2},..,a_{N}. Clearly, w(0) = N. Using Newton's Identities [8,9] the values of w(k) for k = 1,2,... are given by 

_{} 

for k < N, and the recurrence 

_{} 

for all k ³ N. Now let E(k) denote the sum of the kth powers of the values of _{}, i = 1, 2, .., N. Expanding each of the exponential functions as a Taylor's series and combining terms by powers of Dt gives the following expression for E(k): 

_{} 

Letting s_{k} denote (1)^{k} times the kth elementary symmetric function of the quantities _{}, we can now compute s_{0} through s_{N} using the identities 

_{} 

The precision of these computations, including the convergence of (3), can be checked by verifying the equality 

_{} 

The general recurrence for the homogeneous equation in y(t) can now be expressed as: 

_{} 

where S_{A} is the row vector 

_{} 

Figure 1 presents an implementation of the preceding algorithm in BASIC. The inputs to this program are the order N, the stepsize Dt, and the coefficients a_{0}, a_{1 }.., a_{N}. The program computes the values of s_{0}, s_{1},..,s_{N}. The parameter "sumto" determines the number of terms to be evaluated in the summation (3). 

By entering the coefficients b_{k} in place of a_{k} into the preceding algorithm we can determine the recurrence for the homogeneous equation in x(t) 

_{} 

We can then combine the two sides of the recurrence by applying the appropriate scale factors, which are easily deduced from steadystate conditions. The result is 

_{} 

where SS_{A} and SS_{B }denote the sums of the components of S_{A} and S_{B} respectively. 

3.0 Equilateral Recurrences 

If b_{N }= 0 then the characteristic polynomial of the right side of (1) has fewer than N roots. In this case the recurrence (4) would involve fewer "past values" of x(t) than of y(t). Note that the coefficients of x_{n} in (4) were uniquely determined by the requirement that x(t) conform to the true homogeneous behavior when y(t) is constantly zero. However, in practice x(t) is an arbitrary forcing function so we are not really concerned with reproducing the homogeneous "response" of x(t). Instead, we are trying to reconstruct the function x(t) from the discrete samples. Thus it is often considered more appropriate to make use of as many past values of x(t) as is computationally convenient. This typically leads to recurrences that use the same number of past values of x(t) as of y(t), i.e., equilateral recurrences. The usual approach in signal processing is to augment the roots of the righthand characteristic polynomial of (1) with fictitious roots of the form ±pi, and then proceed as before [10]. The effect is to add roots equaling 1 to the characteristic polynomial of the recurrence, which essentially mimics bilinear substitution. 

The rationale for the "method of fictitious roots" just described is based on the frequency response characteristics. In the time domain it is often more appropriate to convert (1) into an equilateral equation by a simple change of variables. If we define the variable z(t) such that y = l(x+z) for some constant l, then (1) can be written 

_{} 

For any choice of l not equal to 0, (b_{N}/a_{N}), or (b_{0}/a_{0}), this equation is equilateral, so we can apply the method of (4) to give the recurrence 

_{} 

where g_{l} is based on the coefficients of the right hand side of (5). Reverting back to the y(t) variable gives 

_{} 

For large values of l the vector l(g_{l}_{ }+ h) is quite insensitive to changes in l, and it converges in the limit as l ® „. In practice, any l greater than about 10 times the largest b_{i} will give very nearly the limiting vector. This "changeofvariables method" enables us to determine equilateral recurrences for nonequilateral differential equations without resorting to fictitious roots. 

4.0 The Expected Value Method 

For comparison we briefly describe a method that explicitly treats x(t) as an independent variable and y(t) as a dependent variable. From this point of view it can be argued that the expected value of x(t) based on the N+1 most recent samples is given by the Nth degree polynomial that passes through those samples. Although the method is a straightforward application of polynomial curve fitting, no explicit formulation of the general Nthorder recurrence formula appears in the literature, so a brief description is given below. 

As in our previous recurrence formulas, we require that the homogeneous response of y(t) be matched exactly, so we use the algorithm shown in Figure 1 to determine S_{A}. Combining this with the particular solution based on the Nth order polynomial curve fit of x(t) leads directly to the total recurrence 

_{} 

where A, B, and T are square matrices such that the elements in the jth columns of the kth rows (where j,k range from 0 to N) are given by 

_{} 

This method gives, arguably, the best possible reconstruction of arbitrary forcing functions x(t), although the matrix inversions make it somewhat laborious when applied to high order equations. 

5.0 Examples 

To illustrate the basic algorithm for determining the rootmatched recurrence relation for a homogeneous equation, consider the differential equation 

_{} 

To simplify the recurrence expressions we will denote y(nDt) by y_{n}. Taking a time interval of 
Dt = 0.1, the algorithm described in Section 2 gives the following recurrence formula for y_{n} 

_{} 

where 
_{} 

(Notice that we've scaled the coefficients to make h_{5} equal to 1.) The accuracy of this recurrence can be illustrated by noting that (6) has the characteristic roots 1, 2+i, 2i, 1+i, and 1i, so one particular analytical solution has the simple form 

_{} 

Taking y_{0} through y_{4} from (8) as initial values, it is easily verified that the recurrence (7) numerically reproduces the exact analytical values of y_{5}, y_{6}, y_{7} ... etc., very precisely. After 35 iterations of the recurrence formula (carried out to 10 significant decimal digits) the value of y_{40} differs from the analytical value by only 0.0000013, this discrepancy being due entirely to accumulated roundoff error. In contrast, using the recurrence relation given by bilinear substitution, the error at y_{40} is 0.0013094. For larger stepsizes the error using bilinear substitution becomes progressively worse, whereas the rootmatched recurrence method is equally applicable to any stepsize. 

We chose a homogeneous equation with simple characteristic roots for the preceding example so that the recurrence values could be easily compared to the analytical solution. Also, the characteristic roots in the preceding example all have negative real parts, so the solutions are "stable". However, our interest is not limited to "stable" systems. For example, problems involving combustion and thermal "runaway" often lead to "unstable" systems. To illustrate the various methods discussed in this note as applied to one of these "unstable" systems, and to include "numerator dynamics" [1], consider the differential equation 

_{} 

Taking a time interval of Dt = 1, the relation between x and y can be simulated using a linear recurrence of the form (2), where the components of h and g are as listed in the table below (normalized for h_{5} = 1). 


All the methods except Bilinear Substitution made use of the BASIC program in Figure 1 to determine the components of h and to assist in determining the components of g. Interestingly, the ChangeofVariables Method gives components of g that closely match those given by the Expected Value Method, even though the former makes no use of a 5th order curve fit of x(t). In contrast, the methods of fictitious roots and bilinear substitution give very different g vectors, and therefore will clearly yield significantly different values of y(t) in response to certain forcing functions x(t). 

6.0 Conclusions 

The technique described in this paper uses Newton's Identities to efficiently determine the rootmatched recurrence relations for numerically simulating linear ordinary differential equations. The algorithm is simple to implement and can easily be incorporated into any numerical simulation program. The main benefit of the method is that it eliminates the need to determine the characteristic roots of the differential equation, which can be difficult for high order equations. This makes the method particularly useful in realtime applications that require frequent computation of recurrence relations corresponding to given continuous transfer functions. The recurrences given by this method (unlike bilinear substitution) are valid for any step size Dt, provided only that the sampling period is small enough to adequately resolve the forcing function x(t). Also, the algorithm requires no special handling of characteristic equations with repeated roots, so it is applicable to any equation of the form (1). We have also shown how the "ChangeofVariables" method can be used to generate robust equilateral recurrences without resorting to fictitious roots. 

References 
1. E. O. Doebelin, System Dynamics, Modeling and Response, 
Charles E Merrill Publishing Company, 1972. 
2. V. Raman, "Explicit Multistep Methods In System Simulation", IEEE International Symposium on Circuits and Systems, 1989, vol 3, p 17921795. 
3. S. Sallam, W. Ameen, "Numerical Solution of General nthOrder Differential Equations Via Splines", Applied Numerical Mathematics, Vol 6, 1989/90, pp 225238. 
4. A. L. Rabenstein, Elementary Differential Equations With Linear Algebra, Academic Press, 1970. 
5. J. M. Smith, Mathematical Modeling and Digital Simulation For Engineers and Scientists, 2nd Ed, John Wiley & Sons, 1987. 
6. T. Miyakodo, "Iterative methods for multiple zeros of a polynomial by clustering", Journal of Computational and Applied Mathematics, vol 28 (1989), p 315326. 
7. V. J. Bucek, Control Systems, Continuous and Discrete, Prentice Hall, 1989. 
8. A. Clark, Elements of Abstract Algebra, Dover Publications, 1984. 
9. I. N. Herstein, Topics In Algebra, 2nd Ed, John Wiley & Sons, 1975. 
10. G. F. Franklin, J. David Powell, M. L. Workman, Digital Control of Dynamic Systems, 2nd Ed, AddisonWesley, 1990. 
