First Order Difference Equations

The questions below are due on Friday April 03, 2026; 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 PRELIMINARY 6.3100 Spring 2026 web site. Please note that information on this page may change throughout the semester, please check back frequently.

First-Order Difference and Differential Equations

In this problem we use first-order differential equations as a simple setting in which to introduce many continuous-time concepts including those you have seen before (e.g. step response and steady-state) as well as new ones (e.g. transfer functions and frequency response).

Discrete Time Reminder

As a reminder, our notation for a general first-order difference equation with input sequence u and output sequence y is

y[n] = \lambda y[n-1] + \gamma \; u[n-1]\;\;,

where \lambda and \gamma as real scalars and n is an integer index. The general solution can be derived by induction as is

y[n] = \lambda^n y[0] + \gamma \sum_{m=0}^{m=n-1} \lambda^{(n-1)-m} u[m],

where y[0] is a given initial condition. We also decomposed y into its zero-input response,

y_{ZIR}[n] = \lambda^n y[0],
and zero-state response (where we assume y[0] = 0),
y_{ZSR}[n] = \gamma \sum_{m=0}^{m=n-1} \lambda^{(n-1)-m} u[m].
We noted that if the system was stable, equivalently if |\lambda| < 1, then
\lim_{n \rightarrow \infty} y_{ZIR}[n] = \lim_{n \rightarrow \infty} \lambda^n y[0] = 0,
and
\lim_{n \rightarrow \infty} y[n] \equiv y[\infty] = \frac{\gamma}{1 -\lambda} u[\infty],
assuming u[n] = u[\infty] for all n. We often refer to y[\infty] as the steady-state.

Aside: The input u does not have to be constant for there to be a steady-state. It is enough for u[n] = u[\infty] for n > n_0 where n_0 is a finite positive integer (u settles to a constant for n large enough).

Continuous Time Solutions

In continuous time, our notation for a general first-order differential equation with input waveform (i.e. function of time) u and output waveform y is

\frac{d}{dt} y(t) = \lambda y(t) + \gamma \; u(t),
where \lambda and \gamma are again real coefficients, and \frac{d}{dt} denotes time-derivative. The general solution is given by
y(t) = e^{\lambda t}y(0) + \gamma \int_0^t e^{\lambda (t - \tau)}u(\tau) d\tau,
where y(0) is a given initial condition. The above solution is easily verified by differentiating the above equation, then simplifying using a property of the exponential,
\frac{d}{dt} e^{\lambda t} = \lambda e^{\lambda t},
and applying integration-by-parts.

Aside: Since the time-derivative of e^{\lambda t} is equal to \lambda e^{\lambda t}, a scaling of the same function function, we often say that e^{\lambda t} is an eigenfunction of the differentiation operator, with eigenvalue \lambda.

We can decompose y(t) into its zero-input response,

y_{ZIR}(t) = e^{\lambda t}y(0),
and zero-state response (where we assume y(0) = 0),
y_{ZSR}(t) = \gamma \int_0^t e^{\lambda (t - \tau)}u(\tau) d\tau.

CT Natural Frequency

Recall that in discrete time, if some non-zero y_{ZIR}[n] satisfies

y_{ZIR}[n] = \lambda y_{ZIR}[n-1]
for all n, then \lambda is the natural frequency of the DT system. If \lambda < 0, then y_{ZIR}[n] oscillates and if |\lambda| < 1, y_{ZIR}[n] decays to zero.

In continuous time, if some non-zero y_{ZIR}(t) satisfies

\frac{d y_{ZIR}(t)}{dt} = \lambda y_{ZIR}(t),
for all t, then we say \lambda is a natural frequency of the CT system. As always, the use of the term "frequency" seems odd, because for any real \lambda, y(t) = e^{\lambda t} y(0) never oscillates, but is a decaying or growing exponential, depending on the sign of \lambda, as in the figure below.

Figure 1: e^{0.02t} (purple) and e^{-0.02t} (yellow) versus t

For higher order difference equations, the DT natural frequencies were roots of higher-order polynomials, and therefore real equations could have complex natural frequencies. The same will be true for CT natural frequencies, they will be the roots of higher-order polynomials and can be complex. But we will return to the higher-order case later, for now, let us consider an example of proportional feedback leading to difference or differential equations.

Steady-State

If the system is stable, equivalently if the real part of \lambda < 0, then

