Appendix B:  Nth Order Recurrence Formulas

 

In Sections 2 and 3 of this document the optimum recurrence formulas for first and second order response were derived based on the assumption that the input function x(t) was the first or second-degree polynomial (respectively) that passes through the given discrete values of x(t).  In this appendix we generalize this method to Nth order response.

 

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

 

 

where the coefficients ak and bk are constants.  We will denote by  r1, r2,..., rN  the roots of the characteristic polynomial of the left hand side of this differential equation.

 

For some instant t0 we are given the discrete values x(t0), x(t0+Dt), x(t0+2Dt), ..., x(t0+NDt) and  y(t0), y(t0+Dt), y(t0+2Dt), ..., y(t0+(N-1)Dt) , where Dt is the execution interval of the simulation.  For brevity we will denote y(t0+kDt) by yk and x(t0+kDt) by xk.  Our objective is to determine the value of yN based on the assumption that the function x(t) over the time interval from t0 to t0+NDt is identical to the Nth degree polynomial that passes through the given discrete values of x(t).

 

Define the N1 column vector Y and the (N+1)1 column vector X as follows:

 

                

 

Then define the square (N+1) (N+1) matrices A, B, and T such that the elements in the jth columns of the ith rows are as follows

 

 

where i and j range from 0 to N.  Then define the 1 (N+1) and 1 N row vectors

 

 

where sk denotes (-1)k times the kth elementary symmetric function of the quantities , i = 1,2,..,N.  For example, if N = 3 we have

 

 

In terms of the matrices defined above, the value of yN can now be expressed as

 

 

Letting gi denote the ith element of the row vector given by STA-1BT-1, we have the explicit recurrence formula

 

 

The values of sk are always purely real for real coefficients ak, even when the roots r1,r2,.., rN are complex.  Also, for reasons discussed in Section 2.2.5, the values of sk are the exact y-term coefficients, regardless of how the x(t) function is interpolated.  Only the gk coefficients would be affected if we changed our assumption as to the form of x(t) (that is, if we identified x(t) with a function different from the Nth order polynomial fit through the discrete values).

 

One way of computing sk is to solve the characteristic equation

 

 

for all the roots rj, j=1 to N, and then directly compute the symmetric functions of the exponentials.  However, finding all the (complex) roots of an Nth degree polynomial can be computationally laborious and involves numerous special cases.  Another approach, one that avoids having to explicitly determine all the roots of the characteristic polynomial, is to use "Newton's sums".  For any non-negative integer k, let w(k) denote the sum of the kth powers of the roots r1,r2,..,rN.  Clearly w(0) = N.  The values of w(k) for k = 1, 2, ... can be found using the identities

 

 

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.  The Taylor series expansion gives

 

 

The values of  sk  can now be computed using the identities

 

 

Return to Table of Contents