Linear Differential Equations

\(
\newcommand{\BE}{\begin{equation}}
\newcommand{\EE}{\end{equation}}
\newcommand{\BA}{\begin{eqnarray}}
\newcommand{\EA}{\end{eqnarray}}
\newcommand\CC{\mathbb{C}}
\newcommand\FF{\mathbb{F}}
\newcommand\NN{\mathbb{N}}
\newcommand\QQ{\mathbb{Q}}
\newcommand\RR{\mathbb{R}}
\newcommand\ZZ{\mathbb{Z}}
\newcommand{\va}{\hat{\mathbf{a}}}
\newcommand{\vb}{\hat{\mathbf{b}}}
\newcommand{\vn}{\hat{\mathbf{n}}}
\newcommand{\vt}{\hat{\mathbf{t}}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\bv}{\mathbf{v}}
\newcommand{\bg}{\mathbf{g}}
\newcommand{\bn}{\mathbf{n}}
\newcommand{\by}{\mathbf{y}}
\)

Linear Differential Equations

Definitions and Terminology

Differential Equations

Ordinary differential equations (ODE)  are equations where the variable we want to solve for is a function of another variable and derivatives of the function with respect to that variable appear.

Let \(x\) be a function of some variable \(t\), and let the dots over the function denote the derivatives with respect to that variable \(t\). Then the equation

\begin{equation}
\ddot x + x = 0
\end{equation}

is a differential equation. Its solution is a function, where the second derivative of the function plus the function itself equals zero, for all values of \(t\).

Such equations arise easily in physics. Imagine a cart on wheels at position x=0, which is attached to a spring with spring constant \(k\). The spring has unstretched length 10cm and is attached to a wall at -10cm in this coordinate system. So there is no force acting on the cart, when it is a x=0.

Now imagine the cart is displaced to the right, and the spring stretches. It will exert a force to the left, which is proportional to the spring constant \(k\) as well as the displacement of the cart, \(x\). The well known equation of motion \(F = ma\) then states that \(-kx = m\ddot x\). Reordering terms yields \(m\ddot x+kx = 0\) , which is exactly the form of the equation above, with additional constant coefficients.

In contrast to ordinary differential equations, which depend only on one variable, partial differential equations (PDE) can depend on several. They contain partial derivatives of multivariable functions with respect to several different variables. They can be found in many areas of physics and engineering, for instance in fluid dynamics and electrodynamics. However, we shall not concern ourselves with them here and rather just focus on ODEs instead.

Order

The order of a differential equation is equal to the highest derivative of the function that appears in the equation. The above example is a second-order differential equation of function \(x(t)\), because the second derivative of \(x\) with respect to \(t\) appears in the equation.

Initial Conditions

The solutions of differential equations by themselves are not unique. In order to get a unique solution, one needs to specify initial conditions. Each additional order requires an extra additional condition. For a first-order differential equation, knowing the value of the function (\x\) at some initial value of the parameter \(t\) is sufficient. This is corresponds to knowing the initial position of the cart above at some initial time \(t_0\).

If it is a second-order differential equation, one must also must know the initial value of the  of the first derivative of the function. In the above cart example, this corresponds to knowing the initial velocity at initial time \(t_0\).

With the cart example, we are done with this, because it is a second-order differential equation. But for higher-order differential equations, we must also know the initial values of higher derivatives to specify a unique solution.

Linear Differential Equations

Linear differential equations are equations, where only multiples of the function and its derivatives appear, i.e. no square terms \(x^2(t)\) or other functional dependencies on \(x(t)\), such as terms like \(e^{x(t)}\) or \(\mathrm{log}(x(t))\). The generic second-order linear differential equation thus takes the form:

\begin{equation}
a(t)\ddot x(t) + b(t)\dot x(t) + c(t) x(t) = f(t) 
\end{equation}

In physics application, the forcing function \(f(t)\) can be thought of as a time-dependent external force independent of the object’s position (unlike the spring force \(-kx(t)\) in the above cart example, which depended on the position of the cart, \(x(t)\)). If \(f(t)=0\) for all \(t\), the differential equation is called homogeneous, and nonhomogeneous otherwise. The coefficients \(a\), \(b\), and \(c\) can generally depend on parameter \(t\), but are often constant in many applications, in which case it is a differential equation with constant coefficients. We shall restrict our discussion to this latter simplifying case.

