Discretization

Discretization of a differential equation¶

Discretization is a key step in simulating a differential equation on a computer. There are many ways to discretize a continuous function:

• foward Euler's method
• backward Euler's method
• bilinear transform
• impulse invariance
• and many more.

Here we'll consider Euler's method and a few variations on impulse invariance.
Consider the following differential equation:
$$\tau\dot{x}=-x+u$$
We need to approximate $x(t)$ and $\dot{x}$ in terms of $x[k]$, where $x[k]$ is the $k$th sample of $x(t)$. That is
$$x[k]=x(t_k)$$
where $t_k = k \Delta t$ for sample period $\Delta t$. We'd like to find an update equation that describes the $x$ over the discrete steps in time. That is, we'd like to find an equation of the form
$$ x[k] = f(x\text{ in the past}) + g(\text{input}) $$

Euler's method¶

Euler's method approximates the differential operator with a difference equation. This is equivalent to linearizing the differential equation about a point. Among the various versions ofthe Euler method, the most common are the forward difference and the backward difference approximations.

\begin{align} \dot{x}(t_k)& \approx \frac{x(t_{k+1})-x(t_{k})}{\Delta t} \quad \text{forward Euler} \\ \dot{x}(t_k)& \approx \frac{x(t_{k})-x(t_{k-1})}{\Delta t} \quad \text{backward Euler} \end{align}

or in discrete notation

\begin{align} \dot{x}(t_k)& \approx \frac{x[k+1]-x[k]}{\Delta t} \quad \text{forward Euler} \\ \dot{x}(t_k)& \approx \frac{x[k]-x[k-1]}{\Delta t} \quad \text{backward Euler} \end{align}

Discretization using the forward Euler's method¶

Substituting the forward Euler approximation into the differential equation under consideration, we have
\begin{align} \tau\frac{x[k+1]-x[k]}{\Delta t} &= -x[k] + u[k] \\ \frac{\tau}{\Delta t}x[k+1]-\frac{\tau}{\Delta t}x[k] &= -x[k] + u[k] \\ \frac{\tau}{\Delta t}x[k+1] &= \left(\frac{\tau}{\Delta t}-1\right)x[k] + u[k] \\ \frac{\tau}{\Delta t}x[k+1] &= \left(\frac{\tau-\Delta t}{\Delta t}\right)x[k] + u[k] \\ x[k+1] &= \left(\frac{\tau-\Delta t}{\tau}\right)x[k] + \frac{\Delta t}{\tau}u[k] \\ x[k+1] &= \alpha x[k] + (1-\alpha) u[k] \qquad \text{where }\alpha\triangleq\left(\frac{\tau-\Delta t}{\tau}\right) \end{align}

Discretization using the backward Euler's method¶

With the backward Euler approximation, we have
\begin{align} \tau\frac{x[k]-x[k-1]}{\Delta t} &= -x[k] + u[k] \\ \frac{\tau}{\Delta t}x[k]-\frac{\tau}{\Delta t}x[k-1] &= -x[k] + u[k] \\ \left(\frac{\tau}{\Delta t}+1\right)x[k] &= \frac{\tau}{\Delta t}x[k-1] + u[k] \\ \left(\frac{\tau+\Delta t}{\Delta t}\right)x[k] &= \frac{\tau}{\Delta t}x[k-1] + u[k] \\ x[k] &= \left(\frac{\tau}{\tau+\Delta t}\right)x[k-1] + \left(\frac{\Delta t}{\tau+\Delta t}\right)u[k] \\ x[k] &= \alpha x[k-1] + (1-\alpha) u[k] \qquad \text{where }\alpha\triangleq\left(\frac{\tau}{\tau+\Delta t}\right) \end{align}

Note that the update equations have the same form in both the forward and backward Euler's method. Both are a weighted sum of the previous state and the input.

Impulse invariance¶

The basic idea of impuse invariance is to discretize a continuous system by sampling the impulse response, finding the transfer function, and then finding the corresponding update equation. That is, for continuous impulse response $h(t)$, the discrete impulse response will be given by
$$h[n]=h(n \Delta t)\Delta t.$$
Note that this formula is not consistent when a causal continuous-time impulse response has a discontinuity at $t=0$ (See here), but for now we'll assume that the correction has a negligible effect (its effect is proportional to the step size) and continue with the uncorrected version for simplicity. We'll look at the corrected version later. For a continuous transfer function with $N$ poles
$$H(s)=\sum_{k=1}^N\frac{A_k}{s-s_k},$$
the discrete transfer function will be given by
$$H(z)=\Delta t\sum_{k=1}^N\frac{A_k}{1-e^{s_k\Delta t}z^{-1}}.$$
Our differential equation $\tau \dot{x}=-x+u$ has impulse response
$$h(t) = \frac{1}{\tau}e^{-t/\tau}u(t)$$
where $u(t)$ is the unit step function, which corresponds to the transfer function

