Recursively computing values of a polynomial using difference equations

In the lecture Introduction to digital control by Peter Corke, he talks about the historical importance of difference equations for computing values of a polynomial. Let's look at this in some more detail.

A first order polynomial

Consider the polynomial $$ p(x) = 4x + 2. $$ The first difference is $$ \Delta p(x) = p(x) - p(x-h) = 4x + 2 - \big( 4(x-h) + 2 \big) = 4h, $$ and the second order difference is zero (as are all higher order differences): $$ \Delta^2 p(x) = \Delta p(x) - \Delta p(x-h) = 4h - 4h = 0. $$

Using the firs order difference, we can also write the second order difference $ \Delta p(x) - \Delta p(x-h) = \Delta^2 p(x) $ as $$ p(x) - p(x-h) - \Delta p(x-h) = \Delta^2p(x) $$ or $$ p(x) = p(x-h) + \Delta p(x-h) + \Delta^2 p(x)$$ which for the first order polynomial above becomes $$ p(x) = p(x-h) + \Delta p(x-h) = p(x-h) + 4h. $$


In [2]:
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
%matplotlib inline

def p1(x): return 4*x + 2 # Our first-order polynomial
# Compute values for x=[0,0.2, 0.4, ... 2] recursively using the difference equation
h = 0.2
x = h*np.arange(11) # Gives the array [0,0.2, 0.4, ... 2]
pd = np.zeros(11)
d1 = 4*h

# Need to compute the first value as the initial value for the difference equation,
pd[0] = p1(x[0])

for k in range(1,10): # Solve difference equation
    pd[k] = pd[k-1] + d1

plt.figure(figsize=(14,6))
plt.plot(x, p1(x), linewidth=2)
plt.plot(x, pd, 'ro')


Out[2]:
[<matplotlib.lines.Line2D at 0x7fc8fab5cac8>]

Second order polynomial

For a second order polynomial $$ p(x) = a_0x^2 + a_1x + a_2 $$ we have $$ p''(x) = 2a_0, $$ and the differences $$ \Delta p(x) = p(x) - p(x-h) = a_0x^2 + a_1x + a_2 - \big( a_0(x-h)^2 + a_1(x-h) + a_2 \big) = h(2a_0x + a_1) -a_0h^2, $$ $$ \Delta^2 p(x) = \Delta p(x) - \Delta p(x-h) = h(2a_0x+a_1) - a_0h^2 - \big( h(2a_0(x-h) + a_1) - a_0 h^2 \big) = h^22a_0 $$

Recall the difference equation using the second order difference $$ p(x) = p(x-h) + \Delta p(x-h) + \Delta^2 p(x)$$ We now get $$ p(x) = p(x-h) + \Delta p(x-h) + \Delta^2 p(x) = p(x-h) + \Delta p(x-h) + h^2 2 a_0,$$ or, using the definition of the first-order difference $\Delta p(x-h)$ $$ p(x) = 2p(x-h) - p(x-2h) + h^2 2 a_0,$$

Consider the second order polynomial $$ p(x) = 2x^2 - 3x + 2, $$ and compute values using the difference equation.


In [38]:
a0 = 2
a1 = -3
a2 = 2
def p2(x): return a0*x**2 + a1*x + a2 # Our second-order polynomial

# Compute values for x=[0,0.2, 0.4, ... 8] recursively using the difference equation
h = 0.2
x = h*np.arange(41) # Gives the array [0,0.2, 0.4, ... 2]
d1 = np.zeros(41) # The first differences
pd = np.zeros(41)
d2 = h**2*2*a0 # The constant, second difference

# Need to compute the first two values to get the initial values for the difference equation,
pd[0] = p2(x[0])
pd[1] = p2(x[1])

for k in range(2,41): # Solve difference equation
    pd[k] = 2*pd[k-1] - pd[k-2] + d2
    
plt.figure(figsize=(14,6))
plt.plot(x, p2(x), linewidth=2) # Evaluating the polynomial
plt.plot(x, pd, 'ro') # The solution using the difference equation


Out[38]:
[<matplotlib.lines.Line2D at 0x7f8e1161cb90>]

Exercise

What order would the difference equation be for computing valuse of a third-order polynomial? What is the difference equation?


In [ ]: