# Non-Homogeneous Second-Order Differential Equations: Variation of Parameters

Depending on $g(t)$ often it is difficult to come up with a suitable trial solution of a given non-homogeneous equation in the method of undetermined coefficient discussed here. In this note, we discuss an alternative method called variation of parameters. First assume that we know the general solution
$$x_{\mathrm{hom}}(t)=c_1x_1(t)+c_2x_2(t)$$
of the homogeneous equation $$\ddot{x}+p(t)\dot{x}+q(t)x=0$$ Let

\label{eq:vp}
x(t)=u_1(t)x_1(t)+u_2(t)x_2(t)

be a solution of the non-homogeneous equation $$\ddot{x}+p(t)\dot{x}+q(t)x=g(t)$$ Differentiating \eqref{eq:vp}
$$\dot{x}(t)=\dot{u}_1x_1+u_1\dot{x}_1+\dot{u}_2x_2+u_2\dot{x}_2$$
Require that

\label{eq:vp2}
\dot{u}_1x_1+\dot{u}_2x_2=0

so that we have

\label{eq:vp3}
\dot{x}(t)=u_1\dot{x}_1+u_2\dot{x}_2

Differentiating \eqref{eq:vp3}

\label{eq:vp4}
\ddot{x}(t)=\dot{u}_1\dot{x}_1+u_1\ddot{x}_1+\dot{u}_2\dot{x}_2+u_2\ddot{x}_2

Substituting $\ddot{x}$, $\dot{x}$, $x$ in the non-homogeneous equation with the corresponding expressions in \eqref{eq:vp}, \eqref{eq:vp3} and \eqref{eq:vp4}, respectively results in the equation

\label{eq:vp5}
\dot{u}_1\dot{x}_1+\dot{u}_2\dot{x}_2=g(t)

Let
$$W(x_1,x_2)(t):=\begin{vmatrix} x_1(t) & x_2(t)\\ \dot{x}_1(t) & \dot{x}_2(t) \end{vmatrix}=x_1(t)\dot{x}_2(t)-x_2(t)\dot{x}_1(t)$$
$W(x_1,x_2)(t)$ is called the Wronskian of $x_1(t)$ and $x_2(t)$. If $x_1(t)$ and $x_2(t)$ are linearly dependent, so are the columns of $W(x_1,x_2)(t)$ hence $W(x_1,x_2)(t)=0$ for all $t$. This means that If $W(x_1,x_2)(t)\ne 0$ for some $t$, $x_1(t)$ and $x_2(t)$ are linearly independent. Also we have the following theorem holds.

Theorem. Let $x_1(t)$ and $x_2(t)$ be solutions of a homogeneous second-order linear differential equation. If $x_1(t)$ and $x_2(t)$ are linearly independent, then $W(x_1,x_2)(t)\ne 0$ for all $t$.

Since $x_1(t)$ and $x_2(t)$ are linearly independent, $W(x_1,x_2)\ne 0$ for all $t$ so by Cramer’s rule the solution of the system of linear equations \eqref{eq:vp2}, \eqref{eq:vp5} in $\dot{u}_1$ and $\dot{u}_2$ is given by

\begin{aligned}
\dot{u}_1(t)&=\frac{\begin{vmatrix}
0 & x_2\\
g(t) & \dot{x}_2
\end{vmatrix}
}{W(x_1,x_2)(t)}=-\frac{g(t)x_2(t)}{W(x_1,x_2)(t)}\\
\dot{u}_2(t)&=\frac{\begin{vmatrix}
x_1 & 0\\
\dot{x}_1 & g(t)
\end{vmatrix}
}{W(x_1,x_2)(t)}=\frac{g(t)x_1}{W(x_1,x_2)(t)}
\end{aligned}\label{eq:vp6}

Integrating \eqref{eq:vp6}, $u_1(t)$ and $u_2(t)$ are determined to be

\label{eq:vp7}
\begin{aligned}
u_1(t)&=-\int\frac{x_2(t)g(t)}{W(x_1,x_2)(t)}dt+c_1\\
u_2(t)&=\int\frac{x_1(t)g(t)}{W(x_1,x_2)(t)}dt+c_2
\end{aligned}

where $c_1$ and $c_2$ are constant.
Therefore

\label{eq:vp8}
\begin{aligned}
x(t)=c_1x_1(t)+&c_2x_2(t)+\left\{-x_1(t)\int\frac{x_2(t)g(t)}{W(x_1,x_2)(t)}dt+\right.\\
&\left.x_2(t)\int\frac{x_1(t)g(t)}{W(x_1,x_2)(t)}dt\right\}
\end{aligned}

is indeed the general solution of the non-homegeneous equation for

\label{eq:vp9}
X(t)=-x_1(t)\int\frac{x_2(t)g(t)}{W(x_1,x_2)(t)}dt+x_2(t)\int\frac{x_1(t)g(t)}{W(x_1,x_2)(t)}dt

being a particular solution of the non-homogeneous equation.

Example. Solve the non-homogeneous equation $\ddot{x}-3\dot{x}-4x=2\sin t$.