\begin{align} \tau sX(s) &= -X(s)+U(s) \\ (\tau s+1)X(s) &= U(s) \\ \frac{X(s)}{U(s)} &= \frac{1}{\tau s+1} \end{align}

\begin{align} H(s) &= \frac{X(s)}{U(s)} \\ &= \frac{1}{\tau s+1} \\ &= \frac{1/\tau}{s+1/\tau} \end{align}

This means that the discrete transfer function will be
$$H(z)=\frac{\Delta t/\tau}{1-e^{-\Delta t/\tau}z^{-1}}.$$
From the transfer function, we can derive the update equation.

\begin{align} H(z) &= \frac{X(z)}{U(z)} \\ \frac{\Delta t/\tau}{1-e^{-\Delta t/\tau}z^{-1}} &= \frac{X(z)}{U(z)} \\ \left(1-e^{-\Delta t/\tau}z^{-1}\right)X(z) &= \frac{\Delta t}{\tau}U(z) \\ X(z) &= e^{-\Delta t/\tau}z^{-1}X(z) + \frac{\Delta t}{\tau}U(z) \\ x[k] &= e^{-\Delta t/\tau}x[k-1]+\frac{\Delta t}{\tau}u[k] \end{align}

See again that the difference equation takes on the same form with the weighted sum of the previous state with the input. In addition, if we use the Taylor series expansion of the exponential we can see the connection between the impulse invariance approximation and the foward Euler method.

\begin{align} e^{\Delta t/\tau} &= 1-\frac{\Delta t}{\tau}+ \frac{1}{2!}\left(\frac{\Delta t}{\tau}\right)^2- \frac{1}{3!}\left(\frac{\Delta t}{\tau}\right)^3\ldots \\ &\approx 1-\frac{\Delta t}{\tau} \\ &= \frac{\tau-\Delta t}{\tau} \end{align}

Substituting the Taylor series approximation, we get $$x[k]=\frac{\tau - \Delta t}{\tau}x[k-1]+\frac{\Delta t}{\tau}u[k]$$ which is the same as the forward Euler method approximation!

Impulse invariance, corrected¶

For a continuous transfer function
$$H(s)=\sum_{k=1}^N\frac{A_k}{s-s_k},$$
the corrected impulse invariance is given by
$$ h[n] = T\left(h(nT)-\frac{1}{2}h(0)\delta[n]\right) $$
which consequently means that
$$ H(z) = \Delta t\sum_{k=1}^{N}\frac{A_k}{1-e^{s_k\Delta t}z^{-1}} - \frac{\Delta t}{2}\sum_{k=1}^{N}A_k. $$
Recall that our differential equation has a transfer function
$$ H(s)=\frac{1/\tau}{s+1/\tau}, $$

so our discrete transfer function will be given by
\begin{align} H(z) &= \frac{\Delta t/\tau}{1-e^{-\Delta t/\tau}z^{-1}} - \frac{\Delta t/\tau}{2} \\ &= \frac{\Delta t}{\tau}\left(\frac{1}{1-e^{-\Delta t/\tau}z^{-1}} - \frac{1}{2}\right) \\ &= \frac{\Delta t}{2\tau}\left(\frac{2-\left(1-e^{-\Delta t/\tau}z^{-1}\right)}{1-e^{-\Delta t/\tau}z^{-1}}\right) \\ &= \frac{\Delta t}{2\tau}\left(\frac{1+e^{-\Delta t/\tau}z^{-1}}{1-e^{-\Delta t/\tau}z^{-1}}\right) \\ \end{align}

Our update equation will then follow as
\begin{align} \frac{X(z)}{U(z)} &= H(z) \\ &= \frac{\Delta t}{2\tau}\left(\frac{1+e^{-\Delta t/\tau}z^{-1}}{1-e^{-\Delta t/\tau}z^{-1}}\right) \\ X(z)\left(1-e^{-\Delta t/\tau}z^{-1}\right) &= \frac{\Delta t}{2\tau}\left(1+e^{-\Delta t/\tau}z^{-1}\right)U(z) \\ &\Updownarrow \\ x[k]-e^{-\Delta t/\tau}x[k-1] &= \frac{\Delta t}{2\tau}u[k]+\frac{\Delta t}{2\tau}e^{-\Delta t/\tau}u[k-1] \\ x[k] &= e^{-\Delta t/\tau}x[k-1] + \frac{\Delta t}{2\tau}u[k]+\frac{\Delta t}{2\tau}e^{-\Delta t/\tau}u[k-1]. \\ \end{align}