We will only need two types of linear differential equations in our courses: first-order and second-order linear differential equations. We shall therefore primarily focus only on these two types in this article. However, we will touch upon higher-order differential equations briefly.

Solutions of Linear Differential Equations

First-Order Linear Differential Equations with Constant Coefficients

Differential equations are solved the easiest with an ansatz, i.e. a good guess of the functional form of the solution, with parameters yet to be determined. Let us start with the generic, homogeneous, first-order linear differential equation with constant coefficients

\begin{equation}
a\dot x + bx = 0,
\end{equation}

where it is understood that \(x\) is a function of parameter \(t\).

Let us guess that the solution will be of the form \(x(t) = A e^{\alpha t}\), with parameters \(A\) and \(\alpha\) yet to be determined. Inserting this ansatz into the equation yields:
\begin{equation}
aA\alpha e^{\alpha t} + bAe^{\alpha t} = 0
\end{equation}

After dividing both sides by \(Ae^{\alpha t}\), we obtain \(a\alpha+b = 0\), which can be solved for \(\alpha\):

\begin{equation}
\alpha = -\frac{b}{a}
\end{equation}

If both \(a\) and \(b\) have the same sign, then \(\alpha\) is negative and the solution \(x(t)=A e^{\alpha t}\) is a exponentially decaying function with time. If \(a\) and \(b\) have different signs, the solution will be exponentially growing. This is an important take-home note for linear first-order differential equations  with constant, real coefficients.

The constant \(A\) can be determined from the initial condition, which is given by specifying that \(x(t)\) shall have the value \(x_0\) for some parameter value \(t_0\), i.e. by specifying \(x(t_0)=x_0\). The parameter \(t\) oftentimes denotes time, which is why this is called an initial condition. In the above example of the cart attached to a spring, we need to know the position \(x_0\) of the cart at some initial time \(t_0\), otherwise we cannot say at all where the cart will be later (in this particular example we will additionally need to know the initial velocity as well, because the differential equation is second order, as we shall see in the next subsection).

Second-Order Linear Differential Equations with Constant Coefficients

The generic, homogeneous second order differential equation with constant coefficients can be written as

\begin{equation}
a\ddot x+b\dot x + cx = 0.
\end{equation}

If we enter this differential equation with the same ansatz as before, we obtain (again after dividing by \(Ae^{\alpha t}\)):

\begin{equation}
a\alpha^2+b\alpha+c = 0
\end{equation}

This is a quadratic equation in \(\alpha\) and is called the characteristic equation of the differential equation (its lefthand side alone being called the characteristic polynomial). This equation is an algebraic equation for the variable \(\alpha\), not containing any derivatives of any functions anymore. It has the two solutions (the characteristic polynomial has two roots, which in some special cases might coincide to be the same number):

\begin{equation}
\alpha = \frac{-b\pm\sqrt{b^2-4ac}}{2a}
\end{equation}

These solutions can either be two real numbers, one single real number, or two complex numbers, which are complex conjugates of each other. The form of the solution of the differential equation will depend on this. (If you happen to be unfamiliar with complex numbers, please follow the above links to Wikipedia articles to catch up on this important topic, before reading on.)

It is easily observed, because the differential equation is linear, that if two functions are solutions of the differential equation, then an arbitrary linear combination of the two is a solution as well. Therefore, the generic solution takes in this case the form:

\begin{equation}
x(t) = A e^{\alpha_1 t} + B e^{\alpha_2 t}
\end{equation}

