User:Johnsy/Advanced Modelling in Biology/Discrete Systems

Discrete Systems

Discrete systems are systems with non-continuous outputs with the result for each time step being determined by the previous time step. The equivalent to differential equations (continuous systems) in discrete system are known as difference equations such as the one shown below. Since the difference equation tells us about the next time step, it is known as a 1st order difference equation.

[math]\displaystyle{ x_{t+1} = rx_t \, }[/math]

This is the discrete equivalent of the differential equation:

[math]\displaystyle{ \frac{dx}{dt} = rx }[/math]

How do we solve these difference equations?

1) By guessing (ansatz - German for "guess")

Let us assume that the solution is in the form

[math]\displaystyle{ x_t = A \beta ^t \, }[/math]

Then substituting into the difference equation, we get

[math]\displaystyle{ A \beta ^{t+1} = rA \beta ^t \, }[/math]

[math]\displaystyle{ r = \beta \, }[/math]

And hence the solution to the difference equation is

[math]\displaystyle{ x_t = x_0 r^t \, }[/math]

2) Z-transform - the discrete equivalent to the Laplace transform, an organized way to find solutions

What are the behaviors of the solution [math]\displaystyle{ x_t = x_0 r^t \, }[/math]?

There are 4 different regimes, depending on the value of r.

If r > 1, then x_t approaches infinity.

If 0 < r < 1, then x_t approaches zero.

If -1 < r < 0, then x_t approaches zero in an oscillatory manner.

If r < -1, then x_t approaches infinity in an oscillatory manner.

Remember that if we have non-linear systems, it is difficult to obtain the global stability analysis of the problem, however, we can see from the above different behaviors that we can generalize it into stable and unstable behaviors as shown in the table below.

Stability Analysis

r States	Stability
[math]\displaystyle{ |r| \lt 1 \, }[/math]	Stable
[math]\displaystyle{ |r| \gt 1 \, }[/math]	Unstable
[math]\displaystyle{ |r| = 1 \, }[/math]	Liapunov Stable

What does Liapunov stability imply? The fixed point is neither stable (goes towards the fixed point) nor unstable (goes away from the fixed point), but the trajectory is bounded. This is similar to a limit cycle or a center in continuous differential equations.

Higher Order Difference Equations

A second order difference equations would tell us about the state of the system two time steps away, for example:

[math]\displaystyle{ x_{t+2} = ax_{t+1} + bx_t \, }[/math]

How do we go about solving this system? We again solve by guessing for the correct solution. We have our initial guess again as [math]\displaystyle{ x_t = Ar^t }[/math]. Substituting into our original equation, we get:

[math]\displaystyle{ Ar^{t+2} = aAr^{t+1} + bAr^t \, }[/math]

The characteristic equation for this system is

[math]\displaystyle{ r^2 - ar - b = 0 \, }[/math]

Solving for r

[math]\displaystyle{ r_{\pm} = \frac{a \pm \sqrt{a^2+4b}}{2} }[/math]

Hence, our general solution will be the sum of all possible solutions, just like in second order differential equations

[math]\displaystyle{ x_t = Ar_{+}^t + Br_{-}^t }[/math]

We can then solve for A and B with our initial conditions.

For the stability of the system, we know that for the system not to blow up to infinity, [math]\displaystyle{ |r_+|, |r_-| \lt 1 }[/math]. For this condition to be met, the following condition must be satisfied (can be solve for by setting r = 1).

[math]\displaystyle{ b \lt 1 - a \, }[/math]

Furthermore, we can generalize this solution for an N-dimensional system:

[math]\displaystyle{ a_0 x_{t+N} + a_1 x_{t+N-1} + \dots + a_N = 0 }[/math]

We again make the guess that our solution is in the form [math]\displaystyle{ x_t = Ar^t }[/math], substitute, and obtain our characteristic polynomial to be:

[math]\displaystyle{ \alpha_0 r^N + \alpha_1 r^{N-1} + \dots + \alpha_N = 0 }[/math]

The roots of our polynomial give the values of r and they must all satisfy [math]\displaystyle{ |r| \lt 1 }[/math] for the entire system to be stable. With the root of the equation being: [math]\displaystyle{ \{r_1^*, \dots, r_N^*\} }[/math], our general solution becomes:

[math]\displaystyle{ x_t = \sum_{i = 1}^N A_i(r_i^*)^t }[/math]

Solving these difference equations is non-trivial and computationally difficult. There are criteria (see Schur or Jury) to check to see if the parameters will lead to a stable conditions (which check to see if the roots are within the unit circle). We can only solve by hand to a maximum of 4 dimensions. With higher order dimensions, a computer becomes necessary to establish the roots of the equation.

Fourier Analysis and Linear Discrete Time Analysis

We can first motivate this by an example. Let us consider a single sinusoid in the time domain.

[math]\displaystyle{ y(t) = Asin(\omega t) \, }[/math]

Is it possible to write this function in terms of a difference equation? Consider the time steps t and t+1 for a sinusoidal function:

[math]\displaystyle{ \begin{alignat}{2} y_t & = A sin(\omega t) \\ y_{t \pm 1} & = A sin(\omega (t \pm 1)) \\ & = A[sin(\omega t)cos(\omega) \pm sin(\omega)cos(\omega t)] \end{alignat} }[/math]

If we now add [math]\displaystyle{ y_{t+1} \, }[/math] and [math]\displaystyle{ y_{t+1} \, }[/math] we obtain:

[math]\displaystyle{ \begin{alignat}{2} y_{t+1} + y_{t-1} & = 2A cos(\omega)sin(\omega t) \\ & = 2A cos(\omega)y_t \end{alignat} }[/math]

Hence we can conclude that:

[math]\displaystyle{ y_{t+1} = 2Acos(\omega)y_t - y_{t-1} \, }[/math]

We can see that a sinusoid function in the time domain can be expressed in terms of a second order difference equation. In general, if we have a sum of N sinusoidal waves making up our function, we can write this as a difference equations with 2N terms. In mathematical terms, this can be expressed as:

[math]\displaystyle{ \begin{alignat}{2} y_t & = \sum_{m=1}^N a_m sin(m \omega t + \phi_m) \\ & \equiv \\ y_{t+2N} & = \sum_{m = 0}^{2N-1} b_m y_{t-m} \end{alignat} }[/math]

Recall that the first equation above is just the Fourier transform of a function in time domain. Hence, for any function in the time domain, this can be written as either a series of sinusoids using Fourier transform or as a difference equation with two terms to each Fourier term. We can also see that the coefficients obtained from performing the Fourier transform are directly related to the coefficients of the difference equation that is generated.

The process of going from the time domain to the difference equation is known as Linear Discrete Time Analysis, since we can see that we obtain a linear difference equation from a sinusoidal function. This is analogous to getting a linear second order differential equation when we are considering the same sinusoid when we do continuous differentiation.

Logistic Map

Let us now consider the logistic map equation as an example of how linear stability analysis is preformed on difference equations/systems.

[math]\displaystyle{ x_{t+1} = rx_t(1-x_t) \, }[/math]

Where 0 < r < 4 are the boundaries set.

First, to find the fixed points, we set [math]\displaystyle{ x_{t+1} = x_t }[/math] to obtain the two fixed points:

[math]\displaystyle{ x_t^* = 0, 1-\frac{1}{r} }[/math]

For linear stability analysis, we need that [math]\displaystyle{ \left | \frac{df}{dx_t} \right | \lt 1 }[/math]

[math]\displaystyle{ f(x_t) = x_{t+1} = rx_t(1-x_t) \, }[/math]

[math]\displaystyle{ \frac{df}{dx_t} = r - 2rx_t }[/math]

Now evaluating at the fixed points

[math]\displaystyle{ f'(x_t^* = 0) = r }[/math]

[math]\displaystyle{ f'(x_t^* = 1 - \frac{1}{r}) = -r + 2 }[/math]

Hence, the stability of the [math]\displaystyle{ x_t^* = 0 }[/math] fixed point is stable if 0 < r < 1 (recall that we are bounded between 0 < r < 4 for our system). The stability of the [math]\displaystyle{ x_t^* = 1 - \frac{1}{r} }[/math] fixed point is stable for 1 < r < 3.

After r = 3, we see that the system is unstable and that there are no fixed points. In fact, we will oscillate between fixed points and observe what is known as period doubling in our system. How do we calculate these fixed points analytically? If we have a period 2 oscillation, that means that

[math]\displaystyle{ x_{t+2} = x_t \, }[/math]

And from our function, we know that

[math]\displaystyle{ x_{t+1} = f(x_t) \, }[/math]

And by substitution

[math]\displaystyle{ x_{t+2} = f(x_{t+1}) = f(f(x_t)) \, }[/math]

In general, to solve for p period solutions, we want that

[math]\displaystyle{ x_{t+p} = x_t \, }[/math]

And we can solve

[math]\displaystyle{ x_{t+p} = f(f...f(f(x_t)))_p \, }[/math]

Returning back to our period-2 oscillations, we substitute our logistic equation into the [math]\displaystyle{ x_{t+2} }[/math] equation

[math]\displaystyle{ x_{t+2} = f(f(x_t)) = r[rx_t(1-x_t)][1-rx_t(1-x_t)]\, }[/math]

We now set this equal to [math]\displaystyle{ x_t }[/math] and solve for the fixed points. As it can be seen that both our initial fixed points are fixed points of this system with two further fixed points being added. The two new fixed points are the points around which the system oscillates.

The bifurcations occur when the looking at the difference plot and seeing when you get more crossings on the graph. As the value of r is increased, the bifurcations occur more closely. After a certain value of r, you can get a point where there are no period solutions exist and this is the point where chaos emerges. Within the chaotic region, numerical analysis has yielded windows of periodicity such as period-3 oscillations, but quickly become chaotic again as r is changed.

With the analysis above, we can easily see that a fixed point is actually a period-1 solution to our system.