\lim_{t \rightarrow +\infty} y_{ZIR}(t) = \lim_{t \rightarrow +\infty} e^{\lambda t}y(0) = 0.
If, in addition to the system being stable, the input is a constant (u(t) = u(\infty) for all t > 0), then the system has a steady-state given by
\lim_{t \rightarrow +\infty} y(t) \equiv y(\infty) = -\frac{\gamma}{\lambda} u(\infty),
which is easily derived by setting \frac{d}{dt} to zero in the general form for the differential equation, or by examining the limit of the ZSR,
\lim_{t \rightarrow +\infty} y_{ZSR}(t) = \lim_{t \rightarrow +\infty} \gamma \int_0^t e^{\lambda (t - \tau)}u(\infty) d\tau.
Simplifying,
\lim_{t \rightarrow +\infty} y_{ZSR}(t) = \gamma u(\infty) \lim_{t \rightarrow +\infty} \int_0^t e^{\lambda (t - \tau)} d\tau = -\frac{\gamma}{\lambda} u(\infty).

Using Linearity and Time-Invariance

First-order differential equation have the same linearity and time-invariance (LTI) properties as their difference-equation counterparts, and we can often use these properties to quickly calculate responses. For example, we are often interested in the zero-state responses of stable systems to step changes in input, such as in

u(t) = 0\;\;\;t < t_0 \;\;\;\;\;\; u(t) = 10\;\;\;t \ge t_0.

We can exploit the time-invariance properties by decomposing u(t) into a constant input, u_B(t), plus a step from a negative value to zero, u_A(t), as shown in the image below.

The ZSR response, given the above u(t), is then the sum of a steady-state (y[\infty] given u[\infty]=10) and a shifted ZIR response. Specifically, y(t) = 0 for t < t_0, and for t \ge t_0 is

y(t) = \frac{-10\gamma}{\lambda}\left(1 - e^{\lambda(t-t_0)}\right).

Suppose we change the input to

