P.5 Review of Ordinary Differential Equations

This appendix provides a refresher on how to solve linear ordinary differential equations with constant coefficients. In particular, we focus on the step response of such systems.

P.5.1 First-Order ODEs

This section reviews how to use an integrating factor to find the solution to a scalar first-order differential equation.

Consider the first order linear constant coefficient ordinary differential equation \[\tag{P.5.1} \frac{dy}{dt}(t) = ay(t) + bu(t),\] with initial condition \(y(0)=y_0\), where \(a\) and \(b\) are constants. Equation P.5.1 can be rearranged as \[\tag{P.5.2} \frac{dy}{dt}(t) - ay(t) = bu(t).\] Multiplying both sides of Equation P.5.2 by the integrating factor \(e^{-at}\) gives \[\tag{P.5.3} e^{-at}\left(\frac{dy}{dt}(t) - ay(t)\right) = e^{-at}bu(t).\] Noting that the left-hand side of Equation P.5.3 is the derivative of \(e^{-at}y(t)\) gives \[\tag{P.5.4} \frac{d}{dt}\left(e^{-at}y(t)\right) = e^{-at}bu(t).\] Integrating both sides of Equation P.5.4 from \(0\) to \(t\) yields \[\int_0^t \frac{d}{d\tau}\left(e^{-a\tau}y(\tau)\right)d\tau = \int_0^te^{-a\tau}bu(\tau)d\tau,\] which implies that \[e^{-at}y(t) - e^{-a\cdot 0}y(0) = \int_0^te^{-a\tau}bu(\tau)d\tau.\] Solving for \(y(t)\) gives \[\tag{P.5.5} y(t) = \underbrace{e^{at}y_0}_{\text{homogeneous response}} + \underbrace{\int_0^t e^{a(t-\tau)}bu(\tau)d\tau}_{\text{forced response}},\] where the first term on the right is the homogeneous response, or the response to initial conditions, and the second term on the right is the forced response, or the response to the input \(u(t)\). Note that the forced response is the convolution of \(u(t)\) with the impulse response \(e^{at}\).

For the discussion to follow, it is important to note that if \(a\) is a complex number, then Equation P.5.5 still holds, but \(y(t)\) will be a complex function of time. To see this, assume that \(a=\sigma+j\omega\). Then \[\begin{aligned} e^{at} &= e^{(\sigma+j\omega)t} \\ &= e^{\sigma t} e^{j\omega t} \\ &= e^{\sigma t}(\cos(\omega t) + j\sin(\omega t)) \\ &= e^{\sigma t}\cos(\omega t) + je^{\sigma t}\sin(\omega t), \end{aligned}\] substituting into Equation P.5.5 yields the complex function \[\begin{gathered} y(t) = \left(e^{\sigma t}\cos(\omega t) y_0 + \int_0^t e^{\sigma (t-\tau)} \cos(\omega \tau) bu(\tau)d\tau\right) \\ + j\left( e^{\sigma t}\sin(\omega t) y_0 + \int_0^t e^{\sigma (t-\tau)} \sin(\omega (t-\tau)) bu(\tau)d\tau \right). \end{gathered}\] Note that if the real part of \(a\) is negative, or in other words if \(a\) is in the left half of the complex plane, then the homogeneous response decays to zero and the forced response is bounded if \(u(t)\) is bounded.

P.5.2 Inverse Laplace Transform

