Linear Constant Coefficient Ordinary Differential Equations

In the most general form, an Nth order ordinary differential equation (ODE) of a single-variable function $y(x)$ can be expressed as

$$f\left(x,y(x),y'(x),y''(x),\cdots,y^{(N)}(x)\right)=0$$ (1)

which can be considered as a special case of a partial differential equation (PDE) for a multi-variable function $y(x_1,\cdots,x_N)$:

$$f\left(x_1,\cdots,x_N,y,\frac{\partial y}{\partial x_1},\cdots,\f... ...rtial x_1},\cdots,\frac{\partial^2 y}{\partial x_1\partial x_n},\cdots\right)=0$$ (2)

The purpose is to obtain a function $y(x)$ that satisfies the given differential equation and the boundary conditions

$$y(x)\vert _{x=a}=y(a)=y_a,\;\;\;\;\;\;\;y(x)\vert _{x=b}=y(b)=y_b$$ (3)

In the following we will only consider the initial value problem (IVP), by assuming the independent variable to be time, $x=t$, and the boundary condition becomes an initial condition $y(t)\vert _{t=0}=y(0)=y_0$.

• An ODE is explicit if it can be written in the following form:

  $$\frac{d^Ny(t)}{dt^N}=y^{(N)}(t) =f(t,y(t),y'(t),y''(t),\cdots,y^{(N-1)}(t))$$ (4)

  otherwise it is implicit.
• An ODE is a autonomous system if it does not explicitly depend on the independent variable $t$, i.e., it can be written in the following form:

  $$\frac{d^Ny(t)}{dt^N}=y^{(N)}(t) =f(y(t),y'(t),y''(t),\cdots,y^{(N-1)}(t))$$ (5)

• An ODE is linear if it can be expressed as:

  $$\sum_{n=0}^N a_i(t) y^{(n)}(t)=\sum_{n=0}^N a_n(t) \frac{d^ny(t)}{dt^n}=x(t)$$ (6)

  Typically this equation describes the behavior of a linear system of which the output $y(t)$ is to be determined as the system's response to a given input $x(t)$ on the right-hand side of the equation, symbolically represented as

  $$y(t)= {\cal T} [ x(t) ]$$ (7)

• A linear ODE is a homogeneous equation if the input $x(t)=0$ on the right-hand side is zero.
• A linear ODE is called a linear constant-coefficient ODE (LCCODE) if all coefficients $a_i(t)=a_i$ are constants:

  $$\sum_{n=0}^N a_i y^{(n)}(t) =\sum_{n=0}^N a_n\frac{d^ny(t)}{dt^n}=x(t)$$ (8)

  Without loss of generality, we typically assume $a_n=1$.

This N-th order LCCODE describes the behavior of a linear system in terms of how its output $y(t)$ is related to the input $x(t)$, symbolically represented by:

We now only consider solving an Nth order LCCODE in the following cases:

• Homogeneous case with $x(t)=0$:

  We assume the output takes the form of a complex exponential $y(t)=c\;e^{st}$ and substitute it together with its derivatives $y^{(n)}(t)=c s^n e^{st}$ into the equation and get

  $$\sum_{n=0}^N a_i y^{(n)}(t) =\left(\sum_{n=0}^N a_n s^n\right) c\,e^{st}=0$$ (9)

  Dividing both sides by $c\,e^{st}$, we get the characteristic equation $\sum_{n=0}^N a_n s^n=0$, which has in general $N$ roots satisfying

  $$\sum_{n=0}^N a_n s^n=\prod_{n=1}^N (s-s_n)=0$$ (10)

  The general homogeneous solution is therefore a linear combination of $N$ such solutions $c_n e^{s_nt}\;(n=1,\cdots,N)$:

  $$y_h(t)=\sum_{n=1}^N c_n e^{s_nt}$$ (11)

  and its kth-order derivatives are

  $$y^{(k)}(t)=\frac{d^k}{dt^k}\left(\sum_{n=1}^N c_n e^{s_nt}\right)... ...}{dt^k}e^{s_nt} =\sum_{n=1}^N c_n s_n^k e^{s_nt} \;\;\;\;\;\;\;(k=0,\cdots,N-1)$$ (12)

  The $N$ coefficients $c_1,\cdots,c_N$ can be determined based on $N$ initial conditions $y^{(k)}(0)=y^{(k)}_0,\;(k=0,\cdots,N-1)$. We evaluate $y_h(t)$ and its derivatives given above at $t=0$ and equate them to the corresponding initial values to get

  $$y^{(k)}(0)=\sum_{n=1}^N c_n s_n^k =y^{(k)}_0\;\;\;\;\;\;\;(k=0,\cdots,N-1)$$ (13)

  which can be written in matrix form:

  $${\bf Sc} =\left[\begin{array}{cccc}1&1&\cdots&1\\ s_1&s_2&\cdots&... ...[\begin{array}{c}y_0\\ y'_0\\ \vdots\\ y_0^{(N-1)}\end{array}\right] ={\bf y}_0$$ (14)

  Solving this system we get the coefficients $c_1,\cdots,c_N$ in ${\bf c}={\bf S}^{-1}{\bf y}_0$.

• Non-Homogeneous case with complex expenential input $x(t)=e^{j\omega t}$:

  We assume the output is $y(t)=c\,e^{j\omega t}$, and substitute its nth-order derivatives

  $$\frac{d^n}{dt^n} y(t)= c (j\omega)^n e^{j\omega t},\;\;\;\;\;(n=0,\cdots,N)$$ (15)

  into the equation, and get

  $$\sum_{n=0}^N a_n c (j\omega)^n e^{j\omega t} = e^{j\omega t}$$ (16)

  Dividing both sides by $e^{j\omega t}\ne 0$ and solving for $c$ we get

  $$c=\frac{1}{\sum_{n=0}^N a_n (j\omega)^n}$$ (17)

  and the solution can be found to be

  $$y(t)={\cal T} [x(t)]={\cal T} [e^{j\omega t}] =c e^{j\omega t}=\frac{e^{j\omega t}}{\sum_{n=0}^N a_n(j\omega)^n}$$ (18)

• Non-Homogeneous case with periodic input $x(t)=x(t+T)$:

  The Fourier series expansion of the periodic input $x(t+T)=x(t)$ is

  $$x(t)=\sum_{k=-\infty}^\infty X_k e^{jk\omega t}$$ (19)

  where $\omega=2\pi/T$ and $X_k$ is the kth Fourier coefficient:

  $$X_k={\cal F}[ x(t) ]=\frac{1}{T}\int_T x(t) e^{-jk\omega t}\;dt$$ (20)

  Based on the previous result, we know that the response to the kth frequency component $e^{jk\omega t}$ of $x(t)$ is

  $${\cal T} [e^{jk\omega t}] =\frac{e^{jk\omega t}}{\sum_{n=0}^N a_n(jk\omega)^n}$$ (21)

  and the total response of this linear system to $x(t)$ is simply the linear combination of these responses:
  

  $$y(t) = {\cal T} [x(t)] ={\cal T}\left[ \sum_{k=-\infty}^\infty X_k e^{jk... ... t} \right] =\sum_{k=-\infty}^\infty X_k {\cal T} \left[ e^{jk\omega t} \right]$$

  $$= \sum_{k=-\infty}^\infty X_k \left[ \frac{e^{jk\omega t}}{\sum_{n=... ...0}^N a_n(jk\omega)^n}e^{jk\omega t} =\sum_{k=-\infty}^\infty Y_k e^{jk\omega t}$$

  
  where $Y_k$ is the kth Fourier coefficient of $y(t)$:

  $$Y_k=\frac{X_k}{\sum_{n=0}^N a_n(jk\omega)^n}$$ (22)

• Non-Homogeneous case with impulse input $x(t)=\delta(t)$:

  The Fourier coefficients of the impulse input are

  $$X_k=\frac{1}{T}\int_T \delta(t) e^{-jk\omega_kt} dt=1\;\;\;\;\;\;(-\infty<k<\infty)$$ (23)

  and the output is the impulse response function:

  $$y(t)=h(t)={\cal T}[\delta(t)] =\sum_{k=-\infty}^\infty \frac{1}{\... ...0}^N a_n(jk\omega)^n}e^{jk\omega t} =\sum_{k=-\infty}^\infty H_k e^{jk\omega t}$$ (24)

  where $H_k$ is the frequency response function:

  $$H_k=H(j\omega)=\frac{1}{\sum_{n=0}^N a_n (jk\omega)^n}$$ (25)

  Now the general response given above can be written as

  $$y(t)=\sum_{k=-\infty}^\infty \frac{X_k}{\sum_{n=0}^N a_n(jk\omega... ...fty}^\infty H_k\;X_k e^{jk\omega t} =\sum_{k=-\infty}^\infty Y_k e^{jk\omega t}$$

  where $Y_k=H_k\;X_k$.

The Nth order LCCODE can also be solved by the method of Laplace transform, similar to the Fourier transform method discussed above. Specifically, we take the Laplace transform on both sides and get

$${\cal L}\left[ \sum_{n=0}^N a_n y^{(n)}(t)\right] = \sum_{n=0}^N a_n {\cal L} \left[ y^{(n)}(t)\right] =\sum_{n=0}^N a_n \left[ s^n Y(s)-\sum_{i=1}^n s^{n-i}y^{(i-1)}(0) \right]$$

$$= Y(s) \sum_{n=0}^N a_n s^n -\sum_{n=1}^N a_n \sum_{i=1}^n s^{n-i}y^{(i-1)}(0) ={\cal L}[ x(t) ]=X(s)$$ (26)

Solving for $Y(s)$ we get

$$Y(s)=\frac{X(s)+\sum_{n=1}^N a_n \sum_{i=1}^n s^{n-i}y^{(i-1)}(0)... ...{X(s)+\sum_{n=1}^N a_n \sum_{i=1}^n s^{n-i}y^{(i-1)}(0)}{\prod_{n=1}^N (s-s_n)}$$ (27)

and the general solution is

$$y(t)={\cal L}^{-1}[ Y(s) ]$$ (28)

Now consider some special cases:

• Homogeneous case:

  Given a zero input $x(t)=0$ with $X(s)=0$, we get the homogeneous solution

  $$y_h(t)={\cal L}^{-1}[ Y_h(s) ] ={\cal L}^{-1}\left[ \frac{\sum_{n=1}^N a_n \sum_{i=1}^n s^{n-i}y^{(i-1)}(0)}{\prod_{n=1}^N (s-s_n)} \right]$$ (29)

  Taking partial fraction expansion of $Y_h(s)$:

  $$Y_h(s)=\frac{\sum_{n=1}^N a_n \sum_{i=1}^n s^{n-i}y^{(i-1)}(0)}{\prod_{n=1}^N (s-s_n)} =\sum_{n=1}^N \frac{c_n}{s-s_n}$$ (30)

  we get

  $$y_h(t)={\cal L}^{-1}\left( \sum_{n=1}^N \frac{c_n}{s-s_n}\right) ... ...n=1}^N {\cal L}^{-1}\left( \frac{c_n}{s-s_n} \right) =\sum_{n=1}^N c_n e^{s_nt}$$ (31)

• Non-homogeneous case with an impulse input:

  Given $x(t)=\delta(t)$ with $X(s)={\cal L}[ \delta(t) ]=1$, we get the output, the impulse response $y(t)=h(t)$ in time domain, and the transfer function $H(s)={\cal L}[ h(t) ]$ in s-domain. For convenience, we assume zero initial conditions $y(0)=y'(0)=\cdots=y^{(N-1)}=0$, so that the transfer function can be written as:

  $$H(s)=Y(s)=\frac{1}{\sum_{n=0}^N a_ns^n}=\prod_{n=1}^N \frac{1}{s-s_n}$$ (32)

  Based on the convolution theorem, we get the impulse response function as the convolution of all $N$ function $e^{s_nt},\;(n=1,\cdots,N)$:

  $$h(t)={\cal L}^{-1}[ H(s) ]={\cal L}^{-1}\left[ \prod_{n=1}^N \frac{1}{s-s_n} \right] =e^{s_1t}*e^{s_2t}*\cdots *e^{s_Nt}$$ (33)

  Alternatively, we can also take partial fraction expansion of $H(s)$:

  $$H(s)=\frac{1}{\prod_{n=1}^N (s-s_n)}=\sum_{n=1}^N\frac{c_n}{s-s_n}$$ (34)

  and get the impulse response function:

  $$h(t)={\cal L}^{-1} [ H(s) ] =\sum_{n=1}^N c_n {\cal L}^{-1}\left[ \frac{1}{s-s_n} \right] =\sum_{n=1}^N c_n e^{s_n t}$$ (35)

• Non-homogeneous case with a general input:

  Given the impulse response $h(t)$, we can get the particular solution corresponding to a general input $x(t)$:

  $$y_p(t)=x(t)\,*\,h(t)=x(t)*e^{s_1t}*e^{s_2t}*\cdots *e^{s_Nt} =\in... ...x(t-\tau) h(\tau)\;d\tau =\sum_{n=1}^N c_n \int_0^t x(t-\tau) e^{s_n\tau} d\tau$$ (36)