Recurrence Relations for Ordinary Differential Equations

Summary: A technique is described for using Newton's identities to compute the coefficients of a multi-step recurrence formulae for simulating the solutions of ordinary linear differential equations of high order. This approach gives the exact root-matched 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

(1)

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 control-system compensation, but in those applications the frequency response is usually emphasized over the exact time-domain 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 multi-step linear recurrence relation of the form

h Y_n = g X_n(2)

where Xn and Yn are the column vectors

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

h = [ h₀ h₁ ... h_N ]g = [ g₀ g₁ ... g_N ]

Given the coefficients ak and bk 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 ai are the roots of the characteristic polynomial of the left side of (1). The usual procedure for determining this "root-matched" recurrence is to solve the characteristic polynomial for the roots ai, 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 root-matched 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₁,a₂,..,a_N. Clearly, w(0) = N. Using Newton's Identities [8,9] the values of w(k) for k = 1,2,... are given by

a_N w(1) + 1a_N-1 = 0

a_N w(2) + a_N-1 w(1) + 2a_N-2 = 0

a_N w(3) + a_N-1 w(2) + a_N-2 w(1) + 3a_N-3 = 0

for k < N, and the recurrence

a_N w(k) + a_N-1 w(k-1) + ... + a₀ w(k-N) = 0

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):

(3)

Letting s_k denote (-1)k times the kth elementary symmetric function of the quantities , we can now compute s₀ through s_N using the identities

E(1) + 1s₁ = 0

E(2) + s₁ E(1) + 2s₂ = 0

E(3) + s₁ E(2) + s₂ E(1) + 3s₃ = 0

etc.

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:

S_AY_n = 0

where SA is the row vector

S_A = [ s_N s_N-1 ... s₀ ]

Figure 1 presents an implementation of the preceding algorithm in BASIC. The inputs to this program are the order N, the step-size Dt, and the coefficients a₀,a₁ ..,a_N. The program computes the values of s₀, s₁,..,s_N. The parameter "sumto" determines the number of terms to be evaluated in the summation (3).

By entering the coefficients bk in place of ak into the preceding algorithm we can determine the recurrence for the homogeneous equation in x(t)

S_BX_n = 0

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

(4)

where S S_A and S S_B denote the sums of the components of SA and SB 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 right-hand characteristic polynomial of (1) with fictitious roots of the form �p i, and then proceed as before [10]. The effect is to add roots equalling -1 to the characteristic polynomial of the recurrence, which essentially mimics bi-linear 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

(5)

For any choice of l not equal to 0, (b_N/a_N), or (b₀/a₀), this equation is equilateral, so we can apply the method of (4) to give the recurrence

h Z_n = gl X_n

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

h Y_n = l (gl + h) X_n

For large values of l the vector l (gl + 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 "change-of-variables method" enables us to determine equilateral recurrences for non-equilateral 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 straight-forward application of polynomial curve fitting, no explicit formulation of the general Nth-order 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 SA. Combining this with the particular solution based on the Nth order polynomial curve fit of x(t) leads directly to the total recurrence

S_AY_n = S_ATA^-1BT^-1X_n

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 root-matched recurrence relation for a homogeneous equation, consider the differential equation

(6)

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

y_n = -h₄ y_n-1 - h₃ y_n-2 - h₂ y_n-3 - h₁ y_n-4 - h₀ y_n-5 (7)

where

h₀ = -0.4965853038

h₁ = 2.8479524294

h₂ = -6.5428490609

h₃ = 7.5263050469

h₄ = -4.3347524368

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

y_n = -8e^-nDt + (-6e^-2nDt + 14e^-nDt) cos(n Dt)(8)

Taking y₀ through y₄ from (8) as initial values, it is easily verified that the recurrence (7) numerically reproduces the exact analytical values of y₅, y₆, y₇ ... etc., very precisely. After 35 iterations of the recurrence formula (carried out to 10 significant decimal digits) the value of y₄₀ differs from the analytical value by only 0.0000013, this discrepancy being due entirely to accumulated round-off error. In contrast, using the recurrence relation given by bilinear substitution, the error at y₄₀ is 0.0013094. For larger step-sizes the error using bilinear substitution becomes progressively worse, whereas the root-matched recurrence method is equally applicable to any step-size.

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 "run-away" 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₅ = 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 Change-of-Variables 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 root-matched recurrence relations for numerically simulating linear ODEs. 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 real-time applications that require frequent computation of recurrence relations corresponding to given continuous transfer functions. The recurrences given by this method (unlike bi-linear 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 "Change-of-Variables" 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 1792-1795.

3.S. Sallam, W. Ameen, "Numerical Solution of General nth-Order Differential Equations Via Splines", Applied Numerical Mathematics, Vol 6, 1989/90, pp 225-238.

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 315-326.

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, Addison-Wesley, 1990.