The time evolution of this function depends on \(\alpha_1\) and \(\alpha_2\). There are three distinct cases:

  • If both \(\alpha_1\) and \(\alpha_2\) are real numbers, the solution is either exponentially decaying (if both are smaller than zero) or exponentially growing (if at least one of them is larger than zero).
  • If \(\alpha_1=\delta + i \omega\) is complex with real part \(\delta=-b/2a\) and imaginary part \(\omega=\sqrt{4ac-b^2}/2a\), and thus \(\alpha_2=\delta – i \omega\), then the solution has the form
    \begin{equation}
    x(t) = A e^{(\delta+i \omega) t} + B e^{(\delta -i \omega)t} = e^{\delta t} (A   e^{(i \omega) t} + B e^{(-i \omega)t}) = e^{\delta t} (C \sin(\omega t) + D \cos(\omega t)),
    \end{equation}
    where in the last step we have used Euler’s formula. This solution is oscillatory with natural frequency \(\omega\) and exponential decay (damping) or growth \(\delta\), depending on the latter parameter being negative or positive, respectively.
  • If \(\alpha_1=\alpha_2=-b/2a\), both being real, then the solution takes the form
    \begin{equation}
    x(t) = (A + Bt)e^{\alpha t}
    \end{equation}
    This is an exponentially growing or decaying solution, depending on whether \(\alpha\) is positive or negative (the exponential always wins over the polynomial term).

Which of these three cases occurs, depends of the discriminant \(\mathcal{D}:=b^2-4ac\). If it is positive, the roots of the characteristic polynomial are both real and different (case 1), if it is zero, the roots are real and the same (case 3), and if \(\mathcal{D}\) is negative, the roots are complex and the solution of the differential equation is oscillatory. We see that in the oscillatory case, the coefficient \(b\) is responsible for the damping, i.e. the coefficient in front of the first-order derivative.

Higher-Order Linear Differential Equations with Constant Coefficients

For an linear differential equation of order \(m\), the characteristic polynomial is an polynomial of degree \(m\), and the above result can be extended analogously. The roots \(\alpha\) of the characteristic polynomial, i.e. the solutions of the characteristic equation, are either real or pairs of complex conjugate numbers. The basis functions, from which a linear combination can be taken to construct a solution, consist of functions of the form e^{\alpha t}, just as we had before. And analogously to before, if a particular root \(\alpha\) occurs \(n\) times (i.e. with multiplicity \(n\)), then the basis functions for that value of \(\alpha\) are 

\begin{equation}
t^{j}e^{\alpha t}
\end{equation}

with \(j=0,\dots,n-1\). For instancee, a root \(\alpha\) occurring four times in the characteristic polynomial would result in basis functions

\begin{equation}
e^{\alpha t}, \hspace{0.5cm} te^{\alpha t}, \hspace{0.5cm}t^2e^{\alpha t}, \hspace{0.5cm} t^3e^{\alpha t}
\end{equation}

for the solution. If the \(\alpha\) is real, the basis function will be exponential, if \(\alpha\) is complex, the basis function will be oscillatory, just as we had for the second-order linear differential equation. This allows the solution of higher-order differential equations to exhibit a variety of oscillatory behaviors as a superposition of sines and cosines of different frequencies (unlike the second-order linear differential equation, which allows only one frequency). This can give the solutions a variety of shapes (as you may know from Fourier transforms, by adding enough sines and cosines with different frequencies, you can obtain essentially an arbitrarily shaped function). 

We shall not dwell on this much here, and would like to refer the interested reader to the Wikipedia article on linear differential equations instead.

Based on the above observations, we shall call a exponentially decaying or increasing solution to be a “first-order response”, and a (damped or diverging) oscillatory solution (with one frequency) to be a “second-order response”, even though the originating system may be a higher-order system.

Harmonic Oscillator

The example with the cart on a spring, which we have encountered at the beginning, is an example of a second-order linear differential equation with constant coefficients, just as discussed above. It is called a simple harmonic oscillator (in the absence of friction). The cart will move back and forth in a sinusoidal motion.

The amplitude will decrease if damping is present due to friction (damped harmonic oscillator). The damping depends on parameter \(b\) in our above generic notation and is thus given by the factor in front of the first-order derivative terms (which was absent in our original formulation of the cart example).

Let us revisit this in a bit more detail. For a damped harmonic oscillator we have in general

\begin{equation}
m\ddot x + c \dot x + kx = 0 
\end{equation} 

