# Where does the $t$ come from in the solution to a critically damped differential equation?

When solving homogeneous second order constant coefficient linear ODEs ($ay''+by'+cy=0$), there are three 'cases' that solutions fall into, based on the roots $r_1$, $r_2$ of the characteristic equation $ar^2+br+c=0$:

• $r_1, r_2 \in \Bbb R,\; r_1 \ne r_2$, aka "overdamped", for which the general solution is of the form $y(t) = c_1 e^{r_1 t} + c_2 e^{r_2 t}$.
• $r_1, r_2 \in \Bbb C,\; r_1=\overline{r_2}$, aka "underdamped", for which the general solution is $y(t) = c_1 e^{\Re(r)\,t} \cos\Im(r)\,t + c_2 e^{\Re(r)\,t} \sin\Im(r)\,t$.
• $r_1, r_2 \in \Bbb R,\; r_1 = r_2$, aka "critically damped", for which the general solution is $y(t) = c_1 e^{rt} + c_2 t e^{rt}$.

I understand the underdamped solution as essentially the same thing as the overdamped solution, equivalent to picking complex $c_1, c_2$ such that $y(t)$ is always real. But the critically damped case doesn't seem to fit neatly as a generalization of the other cases - that extra 't' makes it at least appear like it's something else entirely.

Where does the "$t$" actually come from in the critically damped case? Why don't the other cases require or allow this?

• You gathered some nice answers, so make sure to select the best for you! :) Mar 16, 2017 at 8:50

For a critically damped system, you have a double root so that the usual exponentials $e^{r_1t}$ and $e^{r_2t}$ merge in $e^{rt}$ and no longer form two linearly independent functions. Hence they cannot generate all solutions.

Now if you plug $te^{rt}$ in the equation, you get

$$a(2r+r^2t)e^{rt}+b(1+rt)e^{rt}+cte^{rt}=2ar+b+(ar^2+br+c)t.$$ The $t$ term obviously cancels out (whatever the damping), and so does the constant one because the double root verifies $2ar+b=0$.

A heuristic reasoning can justify the choice of such a term.

Let the damping converge to the critical one, with roots $r$ and $r+h$ getting closer and closer. We consider two independent solutions $e^{rt}$ and $e^{(r+h)t}-e^{rt}$, and amplify the second to avoid cancellation:

$$\lim_{h\to0}\frac{e^{(r+h)t}-e^{rt}}h=e^{rt}\lim_{h\to0}\frac{e^{ht}-1}h=e^{rt}t.$$

The general theory guarantees a two-dimensional solution space, whatever the given coefficients $a\ne0$, $b$, $c$. If the characteristic equation has a double root we obtain just one solution of the form $t\mapsto e^{\lambda t}$. There have to be other solutions, and it so happens that the $t$ thrown in works.

But you want to know why $t$, and not $\log t$, or somethig even more complicated. In this regard consider the following example: $$y'' -\lambda y=0\ .$$ If $\lambda\ne0$ we have two different characteristic values, namely $\lambda$ and $0$, and the general solution is $$y(t)=c_1 e^{\lambda t}+c_2\ .$$ The solution space ${\cal L}$ contains in particular the function $$y_{\rm special}(t):={e^{\lambda t}-1\over\lambda}\ .$$ Now, if the underlying dynamical system is continuously modified (by turning a knob) such that $\lambda\to0$ then this special solution will converge to $y_*(t)=t$. It is then plausible that this $y_*(\cdot)$ will in fact be a solution in the limiting case.

We lose nothing by starting with the equation $$y'' + 2by' + cy = 0,$$ since the coefficient of $y''$ must be nonzero. We can proceed in a couple of different ways:

## Reduction of Order

