The Damped Harmonic Oscillator

On the previous page, we derived the model of the moving-coil galvanometer, and we concluded that it's identical to the damped harmonic oscillator with an external driving force: $$\begin{array}{}\text{(1)}& m\ddot{x}(t)+c\dot{x}(t)+kx(t)=f(t),\end{array}$$ where $x(t)$ is the displacement, $m$ is the mass (cfr. moment of inertia $I$), $c$ is the damping factor, $k$ is the spring constant, and $f(t)$ is an external driving force.

On this page, we'll solve the equation of the damped harmonic oscillator analytically, discussing the different solution regimes, and calculating the key features of the step response, such as overshoot and rise time.

Solution of the step response of the damped harmonic oscillator

To make solving the equation easier, we'll define two constants: $$\begin{array}{}\text{(2)}& \begin{array}{rl}{\omega }_n\triangleq & \sqrt{\frac{k}{m}}\\ \zeta \triangleq & \frac{c}{2\sqrt{km}}\end{array}\end{array}$$ ${\omega }_n$ is called the natural frequency, and $\zeta$ the damping factor. The origin of these names will become clear in the next section.

Equation $\text{(1)}$ then becomes: $$\begin{array}{}\text{(3)}& \ddot{x}(t)+2\zeta {\omega }_n\dot{x}(t)+{\omega }_n^{2}x(t)=\frac{1}{m}f(t)\end{array}$$

Equation $\text{(3)}$ is a linear ordinary differential equation with constant coefficients, so the general solution is the sum of the homogeneous solution and a particular solution: $$x(t)=x_h(t)+x_p(t)$$

In this section, we'll derive the solution of the step response, that is, the solution for $x(t)$ if $f(t)=u(t)$, the Heaviside step function: $$u(t)\triangleq \{\begin{array}{ll}0& t<0\\ 1/2& t=0\\ 1& t>0\end{array}$$

Homogeneous solution

The homogeneous solution $x_h(t)$ is the solution to $$\begin{array}{}\text{(4)}& {\ddot{x}}_h(t)+2\zeta {\omega }_n{\dot{x}}_h(t)+{\omega }_n^2{x}_h(t)=0.\end{array}$$ The characteristic equation is $$\begin{array}{}\text{(5)}& {\lambda }^2+2\zeta {\omega }_n\lambda +{\omega }_n^2=0,\end{array}$$ with the solutions $$\begin{array}{}\text{(6)}& {\lambda }_{1,2}=-\zeta {\omega }_n\pm {\omega }_n\sqrt{{\zeta }^2-1}.\end{array}$$ If ${\lambda }_1\ne {\lambda }_2$, the homogeneous solution is $$\begin{array}{}\text{(7)}& x_h(t)=c_1{e}^{{\lambda }_{1}t}+c_2{e}^{{\lambda }_{2}t},\end{array}$$ where $c_1$ and $c_2$ are two arbitrary real constants of integration that depend on the initial conditions.

If ${\lambda }_1={\lambda }_2=\lambda$, the terms $c_1{e}^{\lambda t}$ and $c_2{e}^{\lambda t}$ are linearly dependent, so they no longer form a basis for the homogeneous solution space. One can form a valid basis by multiplying the second term by $t$: $$\begin{array}{}\text{(8)}& x_h(t)=c_1{e}^{\lambda t}+c_{2}t{e}^{\lambda t}\end{array}$$ It can be verified that these are indeed solutions by simply substituting equations $\text{(7)}$ and $\text{(8)}$ and their derivatives into equation $\text{(4)}$.

Particular solution for the step response

The particular solution $x_p(t)$ is a (preferably simple) solution to $$\begin{array}{}\text{(9)}& {\ddot{x}}_p(t)+2\zeta {\omega }_n{\dot{x}}_p(t)+{\omega }_n^2{x}_p(t)=\frac{1}{m}u(t).\end{array}$$

The Heaviside step function $u(t)$ is flat everywhere, so its derivatives are zero, except in $t=0$, where they are undefined: $$\begin{array}{rl}\mathrm{\forall }t\ne 0:& \\ \dot{u}(t)& \equiv 0\\ \ddot{u}(t)& \equiv 0\end{array}$$