Solution. Recall that $x_1(t)=e^{-t}$, $x_2(t)=e^{4t}$. The Wronskian is
$$W(x_1,x_2)(t)=\begin{vmatrix} e^{-t} & e^{4t}\\ -e^{-t} & 4e^{4t} \end{vmatrix}=5e^{3t}$$
\begin{align*}
\int\frac{x_2(t)g(t)}{W(x_1,x_2)(t)}dt&=\int\frac{e^{4t}(2\sin t)}{5e^{3t}}dt\\
&=-\frac{1}{5}e^t\cos t+\frac{1}{5}e^t\sin t
\end{align*}
and
\begin{align*}
\int\frac{x_1(t)g(t)}{W(x_1,x_2)(t)}dt&=\int\frac{e^{-t}(2\sin t)}{5e^{3t}}dt\\
&=-\frac{2}{85}e^{-4t}\cos t-\frac{8}{85}e^{-4t}\sin t
\end{align*}
The general solution is
\begin{align*}
x(t)&=c_1e^{-t}+c_2e^{4t}\\
&-e^{-t}\left(-\frac{1}{5}e^t\cos t+\frac{1}{5}e^t\sin t\right)+e^{4t}\left(-\frac{2}{85}e^{-4t}\cos t-\frac{8}{85}e^{-4t}\sin t\right)\\
&=c_1e^{-t}+c_2e^{4t}+\frac{3}{17}\cos t-\frac{5}{17}\sin t.
\end{align*}

Example. Solve the non-homogeneous $\ddot{x}+4x=8\tan t$, $-\frac{\pi}{2}<t<\frac{\pi}{2}$

Solution. $x_1(t)=\cos(2t)$, $x_2(t)=\sin(2t)$. The Wronskian is
$$W(x_1,x_2)(t)=\begin{vmatrix} \cos(2t) & \sin(2t)\\ -2\sin(2t) & 2\cos(2t) \end{vmatrix}=2$$
\begin{align*}
\int\frac{x_2(t)g(t)}{W(x_1,x_2)(t)}dt&=\int\frac{\sin(2t)(8\tan t)}{2}dt\\
&=4t-2\sin(2t)
\end{align*}
and
\begin{align*}
\int\frac{x_1(t)g(t)}{W(x_1,x_2)(t)}dt&=\int\frac{\cos(2t)(8\tan t)}{2}dt\\
&=-2\cos(2t)+4\ln(\cos t)
\end{align*}
Thus the general solution is
$$x(t)=c_1\cos(2t)+c_2\sin(2t)-4t\cos(2t)+4\sin(2t)\ln(\cos t)$$

# Non-Homogeneous Second-Order Differential Equations: Method of Undetermined Coefficients

Let us consider the following second-order linear differential equation

\label{eq:nhde}
\ddot{x}+p(t)\dot{x}+q(t)x=g(t)

Let $X_1(t)$ and $X_2(t)$ be two solution of \eqref{eq:nhde}. Then $X_1(t)-X_2(t)$ is a solution of the homogeneous equation

\label{eq:hde}
\ddot{x}+p(t)\dot{x}+q(t)x=0

Thus $X_1(t)-X_2(t)=c_1x_!(t)+c_2x_2(t)$ where $x_1(t)$ and $x_2(t)$ are a fundamental set of solutions of \eqref{eq:hde}. This implies that the general solution of \eqref{eq:nhde} is given by
$$x(t)=c_1x_1(t)+c_2x_2(t)+X(t)$$
where $X(t)$ is a solution of \eqref{eq:nhde}. That is, solving the non-homogeneous equation \eqref{eq:nhde} boils down to finding a solution of \eqref{eq:nhde}. There are two methods of finding a solution of \eqref{eq:nhde}:

1.  Method of Undetermined Coefficients
2. Variation of Parameters

Method of Undetermined Coefficients

This method finds a solution by guessing a particular solution with undetermined coefficients.

Example. Solve the non-homogeneous equation $\ddot{x}-3\dot{x}-4x=3e^{2t}$.

Solution. Let $X(t)=Ae^{2t}$ be a solution of the non-homogeneous equation. Then $\dot{X}=2Ae^{2t}$ and $\ddot{X}=4Ae^{2t}$. By substitution we obtain $\ddot{X}-3\dot{X}-4X=-6Ae^{2t}$ and so $A=-\frac{1}{2}$. Hence the general solution of the non-homogeneous equation is
$$x(t)=c_1e^{-t}+c_2e^{4t}-\frac{1}{2}e^{2t}$$

Example. Solve the non-homogeneous equation $\ddot{x}-3\dot{x}-4x=2\sin t$.