Conservation of area¶

Sometimes we care about preserving the area of the continuous function being discretized. The area, as the integral of the function over time, often has the physical interpretation of energy. Here we look at whether the various discretization techniques preserve the area of an impulse in their impulse responses. An impulse is defined as having unit area. We typically discretize the impulse as a $1$-time step pulse of height $\frac{1}{\Delta t}$ at $u[0]$ so that with a time step of $\Delta t$, the area is given by $\frac{1}{\Delta t}\Delta t=1$.

Forward Euler's method area¶

Recall that the update equation for the forward Euler's method is
$$ x[k+1] = \left(\frac{\tau-\Delta t}{\tau}\right)x[k] + \frac{\Delta t}{\tau}u[k]. $$
Given an impulse of height $\frac{1}{\Delta t}$, the impulse response will then have area
\begin{align} A &= \Delta t\left(\frac{1}{\tau} + \frac{1}{\tau}\left(\frac{\tau-\Delta t}{\tau}\right) + \frac{1}{\tau}\left(\frac{\tau-\Delta t}{\tau}\right)^2 + \ldots\right) \\ &= \frac{\Delta t}{\tau}\sum_{k=0}^{\infty}\left(\frac{\tau-\Delta t}{\tau}\right)^k \\ &= \frac{\Delta t}{\tau}\frac{1}{1-\frac{\tau-\Delta t}{\tau}} \\ &= \frac{\Delta t}{\tau}\frac{\tau}{\tau-(\tau-\Delta t)} \\ &= \frac{\Delta t}{\Delta t} \\ &= 1. \end{align}
The unit area is preserved!

Backward Euler's method area¶

Recall that the update equation for the forward Euler's method is
$$ x[k] = \left(\frac{\tau}{\tau+\Delta t}\right)x[k-1] + \left(\frac{\Delta t}{\tau+\Delta t}\right)u[k] $$
Given an impulse of height $\frac{1}{\Delta t}$, the impulse response will then have area
\begin{align} A &= \Delta t\left(\frac{1}{\tau+\Delta t} + \frac{1}{\tau+\Delta t}\frac{\tau}{\tau+\Delta t} + \frac{1}{\tau+\Delta t}\left(\frac{\tau}{\tau+\Delta t}\right)^2 + \ldots \right) \\ &= \frac{\Delta t}{\tau+\Delta t}\sum_{k=0}^{\infty}\left(\frac{\tau}{\tau+\Delta t}\right)^k \\ &= \frac{\Delta t}{\tau+\Delta t}\frac{1}{1-\frac{\tau}{\tau+\Delta t}} \\ &= \frac{\Delta t}{\tau+\Delta t}\frac{\tau+\Delta t}{\tau+\Delta t-\tau} \\ &= \frac{\Delta t}{\Delta t} \\ &= 1. \end{align}
The unit area is preserved!

Impulse invariance area¶

Recall that the update equation for impulse invariance is
$$ x[k] = e^{-\Delta t/\tau}x[k-1]+\frac{\Delta t}{\tau}u[k] $$
Given an impulse of height $\frac{1}{\Delta t}$, the impulse response will then have area
\begin{align} A &= \Delta t\left(\frac{1}{\tau} + \frac{1}{\tau}e^{-\Delta t/\tau} + \frac{1}{\tau}e^{-2\Delta t/\tau} + \ldots \right) \\ &= \frac{\Delta t}{\tau}\sum_{k=0}^{\infty}\left(e^{-\Delta t/\tau}\right)^k \\ &= \frac{\Delta t}{\tau}\frac{1}{1-e^{-\Delta t/\tau}} \\ &\neq 1. \end{align}
The unit area is not preserved...In fact, taking the Taylor series expansion, we see that the area is increased as
\begin{align} A &= \frac{\Delta t}{\tau}\frac{1}{1-e^{-\Delta t/\tau}} \\ &= \frac{\Delta t}{\tau}\frac{1}{1-\left(1-\frac{\Delta t}{\tau}+ \frac{\Delta t^2}{2\tau^2}- \frac{\Delta t^3}{6\tau^3}+\ldots\right)} \\ &= \frac{\Delta t}{\tau}\frac{1}{\frac{\Delta t}{\tau}- \frac{\Delta t^2}{2\tau^2}+ \frac{\Delta t^3}{6\tau^3}-\ldots} \\ &> 1 \end{align}
because
$$ \frac{\Delta t}{\tau}-\frac{\Delta t^2}{2\tau^2}+ \frac{\Delta t^3}{6\tau^3}-\ldots < \frac{\Delta t}{\tau} $$

