Data Science and Computing with Python for Pilots and Flight Test Engineers
Analytic Solution Example
System Definition
Using what we have learned in the previous four lessons, we will give an application example of all these techniques by finding solutions a system described by the differential equation
$$ \ddot x(t) + 2 \dot x(t) + 17 x(t) = f(t) $$
We will find solutions for the two input functions unit impulse, \(f(t)=\delta(t)\), and unit step, \(f(t)=\theta(t)\).
Equation in Laplace Domain
Performing the Laplace transform on both sides of the equation yields:
$$ s^2 X(s) + 2sX(s) + 17 X(s) = F(s). $$
Here \(X(s)=\mathcal{L}[x(t)]\) is the Laplace transform of \(x(t)\) and likewise for the input function \(F(s)=\mathcal{L}[f(t)]\).
Transfer Function
Correspondingly, the transfer function for this system is:
$$ H(s):=\frac{X(s)}{F(s)} = \frac{1}{s^2+2s+17} $$
System Response \(x(t)\) to Unit Impulse \(f(t)=\delta(t)\)
Finding the Solution \(X(s)\) in the Laplace Domain
For the unit impulse, \(f(t) = \delta(t)\) (Dirac delta function) the Laplace transform is \(F(s)=\delta(s) =1\)
With this, the solution of the differential equation in the Laplace domain is
$$ X(s) = H(s)F(s) = H(s) = \frac{1}{s^2+2s+17} $$
In order to obtain the corresponding solution \(x(t)\) in the time domain, we need to apply the inverse Laplace transform to \(X(s)\): \(x(t)=\mathcal{L}^{-1}[X(s)]\).
Writing \(X(s)\) as a Sum of Know Laplace Transforms
Partial Fraction Decomposition
Performing a partial fraction decomposition is not needed here, because \(H(s)=\frac{1}{s^2+2s+17}\) is already in partial fraction decomposition form.
But we do still need to bring the solution into a suitable form to find the inverse Laplace transform in a table. Recall from the Laplace transform lesson that for the damped sine and cosine functions, respectively
$$\mathcal{L}^{-1}\left[\frac{\omega}{(s+\sigma)^2+\omega^2}\right] = e^{-\sigma t}\sin(\omega t) $$
and
$$\mathcal{L}^{-1}\left[\frac{s+\sigma}{(s+\sigma)^2+\omega^2}\right] = e^{-\sigma t}\cos(\omega t). $$
In our current case, we do not have an \(s\) in the numerator, so we see that we are dealing with the sine version. However, we do need to bring the denominator into a matching form. This is done by completing the square.
Completing the Square
The technique of completing the square is rather straightforward. We take the denominator and try to find values for \(h\) and \(k\), such that
$$s^2+2s+17 = (s+h)^2 + k$$
\(h\) is chosen to be half the middle coefficient on the left-hand side. In this case, this means \(h=1\). Then
$$s^2+2s+17 = (s+h)^2 + k = (s+1)^2 + k $$
\(k\) must be the last term on the left-hand side, minus the last square term of the parenthesis (the one without an \(s\)), i.e. \(k=17-1=16\), yielding
$$s^2+2s+17 = (s+1)^2 + 16 = (s+1)^2 + 4^2 $$
Thus
$$ X(s) = \frac{1}{(s^2+2s+17} = \frac{1}{(s+1)^2 + 16} = \frac{1}{(s+1)^2 + 4^2} = \frac{1}{4}\frac{4}{(s+1)^2 + 4^2} = \frac{1}{4}\frac{\omega}{(s+\sigma)^2 + \omega^2} $$
with \(\sigma=1\) and \(\omega=4\). Notice that we had to determine the prefactor of the whole term as well. Now \(X(s)\) is in a form, which we can find in a table of functions with known Laplace transforms.
Finding the Solution \(x(t)\) in the Time Domain
The inverse Laplace Transform now yields the solution of the differential equation in the time domain
$$ x(t) = \frac{1}{4} e^{-\sigma t} \sin(\omega t) = \frac{1}{4} e^{- 1 t} \sin(4 t) $$
for a unit impulse input \(f(t)=\delta(t)\). We can simply look it up in a table of known Laplace transforms, because we have prepared \(X(s)\) in a suitable form beforehand.
System Response \(x(t)\) to Unit Step \(f(t)=\theta(t)\)
Finding the analytic solution \(x(t)\) for a unit step input \(f(t)=\theta(t)\), where \(\theta\) is the Heaviside step function, works analogously, but is slightly computationally more intensive. We will reuse insights from above and make modifications, where needed.
Finding the Solution \(X(s)\) in the Laplace Domain
First we recall that the Laplace transform of \(f(t)=\theta(t)\) is not equal to 1 this time, but rather \( F(s) = 1/s \). This leads to a somewhat different solution \(X(s)\) of the equation in the Laplace domain,
$$ X(s) = H(s)F(s)= \frac{1}{s(s^2+2s+17)}. $$
Writing \(X(s)\) as a Sum of Know Laplace Transforms and Performing the Transform
Partial Fraction Decomposition
The denominator factors of different solution now: the previous expression without real solution and the solution \(s=0\). Therefore, we must use partial fraction decomposition, before we can complete the square and use the inverse Laplace transform table.
The partial fraction decomposition procedure for this case is
$$ \frac{1}{s(s^2+2s+17)} = \frac{A}{s} + \frac{Bs+C}{s^2+2s+17} $$
Making the denominator the same for all terms yields
$$ \frac{1}{s(s^2+2s+17)} = \frac{A(s^2+2s+17)}{s(s^2+2s+17)} + \frac{Bs^2+Cs}{s(s^2+2s+17)} $$
after which we can forget about the denominator and just equate the numerators:
$$ 1 = A(s^2+2s+17) + Bs^2 + Cs $$
Because we have three variables \(A\), \(B\) and \(C\) to determine, we will need three equations. We decide to obtain them by setting \(s=0\), \(s=1\), and \(s=-1\).
- \(s=0\): \(1 = 17 A\), \(A=1/17\).
- \(s=1\): \(1 = 20 A + B + C\), \(1 = 20/17 + B + C\), \(-3/17 = B + C\), \(B = – 3/17 -C\).
- \(s=-1\): \(1 = 16 A + B – C\), \(1 = 16/17 + B – C \), \(1/17 = B – C \), \(1/17 = -3/17 – C\), \(C = -2/17\), \(B = – 1/17\).
With these values, we obtain
\begin{eqnarray} X(s)&=& \frac{A}{s} + \frac{Bs+C}{s^2+2s+17} = \frac{A}{s} + \frac{Bs}{s^2+2s+17} + \frac{C}{s^2+2s+17} =\\ &=& \frac{1}{17}\frac{1}{s} – \frac{1}{17}\frac{s}{s^2+2s+17} – \frac{2}{17}\frac{1}{s^2+2s+17} \end{eqnarray}
We wish to apply the inverse Laplace transform to all three terms. The first one is recognized to yield the Heaviside step function \(\theta(t)\), while the second term is recognized as a damped cosine, and the third term as a damped sine. However, the last two terms are again not quite yet in the proper form, to be able to apply the Laplace transform tables. Just like in the solution which was the response to an impulse input, we will here also have to complete the square.
Completing the Square
Term 1:
Nothing to do here, other than taking the inverse Laplace transform, since the inverse Laplace transform of \(\frac{1}{s}\) is tabulated and is the Heaviside step function \(\theta(t)\), which is zero for \(t<0\) and equal to 1 for \(t\ge0\). Therefore the part of the solution from the first term is \(x_1(t) = \frac{1}{17} \theta(t)\).
Term 2:
We know that Term 2 will be a damped cosine, because of the polynomial in the denominator does not factor into real solutions, and because the numerator has an \(s\). This is in contrast to the third term, which will be a damped sine, because the numerator does not have an \(s\). Sine and cosine terms often appear together, and it is important in such cases to start working on the cosine first, and process the sine only after.
For the cosine term, we want to determine \(k\), \(\sigma\) and \(\omega\), such that:
$$ -\frac{1}{17}\frac{s}{s^2+2s+17} = k \frac{s+\sigma}{(s+\sigma)^2+\omega^2}$$
We start by completing the square in the denominator. \(\sigma\) needs to be half the coefficient in front of the middle term with the \(s\) in the denominator. Since that coefficient is 2, we set \(\sigma=1\).
$$ s^2+2s+17 = (s+\sigma)^2+\omega^2 = (s+1)^2 – 1 + 17 = (s+1)^2 + 16 = (s+1)^2 + 4^2 $$
where in the third expression we introduced the -1 to undo the \(1^2\) from the preceeding bracket, which was not there in the original expression. From the last expression we see that \(\omega = 4\).
The numerator is also a problem. We must have \(s+\sigma = s+1\) in the numerator, but currently only have \(s\). We can reshape the expression to get:
$$ -\frac{1}{17} \frac{s}{s^2+2s+17} = -\frac{1}{17} \frac{s+1}{(s+1)^2 + 4^2} + \frac{1}{17} \frac{1}{(s+1)^2 + 4^2} $$
The second term in the last expression is a leftover, which we will carry into the computation of the third term. The first term in the last expression becomes a damped cosine after applying the inverse Laplace transform and yields the contribution \(x_2(t) = -\frac{1}{17} e^{-1t}\sin(4t)\) to the final solution.
Term 3:
We are dealing not only with the third term here, but also with the leftover from Term 2. They are of same shape and can be merged:
\begin{eqnarray} & &-\frac{2}{17}\frac{1}{s^2+2s+17} + \frac{1}{17} \frac{1}{(s+1)^2 + 4^2} = -\frac{2}{17}\frac{1}{(s+1)^2 + 4^2} + \frac{1}{17} \frac{1}{(s+1)^2 + 4^2} =\\ & & = – \frac{1}{17} \frac{1}{(s+1)^2 + 4^2} \end{eqnarray}
So far so good, but notice that the sine has an \(\omega\) in the numerator:
$$ – \frac{1}{17} \frac{1}{(s+1)^2 + 4^2} = k \frac{\omega}{(s+\sigma)^2 + \omega^2} $$
The prefactor \(k\) will therefore not be \(-1/17\), but rather \(-1/(17\omega)\) to compensate for the necessary \(\omega=4\) that must appear in the numerator.
$$ – \frac{1}{17} \frac{1}{(s+1)^2 + 4^2} = – \frac{1}{17\cdot 4} \frac{4}{(s+1)^2 + 4^2} = – \frac{1}{68} \frac{4}{(s+1)^2 + 4^2} $$
Inverse Laplace transform yields \(x_3(t) = – \frac{1}{68} e^{-1t}\sin(4t)\) for the third term.
Finding the Solution \(x(t)\) in the Time Domain
The solution \(x(t)\) of the problem in the time domain is the sum over the individual three contributions computed above:
$$ x(t) = x_1(t) + x_2(t) + x_3(t) = \frac{1}{17} \theta(t) – \frac{1}{17} e^{-1t}\cos(4t) – \frac{1}{68} e^{-1t}\sin(4t). $$
Plotting the Solutions
Now that we have found the solutions analytically, we can visualize them by plotting these functions with the code below.
import numpy as np
import matplotlib.pyplot as pltt_impulse = np.linspace(0, 5, 1001)
x_impulse = 1.0/4.0*np.exp(-1.0*t_impulse)*np.sin(4.0*t_impulse)
t_step = np.linspace(0, 5, 1001)
x_step = 1/17.0*1.0 - 1.0/17.0*np.exp(-1.0*t_step)*np.cos(4.0*t_step) - 1/68.0*np.exp(-1.0*t_step)*np.sin(4.0*t_step)
plt.plot(t_impulse, x_impulse, label="Unit Impulse Response")
plt.plot(t_step, x_step, label="Unit Step Response")
plt.xlabel("t [seconds]")
plt.ylabel("x(t)")
plt.legend()
plt.show()Analytic solutions are sometimes desirable and preferable over numerical ones. But it was a lot of work to obtain them. Many times in practice, would it not be much better, if we could just type in, say, the transfer function, which we found easily at the very beginning of this lesson, and then with a couple of lines of code, let the computer find the same solution numerically for us and create the same plot? This – and a lot of associated system analysis to gain further insights – is what we will learn in the rest of the Control Theory Part of this course.
Just to demonstrate that obtaining the numerical solutions is not difficult, we present the code below (it will be explained later).
import numpy as np
import matplotlib.pyplot as plt
from scipy import signalsys = signal.TransferFunction([1], [1, 2, 17])
t_impulse, x_impulse = signal.impulse(sys)
t_step, x_step = signal.step(sys)
plt.plot(t_impulse, x_impulse, label="Unit Impulse Response")
plt.plot(t_step, x_step, label="Unit Step Response")
plt.legend()
plt.show()