u(t) = 0,\;\;t < t_0 \;\;\;\;\;\; u(t) = 10,\;\; t_0 \le t \le 2t_0 \;\;\;\;\;\; u(t) = 0,\;\; t > 2t_0.
Please determine an analytic formula for y(t) for t > 2t_0 (please enter the formula as a legal python expression, using G for \gamma, L for \lambda, t0 for t_0 and e**(L*t) for e^{\lambda t}.

Enter a python expression for y(t) for t > 2t0. y(t) =

Continuously Crane

Figure 2: Crane traveling on a track and measuring

Let us return to the cargo crane traveling on track above a buried beacon, as shown in Figure 2. Suppose we must decide between two control systems, a discrete-time controller that measures distance to the beacon every \Delta T seconds, and a continuous-time controller that measures the distance continuously. One way to compare the two systems is to examine the problem of returning the crane as rapidly as possible to "home" position, that is, back to the beacon at d = 0. And for this comparison, we are going to assume that the crane has an internal feedback controller, one allows us to design a controller that senses position and controls crane velocity (as show in Figure 3 below).

Figure 3: Crane DT proportional controller.

Model

For our velocity-controlled crane with discrete-time sensor,

d[n] = d[n-1] + \Delta T \cdot v[n-1],
where d[n] = d(n \Delta T) is the distance from the beacon at time n \Delta T, and v[n-1] is the crane's constant velocity from from (n-1)\Delta T to n \Delta T.

If we use proportional control to return the robot back to the beacon (d=0), then the DT controller velocity at n \Delta T is given by

v[n] = -K_p d[n],
and the combination of controller and plant leads to a homogenous difference equation for d,
d[n] = (1-\Delta T K_p) d[n-1].

For the continuous-time controller, crane distance and velocity are related by a derivative,

\frac{d}{dt} d(t) = v(t),
where d(t) and v(t) are the robot's distance and velocity at time t. If we use proportional control to return the robot to the beacon,
v(t) = -K_p d(t),
leading to a homogenous differential equation for d,
\frac{d}{dt} d(t) = -K_p d(t).
Note that the solution to the CT homogenous equation is
d(t) = d(0)e^{-K_p t}
where d(0) is the initial distance from the beacon. Note that
\lim_{t->\infty} d(t) = 0
for any K_p > 0. That is, for proportional control on a first-order CT system, THE FEEDBACK SYSTEM IS STABLE FOR ANY POSITIVE GAIN .

Inferring Parameters

Suppose the cranes are being independently monitored with remote LIDAR, which measures the distance between the crane and the beacon every ten seconds. A test was run using the CT controller with an unknown proportional gain K_p, in which two consecutive measurements from the LIDAR system (which we can denote as d(t) and d(t+10)) gave distances of 200 meters and then 10 meters between the beacon and the crane.

What was the third consecutive LIDAR measurement, d(t+20), and to within one percent, what was the controller gain, K_p?

d(t+20)=

Kp=

Complex Exponential Inputs

For the zero-input case, the differential equation

\frac{d}{dt} y(t) - \lambda y(t) = 0,
has a family of solutions (zero-input responses, ZIRs) parameterized by a scalar A, as in
y_{ZIR}(t) = Ae^{\lambda t}.

If the input to the different equation is a complex exponential,

u(t) = Ue^{st}
where both U and s are (possibly complex) scalars, then there is a complex scalar Y such that
y_p(t) = Ye^{st}
where the particular solution, y_p(t), satisfies the differential equation (but not necessarily any initial conditions). To see this, substitute for y and u in the first order differential equation, as in
\frac{d}{dt} Ye^{st} - \lambda Ye^{st} = \gamma Ue^{st},
then divide by e^{st} on both sides and reorganize,
Y = \frac{\gamma}{s-\lambda}U,
and therefore
y_p(t) = \frac{\gamma}{s-\lambda}Ue^{st}.

Since the differential equation is linear,

y_{ZIR}(t) + y_p(t) = Ae^{\lambda t} + \frac{\gamma}{s-\lambda}Ue^{st},
will satisfy the differential equation for any A, so we can adjust A to match any given initial condition y(0). That is,
y(t) =y_{ZIR}(t) + y_p(t) = Ae^{\lambda t} + \frac{\gamma}{s-\lambda}Ue^{st}
where A is chosen so that y(t)|_{t=0} matches y(0), the given initial condition.

As we will see below and more generally in later sections, the function of s that "transfers" the input u to the output y is worth labeling, as in

H(s) \equiv \frac{\gamma}{s-\lambda}.

Step Response and Parameter Estimation

If y(0) = 0 and u(t) = 1 for t \ge 0, then u(t) = 1 e^{st} for s = 0 and

y(t) = Ae^{\lambda t} + \frac{\gamma}{0-\lambda} 1 e^{0t} = \frac{\gamma}{- \lambda} \left(1 - e^{\lambda t} \right),
where the rightmost equality is derived by selecting A to match y(0) = 0.

A plot of the step response for \gamma = 10 and \lambda = -100 is given below. Note that the step response settles to \frac{\gamma}{-\lambda} as t \rightarrow +\infty, and has risen to about two-thirds of the steady-state at t \approx \frac{1}{-\lambda}.

Figure 4: Step Response for \gamma = 10 and \lambda = -100 , halfway point circled in red.

When calibrating a model using experimental data, one should design the experiments so that the needed data is both easy to measure and informative. For example, it is easy to measure the time for a step response to rise half way to steady-state (the half way point in the above figure is circled in red, and the time to rise to that halfway point is labelled Tau). Assuming the step response is well-described by the above first-order model, the half way point is given by

y_{Halfway} = 0.5 \frac{\gamma}{- \lambda} = \frac{\gamma}{- \lambda} \left(1 - e^{\lambda \tau} \right),
where we denote the time to rise half way as \tau.

What model parameter can you estimate from the "half-way rise time" \tau, and what is the analytic formula for that parameter once you are given \tau? (please enter your answers as legal python expressions, using G for \gamma, L for \lambda, Tau for \tau, e**(L*t) for e^{\lambda t} and log(a) for the natural log of a (that is, log(e^a) = a).

Which model parameter can be estimated from half-way risetime?

Give an expression for the estimated parameter as a function of Tau

Crane Control with Step Input

If we want to move the crane to positions other than the beacon position at d=0, we will need to specify a desired position along the track. The desired position, which we will denote as d_{in}(t), will become an input to our position controller, which will try to keep the drane position, d(t), close to the desired position d_{in}(t).

A proportional feedback approach to controlling the crane's position would be to make the crane velocity proportional to the difference between the desired and measured crane positions, as diagrammed for the DT case in the figure above. Then the velocity is given by

v(t) = K_p(d_{in}(t) - d(t)).
where K_p is the proportional gain and d_{in}(t) is the input, or desired, position at time t.

In order to determine how the performance of this proportional-feedback controller changes as we increase or decrease gain, we substitute the control formula in to the distance-velocity relation for the crane,

\frac{d}{dt} d(t) = K_p \left( d_{in}(t) - d(t) \right),
and then reorganize as
\frac{d}{dt} d(t) = - K_p d(t) + K_p d_{in}(t).

Since d_{in}(t) = 1 = e^{0t}, and d(0) = 5, we can determine the analytic solution for d(t) as a function of K_p and t (Hint: specialize the general result from the previous subsection). Use a python expression, use Kp to for K_p and use, for example, e**(-a*t) to represent e^{-at}.

Enter a python expression for d(t) as a function of Kp and t

Figure 5: d(t) versus t for four different values of K_p

Above are four plots of d(t) versus time for the case d(0) = 5 and d_{in}(t) = 1. There are four plots, blue, orange,yellow and purple, and each plot is rassociated with a different value of K_p. Assuming that the four K_p values are integers in the range \{1,...20\}, what are the four gains in blue, orange, yellow, and purple (hint, e^{-1} \approx \frac{1}{3}). Enter your solution as a python list, e.g. [5,10,15,20].

K_p for [blue,orange,yellow,purple] =

Frequency Response

As you have already seen in lab, we can learn many characteristics of a stable system by applying sine waves to the system, waiting for the system to achieve sinusoidal steady-state, and then measuring that steady-state response. We show a diagram of that process below (note that in 6.3100 the "oscillator" and "oscilloscope" are virtual. That is, they are implemented using a combination of software, the Teensy controller board, and your laptop).

For any LTI system in sinusoidal-steady-state, the output frequency always matches the input frequency. It is the frequency-dependence of the input-output amplitude ratio and input-output phase difference, collectively referred to as the frequency response, that can reveal the specifics of an LTI system. As we will see below for the first-order differential equation case (and later more generally), we can connect the frequency response to the transfer function, H(s), mentioned above.

Sums of Complex Exponentials

We noted above that if the input, u(t), to our general first-order differential equation,

\frac{d}{dt} y(t) = \lambda y(t) + \gamma \; u(t),
is of the form
u(t) = Ue^{st},
then (using the transfer function H(s) \equiv \frac{\gamma}{s-\lambda} mentioned above),
y_p(t) = \frac{\gamma}{s-\lambda}Ue^{st} = H(s)Ue^{st}
satisfies the first-order differential equation and
y(t) = Ae^{\lambda t} + H(s)Ue^{st}
satisfies the differential equation and will match the initial condition if A is chosen correctly.

Linearity insures that we can generalize the above to inputs which are sums of complex exponentials. That is, if

u(t) = \sum_{k=1}^{k=K} U_ke^{s_k t},
the particular solution for the sum input would be
y_p(t) = \sum_{k=1}^{k=K} \frac{\gamma}{s_k-\lambda}U_ke^{s_k t} = \sum_{k=1}^{k=K} H(s_k)U_ke^{s_k t},
and
y(t) = Ae^{\lambda t} + \sum_{k=1}^{k=K}H(s_k) U_ke^{s_k t},
will satisfy the differential equation and, with the correct A, match the initial condition.

If our system is stable, that is \operatorname{Re}(\lambda) < 0, and the ZIR decays faster than any of the inputs, that is \operatorname{Re}(\lambda) < \operatorname{Re}(s_k) for all s_k, then if t is large enough

y(t) \approx \sum_{k=1}^{k=K} H(s_k) U_ke^{s_k t}.
That is, if we wait long enough, and our INPUTS DECAY SLOWLY ENOUGH, the response due to initial conditions and natural frequencies will decay away, and y(t) \approx y_p(t).

There is an important class of inputs that decay more slowly than the natural response of any stable system. Sines and cosines NEVER decay, which is why they are used for system characterization.

Sines, Cosines, Magnitude and Phase.

If u(t) = Ue^{st} with U real and s = j\omega t, then

u(t) = Ue^{st} = Ue^{j \omega t} = U\left(\cos{\omega t} + j\sin{\omega t}\right),
a complex (part real and part imaginary) function of time.

We can construct a real input, u_r(t), while retaining the exponential form by using the sum of two complex exponentials,

u_r(t) = {Ue^{j\omega t} + Ue^{-j\omega t} \over 2} = U\cos(\omega t).
Note that u_r(t) is REAL, and if the input to a real differential equation is real, the results must be real.

If our first-order differential equation is stable (\operatorname{Re}(\lambda) < 0) and the input is u_r(t), then for large enough t,

y(t) \approx {H(j\omega) Ue^{j\omega t} + H(-j\omega)Ue^{-j\omega t} \over 2}.
Rewriting H(j\omega) and H(-j\omega) in polar form
H(j\omega) = \left| H(j\omega)\right| e^{j \angle H(j\omega)} \;\;\;\; H(-j\omega) = \left| H(j\omega)\right| e^{-j \angle H(j\omega)},
and substituting,
y(t) \approx \frac{1}{2}\left|H(j\omega)\right|Ue^{j(\omega t+\angle H(j\omega))} + \frac{1}{2}\left|H(j\omega)\right|Ue^{-j(\omega t+\angle H(j\omega))} = \left|H(j\omega)\right|U\cos(\omega t + \angle H(j\omega)).

As the above equation makes clear, magnitude and angle (or phase) are the key characteristics of frequency response.

Measuring and Plotting Frequency Responses

Setting \gamma = 1 and \lambda = -1 in our first-order differential equation

\frac{d}{dt} y(t) = \lambda y(t) + \gamma \; u(t) = -y(t) + u(t),
for which the associated H(s) is
H(s) = {1 \over s + 1}.
or for s = j \omega
H(j\omega) = {1 \over j\omega + 1}.

For t >> 1 (so that e^{\lambda t} = e^{-1t} \approx 0) the sinusoidal steady state approximation is valid and

y(t) \approx |H(j\omega)|\cos(\omega t +\angle H(j\omega)) = {1 \over \sqrt{\omega^2 + 1}}\cos(\omega t - \tan^{-1}\omega).

The process of measuring the frequency response, applying sinusoids of progressively higher frequency and measuring the associated outputs, is best visualized as a sequence of experiments, as shown below.

The Bode plot format,in which magnitude and phase are plotted separately and stacked, is the most common approach for plotting frequency response, or H(j\omega) versus \omega. The top plot in the Bode format is a log-log plot of |H(j\omega)| versus \omega and the bottom plot is a log-linear plot of \angle H(j\omega) (in degrees) versus \omega. Below we show Bode plots for several difference transfer functions.

Frequency Response and Parameter Estimation

Measuring the halfway-rise-time of a system's step response is helpful for model calibration because it is easy to measure, informative, and a single step-response experiment will generate the needed data. As the above diagram makes clear, measuring frequency response involves many experiments, so is likely to be more time-consuming (but is also more reliable, as you may have already seen in lab, and more generalizable, as we will see later).

Nevertheless, one might hope that the sinusoidal-steady-state response to a sinusoidal input at a well-chosen frequency might simplify the parameter estimation problem. To see an example of this, consider the graph below and please answer the questions that follow.

The graph below shows two plots, an input u(t) (in pink) and output y(t) (in blue) for a system that is accurately modeled by our general first-order differential equation model,

\frac{d}{dt} y(t) = \lambda y(t) + \gamma \; u(t).

Use the graph, along with your insights about H(j\omega), to determine numerical values for \gamma and \lambda in the above model (be careful about 2*pi)!

\lambda =

\gamma =

Crane Control Frequency Response

The differential equation that describes the crane control system is given by

\frac{d}{dt} d(t) = - K_p d(t) + K_p d_{in}(t).
Since the system is stable for all K_p > 0, we can examine its frequency respoonse
H(j\omega) = H(j 2\pi f) = \frac{K_p}{j2\pi f + K_p}.
where \omega is the frequency in radians/second and f is the frequency in cycles/second.

Below is a Bode-ish plot of H(j\omega), or more precisely, a log-log plot of |H(j2\pi f)| versus f, in cycles/second, and log-linear plot of \angle H(j2\pi f) (in degrees) versus f, again in cycles/second. As you can see, we have removed the axis labels, so please determine magnitude A, angle C (assume angle B=0), and frequency D (in cycles/second NOT radians/second) as functions of K_p. (Hint: only one of the labels is dependent on K_p).

In answering the questions below, use Kp for K_p and pi for \pi.

Enter a python expression for magnitude A as a function of Kp

Enter a python expression for angle C as a function of Kp

Enter a python expression for frequency (in cycles/second) D as a function of Kp