A guess for the particular solution could be $$\begin{array}{}\text{(10)}& x_p(t)=c_{3}u(t).\end{array}$$ You can substitute this into equation $\text{(9)}$ to verify that this is indeed a solution, if $c_3=\frac{1}{m{\omega }_n^2}=\frac{1}{k}$: $$\cancel{c_3\ddot{u}(t)}+\cancel{2c_3\zeta {\omega }_n\dot{u}(t)}+c_3{\omega }_n^{2}u(t)=\frac{1}{m}u(t)$$

General solution for the step response

The general solution is the sum of the homogeneous and the particular solution: $$\begin{array}{}\text{(11)}& x(t)=c_1{e}^{{\lambda }_{1}t}+c_2{e}^{{\lambda }_{2}t}+\frac{1}{k}u(t)\end{array}$$ In the case where ${\lambda }_1={\lambda }_2=\lambda$, this becomes: $$\begin{array}{}\text{(12)}& x(t)=c_1{e}^{\lambda t}+c_{2}t{e}^{\lambda t}+\frac{1}{k}u(t)\end{array}$$

Determining the constants of integration

To find $c_1$ and $c_2$, we use the following initial conditions: $$\{\begin{array}{l}x(0^{+})=0\\ \dot{x}(0^{+})=0\end{array}$$ They are specified at $t=0^{+}$, which should be understood as a limit for $t$ approaching zero from the right, because the derivatives are undefined at $t=0$. $$\dot{x}(0^{+})=0\iff \underset{t\to 0^{+}}{lim}\dot{x}(t)=0$$ For the practical computation of $c_1$ and $c_2$, this doesn't matter, but the discontinuity of $\dot{x}$ does imply that we cannot extend the solution for $t<0$.

Case 1:   ${\lambda }_1\ne {\lambda }_2$

In order to be able to evaluate the second boundary condition, we need an expression for $\dot{x}(t)$. Luckily, $x(t)$ consists of just exponentials and a flat Heaviside step (recall that $\dot{u}(t)\equiv 0$). $$\dot{x}(t)=c_1{\lambda }_1{e}^{{\lambda }_{1}t}+c_2{\lambda }_2{e}^{{\lambda }_{2}t}$$ Evaluating for $t\to 0^{+}$: $$\{\begin{array}{ll}x(0^{+})=c_1+c_2+{\displaystyle \frac{1}{k}}& =0\\ \dot{x}(0^{+})=c_1{\lambda }_1+c_2{\lambda }_2& =0\end{array}$$ This is simply a system of two equations and two unknowns, $c_1$ and $c_2$. The solutions are: $$\{\begin{array}{l}{c}_1={\displaystyle \frac{1}{k({\lambda }_1/{\lambda }_2-1)}}\\ c_2={\displaystyle \frac{1}{k({\lambda }_2/{\lambda }_1-1)}}\end{array}$$

Case 2:   ${\lambda }_1={\lambda }_2=\lambda$

In this case, the derivative is a bit more complicated: $$\dot{x}(t)=c_1\lambda e^{\lambda t}+c_2{e}^{\lambda t}+c_2\lambda t{e}^{\lambda t}$$ Evaluating for $t\to 0^{+}$: $$\{\begin{array}{ll}x(0^{+})=c_1+0+{\displaystyle \frac{1}{k}}& =0\\ \dot{x}(0^{+})=c_1\lambda +c_2& =0\end{array}$$ The solutions to this system of equations are: $$\{\begin{array}{l}{c}_1=-{\displaystyle \frac{1}{k}}\\ c_2=\phantom{-}{\displaystyle \frac{\lambda }{k}}\end{array}$$

Solution regimes of the damped harmonic oscillator

Depending on the parameters $k$, $m$ and $d$, the solution $x(t)$ can look qualitatively very different. If the discriminant of equation $\text{(5)}$ is positive, its roots ${\lambda }_{1,2}$ will be real. If the discriminant is negative, ${\lambda }_{1,2}$ will be complex. If the roots are real, we say that the oscillator is in the overdamped regime, and if the roots have a nonzero imaginary part, we say that the oscillator is in the underdamped regime. On the boundary between these two regimes, the discriminant of $\text{(5)}$ is exactly zero, and ${\lambda }_1={\lambda }_2$. This is called the critically damped regime.

Looking at equation $\text{(6)}$, you can see that the roots will be real if ${\zeta }^2>1$, and complex if ${\zeta }^2<1$. In conclusion:

Regime	Damping factor	Physical constants
Underdamped	$\zeta <1$	$c<2\sqrt{km}$
Critically damped	$\zeta =1$	$c=2\sqrt{km}$
Overdamped	$\zeta >1$	$c>2\sqrt{km}$

The figure below shows the qualitative differences between the different regimes.

Image source code

Underdamped regime

We know that in this case, ${\zeta }^2-1<0$, so we can revise the formula for ${\lambda }_{1,2}$: $$\begin{array}{}\text{(13)}& \begin{array}{rl}{\lambda }_{1,2}& =-\zeta {\omega }_n\pm {\omega }_n\sqrt{{\zeta }^2-1}\\ & =-\zeta {\omega }_n\pm j{\omega }_n\sqrt{1-{\zeta }^2}\\ & =-\alpha \pm j\beta ,\end{array}\end{array}$$ where $j$ is the imaginary unit $(j^2=-1)$, $\alpha$ and $\beta$ are real, positive constants defined as $$\begin{array}{}\text{(14)}& \begin{array}{rl}\alpha & \triangleq \zeta {\omega }_n\\ \beta & \triangleq {\omega }_n\sqrt{1-{\zeta }^2}.\end{array}\end{array}$$

Substituting this into the solution (equation $\text{(11)}$), and using Euler's formula $(e^{j\theta }=\cos \theta +j\sin \theta )$: $$\begin{array}{}\text{(15)}& \begin{array}{rl}x(t)& =c_1{e}^{(-\alpha +j\beta )t}+c_2{e}^{(-\alpha -j\beta )t}+\frac{1}{k}u(t)\\ & =e^{-\alpha t}(c_1{e}^{j\beta t}+c_2{e}^{-j\beta t})+\frac{1}{k}u(t)\\ & =e^{-\alpha t}(c_1(\cos \beta t+j\sin \beta t)+c_2(\cos \beta t-j\sin \beta t))+\frac{1}{k}u(t)\\ & =e^{-\alpha t}((c_1+c_2)\cos \beta t+j(c_1-c_2)\sin \beta t)+\frac{1}{k}u(t)\\ & =e^{-\alpha t}({\tilde{c}}_1\cos \beta t+{\tilde{c}}_2\sin \beta t)+\frac{1}{k}u(t)\end{array}\end{array}$$ You could compute the new constants ${\tilde{c}}_1$ and ${\tilde{c}}_2$ using their definitions $$\begin{array}{}\text{(16)}& \{\begin{array}{l}{\tilde{c}}_1=c_1+c_2\\ {\tilde{c}}_2=j(c_1-c_2)\end{array}\end{array}$$ but it's probably easiest to just solve the system of initial conditions again: $$\begin{array}{}\text{(17)}& \begin{array}{rl}\dot{x}(t)=& -\alpha e^{-\alpha t}({\tilde{c}}_1\cos \beta t+{\tilde{c}}_2\sin \beta t)\\ & +e^{-\alpha t}(-{\tilde{c}}_1\beta \sin \beta t+{\tilde{c}}_2\beta \cos \beta t)\\ =& -e^{-\alpha t}(({\tilde{c}}_1\alpha -{\tilde{c}}_2\beta )\cos \beta t+({\tilde{c}}_1\beta +{\tilde{c}}_2\alpha )\sin \beta t)\end{array}\end{array}$$ $$\{\begin{array}{ll}x(0^{+})={\tilde{c}}_1+0+{\displaystyle \frac{1}{k}}& =0\\ \dot{x}(0^{+})=-{\tilde{c}}_1\alpha +{\tilde{c}}_2\beta & =0\end{array}$$ $$\begin{array}{}\text{(18)}& \{\begin{array}{l}{\tilde{c}}_1=-{\displaystyle \frac{1}{k}}\\ {\tilde{c}}_2=-{\displaystyle \frac{\alpha }{\beta k}}\end{array}\end{array}$$

Peak time

Now that we have a nice sinusoidal expression for the solution, we can determine the position of the first peak, the global maximum of the step response (see the previous figure).

The factor $e^{-\alpha t}$ in equation $\text{(15)}$ is monotonically decreasing, and the other factor is sinusoidal with a constant amplitude, so the global maximum of $x(t)$ occurs at the same time as the first local maximum of ${\tilde{c}}_1\cos \beta t+{\tilde{c}}_2\sin \beta t$. To find the value for $t$ at this maximum, $t_p$, we look at the roots of the derivative $\dot{x}$ (equation $\text{(17)}$): $$\begin{array}{}\text{(19)}& \begin{array}{rlrl}& & \dot{x}(t_p)& =0\\ \iff & & -e^{-\alpha t_p}(({\tilde{c}}_1\alpha -{\tilde{c}}_2\beta )\cos \beta t_p+({\tilde{c}}_1\beta +{\tilde{c}}_2\alpha )\sin \beta t_p)& =0\\ \iff & & \cancel{({\tilde{c}}_1\alpha -{\tilde{c}}_2\beta )}\cos \beta t_p+({\tilde{c}}_1\beta +{\tilde{c}}_2\alpha )\sin \beta t_p& =0\\ \iff & & \sin \beta t_p& =0\end{array}\end{array}$$ The peak time is the first solution of $\sin \beta t_p=0$ for $t>0$: $$\begin{array}{}\text{(20)}& t_p=\frac{\pi }{\beta }\end{array}$$

Overshoot

By plugging in this value $t_p$ into the expression for $x(t)$, we can determine the amplitude of peak, and by how much it overshoots the limit $x(\mathrm{\infty })=\underset{t\to \mathrm{\infty }}{lim}x(t)=1/k$. $$\begin{array}{rl}x(t_p)& =e^{-\alpha \pi /\beta }({\tilde{c}}_1\cos \frac{\beta \pi }{\beta }+{\tilde{c}}_2\sin \frac{\beta \pi }{\beta })+\frac{1}{k}\\ & =-e^{-\alpha \pi /\beta }{\tilde{c}}_1+\frac{1}{k}\\ & =e^{-\alpha \pi /\beta }\frac{1}{k}+\frac{1}{k}\end{array}$$ Expressed as a percentage, the overshoot is $$\begin{array}{}\text{(21)}& \begin{array}{rl}\text{P.O.}& =100\cdot \frac{x(t_p)-x(\mathrm{\infty })}{x(\mathrm{\infty })}\mathrm{\%}\\ & =100\cdot \exp (-\frac{\alpha \pi }{\beta })\mathrm{\%}\\ & =100\cdot \exp (-\frac{\zeta \pi }{\sqrt{1-{\zeta }^2}})\mathrm{\%}.\end{array}\end{array}$$

Image source code

Rise time

The rise time $t_r$ is the time it takes for the step response to reach $x(\mathrm{\infty })$ for the first time. $$\begin{array}{rlrl}& & x(t_r)& =x(\mathrm{\infty })\\ \iff & & e^{-\alpha t_r}({\tilde{c}}_1\cos \beta t_r+{\tilde{c}}_2\sin \beta t_r)+\frac{1}{k}& =\frac{1}{k}\\ \iff & & {\tilde{c}}_1\cos \beta t_r+{\tilde{c}}_2\sin \beta t_r& =0\\ \iff & & \tan \beta t_r& =-\frac{{\tilde{c}}_1}{{\tilde{c}}_2}\\ \iff & & \tan \beta t_r& =-\frac{\beta }{\alpha }\end{array}$$ The first solution gives a negative value for $t_r$, so $\pi$ is added to find the second solution, which is the first solution for $t_r>0$. $$\begin{array}{}\text{(22)}& \begin{array}{rl}{t}_r& =\frac{1}{\beta }(\pi -\mathrm{atan}\frac{\beta }{\alpha })\\ & =\frac{1}{{\omega }_n\sqrt{1-{\zeta }^2}}(\pi -\mathrm{atan}\frac{\sqrt{1-{\zeta }^2}}{\zeta })\end{array}\end{array}$$

Since the rise time is inversely proportional to the natural frequency of the system, it's common to eliminate it from the equation: $$\begin{array}{}\text{(23)}& t_r{\omega }_n=\frac{1}{\sqrt{1-{\zeta }^2}}(\pi -\mathrm{atan}\frac{\sqrt{1-{\zeta }^2}}{\zeta })\end{array}$$

Image source code

The lower the damping ratio, the faster the output rises, so the smaller the rise time. As the damping ratio increases, the system reacts more slowly, and the rise time increases. For $\zeta =1$, the system becomes critically damped, and the output will get asymptotically closer to $x(\mathrm{\infty })$, but never reach it, so the rise time is infinite.

Settling time

The settling time $t_s$ is the time after which the step response settles in a small error band around $x(\mathrm{\infty })$, never to leave it again. Usually, the width of the error band is chosen to be $2\mathrm{\%}$ of $x(\mathrm{\infty })$: $$\frac{|x(t_s)-x(\mathrm{\infty })|}{x(\mathrm{\infty })}<1\mathrm{\%}$$

In equation $\text{(19)}$, we derived that the extrema of the step response occur when $\sin \beta t=0$, or $t_{x,n}=n\pi /\beta$. Combining this with equations $\text{(15)}$ and $\text{(18)}$, we find that the value of $x(t)$ at these extrema is: $$\begin{array}{}\text{(24)}& \begin{array}{rl}x(t_{x,n})& =e^{-\alpha t_{x,n}}{\tilde{c}}_1(-1{)}^n+\frac{1}{k}\\ & =\frac{1}{k}(1-(-1{)}^n{e}^{-\alpha t_{x,n}})\\ & =x(\mathrm{\infty })(1-(-1{)}^n{e}^{-\alpha t_{x,n}})\end{array}\end{array}$$

The amplitude of these extrema can be used to approximate the settling time: $$\begin{array}{}\text{(25)}& \begin{array}{rl}{e}^{-\alpha t_s}& \approx 1\mathrm{\%}\\ t_s& \approx -\frac{\log (0.01)}{\alpha }\\ & \approx \frac{4.61}{\alpha }\\ & =\frac{4.61}{\zeta {\omega }_n}\end{array}\end{array}$$

The following image visualizes the exponential decay of the extrema, as derived in equation $\text{(24)}$.

Image source code

Damped frequency

The period of the sinusoidal factor of the step response (equation $\text{(15)}$) is $$T_d=\frac{2\pi }{\beta }.$$ From equation $\text{(19)}$, we also know that this is equal to distance between two successive maxima or minima of the step response.

The corresponding frequency $T_d^{-1}$ is therefore called the damped frequency of the oscillator. It's the frequency of the oscillations of the step response of the damped harmonic oscillator. Usually, the damped frequency is expressed as an angular frequency: $$\begin{array}{}\text{(26)}& \begin{array}{rl}{\omega }_d& =\frac{2\pi }{T_d}\\ & =\beta \\ & ={\omega }_n\sqrt{1-{\zeta }^2}\end{array}\end{array}$$

Image source code

Natural frequency

Looking at equation $\text{(26)}$ and the previous figure, you can see that as the damping ratio decreases to zero (no damping), the frequency of the oscillations becomes equal to ${\omega }_n$. This explains the choice of parameters at the beginning of this page: the natural frequency ${\omega }_n$ is the frequency at which the oscillator would oscillate if all damping were removed.

Overdamped regime

In this case, ${\lambda }_1$ and ${\lambda }_2$ will be real, and the step response (equation $\text{(11)}$) is made up of two exponential terms and the step function: $$x(t)=c_1{e}^{{\lambda }_{1}t}+c_2{e}^{{\lambda }_{2}t}+\frac{1}{k}u(t)$$ The following image clearly shows the interaction between the two exponential terms. The step response is the blue curve, it is the sum of the orange and green curves.

Image source code

At $t=0^{+}$, the two terms cancel out (this was imposed by the initial condition $x(0^{+})=0$). The green curve decays to zero relatively quickly, so the actual step response is pretty close to the orange curve, which approaches the final value more slowly. The green curve does have an important contribution at the beginning: it ensures that the slope of the step response gradually increases. This was the second initial condition: $\dot{x}(0^{+})=0$.

Rise time

Unlike in the underdamped case, the step response of an overdamped oscillator never reaches its final value $x(\mathrm{\infty })$, it approaches it asymptotically. When using the same definition of rise time as in the underdamped case, you could only conclude that $t_r=\mathrm{\infty }$ for any overdamped oscillator, which is of course meaningless. As a solution, we define the rise time of an overdamped system as the time it takes to go from $10\mathrm{\%}$ of $x(\mathrm{\infty })$ to $90\mathrm{\%}$ of $x(\mathrm{\infty })$.

To compute the rise time, we need to find $t_1$ and $t_2$ such that: $$\begin{array}{}\text{(27)}& \begin{array}{rl}x(t_1)& =0.1x(\mathrm{\infty })\\ x(t_2)& =0.9x(\mathrm{\infty })\\ t_r& =t_2-t_1\end{array}\end{array}$$ This results in equations of the form $$\begin{array}{}\text{(28)}& \frac{1}{{\lambda }_1/{\lambda }_2-1}{e}^{{\lambda }_{1}t}+\frac{1}{{\lambda }_2/{\lambda }_1-1}{e}^{{\lambda }_{2}t}=\gamma -1,\end{array}$$ where $\gamma$ is either $0.1$ or $0.9$.
Note that the coefficients $\frac{1}{{\lambda }_1/{\lambda }_2-1}$ and $\frac{1}{{\lambda }_2/{\lambda }_1-1}$ do not depend on ${\omega }_n$, only on the damping factor $\zeta$.

Unfortunately, this type of equation doesn't have a general analytical solution, so we'll have to settle for an approximated or numerical solution.

A first observation is that as the damping factor $\zeta$ increases, ${\lambda }_1$ decreases in magnitude, while ${\lambda }_2$ increases. The inverses of ${\lambda }_{1,2}$ are the time constants ${\tau }_{1,2}$ of the exponential terms in $x(t)$. As ${\lambda }_2$ increases, the corresponding time constant ${\tau }_2=1/{\lambda }_2$ becomes very small. This means that the second term $c_2{e}^{{\lambda }_{2}t}$ decays very quickly compared to the first term. This can be seen in the previous image, where the green curve decays much faster than the orange one.

As a consequence of ${\lambda }_1$ decreasing and ${\lambda }_2$ increasing, $c_1$ approaches $-1$, and $c_2$ goes to zero, as shown in the following figure.

Image source code

This means that if $\zeta$ is large enough, the second term of $x(t)$ can be ignored: $$x(t)\approx c_1{e}^{{\lambda }_{1}t}+c_3$$ An approximation of $t_1$ and $t_2$ would therefore be the solution of $$\frac{1}{{\lambda }_1/{\lambda }_2-1}{e}^{{\lambda }_1\tilde{t}}=\gamma -1.$$ Moving the denominator to the right-hand side and taking the logarithm from both sides: $$\begin{array}{}\text{(29)}& \begin{array}{rl}{\tilde{t}}_{1,2}& =\frac{\log ((\gamma -1)(\frac{{\lambda }_1}{{\lambda }_2}-1))}{{\lambda }_1}\\ & =\frac{\log ((\gamma -1)(\frac{{\lambda }_1}{{\lambda }_2}-1))}{-\zeta {\omega }_n+{\omega }_n\sqrt{{\zeta }^2-1}}\\ & =\frac{1}{{\omega }_n}\frac{\log ((\gamma -1)(\frac{{\lambda }_1}{{\lambda }_2}-1))}{-\zeta +\sqrt{{\zeta }^2-1}}.\end{array}\end{array}$$ The corresponding rise time is then $$\begin{array}{}\text{(30)}& \begin{array}{rl}{\tilde{t}}_r& ={\tilde{t}}_2-{\tilde{t}}_1\\ & =\frac{1}{{\omega }_n}\frac{\log ((0.9-1)(\frac{{\lambda }_1}{{\lambda }_2}-1))}{-\zeta +\sqrt{{\zeta }^2-1}}-\frac{1}{{\omega }_n}\frac{\log ((0.1-1)(\frac{{\lambda }_1}{{\lambda }_2}-1))}{-\zeta +\sqrt{{\zeta }^2-1}}\\ & =\frac{1}{{\omega }_n}\frac{\log \left(\frac{1-0.9}{1-0.1}\right)}{-\zeta +\sqrt{{\zeta }^2-1}}\\ & =\frac{1}{{\omega }_n}\frac{\log \left(9\right)}{\zeta -\sqrt{{\zeta }^2-1}}.\end{array}\end{array}$$

We can simplify the denominator of equation $\text{(30)}$: $$\begin{array}{rl}\frac{1}{\zeta -\sqrt{{\zeta }^2-1}}& =\frac{\zeta +\sqrt{{\zeta }^2-1}}{(\zeta -\sqrt{{\zeta }^2-1})\left(\zeta +\sqrt{{\zeta }^2-1}\right)}\\ & =\frac{\zeta +\sqrt{{\zeta }^2-1}}{{\zeta }^2-\zeta \sqrt{{\zeta }^2-1}+\zeta \sqrt{{\zeta }^2-1}-{\zeta }^2+1}\\ & =\zeta +\sqrt{{\zeta }^2-1}\\ & \approx 2\zeta & (\zeta \gg 1)\end{array}$$ In conclusion, the approximation of rise time by ignoring the second exponential is $$\begin{array}{}\text{(31)}& {\tilde{t}}_r=\frac{1}{{\omega }_n}\log \left(9\right)(\zeta +\sqrt{{\zeta }^2-1}).\end{array}$$ When $\zeta$ is large, this simplifies to $$\begin{array}{}\text{(32)}& \begin{array}{rlr}{\tilde{t}}_r& \approx \frac{2\log \left(9\right)}{{\omega }_n}\zeta & (\zeta \gg 1)\\ & \approx \frac{4.39}{{\omega }_n}\zeta .\end{array}\end{array}$$

To solve equation $\text{(28)}$, you could use Newton's or Halley's method, since you can easily write out $\dot{x}(t)$ and $\ddot{x}(t)$. However, it's even easier to use fixed-point iteration to solve: $$t=\frac{1}{{\lambda }_1}\log \left(\frac{\gamma -1-\frac{1}{{\lambda }_2/{\lambda }_1-1}{e}^{{\lambda }_{2}t}}{{\lambda }_1/{\lambda }_2-1}\right)$$

The following image shows the numerical solutions for the true rise time, the approximation by ignoring the second exponential, and the linear approximation, i.e. the asymptote for $\zeta \to \mathrm{\infty }$.

Image source code

As expected, the approximation that ignores the second exponential gives poor results for $\zeta$ close to $1$, but matches the true solution very nicely for large $\zeta$. Eventually, both curves approach the asymptote.

Settling time

In the overdamped case, the rise time and the settling time are tightly coupled. The settling time gives requires solving the same kind of transcendental equation of the same form as $t_{1,2}$ (equation $\text{28}$). To get the $1\mathrm{\%}$ settling time, $\gamma =0.99$.

The same kind of approximation is used as well, by ignoring the second exponential term. $$\begin{array}{}\text{(33)}& \begin{array}{rl}{\tilde{t}}_s& =\frac{1}{{\omega }_n}\log \left(\frac{100}{1-{\lambda }_1/{\lambda }_2}\right)(\zeta +\sqrt{{\zeta }^2-1})\\ & \approx \frac{2\log (100)}{{\omega }_n}\zeta & (\zeta \gg 1)\\ & \approx \frac{9.21}{{\omega }_n}\zeta \end{array}\end{array}$$

Image source code

As expected, the approximation is more accurate for the settling time than for the rise time, because the second exponential has more time to decay.

Critically damped regime

The critically damped regime is the regime with the fastest rise time while not having any overshoot or oscillations.

Image source code

Rise time

Like the overdamped regime, the $100\mathrm{\%}$ rise time of the critically damped regime is undefined, because the step response approaches its final value asymptotically, it never reaches or crosses it. The $10\mathrm{\%}-90\mathrm{\%}$ rise time is defined in the same way as for the overdamped case (equation $\text{(27)}$), and leads to similar transcendental equations: $$\begin{array}{}\text{(34)}& -e^{\lambda t}+\lambda t{e}^{\lambda t}=\gamma -1\end{array}$$ Rearranging this equation for fixed-point iteration: $$\begin{array}{}\text{(35)}& t=\frac{1}{\lambda }\log \left(\frac{\gamma -1}{\lambda t-1}\right)\end{array}$$ Solving this equation numerically yields $$\begin{array}{}\text{(36)}& t_r\approx \frac{3.35790856147781}{{\omega }_n}.\end{array}$$

Settling time

The formulation of the settling time is of same the form of equation $\text{(34)}$. The numerical result for the $1\mathrm{\%}$ settling time ($\gamma =0.99$) is $$\begin{array}{}\text{(37)}& t_s\approx \frac{6.638352067993811}{{\omega }_n}.\end{array}$$