Convolution

In this notebook, we show that infinetly many impulses can approximate an arbitrary input, and thus justify the convolution integral.

In [1]:
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

from matplotlib import animation, rc
from IPython.display import HTML

import sympy
sympy.init_printing()

The convolution integral can be viewed as the summation of the impulse responses of a linear time invariant system.

$y(t) \approx \sum\limits_{k=0}^{t/\Delta} g(t - k \Delta ) u(k \Delta) \Delta$

As $\Delta \rightarrow 0$, this approaches the integral:

$y(t) = \int_0^t g(t-\sigma) u(\sigma) d\sigma$

where: $\sigma = k \Delta$, and $d\sigma = \Delta$.

Convolution

$y(t) = (g \otimes u)(t) \equiv \int_0^t g(t-\sigma) u(\sigma) d\sigma$

Example

The Differential Equation

$\dot{y}(t) + y(t) = u(t)$

Laplace Transform

$sY(s) + Y(s) = U(s)$

$(s + 1)Y(s) = U(s)$

$G(s) = Y(s)/U(s) = 1/(s + 1)$

Impulse Response Function

$g(t) = \mathcal{L}^{-1}[G(s)] = e^{-t}$

Input Function

$u(t) = \sin(t)$

Input Laplace Transform

$U(s) = 1/(s^2 +1)$

Output

$Y(s) = 1/(s^2+1) * 1/(s+1)$

$ \mathcal{L}^{-1}[Y(s)] = \frac{1}{2} \left(e^{-t} + sin(t) - cos(t)\right)$

In [2]:
def g(t):
    """The impulse response of the differential equation."""
    return np.exp(-t)

def u(t):
    return np.sin(t)

def delayed_impulse_response(g, t, t0):
    return g(t-t0)*(t>t0)

def ydot(y, u):
    """The differential equation."""
    return -y + u
In [3]:
def sim(dt_impulse):
    tf = 10
    dt = 0.001

    t = np.arange(0, tf + dt, dt)
    y = t*0
    hits = t*0
    y_analytic = np.exp(-t)/2 + (np.sin(t) - np.cos(t))/2
    y_euler = t*0
    yi_euler = 0
    t_last_pulse = 0
    u_true = u(t)

    for i, t0 in enumerate(t):
        if (t0 - t_last_pulse) >= dt_impulse:
            t_last_pulse = t0
            A = u_true[i]*dt_impulse
            hits[i] = A
            y += A*delayed_impulse_response(g, t, t0)
        yi_euler += dt*ydot(yi_euler, u_true[i])
        y_euler[i] = yi_euler
    
    return {
        't': t,
        'hits': hits,
        'u_true': u_true,
        'y': y,
        'y_analytic': y_analytic,
        'y_euler': y_euler
    }
In [4]:
def do_plot():
    
    data = sim(1)
    
    plt.figure()

    plt.subplot(211)
    plt.title('input')
    plt.plot(data['t'], data['hits'], label='impulse approx.')
    plt.plot(data['t'], data['u_true'], label='input')
    plt.legend()

    plt.subplot(212)
    plt.title('output')
    plt.plot(data['t'], data['y'], label='impulse response approx')
    plt.plot(data['t'], data['y_analytic'], label='analytic')
    plt.plot(data['t'], data['y_euler'], label='euler')
    plt.legend()
    
do_plot()
In [5]:
def make_movie():
    data = sim(1)

    fig = plt.figure(figsize=(10, 10))

    ax = plt.subplot(211)
    plt.title('input')
    h_hits, = plt.plot(data['t'], data['hits'], label='impulse approx.')
    plt.plot(data['t'], data['u_true'], label='input')
    plt.legend()

    ax = plt.subplot(212)
    plt.title('output')
    h_y, = plt.plot(data['t'], data['y'], label='impulse response approx')
    plt.plot(data['t'], data['y_analytic'], label='analytic')
    plt.plot(data['t'], data['y_euler'], label='euler')
    plt.legend()

    # initialization function: plot the background of each frame
    def init():
        h_hits.set_data([], [])
        h_y.set_data([], [])
        return (h_hits, h_y,)

    # animation function. This is called sequentially
    def animate(i):
        dt_impulse = 1 - 0.01*i
        if dt_impulse < 0.01:
            dt_impulse = 0.01
        data = sim(dt_impulse)
        h_hits.set_data(data['t'], data['hits'])
        h_y.set_data(data['t'], data['y'])
        return (h_hits, h_y,)

    # call the animator. blit=True means only re-draw the parts that have changed.
    anim = animation.FuncAnimation(fig, animate, init_func=init,
                                   frames=100, interval=100, blit=True)
    plt.close()
    return anim

anim = make_movie()
anim.save('convolution.mp4')
HTML(anim.to_html5_video())
Out[5]:
In [6]:
def find_inv_laplace():
    t = sympy.symbols('t', real=True)
    s= sympy.symbols('s', complex=True)
    return sympy.inverse_laplace_transform(1/(s+1)*1/(s**2+1), s, t)

find_inv_laplace()
Out[6]:
$\displaystyle \frac{\sin{\left(t \right)} \theta\left(t\right)}{2} - \frac{\cos{\left(t \right)} \theta\left(t\right)}{2} + \frac{e^{- t} \theta\left(t\right)}{2}$
In [7]:
def find_laplace():
    t = sympy.symbols('t', real=True)
    s= sympy.symbols('s', complex=True)
    return sympy.laplace_transform(sympy.sin(t), t, s)

find_laplace()
Out[7]:
$\displaystyle \left( \frac{1}{s^{2} + 1}, \ 0, \ \text{True}\right)$