Trapezoidal Rule for Numerical Integration

The velocity $v(t)$ can be written as

$$\begin{eqnarray*} v(t)=v(0)+\left(\int_{0}^{t}\dot v(\tau)d\tau\right) \end{eqnarray*}$$

In particular,

$$\begin{eqnarray*} v(nT) % &=& v(0) + \int_{0}^{nT}\vt(\tau)d\tau \\ [10pt] &=& v(0) + \int_{0}^{(n-1)T}\dot v(\tau)d\tau + \int_{(n-1)T}^{nT}\dot v(\tau)d\tau \\ [10pt] &=& v[(n-1)T]+\int_{(n-1)T}^{nT}\dot v(\tau)d\tau \\ [10pt] &\approx& v[(n-1)T] +T\frac{\dot v[(n-1)T] + \dot v(nT)}{2} \end{eqnarray*}$$

• This approximation replaces a one-sample integral by the area under the trapezoid having vertices $(n-1,0), (n-1,\dot v_{n-1}),(n,0), (n,\dot v_n)$
• In other words, $\dot v(t)$ is approximated by a straight line between time $n-1$ and $n$
• This is a first-order approximation of $\dot v(t)$ in contrast to the zero-order approximation used by forward and backward Euler schemes
• We will see that the commonly used bilinear transform is equivalent
• Model is exact if driving force is piecewise linear, having a constant slope over each sampling interval
• (Backward Euler is similarly exact for a piecewise-constant driving force)

Bilinear Transform as Compensated BE/FE

In Newton's law $f=m\dot v$ , look at the Backward Euler (BE) approximation of the time-derivative:

$$f(t) \;=\;m\,\dot v\,\;\approx\;\,m\, \frac{v(t) - v(t-T)}{T}$$

We see there is a 1/2 sample delay in the first-order difference on the right. This misaligns the force $f(t)$ and subsequent velocity by half a sample. A very simple delay compensation is to use a two-point average on the left:

$$\frac{f(n)+f(n-1)}{2} \,\;\approx\;\,m\, \frac{v(n) - v(n-1)}{T}$$

The extra attenuation at high frequencies due to the two-point average actually helps. Taking the $z$ transform:

$$\frac{1+z^{-1}}{2}F(z) \,\;\approx\;\,m\, \frac{1-z^{-1}}{T} V(z)$$

or

$$F(z) \,\;\approx\;\,m\, \left(\frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}}\right) V(z)$$

which is the bilinear transform of $F(s) = ms\,V(s)$ :

$$\zbox{s \mapsto \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}}$$

Frequency Warping is the Only Error

We have

$$F(z) \,\;\approx\;\,m\, \left(\frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}}\right) V(z)$$

using the bilinear transform (trapezoidal integration in the time domain)

Let's look along the unit circle in the $z$ plane:

$$\frac{F(e^{j\omega T})}{V(e^{j\omega T})} \,\;\approx\;\,m\, \left(\frac{2}{T} \frac{1-e^{-j\omega T}}{1+e^{-j\omega T}}\right) \;=\;m\, j \left(\frac{2}{T} \tan\left(\frac{\omega T}{2}\right)\right)$$

Since the exact formula is $F(e^{j\omega T})/V(e^{j\omega T}) = m\, j\omega$ , we can push all of the error into a frequency warping:

$$\omega_d \mathrel{\stackrel{\mathrm{\Delta}}{=}}\frac{2}{T} \tan\left(\frac{\omega_a T}{2}\right)$$

• Frequency-warping is the only error over the unit circle when using the bilinear transform
• What started out as different gain errors on the left and right became the correct gains at warped frequency locations
• Frequency-warping implications should also be considered in the time domain