where \(m\) is the mass, \(c\) is the damping coefficient, and \(k\) is the spring constant. The friction force has been modeled above by first-order term \(c\dot x\), which makes the friction force proportional to velocity. This is often a good model for many vibrating systems. However, it has its limitations and one needs to be aware that if the cart was sliding on a surface, friction force would depend only on the direction of motion (and the carts normal force pressing it again the surface) and not the cart’s velocity (Coulomb damping). And aerodynamic drag force (air resistance) goes as velocity squared. Nonetheless, we shall discuss the damped harmonic oscillator as written above with friction force modeled proportional to velocity.

The critical damping coefficient for the above equation of motion is defined as:

\begin{equation}
c_c: = \sqrt{km}=2m\omega_n .
\end{equation} 

where in the last expression we have written the critical damping in terms of the natural frequency \(\omega_n:=\sqrt{k/m}\) of the system, i.e. its frequency of oscillation in the absence of any damping.

For the actual damping \(c\) of the system, four distinct cases exist, resulting in different behavior over time:

  • A system with no damping at all, \(c=0\) is called undamped. It continues to oscillate with the same amplitude at its natural frequency \(\omega_n = \sqrt{k/m}\) indefinitely.
  • If \(c<c_c\) the motion of the object will be oscillatory with decreasing amplitude. The system is called underdamped. The damping is weak and while it decreases the amplitude of the motion over time, the object overshoots its equilibrium position (usually several times) before slowly converging to it.
  • If \(c>c_c\), then the system will converge directly to the final equilibrium position without performing oscillations (i.e. without overshooting it). Such a system is called overdamped. The damping is so strong in this case that the object does not have enough kinetic energy to overshoot the final equilibrium. However, a heavily overdamped system may take a long time to converge to the final equilibrium, because the damping also slows down the motion towards the equilibrium in general.
  • When the system is critically damped, i.e. \(c=c_c\) the system converges to the final equilibrium the fastest without oscillating.

Written in terms of natural frequency \(\omega_n:=\sqrt{k/m}\) and damping ratio \(\zeta:=c/c_c\) the damped harmonic oscillator equation can be recast in the often used form

\begin{equation}
\ddot x + 2\zeta\omega_n\dot x + \omega_n^2 x = 0.
\end{equation}

The characteristic polynomial then has the roots (in the oscillatory case, when they are complex)

\begin{equation}
\alpha = -\zeta\omega_n \pm i \omega_n\sqrt{1-\zeta^2}.
\end{equation}

Occurrences in Our Courses

We will encounter differential equations on several occasions in our courses. The dynamic equations of motion are first-order nonlinear differential equations in the translational and angular velocities. While we will derive them, we will not solve them, so we do not have to worry about them being nonlinear. In order to study aircraft stability, we will linearize them first, assuming that the perturbations around the equilibrium condition are small. They then become a system of linear differential equations, to which the above insights can be applied.

The two dynamic equations (actually six, because each of them is a vector-valued equation containing three scalar ones) are supplemented by two first-order kinematic equations to form a twelve-equation system. However, just because the equations are first-order does not mean that the solution of the whole system is a first-order response (see next section below, where we discuss systems of differential equations). In fact, the solution is generally higher order, but oftentimes (decoupled) modes can be identified which correspond approximately to first- or second- order responses.

