Root-Matched Recurrences For DiffEQs

Here is a method for determining the "root-matched" recurrence 
corresponding to a given ordinary differential equation without 
having to explicitly solve for the roots of the characteristic 
polynomial.  More generally, this method enables us to find the 
exact root-matched recurrence corresponding to any continuous 
linear system without having to solve for the eigen values of 
the system.  The following description covers the just simple 
homogeneous case.
 
Suppose you have an ordinary linear differential equation with constant
coefficients:

  a5 y''''' + a4 y'''' + a3 y''' + a2 y'' + a1 y' + a0 y  =  0      (1)

where the (') denotes derivatives with respect to t.  The response of
y(t) can be reproduced discretely by a recurrence relation of the form

c5 y(5dt) + c4 y(4dt) + c3 y(3dt) + c2 y(2dt) + c1 y(dt) + c0 y(0) = 0

where "dt" is the time increment of the simulation.  The question is,
given the coefficients ai of the differential equation, how can you
compute the coefficients ci of the corresponding recurrence relation.

The coefficients ci must be proportional to the elementary symmetric
functions (with alternating signs) of the quantities e^(ri dt), where 
ri (i=1,2,..,5) are the roots of the characteristic polynomial of (1).  
Thus, a straightforward method of computing the ci is to solve the 
characteristic polynomial of (1) for the roots ri, and then compute 
the elementary symmetric functions of these roots.  However, determining 
all the roots of a high order polynomial can be computationally tricky, 
particularly for embedded code, where the use of iterative processes 
can cause problems.

Fortunately, there is a way to compute the ci without having to solve 
for the characteristic roots of the equation.  The technique makes use 
of "Newton's Identities" and the Taylor's series expansion of the 
exponential function.

Assume our basic ODE is of order N.  For any integer k, let w(k) 
denote the sum of the kth powers of the roots r1,r2,...,rN.  
Clearly, w(0)=N.  Using Newton's Identities, the values of w(k) 
for k=1,2,... are given by

      w(1)a(N) + 1a(N-1)  =  0

      w(2)a(N) + w(1)a(N-1) + 2a(N-2)  =  0

      w(3)a(N) + w(2)a(N-1) + w(1)a(N-2) + 3a(N-3)  =  0

                            etc

for k=N.  Now let E(k) denote the sum of the kth powers of the
values e^(ri dt), i=1,2,..,N.  Expanding each exponential function as 
a Taylor series and combining terms by powers of dt gives the 
following expression for E(k):

      E(k) = SUM from j=0 to inf of { ((k dt)^j) w(j) / j! }

for k=1,2,..,N.  The coefficients of the recurrence relation can then 
be computed using the identities

                                              c(N)  =  1

                                   E(1) + 1 c(N-1)  =  0

                     E(2) + E(1) c(N-1) + 2 c(N-2)  =  0

       E(3) + E(2) c(N-1) + E(1) c(N-2) + 3 c(N-3)  =  0

This gives the exact root-matched recurrence coefficients without 
ever having to determine the roots.  It also has the advantage that 
no special provisions are needed to cover the case of repeated roots.