Taking the Laplace transform of Equation P.5.1 gives \[(sY(s) - y_0) = aY(s) + bU(s).\] Solving for \(Y(s)\) yields \[Y(s) = \left(\frac{1}{s-a}\right)y_0 + \left(\frac{b}{s-a}\right)U(s).\] Taking the inverse Laplace transform gives \[\tag{P.5.6} y(t) = \mathcal{L}^{-1}\left(\frac{1}{s-a}\right) y_0 + \mathcal{L}^{-1}\left(\frac{b}{s-a}\right)*\mathcal{L}^{-1}(U(s)),\] where \(*\) denotes the convolution operator. Comparing Equation P.5.5 with Equation P.5.6 implies that \[\tag{P.5.7} \mathcal{L}^{-1}\left(\frac{b}{s-a}\right) = be^{at}, \quad t\geq 0.\] In the limit, as \(a\to 0\) we get \[\tag{P.5.8} \mathcal{L}^{-1}\left(\frac{b}{s}\right) = b, \quad t\geq 0.\] Again we note that Equations P.5.7 and P.5.8 are true even when \(a\) and \(b\) are complex numbers.

P.5.3 Partial Fraction Expansion

The inverse Laplace transform, together with partial fraction expansion, can be used to find the response of first-order and higher-order systems to step inputs. For example, consider the task of finding the step response for the first-order system \[\tag{P.5.9} \dot{y} = -ay + bu,\] where \(a\) and \(b\) are real positive constants, where \(u\) is a step input of size \(A\), i.e., \[u(t) = \begin{cases} A & t\geq 0 \\ 0 & \text{~otherwise} \end{cases}.\] The Laplace transform of \(u(t)\) is \(U(s) = \frac{A}{s}\). Similarly, taking the Laplace transform of Equation P.5.9 and solving for \(Y(s)\) gives \[\tag{P.5.10} Y(s) = \left(\frac{b}{s+a}\right)U(s) = \frac{bA}{s(s+a)}.\] The partial fraction expansion of \(Y(s)\) expresses \(Y(s)\) as a linear combination of each factor in the denominator, i.e., \[\tag{P.5.11} Y(s) = \frac{c_1}{s} + \frac{c_2}{s+a}.\] To find the coefficients \(c_1\) and \(c_2\), set Equation P.5.11 equal to Equation P.5.10 to obtain \[\tag{P.5.12} \frac{c_1}{s} + \frac{c_2}{s+a} = \frac{bA}{s(s+a)}.\] To find \(c_1\), multiply both sides of Equation P.5.12 by \(s\) to obtain \[\tag{P.5.13} c_1 + \frac{s c_2}{s+a} = \frac{bA}{s+a}.\] To remove the term involving \(c_2\) in Equation P.5.13, set \(s=0\) to obtain \[c_1 = \frac{bA}{a}.\] To find \(c_2\), multiply both sides of Equation P.5.12 by \(s+a\) to obtain \[\tag{P.5.14} \frac{c_1(s+a)}{s} + c_2 = \frac{bA}{s}.\] To remove the term involving \(c_1\) in Equation P.5.14, set \(s=-a\) to obtain \[c_2 = -\frac{bA}{a}.\] Therefore, \[\tag{P.5.15} Y(s) = \frac{\frac{bA}{a}}{s} - \frac{\frac{bA}{a}}{s+a}.\] Since the inverse Laplace transform is a linear operator, we get \[\begin{aligned} y(t) &= \mathcal{L}^{-1}\left(\frac{\frac{bA}{a}}{s}\right) - \mathcal{L}^{-1}\left(\frac{\frac{bA}{a}}{s+a}\right) \\ &= \frac{bA}{a} - \frac{bA}{a}e^{-at}, \qquad t\geq 0 \\ &= \frac{bA}{a}(1-e^{-at}), \qquad t\geq 0. \end{aligned}\]

