Root-Matched Recurrences For DiffEQs

The "root-matched" recurrence formula corresponding to a given ordinary 
differential equation can be determined without explicitly solving for the 
roots of the characteristic polynomial.  This enables us to find the exact 
root-matched recurrence formula 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 we 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 we
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 difficult, 
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.