Solution. Let $X(t)=A\sin t$. Then $\dot{X}(t)=A\cos t$ and $\ddot{X}=-A\sin t$. By substitution we obtain $\ddot{X}-3\dot{X}-4X=-5A\sin t -3A\cos t$. This has to be the same as $2\sin t$ which leads to the equation $(5A+2)\sin t+3A\cos t$. Since $\sin t$ and $\cos t$ are linearly independent, $5A+2=0$ and $3A=0$, a contradiction! This time we assume that $X(t)=A\sin t+B\cos t$. Then $\dot{X}(t)=A\cos t-B\sin t$, $\ddot{X}(t)=-A\sin t-B\cos t$ and
$$\ddot{X}-3\dot{X}-4X=(-5A+3B)\sin t+(-3A-5B)\cos t$$
Comparing this with $2\sin t$ we get the system of linear equations
\left\{\begin{aligned}-5A+3B&=2\\ -3A-5B&=0 \end{aligned}\right.
of which solution is $A=-\frac{5}{17}$ and $B=\frac{3}{17}$. Therefore, the general solution is
$$x(t)=c_1e^{-t}+c_2e^{4t}-\frac{5}{17}\sin t+\frac{3}{17}\cos t$$

Example. Solve the non-homogeneous equation $\ddot{x}-3\dot{x}-4x=4t^2-1$.

Solution. Let $X(t)=At^2+Bt+C$ be a solution of the non-homogeneous equation. Then $A=-1$, $B=\frac{3}{2}$ and $C=-\frac{11}{8}$. Hence the general solution is
$$x(t)=c_1e^{-t}+c_2e^{4t}-t^2+\frac{3}{2}t-\frac{11}{8}$$

Example. Solve the non-homogeneous equation $\ddot{x}-3\dot{x}-4x=-8e^t\sin 2t$.

Solution.
\begin{align*}
\dot{X}(t)&=(-2A+B)e^t\sin 2t+(A+2B)e^t\cos 2t\\
\ddot{X}(t)&=(-3A+4B)e^t\cos 2t+(-4A-3B)e^t\sin 2t
\end{align*}
and
$$\ddot{X}-3\dot{X}-4X=(-10A-2B)e^t\cos 2t+(2A-10B)e^t\cos 2t$$
Comparing this with $-8e^t\cos 2t$ we get the system of linear equations
\left\{\begin{aligned} -10A-2B&=-8\\ 2A-10B&=0 \end{aligned}\right.
of which solution is $A=\frac{10}{13}$ and $B=\frac{1}{13}$.
Hence the general solution is
$$x(t)=c_1e^{-t}+c_2e^{4t}+\frac{10}{13}e^t\cos 2t+\frac{2}{13}e^t\sin t$$

# Second-Order Linear Differential Equations and Linear Algebra

When we studied second-order linear differential equations through undamped and damped harmonic motion, we made a hand waving argument that the general solution of a second-order linear differential equation is the linear combination of two distinct solutions that are linearly independent in the same sense as for vectors i.e. one solution is not a constant multiple of another solution. These solutions are called fundamental solutions. In this note, we discuss an important and intriguing relationship between second-order linear differential equations and linear algebra and explain why the general solution is given by the linear combination of two fundamental solutions. For the sake of simplicity, we limit our discussion to the case that characteristic equation has two distinct real solutions. In order to maintain this lecture note as much self-contained as possible, we include some of the basic concepts in linear algebra that we need for our discussion below.

Let us consider a second-order linear differential equation $$\label{eq:2lde}\ddot{x}+p\dot{x}+rx=0$$ \eqref{eq:2lde} can be written as a system of two first-order linear differential equations \left\{\begin{aligned}\frac{dx}{dt}&=s\\\frac{ds}{dt}&=-ps-rx\end{aligned}\right.\label{eq:ldesys} Let $X=\begin{pmatrix}x\\s\end{pmatrix}$ and $A=\begin{pmatrix}0 & 1\\-r & -p\end{pmatrix}$. Then \eqref{eq:ldesys} can be written as the matrix differential equation $$\frac{dX}{dt}=AX$$

Definition. Let $A$ be a $2\times 2$ real matrix or equivalently a linear map $A: \mathbb{R}^2\longrightarrow\mathbb{R}^2$. A vector $v\in\mathbb{R}^2$ is called an eigenvector of $A$ if there exists a number $q\in \mathbb{R}$ such that $Av=qv$. The number $q$ is called an eigenvalue of $A$ belonging to the eigenvector $v$. We also say $v$ is an eigenvector associated with the eigenvalue $q$. Eigen is a German word and it means own or self. As you will see below, given an eigenvalue there are infinitely many eigenvectors that are associated with the eigenvalue but they are linearly dependent i.e one eigenvector is a scalar multiple of another. So the name makes sense.

How do we find eigenvalues of a matrix $A$? The equation $Av=qv$ is the same as $(A-qI)v=0$. In order for this equation to have a non-trivial solution ($v\ne 0$) it must be that $$\label{eq:cheq}\det(A-qI)=0$$ The equation \eqref{eq:cheq} is called the characteristic equation. You heard the name characteristic equation before when we discussed harmonic motion. While you must not see any resemblance, that characteristic equation and \eqref{eq:cheq} are the same thing, hence the name characteristic equation. For example, if $A=\begin{pmatrix}0 & 1\\-r & -p\end{pmatrix}$, $\det(A-qI)=q^2+pq+r$, so we see that $\eqref{eq:cheq}$ is the same as the characteristic equation $$q^2+pq+r=0$$ of the second-order linear differential equation \eqref{eq:2lde}. Again for the sake of simplicity, the rest of the discussion will be done using a simple example but the same idea applies to the general case. Let us now consider the second-order differential equation $$\frac{d^2}x{dt^2}+5\frac{dx}{dt}+6x=0$$ The matrix $A$ is $A=\begin{pmatrix}0 & 1\\-r & -p\end{pmatrix}=\begin{pmatrix}0 & 1\\-6 & -5\end{pmatrix}$ and $\det(A-qI)=q^2+5q+6=0$ has two distinct real solutions $q=-3, -2$. These are the eigenvalues of $A$. Now we find eigenvectors. For $q_1=-3$, $Av_1=q_1v_1$ with $v_1=\begin{pmatrix}a\\b\end{pmatrix}$ leads to the equation $b=-3a$. So we may choose $v_1=\begin{pmatrix}1\\-3\end{pmatrix}$. Similarly for $q_2=-2$, we find an eigenvector $v_2=\begin{pmatrix}1\\-2\end{pmatrix}$. These eigenvectors can be used to find solutions of $\frac{dX}{dt}=AX$. To see this let $Av=qv$ and $X(t)=f(t)v$. Suppose that $\frac{dX}{dt}=AX$. Then the function $f(t)$ is to be determined. \begin{align*}\frac{df(t)}{dt}v&=A(f(t)v)\\&=f(t)Av\\&=f(t)qv\end{align*} This implies that $$\frac{df(t)}{dt}=f(t)q$$ whose solution is $f(t)=Ae^{qt}$ where $A$ is a constant. So we see that $$X_1(t)=A_1e^{-3t}\begin{pmatrix}1\\-3\end{pmatrix}$$ and $$X_2(t)=A_2e^{-2t}\begin{pmatrix}1\\-2\end{pmatrix}$$ are solutions of $\frac{dX}{dt}=AX$. Since the equation is linear, their sum \begin{aligned}X_1(t)+X_2(t)&=A_1e^{-3t}\begin{pmatrix}1\\-3\end{pmatrix}+A_2e^{-2t}\begin{pmatrix}1\\-2\end{pmatrix}\\&=\begin{pmatrix}A_1e^{-3t}+A_2e^{-2t}\\-3A_1e^{-3t}-2A_2e^{-2t}\end{pmatrix}\end{aligned}\label{eq:ldesyssol} is also a solution. It turns out that \eqref{eq:ldesyssol} is the most general solution of $\frac{dX}{dt}=AX$ meaning any solution would be in the form of \eqref{eq:ldesyssol}. To understand why this is the case let us first suppose that $A$ is a diagonal matrix $$A=\begin{pmatrix}q_1 & 0\\0 & q_2\end{pmatrix}$$ with $q_i\ne 0$, $i=1,2$. Let $X(t)=\begin{pmatrix}x_1(t)\\x_2(t)\end{pmatrix}$ be a solution of $\frac{dX}{dt}=AX$. Then $$\frac{dx_1(t)}{dt}=q_1x_1(t),\ \frac{dx_2(t)}{dt}=q_2x_2(t)$$ of which solutions are $$x_1(t)=A_1e^{q_1t},\ x_2(t)=A_2e^{q_2t}$$ Now $X(t)$ can be written as $$X(t)=\begin{pmatrix}A_1e^{q_1t}\\A_2e^{q_2t}\end{pmatrix}=A_1e^{q_1t}\begin{pmatrix}1\\0\end{pmatrix}+A_2e^{q_2t}\begin{pmatrix}0\\1\end{pmatrix}$$ Note that $\begin{pmatrix}1\\0\end{pmatrix}$ and $\begin{pmatrix}0\\1\end{pmatrix}$ are the eigenvectors of $A$ associated with the eigenvalues $q_1$ and $q_2$, respectively. Conversely, any matrix-valued function $X(t)$ of the form $X(t)=\begin{pmatrix}A_1e^{q_1t}\\A_2e^{q_2t}\end{pmatrix}$ satisfies the differential equation $\frac{dX}{dt}=AX$. Let $V$ be the set of all solutions of $\frac{dX}{dt}=AX$. Then $V$ is a vector space over $\mathbb{R}$. (The verification is straightforward and is left for readers.) The above argument shows that the linearly independent solutions $e^{q_1t}\begin{pmatrix}1\\0\end{pmatrix}$ and $e^{q_2t}\begin{pmatrix}0\\1\end{pmatrix}$ form a basis for $V$. So the dimension of $V$ is 2.

Remark. Let $V$ be the set of all infinitely differentiable functions. For $f,g\in V$ and $c\in\mathbb{R}$, define $f+g$ and $cf$ by \begin{align*}(f+g)(t)&=f(t)+g(t)\\(cf)(t)&=cf(t)\end{align*} Then $V$ forms a vector space over $\mathbb{R}$. So infinitely differentiale functions can be considered as vectors. The derivative $\frac{d}{dt}$ is a map $$\frac{d}{dt}: V\longrightarrow V;\ f\longmapsto \frac{df}{dt}$$ The well-known properties of the derivative: \begin{align*}\frac{d(f+g)}{dt}&=\frac{df}{dt}+\frac{dg}{dt}\\\frac{d(cf)}{dt}&=c\frac{df}{dt}\end{align*} ensure that $\frac{d}{dt}:V\longrightarrow V$ is indeed a linear map.  Let $\lambda\in\mathbb{R}$. Then $f(t)=e^{\lambda t}$ is an eigenvector of $\frac{d}{dt}$ associated with the eigenvalue $\lambda$ because $\frac{de^{\lambda t}}{dt}=\lambda e^{\lambda t}$.

In our case, $A=\begin{pmatrix}0 & 1\\-r & -p\end{pmatrix}$ is not a diagonal matrix so the previous argument does not apply straightforwardly. However if $A$ has two distinct eigenvalues. then it is diagonalizable namely there is an invertible matrix $M$ such that $MAM^{-1}$ is a diagonal matrix. Such a matrix $M$ is called change of basis matrix. Note that $MAM^{-1}$ has exactly the same eigenvalues as those of $A$: \begin{align*}\det(MAM^{-1}-qI)&=\det(MAM^{-1}-M(qI)M^{-1})\\&=\det[M(A-qI)M^{-1}]\\&=\det(M)\det(A-qI)\det(M^{-1})\\&=\det(A_qI)\end{align*} If $v$ is the eigenvector of $A$ associated with the eigenvalue $q$, then $Mv$ is the eigenvector of $MAM^{-1}$ associated with the same eigenvalue $q$. To see this let $Av=qv$. Then \begin{align*}(MAM^{-1})(Mv)&=(MA)v\\&=M(Av)\\&=M(qv)\\&=q(Mv)\end{align*} Let $V$ be the ssolution space of $\frac{dX}{dt}=AX$ and $W$ the solutions apce of $\frac{dY}{dt}=(MAM^{-1})Y$. Since $M$ is invertible, $M: V\longrightarrow W$ is an isomorphism. For example, let $A=\begin{pmatrix}0 & 1\\-6 & -5\end{pmatrix}$ and $M$ change of basis matrix such that $MAM^{-1}=\begin{pmatrix}-3 & 0\\0 & -2\end{pmatrix}$. Set \begin{aligned}M\begin{pmatrix}1\\-3\end{pmatrix}&=\begin{pmatrix}1\\0\end{pmatrix}\\M\begin{pmatrix}1\\-2\end{pmatrix}&=\begin{pmatrix}0\\1\end{pmatrix}\end{aligned}\label{eq:isobasis} Let $M=\begin{pmatrix}a & b\\c & d\end{pmatrix}$. Then \eqref{eq:isobasis} results in the systems of linear equations \left\{\begin{aligned}a-3b&=1\\a-2b&=0\end{aligned}\right. and \left\{\begin{aligned}c-3d&=0\\c-2d&=1\end{aligned}\right. whose solutions are $a=-2$, $b=-1$, $c=1$, and $d=3$. That is, $M=\begin{pmatrix}-2 & -1\\3 & 1\end{pmatrix}$ and we have $$MAM^{-1}=\begin{pmatrix}-2 & -1\\3 & 1\end{pmatrix}\begin{pmatrix}0 & 1\\-6 & -5\end{pmatrix}\begin{pmatrix}1 & 1\\-3 & 2\end{pmatrix}=\begin{pmatrix}-3 & 0\\0 & -2\end{pmatrix}$$ as expected. Recall from linear algebra that an isomorphism $M: V\longrightarrow W$ maps a basis of $V$ to a basis of $W$. So $e^{-3t}\begin{pmatrix}1\\-3\end{pmatrix}$ and $e^{-2t}\begin{pmatrix}1\\-2\end{pmatrix}$ form basis for $V$. Therefore $X(t)=A_1e^{-3t}\begin{pmatrix}1\\-3\end{pmatrix}+A_2​e^{-2t}\begin{pmatrix}1\\-2\end{pmatrix}$ is the general solution of $\frac{dX}{dt}=\begin{pmatrix}0 & 1\\-6 & -5\end{pmatrix}X$ and consequently $x(t)=A_1e^{-3t}+A_2e^{-2t}$ is the general solution of the second-order linear differential equation $\ddot{x}+5\dot{x}+6=0$.

# Harmonic Motion: Damped

The harmonic motion we discussed here is not physically realistic because it does not take damping due to friction into account. In this note, we discuss damped harmonic motion which is physically more realistic. If $x$ is the displacement from the equilibrium position, then the restoring force exerted by a spring is $-kx$ and the retarding force is $-cv$ (friction is proportional to the velocity $v$). Consequently we have the equation $$F=-kx-cv=-kx-c\dot{x}$$ or equivalently the second-order linear equation $$\label{eq:damped}m\ddot{x}+c\dot{x}+kx=0$$ To solve the equation \eqref{eq:damped} we still attempt to use the trial solution $x=e^{qt}$. As a result we obtain $$\label{eq:chareq}mq^2+cq+k=0$$ \eqref{eq:chareq} is called the auxiliary equation or the characteristic equation. Its solution is given by $$q=\frac{-c\pm\sqrt{c^2-4mk}}{2m}$$ There are three physically distinct cases:

1. $c^2>4mk$: overdamping
2. $c^2=4mk$: critically damping
3. $c^2<4mk$: underdamping

Case 1. Let $-\gamma_1<0$ and $-\gamma_2<0$ be two real values of $q$. Then the general solution (we will discuss, using linear algebra, why this is indeed the general solution later) is $$x=A_1e^{-\gamma_1 t}+A_2e^{-\gamma_2 t}$$ The motion is nonoscillatory and the displacement $x$ decays to 0 in an exponential matter.

Example. The second-order linear differential equation $\ddot{x}+3\dot{x}+2x=0$ has the general solution $$x(t)=A_1e^{-2t}+A_2e^{-t}$$

x(t)=exp(-2t)+2exp(-t), t=0..8

Case 2. $q$ has one real value $q=-\gamma$ where $\gamma=\frac{c}{2m}$. The equation \eqref{eq:damped} can be written as $$\label{eq:criticdamped}\left(\frac{d}{dt}+\gamma\right)^2x=\left(\frac{d}{dt}+\gamma\right)\left(\frac{d}{dt}+\gamma\right)x=0$$ Let $u=\frac{d}{dt}+\gamma$. Then \eqref{eq:criticdamped} reduces to  a first-order differential equation $$\left(\frac{d}{dt}+\gamma\right)u=\frac{du}{dt}+\gamma u=0$$ which is separable. It’s solution is $u=A_1e^{-\gamma t}$. Now we have a first-order linear differential equation $$\frac{dx}{dt}+\gamma x=A_1e^{-\gamma t}$$ whose solution is given by $$x(t)=e^{-\gamma t}(A_1 t+A_2)$$

Example. The second-order linear differential equation $$\ddot{x}+2\dot{x}+x=0$$ has the general solution $$x(t)=e^{-t}(A_1t+A_2)$$

x(t)=exp(-t)(t+2), t=0..10

This also represents a nonoscillatory motion and the displacement $x$ decays to zero asymptotically. Critical damping produces an optimal return to the equilibrium position, so it is used, for example, for galvanometer suspensions.

Case 3. Suppose that $c$ is small enough so that $c^2-4mk<0$. In this case, $q$ are two complex numbers $-\gamma\pm i\omega_1$ where $\gamma=\frac{c}{2m}$, $\omega_1=\sqrt{\frac{k}{m}-\frac{c^2}{4m^2}}=\sqrt{\omega_0^2-\gamma^2}$. So $e^{(-\gamma+i\omega_1)t}$, $e^{(-\gamma-i\omega_1)t}$ are solutions of \eqref{eq:damped}. Due to the linearity of \eqref{eq:damped}, the real part $e^{-\gamma t}\cos\omega_1 t$and the imaginary part $e^{-\gamma t}\sin\omega_1 t$ of $e^{(-\gamma+i\omega_1)t}$ are also solutions of \eqref{eq:damped}. Hence, $$x(t)=ae^{-\gamma t}\cos\omega_1 t+be^{-\gamma t}\sin\omega_1 t$$ is the general solution. This can be written as $$x(t)=Ae^{-\gamma t}\cos(\omega_1 t-\theta_0)$$ where $A=\sqrt{a^2+b^2}$ and $\theta_0=\tan^{-1}\left(\frac{b}{a}\right)$. Note that the angular frequency (or natural frequency) $\omega_1$ is smaller than that of undamped harmonic oscillator $\omega_0$. $$\omega_1=\sqrt{\omega_0^2-\gamma^2}=\omega_0\sqrt{1-\left(\frac{\gamma}{\omega_0}\right)^2}$$ If $\frac{\gamma}{\omega_0}<1$ then $\sqrt{1-\left(\frac{\gamma}{\omega_0}\right)^2}\approx 1-\frac{1}{2}\frac{\gamma^2}{\omega_0^2}$. So $$\omega_1\approx \omega_0-\frac{\gamma^2}{2\omega_0}$$

Example. Let $m=1$, $c=4$ and $k=404$. Then the resulting equation of damped harmonic oscillator is $$\ddot{x}+4\dot{x}+404=0$$ Let us solve this equation with $x(0)=1$ and $\dot{x}(0)=0$. The characteristic equation $$q^2+4q+404=0$$ has two complex solutions $q=-2\pm 20i$. Thus $$x(t)=ae^{-2t}\cos(20t)+be^{-2t}\sin(20t)$$ $x(0)=1$ results in $a=1$. To determine $b$ we need $\dot{x}$. $$\dot{x}=-2e^{-2t}\cos(20t)-20e^{-2t}\sin(20t)-2be^{-2t}\sin(20t)+20be^{-2t}\cos(20t)$$ The initial condition $\dot{x}(0)=0$ results in $b=\frac{1}{10}$. Hence $x(t)$ is given by $$x(t)=e^{-2t}\cos(20t)+\frac{1}{10}e^{-2t}\sin(20t)$$ Since $A=\sqrt{a^2+b^2}=\frac{\sqrt{101}}{10}$ and $\theta_0=\tan^{-1}\left(\frac{b}{a}\right)=\tan^{-1}\left(\frac{1}{10}\right)\approx 0.099669$, $x(t)$ also can be written as $$x(t)=\frac{\sqrt{101}}{10}e^{-2t}\cos(20t-0.099669)$$

x(t)=1.00499exp(-2t)cos(20t-0.099669), t=0..pi

Differentiate the total energy $$E=\frac{1}{2}m\dot{x}^2+\frac{1}{2}kx^2$$ with respect to $t$. \begin{align*}\frac{dE}{dt}&=m\ddot{x}\dot{x}+k\dot{x}x\\&=(m\ddot{x}+kx)\dot{x}\\&=(-c\dot{x})\dot{x}\\&=-c(\dot{x})^2<0\end{align*} This is the rate at which energy is dissipated into heat by friction.

# Harmonic Motion: Undamped

A force exerted by an elastic cord or by a spring obeys Hooke’s law $F=-kx$ where $x$ is the displacement of the equilibrium position.

From Newton’s second law of motion we have $F=ma=m\ddot{x}$, so we obtain the second order linear differential equation $$\label{eq:undamped}m\ddot{x}+kx=0$$ Solving \eqref{eq:undamped} for $\ddot{x}$ we obtain $\ddot{x}=-\frac{k}{m}x\sim -x$. For a trial solution, $x=e^{qt}$ is a candidate. To see if this trial solution works, plug it back into \eqref{eq:undamped}. $$m\ddot{x}+kx=mq^2e^{qt}+ke^{qt}=0$$ i.e. $$\label{eq:auxeq}mq^2+k=0$$ whose solutions are $$q=\pm i\sqrt{\frac{k}{m}}=\pm i\omega_0$$ where $\omega_0=\sqrt{\frac{k}{m}}$. The equation \eqref{eq:auxeq} is called the auxiliary equation or the characteristic equation. So, $x_1=e^{i\omega_0 t}$ and $x_2=e^{-i\omega_0 t}$ are solutions of \eqref{eq:undamped}. It can be easily shown that their linear combination $$\label{eq:undamped2}x=A_1e^{i\omega_0 t}+A_2e^{-i\omega_0 t}$$ is also a solution of \eqref{eq:undamped}. \eqref{eq:undamped2} is a complex solution so it is not suitable for the physical analysis of a motion governed by \eqref{eq:undamped}. What we need is a real solution. It turns out that the real part and the imaginary part of $e^{i\omega_0 t}$ also, respectively, satisfy \eqref{eq:undamped}. This is due to the linearity of \eqref{eq:undamped}. Hence, an alternative form of the solution which is real is $$\label{eq:undamped3}x=a\cos\omega_0 t+b\sin\omega_0 t$$ Using a trigonometric identity \eqref{eq:undamped3} can be written as $$\label{eq:undamped4}x=\sqrt{a^2+b^2}\cos(\omega_0 t-\theta_0)$$ where $\theta_0=\tan^{-1}\frac{b}{a}$ or as $$\label{eq:unddamped4a}x=\sqrt{a^2+b^2}\sin(\omega_0 t+\phi_0)$$ where $\phi_0=\tan^{-1}\frac{a}{b}$. The angles $\theta_0$ and $\phi_0$ are called the phase.

Some Terminologies

The equation \eqref{eq:undamped} is called the differential equation of the harmonic oscillator. $\sqrt{a^2+b^2}$ is the amplitude (the maximum value of $x$) The period $T_0$ of the oscillation is the time required for one complete cycle. $$T_0=\frac{2\pi}{\omega_0}=2\pi\sqrt{\frac{m}{k}}$$ The linear frequency of oscillation $f_0$ is the number of cycles in unit time. $$f_0=\frac{1}{T_0}=\frac{\omega_0}{2\pi}=\frac{1}{2\pi}\sqrt{\frac{k}{m}}$$ $\omega_0=2\pi f_0$ is called the angular frequency and is also called the natural frequency.

Example. Suppose that a mass weighing 10 lb stretches a spring 2 in. If the mass is displaced an additional 2 in and is then set in motion with an initial upward velocity of 1 ft/s, determine the position of the mass at any later time. Also determine the period, amplitude, and phase of the motion.

Solution. $k=10\mathrm{lb}/2\mathrm{in}=60\mathrm{lb}/\mathrm{ft}$ and $m=\mathrm{weight}/g=10\mathrm{lb}/32\mathrm{ft}/\mathrm{s}^2$. Thus $$m\ddot{x}+kx=\frac{10}{32}\ddot{x}+60x=0$$ i.e. $$\ddot{x}+192x=0$$ The auxiliary equation is $q^2+192=0$ and $q=\pm i8\sqrt{3}$. The solution $x(t)$ is then given by $$x(t)=a\cos(8\sqrt{3}t)+b\sin(8\sqrt{3}t)$$ From the initial conditions $x(0)=\frac{1}{6}$ ft and $\dot{x}(0)=-1$ ft/s, we find $a=\frac{1}{6}$ and $b=-\frac{1}{8\sqrt{3}}$. The natural frequency is $\omega_0=8\sqrt{3}\approx 13.856$ rad/s. The period is $T_0=\frac{2\pi}{\omega_0}=\frac{2\pi}{8\sqrt{3}}\approx 0.453$ sec. The amplitude is $\sqrt{a^2+b^2}=\sqrt{\frac{19}{576}}\approx 0.182$ ft. The phase is $\theta_0=\tan^{-1}\left(\frac{b}{a}\right)=-\frac{\sqrt{3}}{4}\approx -0.408645$ rad. Therefore, $x(t)$ can be written as $$x(t)=0.182\cos(13.856t+0.40864)$$ Figure 1 shows a complete cycle.

Figure 1. A complete cycle of x(t)=0.182cos(13.856t+0.40864)

Figure 2. shows an animation of $x(t)$.

Figure 2. An animation of x(t)=0.182cos(13.856t+6.283185308k+0.40864) with k=1..5, t=1..100, and FPS=10

The Conservation of Total Energy

Let us calculate the work done by an external force $F_{\mathrm{ext}}$ in moving the mass from the equilibrium position ($x=0$) to some position $x$. $F_{\mathrm{ext}}=-F=kx$ and $$W=\int F_{\mathrm{ext}} dx=\int_0^x kx dx=\frac{1}{2}kx^2$$ The work $W$ is stored in the spring as potential energy $V(x)=W=\frac{1}{2}kx^2$. The potential energy $V(x)$ for a force $F(x)$ is defined by $$\label{eq:conservative}F=-\frac{dV}{dx}$$ In our case, $F=-\frac{dV}{dx}=-kx$. For any force $F$ satisfying \eqref{eq:conservative}, the total energy $$E=T+V,$$ the sum of the kinetic energy $T=\frac{1}{2}m\dot{x}^2$ and the potential energy $V$, is  constant. For this reason, a force satisfying \eqref{eq:conservative} is called a conservative force.  Since $$\ddot{x}=\frac{d\dot{x}}{dt}=\frac{d\dot{x}}{dx}\frac{dx}{dt}=v\frac{dv}{dx},$$ $$F(x)=m\ddot{x}=mv\frac{dv}{dx}=\frac{1}{2}m\frac{dv^2}{dx}=\frac{dT}{dx}$$ The work done on the particle by impress force $F(x)$ is $$\int F(x)dx=\int dT=\frac{1}{2}m\dot{x}^2+C_1$$ where $C_1$ is a constant. On the other hand, from \eqref{eq:conservative} we also have $$\int F(x)dx=-\int dV=-V(x)+C_2$$ where $C_2$ is a constant. Therefore, $$T+V=\frac{1}{2}m\dot{x}^2+V(x)=E$$ is a constant. In our case the conservation of total energy $E$ $$\label{eq:totalenergy}E=\frac{1}{2}m\dot{x}^2+\frac{1}{2}kx^2$$ can be more directly shown. Differentiating the total energy \eqref{eq:totalenergy} with respect to $t$ \begin{align*}\frac{dE}{dt}&=m\ddot{x}\dot{x}+kx\dot{x}\\&=(m\ddot{x}+kx)\dot{x}\\&=0\end{align*} So the restoring force $F=-kx$ for an undamped harmonic motion is conservative.

The solution \eqref{eq:undamped4} can be also obtained by considering energy. Solving \eqref{eq:totalenergy} for $\dot{x}$ $$\label{eq:velocityundamped}\dot{x}=\pm\sqrt{\frac{2E-kx^2}{m}}$$ This is a separable equation. For $\dot{x}=\sqrt{\frac{2E-kx^2}{m}}$ it’s solution is  $$x=\sqrt{\frac{2E}{k}}\sin\left(\sqrt{\frac{k}{m}}t+\phi_0\right)$$ and for $\dot{x}=-\sqrt{\frac{2E-kx^2}{m}}$, it’s solution is $$x=\sqrt{\frac{2E}{k}}\cos\left(\sqrt{\frac{k}{m}}t-\theta_0\right)$$ So the amplitude is $A=\sqrt{\frac{2E}{k}}$ and the angular frequency is $\omega_0=\sqrt{\frac{k}{m}}$. In order for \eqref{eq:velocityundamped} to make sense $\frac{2E-kx^2}{m}\geq 0$ i.e. the total energy is greater than or equal to the potential energy $V(x)=\frac{1}{2}kx^2$. This means that the particle is confined to the region $-\sqrt{\frac{2E}{k}}\leq x\leq\sqrt{\frac{2E}{k}}$. Also the speed becomes zero when $V(x)=E$. This means that the particle must come to rest and reverse its motion at $x=\pm\sqrt{\frac{2E}{k}}$. The points $x=\pm\sqrt{\frac{2E}{k}}$ are called the turning points of the motion. The maximum value of $\dot{x}$ occurs at $x=0$. $v_{\max}=\frac{k}{m}A=\omega_0A$ and $E=\frac{1}{2}mv_{\max}^2=\frac{1}{2}kA^2$.