The same technique works for finding the step response for higher-order systems. As an example, find the unit step response of the system described by \[y^{(4)} + 10y^{(3)} + 35\ddot{y} + 50\dot{y} + 24 y = u.\] Taking the Laplace transform and solving for \(Y(s)\), and using the fact that \(U(s)=\frac{1}{s}\) gives \begin{align} Y(s) &= \left(\frac{1}{s^4+10s^3+35s^2+50s+24}\right)U(s) \notag \\ &= \left(\frac{1}{s^4+10s^3+35s^2+50s+24}\right)\frac{1}{s} \notag \\ & = \frac{1}{s}\left(\frac{1}{s+1}\right)\left(\frac{1}{s+2}\right)\left(\frac{1}{s+3}\right)\left(\frac{1}{s+4}\right). \tag{P.5.16} \end{align} Using partial fraction expansion we have that \[\tag{P.5.17} Y(s) = \frac{c_1}{s} + \frac{c_2}{s+1} + \frac{c_3}{s+2} + \frac{c_4}{s+3} + \frac{c_5}{s+4}.\] Setting Equation P.5.17 equal to Equation P.5.16 gives \begin{gather} \tag{P.5.18} \frac{c_1}{s} + \frac{c_2}{s+1} + \frac{c_3}{s+2} + \frac{c_4}{s+3} + \frac{c_5}{s+4} \\ = \frac{1}{s}\left(\frac{1}{s+1}\right)\left(\frac{1}{s+2}\right)\left(\frac{1}{s+3}\right)\left(\frac{1}{s+4}\right). \end{gather} The coefficient \(c_1\) is found by multiplying both sides of Equation P.5.18 by \(s\) and then setting \(s=0\) to remove the coefficients \(c_2\)–\(c_5\) to obtain \(c_1 = 1/24\). Following a similar procedure for each coefficient yields \(c_2=-1/6\), \(c_3=1/4\), \(c_4=-1/6\), and \(c_5=1/24\). The step response is therefore given by \[\begin{aligned} y(t) &= \mathcal{L}^{-1}\left(\frac{1/24}{s}\right) + \mathcal{L}^{-1}\left(\frac{-1/6}{s+1}\right) + \mathcal{L}^{-1}\left(\frac{1/4}{s+2}\right) \\ &\qquad + \mathcal{L}^{-1}\left(\frac{-1/6}{s+3}\right) + \mathcal{L}^{-1}\left(\frac{1/24}{s+4}\right) \\ &= \frac{1}{24} -\frac{1}{6}e^{-t} + \frac{1}{4}e^{-2t} - \frac{1}{6}e^{-3t} + \frac{1}{24}e^{-4t}, \qquad t\geq 0 \\ &= \frac{1}{24}\left(1 -4e^{-t} + 6e^{-2t} - 4e^{-3t} + e^{-4t}\right), \qquad t\geq 0. \end{aligned}\]

In MATLAB, the coefficients can be found by using the residue command. For the problem above we could use the matlab commands

num = 1;  %
%
den = poly([0,-1,-2,-3,-4]); 
%
[R,P,K] = residue(num,den);  

The vector R will contain the coefficients \(c_1\)–\(c_5\), the vector P contains the associated roots, and K will be the remainder polynomial if the numerator is a higher order polynomial than the denominator.