The simplest way is to suppose that $\alpha$ is a root of the characteristic equation $k^2+2bk+c=0$, and substitute $y=ue^{\alpha x}$. This gives $$0 = y'' + 2by' + cy = e^{\alpha x} \left( (\alpha^2+2b\alpha+c)u + u''+ 2(\alpha+b) u' \right).$$ The first term in the bracket vanishes since $\alpha$ is a root, so we end up having to solve $$u'' + 2(\alpha+b)u' = 0,$$ a first-order equation. Now, suppose that $\beta$ is the other root of the characteristic equation. We have $$k^2+2bk+c = (k-\alpha)(k-\beta) = k^2 - (\alpha+\beta)k+\alpha\beta,$$ so $2b = -\alpha-\beta$. Hence the new differential equation can be written as $$u'' + (\alpha-\beta)u' = 0.$$ There are now two possibilities:

• if $\alpha=\beta$, this is $u''=0$, which has solution $At+B$
• on the other hand, if $\alpha \neq \beta$, this equation can be solved by using the integrating factor $e^{(\alpha-\beta)x}$, so $$(e^{(\alpha-\beta)x}u')' = 0 \implies u' = Ae^{(\beta-\alpha)x},$$ and hence $$u = A\frac{e^{(\beta-\alpha)x}}{\beta-\alpha} +B,$$ which obviously only works if $\alpha-\beta \neq 0$.

## Limit of the Initial Value Problem

A differential equation normally comes with boundary conditions of some kind; the type that are easiest to use here for understanding are initial conditions, say $y(0) = q$, $y'(0) = p$. We can write down a general solution with unequal roots to this, $$y = Ae^{-bx}\cos{\omega t} + B e^{-bx} \frac{\sin{\omega t}}{\omega},$$ where $\omega = \sqrt{c-b^2}$ is the square root of the discriminant. (This solution works for both real and imaginary $\omega$, and is always real, because $(\sin{\omega t})/\omega$ is real in both cases.) Note that this works because the quadratic formula gives the roots as $-b \pm i\omega$.

Now, the solution with the right initial conditions can be found to be $$y = q e^{-bx}\cos{\omega t} + (p+bq) e^{-bx} \frac{\sin{\omega t}}{\omega}$$ by solving the equations simultaneously, and this is a good reason for writing the solution in the form we did in the first place. Now, we can rewrite this in terms of the roots: $$b = -\frac{\alpha+\beta}{2}, i\omega = \frac{\alpha-\beta}{2},$$ where we chose which is $\alpha$ and which $\beta$ so that the sign makes sense. (Since $\cos{\omega t}$ and $(\sin{\omega t})/\omega$ are both even, it actually makes no difference.)

The upshot of all this is that taking $\beta \to \alpha$ is the same as taking $\omega$ to $0$. We have the limits $$\cos{\omega t} \to 1, \quad \frac{\sin{\omega t}}{\omega} \to t$$ as $\omega \to 0$, so the solution assumes the form we wanted.

The moral of the latter story is that $(\sin{\omega t})/\omega$ is the important function to consider in initial value problems, and has a limit that gives the answer we expect. (This function is also important in Fourier theory, and geometry of surfaces with constant curvature, as it happens.)

• The part with the limit of the underdamped equation as $\omega\to 0$ makes sense, but is there a way to take the limit of the overdamped form to get the critical form? I tried writing it out in terms of the fundamental solution set, but all the terms appear to diverge when taking the limit. Mar 15, 2017 at 17:34
• Or for example taking $\lim_{a\to b} c_1 e^{at} + c_2 e^{bt}$, with $c_1=\frac{by_0-y_0'}{b-a}$ and $c_2=\frac{y_0-ay_0'}{b-a}$. Mar 15, 2017 at 17:41
• You can write the overdamped solution in the same way, as $$Ae^{-bt}\cosh{\xi t} + Be^{-bt} \frac{\sinh{\xi t}}{\xi},$$ where $\xi = \sqrt{b^2-c} = \frac{\alpha-\beta}{2}$, and the limit works in the same way (indeed, this is what I meant when I said that it held for real and imaginary $\omega$: the two forms are the same with $i\omega \leftrightarrow \xi$). Mar 15, 2017 at 18:12
• For the other form, it should work if you write $b = a+\epsilon$ and take $\epsilon \to 0$: then $$c_1 e^{at}+c_2 e^{bt} = e^{at}\left( \frac{(a+\epsilon)y_0-y_0'}{\epsilon} + (y'_0-ay_0)\frac{e^{\epsilon t}}{\epsilon} \right) = e^{at} \left( y_0 + (y_0'-ay_0)\frac{e^{\epsilon t}-1}{\epsilon} \right),$$ and the second term tends to $t$ as $\epsilon \to 0$. (Note that your $c_2$ has $y_0$ and $y_0'$ the wrong way round: should be $(y'_0-ay_0)/(b-a)$, if only on dimensional grounds.) Mar 15, 2017 at 18:23
• Ok, thanks! That worked perfectly to get me to what I needed: $$y(t) = \left(1 + a\frac{e^{\epsilon t} - 1}{\epsilon}\right)e^{at}y_0 + \frac{e^{\epsilon t} - 1}{\epsilon}e^{at}y_0' = (1+at)e^{at}y_0+te^{at}y_0' = \left(\cos \omega t - a \frac{\sin \omega t}{\omega}\right)e^{at}y_0+\frac{\sin \omega t}{\omega}e^{at}y_0'$$ Mar 15, 2017 at 18:59

By dividing through by the coefficient $a$, we may as well assume our o.d.e. has the form $$y'' + \beta y' + \gamma y = 0 .$$ As you say, critical damping means that the characteristic polynomial, $$r^2 + \beta r + \gamma,$$ has a double root, say, $r_0$, so that $r^2 + \beta r + \gamma = (r - r_0)^2$; in particular, $\beta = - 2 r_0$ and $\gamma = r_0^2$. Now, expanding the l.h.s. of the differential equation $$(y e^{-r_0 t})'' = 0$$ by applying the Leibniz Rule twice yields $$e^{-r_0 t} (y'' - 2 r_0 y' + r_0^2 y) = 0 ,$$ and multiplying through by the exponential factor and the above expressions for $\beta, \gamma$ shows that this differential equation is equivalent to our original one. Integrating twice both sides of $$(y e^{-r_0 t})'' = 0$$ gives that $y e^{-r_0 t}$ is some affine function $a t + b$, and multiplying through gives the general solution $$y = e^{r_0 t}(a t + b)$$ as desired.

We have $a \ddot y + b \dot y + c y = 0$. Let $x_1 := y$ and $x_2 := \dot y$. Hence,

$$\begin{bmatrix} \dot x_1\\ \dot x_2\end{bmatrix} = \begin{bmatrix} 0 & 1\\ -\frac ca & -\frac ba\end{bmatrix} \begin{bmatrix} x_1\\ x_2\end{bmatrix}$$

More succinctly,

$$\dot{\mathrm x} = \mathrm A \mathrm x$$

Let a Jordan decomposition of $\mathrm A$ be $\mathrm A = \mathrm P \mathrm J \mathrm P^{-1}$. Let $\eta := \mathrm P^{-1} x$. Hence, we have

$$\dot{\eta} = \mathrm J \eta$$

Integrating, we obtain

$$\eta (t) = \exp(t \mathrm J) \, \eta_0$$

where $\eta_0$ is a vector of initial conditions. In the critically damped case, $\mathrm A$ is non-diagonalizable and, thus, $\mathrm J$ has a single $2 \times 2$ Jordan block

$$\mathrm J = \begin{bmatrix} r & 1\\ 0 & r\end{bmatrix}$$

and, thus,

$$\exp(t \mathrm J) = \begin{bmatrix} e^{r t} & t \, e^{r t}\\ 0 & e^{r t}\end{bmatrix}$$

If $r < 0$ and assuming nonzero initial conditions, then cascading two leaky integrators with the same time constant leads to (non-oscillatory) resonance. Hence, the $t$ multiplying $e^{r t}$ in $t \, e^{r t}$.