We will encounter linear differential equations and their solutions primarily during our study of dynamic stability of airplanes, which studies the response of the aircraft over time. We will make the following observations.

  • The short period and phugoid modes of longitudinal stability are oscillatory and therefore definitely at least second-order responses. They arise from the six linearized equation of longitudinal stability for the variables \(u, w, \alpha, \theta, x, z\) (sometimes \(u\) and \(w\) are replaced by airspeed \(V\) and flight path angle \(gamma\). The whole system is fourth-order in the first four of the variables, because \(x\) and\(z\) do not couple back in. Among the first four variables, the cross-coupling between the phugoid and the short period mode in usually small, resulting to a good approximation in two separate second-order responses for the phugoid and the short-period mode.
  • As we will see in a root locus analysis, the characteristic equation of the linearized coupled equations of motion of lateral-directional stability has four roots: two roots are complex conjugates of each other, giving rise to the second-order Dutch roll mode. The two other roots are real, without imaginary parts, and therefore show non-oscillatory, exponential (“first-order-like”) behavior. The first real root leads to the roll mode, which looks like a first-order response by virtue of heavy damping. The second one leads to the spiral mode. The spiral mode is a good example of a mode which can be either convergent (exponentially decaying) or divergent (exponentially increasing), as we have observed above by having exponentially growing or decaying solutions of first- and second-order linear differential equations. It is to be noted that for a few aircraft, the two real roots come together on the real axis and split off into a second complex pair in the complex plane, leading to a combined, oscillatory “coupled roll spiral mode” instead, which can be difficult to control for the pilot. We will discuss all of this in our test pilot course series and, in much less detail, in some of our somewhat less technical courses (e.g. in upset prevention and recovery training).

Systems of Linear Differential Equations

A proper treatment of solving systems of linear differential equations is beyond the current state of this article. But we want to make one comment, which may address an astute reader’s concern about how many different roots we alluded to in the above lateral-directional case (more than two).

Consider the following general system of two first-order homogeneous linear differential equations with constant coefficients:

\begin{eqnarray}
\dot x &=& ax + by\\
\dot y &=& cx + dy
\end{eqnarray}

Naively, one would think that the solution would be exponentially growing or decaying functions for \(x(t)\) and for \(y(t)\), and no oscillatory solutions, because we are dealing with first-order linear differential equations.

But this is not the case in general, as we show by picking a specific example of the above system. A special case of the above system of equations is if \(c=1\) and \(d=0\). Then the second equation reads simply \(\dot y = x\).

In order to solve this system, we can substitute the second equation into the first. To accomplish this, we take the time derivative of the first equation and obtain \(\ddot x = a\dot x + b\dot y\).

Eliminating \(\dot y\) in this equation by substituting in \(\dot y = x\) from the second equation obtained previously, yields the second-order differential equation for \(x(t)\):

\begin{equation}
\ddot x = a \dot x + b x.
\end{equation}

Since this is a second-order differential equation, it may very well have an oscillatory solution for \(x(t)\), depending on the values of \(a\) and \(b\). And for \(y(t)\) we may obtain an oscillatory solution as well. For instance, if \(x(t)\) turns out to be a cosine function, then \(y(t)\) would be a sine function (because its derivative is a cosine function). This is completely contrary to our naive initial expectation that the solutions for \(x(t)\) and \(y(t)\) would have to be first-order exponentially growing or decaying functions.

The key to solving this problem in a general case is to find the eigenvalues of the matrix on the righthand side. These eigenvalues later become the exponential coefficients in the solution (matrix exponentials of diagonal matrices are easy to compute). Some real matrices have complex eigenvalues, which can then lead to oscillatory behavior of the functions.

This was just a very quick illustration. More on solving systems of linear differential equations can be found on the Wikipedia pages for matrix differential equations and matrix exponentials, including some calculated examples.

Laplace Transform

The Laplace transform \(X(s)\) of a function x(t) is defined as

\begin{equation}
X(s):=L[x(t)]:=\int_0^\infty x(t)e^{-st}\,dt
\end{equation}

The Laplace transform itself is denoted by the letter L (and applied to a whole function, \(x(t)\) in this case). We shall use the uppercase letter (“X” in this case) of the original function letter (“x” in this case) to denote the Laplace transformed function, after the Laplace transform was applied. Note that the variable has changed from \(t\) to \(s\) after the transform. (This is similar in style to the frequently used Fourier transform, where the Fourier transformed function depends on the new variable \(k\) in the frequency domain. The Fourier transform uses a different exponent and turns a real function into a complex one in the frequency domain. The Laplace transform is somewhat different and in general turns a complex function into a complex function.)

What \(X(s)\) looks like depends on function \(x(t)\) and has to be determined for every new function \(x(t)\) individually (extensive lookup tables have been compiled in the mathematical and engineering literature for many common functions). However, without even worrying what \(X(s)\) might look like, we can make the following generally true statements which hold regardless of what kind of function \(x(t)\) is. If \(X(s)\) denotes \(L[x(t)]\), i.e. \(L[x(t)]=X(s)\), then we also have (after some calculation steps which we have omitted here and left as an exercise for the reader)

\begin{eqnarray}
L[cx(t)] &=& cX(s)\\
L[b\dot x(t) &=& b\left(sX(s)-x(0)\right)\\
L[a\ddot x(t)] &=& a \left(s^2X(s)-sx(0)-\dot x(0)\right)
\end{eqnarray}

The \(x(0)\) and \(\dot x(0)\) denote the initial conditions, i.e. the value of the original function \(x(t)\) at \(t=0\). They are fixed numbers. We obtain for higher derivatives

\begin{equation}
L[\frac{d^k x}{dt}(t)] = s^k X(s) – p_0(s)
\end{equation}

where \(p_0(s)\) stands for the \(k-1\)-degree polynomial in \(s\) consisting of various terms containing initial conditions \(x(0)\), \(\dot x(0)\), \(\ddot x(0)\), etc. For zero initial conditions, \(p_0(s)\) vanishes. The above expression is extremely convenient, because we realize that the Laplace transform can turn derivatives into algebraic expressions without derivatives.From the above we can now see that the Laplace transform, when applied to both sides of a whole linear differential equation, can turn the differential equation for \(x(t)\) into an algebraic equation for function \(X(s)\) (i.e. there are no derivatives of the function \(X(s)\) with respect to \(s\), only terms with factors of the function \(X(s)\) itself\). Indeed, the second-order linear differential equation

\begin{equation}
a\ddot x(t) + b\dot x(t) + cx(t) = f(t)
\end{equation}

becomes

\begin{equation}
a \left(s^2X(s)-sx(0)-\dot x(0)\right) + b\left(sX(s)-x(0)\right) + cX(s) = F(s)
\end{equation}

which can be regrouped into

\begin{equation}
(as^2+bs + c)X(s) – \underbrace{(\dot x(0) +(as+b) x(0))}_{=:p_0(s)} = F(s)
\end{equation}

and easily solved for \(X(s)\)

\begin{equation}
X(s) = \frac{F(s) + \overbrace{\dot x(0)+(as+b) x(0)}^{=p_0(s)}}{(as^2+bs + c)}
\end{equation}

Just like our ansatz approach earlier, where we guessed the functional form of the solution and were able to turn the differential equation into a characteristic equation without derivative, the Laplace transform approach gets rid of the derivatives as well, and we will encounter it in a few places in some of our courses. The advantage of this approach compared to the ansatz approach above is that it can be applied to nonhomogeneous linear differential equations as well (in general, the solution of a nonhomogeneous linear differential equation is the transient solution of the homogeneous equation plus a particular solution including the forcing function). All one needs to do next is to compute the Laplace transform \(F(s)\) of function \(f(t)\) and then use the inverse Laplace transform to transform the obtained \(X(s)\) back into an \(x(t)\).

We see that the derivative terms in the original equation result in a polynomial prefactor in front of \(X(s)\), which ends up in the denominator eventually. The polynomial \(p_0(s)\) on the other hand, ends up in the numerator. A typical approach to finding a solution is, therefore, to do a partial fraction decomposition of the fraction of two polynomials on the right hand side and then use a lookup table for the back transform of the partial fraction terms, which are simpler function that are often tabulated in the literature. A few notable Laplace back-transforms of functions are:

\begin{eqnarray}
L^{-1}[\frac{n!}{s^{n+1}}] &=& t^n\\
L^{-1}[\frac{n!}{(s+a)^{n+1}}] &=& t^ne^{-at}\\
L^{-1}[\frac{1}{(s+a)(s+b)}] &=& \frac{1}{b-a}(e^{-at}-e^{-bt})\\
L^{-1}[\frac{a}{s^2+a^2}]&=&\sin(at)\\
L^{-1}[\frac{s}{s^2+a^2}]&=&\cos(at)\\
\end{eqnarray}

for \(n=1,2,3,\dots\) and \(a\not= b\).