Impulse invariance, corrected, area¶

Recall that the update equation for corrected impulse invariance is
$$ x[k] = e^{-\Delta t/\tau}x[k-1] + \frac{\Delta t}{2\tau}u[k]+\frac{\Delta t}{2\tau}e^{-\Delta t/\tau}u[k-1] $$

Given an impulse of height $\frac{1}{\Delta t}$, the impulse response will then have area

\begin{align} A &= \Delta t\left(\overbrace{\frac{1}{2\tau}}^{k=0} + \overbrace{\frac{1}{2\tau}e^{-\Delta t/\tau} + \frac{1}{2\tau}e^{-\Delta t/\tau}}^{k=1} + \overbrace{\left(\frac{1}{2\tau}e^{-\Delta t/\tau} + \frac{1}{2\tau}e^{-\Delta t/\tau}\right)e^{-\Delta t/\tau}}^{k=2} + \ldots \right) \\ &= \Delta t\left(\frac{1}{2\tau} + \frac{1}{\tau}e^{-\Delta t/\tau} + \frac{1}{\tau}e^{-2\Delta t/\tau} + \ldots \right) \\ &= \frac{\Delta t}{2\tau}+\frac{\Delta t}{\tau}\sum_{k=1}^{\infty}\left(e^{-\Delta t/\tau}\right)^k \\ &= \frac{\Delta t}{\tau}\left(\frac{1}{2}+\sum_{k=1}^{\infty}\left(e^{-\Delta t/\tau}\right)^k\right) \\ &= \frac{\Delta t}{\tau}\left(\frac{1}{2}+\frac{e^{-\Delta t/\tau}}{1-e^{-\Delta t/\tau}}\right) \\ &= \frac{\Delta t}{\tau}\left(\frac{1-e^{-\Delta t/\tau}+2e^{-\Delta t/\tau}}{2\left(1-e^{-\Delta t/\tau}\right)}\right) \\ &= \frac{\Delta t}{2\tau}\frac{1+e^{-\Delta t/\tau}}{1-e^{-\Delta t/\tau}} \\ &\neq 1 \end{align}
Again, the unit area is not preserved...

Area preserving impulse response discretization¶

There is a modification to the impulse response update equation that will preserve area. If we let
$$ x[k] = e^{-\Delta t/\tau}x[k-1]+\left(1-e^{-\Delta t/\tau}\right)u[k] $$
then the unit area of an impulse will be preserved. To prove this,
\begin{align} A &= \Delta t\left(\frac{1-e^{-\Delta t/\tau}}{\Delta t} + \frac{1-e^{-\Delta t/\tau}}{\Delta t}e^{-\Delta t/\tau} + \frac{1-e^{-\Delta t/\tau}}{\Delta t}e^{-2\Delta t/\tau} + \ldots \right) \\ &= \left(1-e^{-\Delta t/\tau}\right)\sum_{k=0}^{\infty}\left(e^{-\Delta t/\tau}\right)^k \\ &= \left(1-e^{-\Delta t/\tau}\right)\frac{1}{1-e^{-\Delta t/\tau}} \\ &= 1. \end{align}

Back to top

Noisy Dynamical Systems

Today we're going to play with a nonlinear dynamical system and noise. Code for singing along is posted on github. Our system is described, in general, by

\[\dot{x} = f(x),\]

which says that the rate of change of \(x\) is just some arbitrary function \(f(\cdot)\) of \(x\). In block diagram format, the system looks like this:

System block diagram. \(1/s\) is an integrator.  The derivative of \(x\) is its input, and \(x\) is its output.

Specifically, let's let our nonlinear function take the form 

\[f(x) = e^{-x}-1.\]

To start analyzing this system, we'll plot the phase portrait, which shows us how \(\dot{x}\) depends on \(x\).

The phase portrait for our nonlinear system tells us how the system changes. We'll be running this system from several different initial conditions indicated with the black to red  dots. Their positions along the x-axis are their initial condition, and their positions along the y-axis are their initial derivatives. 

From the phase portrait, we can tell that the system's state, or value of \(x\), increases when the state is negative and decreases when the state is positive. Therefore, we can predict that the system will stabilize at 0.

