Second Order Homogenous Linear Difference Equations

The questions below are due on Friday February 21, 2025; 10:00:00 AM.
 
You are not logged in.

Please Log In for full access to the web site.
Note that this link will take you to an external site (https://shimmer.mit.edu) to authenticate, and then you will be redirected back to this page.
You are viewing the 6.3100 Spring 2025 web site. Please note that information on this page may change throughout the semester, please check back frequently.

Crane Mutiny Example

Crane Seeking a Beacon

A heavy cargo crane travels on track above a buried beacon, and measures its signed distance to the beacon 100 times a second (by signed distance, we mean that the distance is positive if the crane has not yet reached the beacon and negative if the crane has traveled past it). A human operator usually drives the crane and parks it just above the beacon. The crane then keeps itself in place automatically by always accelerating towards the beacon at a rate proportional to its distance from the beacon. On one unfortunate day, the operator noticed that the crane seemed to be "rocking" back and forth over the beacon, and tried to fix the problem by flipping the sign of the distance sensor. The crane almost immediately sped off down the track, halted only by the quick-thinking operator, who shut down the crane. Why did this happen?

In the previous prelab, we assumed that we could control the crane by sending it velocity commands, which led to a first-order difference equation model. A more realistic model would be to assume that we can send the Crane acceleration commands, and that will lead to a more complicated higher-order difference equation model. We will start in the same way, and will represent the crane's measured distance to the beacon as a sequence, d, where d[n] is the n^{th} measured distance. Then d[n] can be related to d[n-1] by

d[n] \approx d[n-1] + \Delta T \; v[n-1]
where \Delta T is the time between distance measurements, v[n-1] is the crane velocity at the time of the (n-1)^{th} distance measurement, and we note that the relationship is approximate because the velocity of an accelerating crane varies continuously. Since crane acceleration is proportional to measured distance, it does not vary between samples, so we have an exact relation between velocity and acceleration. That is,
v[n] = v[n-1] + \Delta T \; a[n-1] = v[n-1] - \Delta T \; K_p d[n-1]
where a[n-1] is the crane acceleration at the time of the (n-1)^{th} distance measurement, and K_p is the proportional gain relating acceleration to measured distance.

There are now two sequences, v and d, and two coupled linear difference equations. As we will see below, we can use vectors and matrices to generalize many of the results from the first prelab, and use them to analyze dour Crane problem.

Matrix Difference Equations

The generic form for a vector first-order homogenous difference equation is

x[n] = \textbf{A} \; x[n-1],
where x is a vector sequence whose n^{th} entry, x[n], is a vector of length N, and \textbf{A} is an N \times N matrix. The solution, x[n] as a function of an initial condition, x[0], is
x[n] = \textbf{A}^n x[0].
The solution to the scalar first-order homogenous difference equation, y[n] = \lambda y[n-1], has the same form,
y[n] = \lambda^n y[0],
and we already know that the behavior of the sequence y depends on properties of \lambda. For example, if | \lambda | < 1, then \lim_{n\rightarrow \infty} y[n]=0, and if \lambda < 0, then y oscillates.

So what properties of \textbf{A} determine the properties of the sequence x? As we will see below, for the vector case, sequence properties are determined by the eigenvalues of \textbf{A} (assuming diagonalizability). For example, if all N eigenvalues of \textbf{A}, |\lambda_i| i \in {1,...N}, are less than one in magnitude (strictly inside the unit circle), then \lim_{n\rightarrow \infty} x[n] = 0.

To see the above result, recall that if \lambda_i is the i^{th} eigenvalue of an N \times N matrix \textbf{A}, and \textbf{V}_i is the associated N-length eigenvector, then

\textbf{A} \textbf{V}_i = \lambda_i \textbf{V}_i.
If \textbf{A} is diagonalizable (as almost all matrices are), then there exists an N \times N matrix \textbf{V} whose columns are a linearly independent set of eigenvectors of \textbf{A}. That is,
\textbf{A} \begin{bmatrix} \textbf{V}_1 \textbf{V}_2 ... \textbf{V}_N \end{bmatrix} = \begin{bmatrix} \textbf{V}_1 \textbf{V}_2 ... \textbf{V}_N \end{bmatrix} \begin{bmatrix} \lambda_1 & 0 & ... & ... & 0\\ 0 & \lambda_2 & 0 & ... & 0 \\ ... & ... & ... & ... & ...\\ 0 & ... & ... & 0 & \lambda_N \end{bmatrix},
or
\textbf{A} \textbf{V} = \textbf{V} \Lambda.
And since the columns of \textbf{V} are linearly independent, \textbf{V} is invertible and
\textbf{A} = \textbf{V} \Lambda \textbf{V}^{-1},
which is sometimes referred to as the eigendecomposition of \textbf{A}.

We can use the eigendecomposition to determine a formula for powers of \textbf{A}. For example, \textbf{A}^2 is

\textbf{A}^2 = \textbf{V} \Lambda \textbf{V}^{-1} \textbf{V} \Lambda \textbf{V}^{-1} = \textbf{V} \Lambda \Lambda \textbf{V}^{-1} = \textbf{V} \Lambda^2 \textbf{V}^{-1}.
and by induction
\textbf{A}^n = \textbf{V} \Lambda^n \textbf{V}^{-1}.
Plugging into the matrix difference equation,
x[n] = \textbf{A}^n x[0] = \textbf{V} \Lambda^n \textbf{V}^{-1} x[0]= \textbf{V} \begin{bmatrix} \lambda_1^n & 0 & ... & ... & 0\\ 0 & \lambda_2^n & 0 & ... & 0 \\ ... & ... & ... & ... & ...\\ 0 & ... & ... & 0 & \lambda_N^n \end{bmatrix} \textbf{V}^{-1} x[0].

As the above equation makes clear, the behavior of the vector sequence x is determined by the eigenvalues of \textbf{A} just like the behavior of the scalar sequence y is determined by \lambda. For example, if \max_{i} |\lambda_i| < 1 then

\lim_{n \rightarrow \infty} x[n] = 0.
The above equation also implies that
x_i[n] = c_{i,1} \lambda_1^n + c_{i,2} \lambda_2^n + ... + c_{i,N} \lambda_N^n,
where the constants c_{i,j} can be determined from the initial condition x[0], \textbf{V}, and \textbf{V}^{-1}. Alternatively, one can solve for c_{i,1}, c_{i,2}, ... c_{i,N} given N values of x_i[n]. For example, given x_i[0],...,x_i[N-1], solve for the constants from
x_i[0] = c_{i,1} + c_{i,2} + ... + c_{i,N},
x_i[1] = c_{i,1} \lambda_1 + c_{i,2} \lambda_2 + ... + c_{i,N} \lambda_N,
...
x_i[N-1] = c_{i,1} \lambda_1^{N-1} + c_{i,2} \lambda_2^{N-1} + ... + c_{i,N} \lambda_N^{N-1}.

Crane Example

We can write the equations for the crane as as a first-order vector equation

\begin{bmatrix}d[n] \\ v[n]\end{bmatrix} = \begin{bmatrix} 1 & \Delta T\\ -\Delta T K_p & 1\end{bmatrix} \begin{bmatrix}d[n-1] \\ v[n-1]\end{bmatrix},
or more generically,
x[n] = \textbf{A}x[n-1],
where x[n] is a vector (of length 2 in this case) and \textbf{A} is a matrix (of size 2 \times 2 in the above case).

x[n] \equiv \begin{bmatrix}d[n] \\ v[n]\end{bmatrix}
and
\textbf{A} \equiv \begin{bmatrix} 1 & \Delta T\\ -\Delta T K_p & 1\end{bmatrix}.

For a two-by-two system, we can determine the two eigenvalues using determinants, in which case the eigenvalues are the roots of the quadratic equation given by

det(\lambda \textbf{I}- \textbf{A}) = 0 = (\lambda-1)^2 +\left(\Delta T\right)^2 K_p = \lambda^2 - 2 \lambda + 1 + \left(\Delta T\right)^2 K_p,
which are (assuming K_p > 0)
\lambda_1,\lambda_2 = 1+j\Delta T \sqrt{K_p}, 1-j\Delta T \sqrt{K_p},
and therefore
d[n] = c_1 \left( 1+j\Delta T \sqrt{K_p} \right) ^n + c_2 \left( 1-j\Delta T \sqrt{K_p} \right) ^n.
We could then determine the values of c_1 and c_2 given d[0] and
d[1] = d[0] + \Delta T \; v[0],
which follows from the orginal difference equation.

The Crane example has complex natural frequencies, which we will return to in a later section. For now, we will stick to examples with real natural frequencies.

Fibonacci as LDE

Consider the Fibonacci sequence, in which the next element in the sequence is the sum of the previous two. The sequence has a second-order LDE description,

y[n] = y[n-1] + y[n-2].
can be written as a first-order vector difference equation
\begin{bmatrix}y[n] \\ y[n-1]\end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0\end{bmatrix} \begin{bmatrix}y[n-1] \\ y[n-2]\end{bmatrix}.

If we specify y[0] and y[1], which we often refer to as initial conditions, then we can use the vector difference equation to iteratively compute y[n] for all n \ge 2.

For example, to generate the classic Fibonacci numbers, set y[0] = 0 and y[1] = 1, and repeatedly apply the above formula with n =2 , then n=3, n=4, and so on, to determine y[2]=1, y[3]=2, y[4]=3, y[5]=5, ... etc. In general, given an LDE and enough initial conditions, it is easy to explicitly enumerate the rest of the sequence values.

These questions refer to the Fibonacci difference equation

y[n] = y[n-1] + y[n-2].

  1. If you pick initial conditions y[0] = 8 and y[1] = 13 for the above difference equation, what will be true about the generated sequence y[2], y[3], y[4], y[5], ... ? Please list all that apply below (separated by commas).
  • A: The sequence entries will grow without bound.
  • B: The sequence entries match the values of a shifted Fibonacci sequence.
  • C: \lim_{n \rightarrow \infty} \frac{y[n]}{y[n-1]} > 2.

Suppose the initial conditions for the difference equation are given by two constants, A and B , where y[0] = A and y[1] = B . What will be true about the generated sequence y[2], y[3], y[4], y[5], ... ? Please check all that apply.

  • A: If A is negative and B is positive, and |A| \gt B , then it is guaranteed that y[n] \lt 0 for all n \gt 2 .
  • B: If A is negative and B is positive, and |A| \gt 2\cdot B , then it is guaranteed that y[n] \lt 0 for all n \gt 2 .
  • C: If A is positive and B is negative, and A \lt |B| , then it is guaranteed that y[n] \lt 0 for all n \gt 2 .
  • D: If A is positive and B is negative, and |A| \gt |2\cdot B |, then it is guaranteed that y[n] \lt 0 for all n \gt 2 .

"Generalized" Fibonacci

Sequences described by second-order homogeneous LDEs, like the Fibonacci sequence, have the property that y[n] is given by a weighted combination of y[n-1] and y[n-2]. Their general form is

y[n] = a_1 y[n-1] + a_2 y[n-2],
and have solutions of the form
y[n] = c_1 \lambda _1^n + c_2 \lambda_2^n
where \lambda _1 and \lambda_2 are the eigenvalues of the matrix associated with the vector form of the above difference equation, and c_1 and c_2 can be determined from the initial conditions y[0] and y[1].

Suppose we do not know a_1 and a_2, but we do know the solution has the form

y[n] = c_1 \cdot 3^n + c_2 \left(\frac{-1}{2}\right)^n.

If y[n] goes to zero as n \rightarrow \infty , what is c_1 ?

If y[n] \ge 0 for all n \ge 0 and c_2 = 3 , then what is the smallest possible value for c_1 (give answer as a decimal number)?

For what values of a_1 and a_2 will y[n] = a_1 y[n-1] + a_2 y[n-2] have as solution of thee form

y[n] = c_1 \cdot 3^n + c_2 \left(\frac{-1}{2}\right)^n.
(You will have to form the two-by-two matrix and compute is eigenvalues!)

What is the value of a_1:

...and a_2:?

For the following problems, consider the second-order linear difference equation

y[n] = a_1 y[n-1] + a_2 y[n-2],

with the solution

y[n] = 10(-2)^ n + 5(-1)^n.

Determine a_1, a_2, y[0], and y[1].

a_1=

a_2=

y[0]=

y[1]=

Matching Initial Conditions

Now that we have the two natural frequencies, there is the matter of initial conditions. How to incorporate the initial conditions is best seen by example, so we will examine solving the second-order homogeneous LDE for the Fibonacci sequence. Since the next value in the Fibonacci sequence is the sum of the previous two,

y[n] = y[n-1] + y[n-2],

or a_1 = a_2 = 1 in our general second-order homogeneous LDE.

From the eigenvalues of the matrix \textbf{A} of the vector form of the above difference equation, we know that the Fibonacci LDE has two distinct natural frequencies,

(\lambda _1, \lambda _2) = \frac{1 \pm \sqrt {5}}{2},

therefore any sequence of the form

y[n] = c_1 \left(\frac{1 + \sqrt {5}}{2}\right)^ n + c_2 \left(\frac{1 - \sqrt {5}}{2}\right)^ n,

will satisfy the Fibonacci LDE.

In order to produce the classical Fibonacci sequence, we need to find values for c_1 and c_2 so that y[0] = 0 and y[1]= 1. To start, consider how y[0] = 0 constrains c_1 and c_2,

y[0] = c_1 \left(\frac{1 + \sqrt {5}}{2}\right)^0 + c_2 \left(\frac{1 - \sqrt {5}}{2}\right)^0=0
c_1 + c_2=0
c_1=-c_2

Including the constraint from y[1] = 1,

y[1] = c_1 \left(\frac{1 + \sqrt {5}}{2}\right)^1 + c_2 \left(\frac{1 - \sqrt {5}}{2}\right)^1=1

c_1 \frac{1 + \sqrt {5}}{2} - c_1 \frac{1 - \sqrt {5}}{2}=1
c_1 ((1 + \sqrt {5}) - (1 - \sqrt {5}))=2
c_1 2 \sqrt {5}=2
c_1 = \frac{1}{\sqrt {5}}

So, the final equation for the Fibonacci sequence is

y[n] = \frac{1}{\sqrt {5}} \left(\frac{1 + \sqrt {5}}{2}\right)^ n - \frac{1}{\sqrt {5}} \left(\frac{1 - \sqrt {5}}{2}\right)^ n\; \; .

Evaluating the numeric constants (approximately)

y[n] = 0.44721 \cdot 1.618^ n - 0.44721 \cdot (-0.618)^ n

Plots of the two components of the Fibonacci sequence, and their sum. The component in the leftmost graph is due to the negative natural frequency, and we can clearly see its oscillation, but note that its magnitude is very small. The lower sequence is the sum of the first two, and we can clearly see that the exponential growth associated with the large positive natural frequency dominates, and reveals the familiar large n approximation for the Fibonacci sequence.

Are there values of y[0] and y[1] for which the sequence given by

y[n] = c_1 \left(\frac{1 + \sqrt {5}}{2}\right)^ n + c_2 \left(\frac{1 - \sqrt {5}}{2}\right)^ n

decays to 0 as n approaches \infty?

Select the most accurate response below:

Consider a system modeled by the difference equation

y[n]= -2y[n-1]-y[n-2]

with initial conditions y[0]=1 and y[1]=-1.

Which of the following most accurately characterizes the eventually (large n) behavior of the sequence y?