The method of partial fraction expansion is also valid if the poles of the system are complex. For example, consider finding the response of the second-order differential equation \[\ddot{y} + 2\zeta\omega_n\dot{y} + \omega_n^2y = \omega_n^2u,\] to a step of size \(A\). Taking the Laplace transform and letting \(U(s)=A/s\) gives \begin{align} Y(s) &= \left(\frac{\omega_n^2}{s^2+2\zeta\omega_ns+2\omega_n^2}\right)\frac{A}{s} \notag \\ & = \frac{A\omega_n^2}{s}\left(\frac{1}{s+\zeta\omega_n+j\omega_n\sqrt{1-\zeta^2}}\right) \\ &\qquad \cdot \left(\frac{1}{s+\zeta\omega_n-j\omega_n\sqrt{1-\zeta^2}}\right). \tag{P.5.20} \end{align} Expanding \(Y(s)\) as a sum of partial fractions yields \begin{gather} \tag{P.5.21} Y(s) = \frac{c_1}{s} + \left(\frac{c_2}{s+\zeta\omega_n+j\omega_n\sqrt{1-\zeta^2}}\right) \\ + \left(\frac{c_3}{s+\zeta\omega_n-j\omega_n\sqrt{1-\zeta^2}}\right). \end{gather} Setting Equation P.5.21 equal to Equation P.5.20 gives \[\begin{gathered} \frac{c_1}{s} + \left(\frac{c_2}{s+\zeta\omega_n+j\omega_n\sqrt{1-\zeta^2}}\right) + \left(\frac{c_3}{s+\zeta\omega_n-j\omega_n\sqrt{1-\zeta^2}}\right) \\ = \frac{A\omega_n^2}{s}\left(\frac{1}{s+\zeta\omega_n+j\omega_n\sqrt{1-\zeta^2}}\right)\left(\frac{1}{s+\zeta\omega_n-j\omega_n\sqrt{1-\zeta^2}}\right). \end{gathered}\] Solving for \(c_1\), \(c_2\), and \(c_3\) gives \[\begin{aligned} c_1 &= A \\ c_2 &= -\frac{A}{2}\left(1+j\frac{\zeta}{\sqrt{1+\zeta^2}}\right)\\ c_3 &= -\frac{A}{2}\left(1-j\frac{\zeta}{\sqrt{1-\zeta^2}}\right). \end{aligned}\] Therefore, \[Y(s) = \frac{A}{s} + \frac{-\frac{A}{2}\left(1+j\frac{\zeta}{\sqrt{1+\zeta^2}}\right)}{s+\zeta\omega_n+j\omega_n\sqrt{1-\zeta^2}} + \frac{-\frac{A}{2}\left(1-j\frac{\zeta}{\sqrt{1-\zeta^2}}\right)}{s+\zeta\omega_n-j\omega_n\sqrt{1-\zeta^2}}.\] Taking the inverse Laplace transform gives \[\begin{aligned} y(t) &= A - \frac{A}{2}\left(1+j\frac{\zeta}{\sqrt{1+\zeta^2}}\right)e^{-\zeta\omega_nt-j\omega_n\sqrt{1-\zeta^2}t} \\ &\qquad\qquad - \frac{A}{2}\left(1-j\frac{\zeta}{\sqrt{1-\zeta^2}}\right)e^{-\zeta\omega_nt+j\omega_n\sqrt{1-\zeta^2}t} \\ &= A\Bigg(1 - \frac{1}{2}\left|1+j\frac{\zeta}{\sqrt{1+\zeta^2}}\right|e^{-\zeta\omega_nt} \\ &\qquad\qquad \cdot \left(e^{-j\omega_n\sqrt{1-\zeta^2}t}e^{j\angle\left(1+j\frac{\zeta}{\sqrt{1+\zeta^2}}\right)} + e^{+j\omega_n\sqrt{1-\zeta^2}t}e^{j\angle\left(1-j\frac{\zeta}{\sqrt{1+\zeta^2}}\right)} \right)\Bigg) \\ &= A\Bigg(1 - \frac{1}{\sqrt{1-\zeta^2}}e^{-\zeta\omega_nt} \\ &\qquad\qquad \cdot \left(\frac{e^{-j\left(\omega_n\sqrt{1-\zeta^2}t-\tan^{-1}\left(\frac{\zeta}{\sqrt{1-\zeta^2}}\right)\right)} + e^{j\left(\omega_n\sqrt{1-\zeta^2}t-\tan^{-1}\left(\frac{\zeta}{\sqrt{1-\zeta^2}}\right)\right)}}{2} \right)\Bigg) \\ &= A\left(1 - \frac{1}{\sqrt{1-\zeta^2}}e^{-\zeta\omega_nt} \cos\left(\omega_n\sqrt{1-\zeta^2}t-\tan^{-1}\left(\frac{\zeta}{\sqrt{1-\zeta^2}}\right)\right)\right), \end{aligned}\] for \(t\geq 0\). In general, complex conjugate poles will always lead to complex conjugate coefficients that end up canceling to produce a decaying exponential multiplied by a cosine term with a phase shift. If the coefficients of the original ordinary differential equation are real, then the response will always be real.

second-order