Do we see this stable point?

System state trajectory over time for initial states scattered between -1.5 and 1.5 .

Yes! Great!

What if we add some noise? 

\[\dot{x} = f(x) + u \]

\[u \text{ is } \mathcal{N}(0, \sigma)\]

That is, we've added Gaussian white noise of variance \(\sigma^{2}\) to our system.

Our system block diagram now looks like this: 

Block diagram of system with noise added to input of integrator.

Do we still see the system go to same stable point?

System simulation with noise added to the input of the integrator.  Noise is  Gaussian white noise with variance indicated by \(\sigma\) above each subplot.

It looks like the system goes to the same stable point at low levels of noise, but system get's pretty swamped with noise at higher variances, which makes it hard to tell if the system settles around the predicted stable point. Running the simulation many times and averaging across trials helps us see the long term behavior more clearly:

Each line is the average of 25 trials.

That's better. It looks like the system settles, on average, at the predicted state for low noise levels, but at higher levels of noise, the system seems to be averaging at a higher state value. If we do the same averaging over trials but with much larger noise, we see this effect more clearly:

The system no longer averages out anywhere close to the predicted stable point without noise. The system doesn't even look like it reaches its average stable value by the end of the simulation for very high levels of noise. Why does the stable point shift?

Does the average stable point change because we are now inputting a low pass filtered version of the noise into the nonlinearity? I'm not sure of an intuitive way of understanding this yet.

Now let's tweak the system a little:

\[\dot{x} = f(x+u)\]

Now, noise is passed through the nonlinearity before integration.

Show system diagram.

Block diagram of system with noise added to the input to the nonlinearity.

Running the simulation for this system yields:

System simulation with noise added to the input to the nonlinearity.  Each line is a single trial, not averaged. 

Look carefully at the value that the system settles at for each level of noise. The system is less noisy, but fixed point position is far more sensitive to the level of noise!

Think about it this way. If you are the system, and you are at a state x, what do you see? Because of the noise you see the everything around you. So if we take the noise distribution and average each point of the phase portrait with our noise (i.e. convolve our system function \(f(\cdot)\) with the pdf (probability density function) of our noise), we get a better sense of the fixed points of our system:

Phase portrait of the the system after convolving the system function, \(f\), with the noise distribution pdf.

Another way of analyzing the effect of noise is via the Taylor series expansion of the \[f(\cdot)\]. With noise, \(\dot{x}=f(x+u)\), and

\[f(x+u)=f(x)+f'(x)u+\frac{1}{2}f''(x)u^{2}+\frac{1}{3!}f^{(3)}(x)u^3+...\]

Taking the expected value,

\[E[f(x+u)]=f(x)+f'(x)E[u]+\frac{1}{2}f''(x)E[u^{2}]+\frac{1}{3!}f^{(3)}(x)E[u^3]+...\]

Since we're assuming \(u\) is symmetric about 0, the expected value of the odd powers of \(u\) are also 0,  so

\[E[f(x+u)]=f(x)+\frac{1}{2}f''(x)E[u^{2}]+\frac{1}{4!}f^{(4)}(x)E[u^4]+...\].

Before moving forward, we can see here why the phase portrait changes with noise. When noise entering the nonlinearity, \(\dot{x}\) is no longer just \(f(\cdot)\), but \(f(\cdot)\) summed with its derivatives weighted by the higher order moments of \(u\). In our example case,

\[f(x)=e^{-x}-1\]

\[f'(x)=-e^{-x}\]
\[f''(x)=e^{-x}\]

\[f^{(3)}(x)=-e^{-x}\]

\[f^{(4)}(x)=e^{-x}\]

\[\vdots\]

Notice a pattern? All of the odd ordered derivatives are the same, and all of the even ordered derivatives are the same! So for our example,

\[E[f(x+u)]=e^{-x}-1+\frac{1}{2}e^{-x}E[u^{2}]+\frac{1}{4!}e^{-x}E[u^4]+...\]

\[E[f(x+u)]=e^{-x}-1+e^{-x}\left(\frac{\sigma^{2}}{2}+\frac{3\sigma^{4}}{4!}+...\right)\]

This is neat! We can analytically see that noise changes the phase portrait.


What's different between the two noise cases?
Both noise cases change the dynamics of the system, but the system is more sensitive to noise input to the nonlinearity.

Keep in mind that in reality we can really see both kinds of noise:.

Block diagram of system with noise added to both 

If you can control its distribution, noise can even become another knob to tune your system!