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: (1)mx¨(t)+cx˙(t)+kx(t)=f(t), 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: (2)ωnkmζc2km ωn is called the natural frequency, and ζ the damping factor. The origin of these names will become clear in the next section.

Equation (1) then becomes: (3)x¨(t)+2ζωnx˙(t)+ωn2x(t)=1mf(t)

Equation (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)=xh(t)+xp(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){0t<01/2t=01t>0

Homogeneous solution

The homogeneous solution xh(t) is the solution to (4)x¨h(t)+2ζωnx˙h(t)+ωn2xh(t)=0. The characteristic equation is (5)λ2+2ζωnλ+ωn2=0, with the solutions (6)λ1,2=ζωn±ωnζ21. If λ1λ2, the homogeneous solution is (7)xh(t)=c1eλ1t+c2eλ2t, where c1 and c2 are two arbitrary real constants of integration.


If λ1=λ2=λ, the terms c1eλt and c2eλ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: (8)xh(t)=c1eλt+c2teλt Compute the first two derivatives of equation (8) and substitute them into equation (4) to convince yourself that these are indeed solutions.

Particular solution for the step response

The particular solution xp(t) is a (preferably simple) solution to (9)x¨p(t)+2ζωnx˙p(t)+ωn2xp(t)=1mu(t).

The Heaviside step function u(t) is flat everywhere, so its derivatives are zero, except in t=0, where they are undefined: t0:u˙(t)0u¨(t)0

A guess for the particular solution could be (10)xp(t)=c3u(t). You can substitute this into equation (9) to verify that this is indeed a solution, if c3=1mωn2=1k: c3u¨(t)+2c3ζωnu˙(t)+c3ωn2u(t)=1mu(t)

General solution for the step response

The general solution is the sum of the homogeneous and the particular solution: (11)x(t)=c1eλ1t+c2eλ2t+1ku(t) In the case where λ1=λ2=λ, this becomes: (12)x(t)=c1eλt+c2teλt+1ku(t)

Determining the constants of integration

To find c1 and c2, we use the following initial conditions: {x(0+)=0x˙(0+)=0 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. x˙(0+)=0limt0+x˙(t)=0 For the practical computation of c1 and c2, this doesn't matter, but the discontinuity of x˙ does imply that we cannot extend the solution for t<0.

Case 1:   λ1λ2

In order to be able to evaluate the second boundary condition, we need an expression for x˙(t). Luckily, x(t) consists of just exponentials and a flat Heaviside step (recall that u˙(t)0). x˙(t)=c1λ1eλ1t+c2λ2eλ2t Evaluating for t0+: {x(0+)=c1+c2+1k=0x˙(0+)=c1λ1+c2λ2=0 This is simply a system of two equations and two unknowns, c1 and c2. The solutions are: {c1=1k(λ1/λ21)c2=1k(λ2/λ11)

Case 2:   λ1=λ2=λ

In this case, the derivative is a bit more complicated: x˙(t)=c1λeλt+c2eλt+c2λteλt Evaluating for t0+: {x(0+)=c1+0+1k=0x˙(0+)=c1λ+c2=0 The solutions to this system of equations are: {c1=1kc2=λk

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 (5) is positive, its roots λ1,2 will be real. If the discriminant is negative, λ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 (5) is exactly zero, and λ1=λ2. This is called the critically damped regime.

Looking at equation (6), you can see that the roots will be real if ζ2>1, and complex if ζ2 <1. In conclusion:

Regime Damping factor Physical constants
Underdamped ζ<1 c<2km
Critically damped ζ=1 c=2km
Overdamped ζ>1 c>2km

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

Different damped harmonic oscillator regimesImage source code

Underdamped regime

We know that in this case, ζ21<0, so we can revise the formula for λ1,2: (13)λ1,2=ζωn±ωnζ21=ζωn±jωn1ζ2=α±jβ, where j is the imaginary unit (j2=1), α and β are real, positive constants defined as (14)αζωnβωn1ζ2.

Substituting this into the solution (equation (11)), and using Euler's formula (ejθ=cosθ+jsinθ): (15)x(t)=c1e(α+jβ)t+c2e(αjβ)t+1ku(t)=eαt(c1ejβt+c2ejβt)+1ku(t)=eαt(c1(cosβt+jsinβt)+c2(cosβtjsinβt))+1ku(t)=eαt((c1+c2)cosβt+j(c1c2)sinβt)+1ku(t)=eαt(c~1cosβt+c~2sinβt)+1ku(t) You could compute the new constants c~1 and c~2 using their definitions (16){c~1=c1+c2c~2=j(c1c2) but it's probably easiest to just solve the system of initial conditions again: (17)x˙(t)=αeαt(c~1cosβt+c~2sinβt)+eαt(c~1βsinβt+c~2βcosβt)=eαt((c~1αc~2β)cosβt+(c~1β+c~2α)sinβt) {x(0+)=c~1+0+1k=0x˙(0+)=c~1α+c~2β=0 (18){c~1=1kc~2=αβk

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αt in equation (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 c~1cosβt+c~2sinβt. To find the value for t at this maximum, tp, we look at the roots of the derivative x˙ (equation (17)): (19)x˙(tp)=0eαtp((c~1αc~2β)cosβtp+(c~1β+c~2α)sinβtp)=0(c~1αc~2β)cosβtp+(c~1β+c~2α)sinβtp=0sinβtp=0 The peak time is the first solution of sinβtp=0 for t>0: (20)tp=πβ

Overshoot

By plugging in this value tp into the expression for x(t), we can determine the amplitude of peak, and by how much it overshoots the limit x()=limtx(t)=1/k. x(tp)=eαπ/β(c~1cosβπβ+c~2sinβπβ)+1k=eαπ/β c~1+1k=eαπ/β 1k+1k Expressed as a percentage, the overshoot is (21)P.O.=100x(tp)x()x()%=100exp(απβ)%=100exp(ζπ1ζ2)%.

images/overshoot.svgImage source code

Rise time

The rise time tr is the time it takes for the step response to reach x() for the first time. x(tr)=x()eαtr(c~1cosβtr+c~2sinβtr)+1k=1kc~1cosβtr+c~2sinβtr=0tanβtr=c~1c~2tanβtr=βα The first solution gives a negative value for tr, so π is added to find the second solution, which is the first solution for tr>0. (22)tr=1β(πatanβα)=1ωn1ζ2(πatan1ζ2ζ)

Since the rise time is inversely proportional to the natural frequency of the system, it's common to eliminate it from the equation: (23)trωn=11ζ2(πatan1ζ2ζ)

images/risetime.svgImage source code

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

Settling time

The settling time ts is the time after which the step response settles in a small error band around x(), never to leave it again. Usually, the width of the error band is chosen to be 2% of x(): |x(ts)x()|x()<1%

In equation (19), we derived that the extrema of the step response occur when sinβt=0, or tx,n=nπ/β. Combining this with equations (15) and (18), we find that the value of x(t) at these extrema is: (24)x(tx,n)=eαtx,n c~1(1)n+1k=1k(1(1)n eαtx,n)=x()(1(1)n eαtx,n)

The amplitude of these extrema can be used to approximate the settling time: (25)eαts1%tslog(0.01)α4.61α=4.61ζωn

The following image visualizes the exponential decay of the extrema, as derived in equation (24).

images/underdamped.svgImage source code

Damped frequency

The period of the sinusoidal factor of the step response (equation (15)) is Td=2πβ. From equation (19), we also know that this is equal to distance between two successive maxima or minima of the step response.

The corresponding frequency Td1 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: (26)ωd=2πTd=β=ωn1ζ2

The damped frequency ω_d in function of the damping factor ζImage source code

Natural frequency

Looking at equation (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 ωn. This explains the choice of parameters at the beginning of this page: the natural frequency ωn is the frequency at which the oscillator would oscillate if all damping were removed.

Overdamped regime

In this case, λ1 and λ2 will be real, and the step response (equation (11)) is made up of two exponential terms and the step function: x(t)=c1eλ1t+c2eλ2t+1ku(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.

images/overdamped.svgImage 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: x˙(0+)=0.

Rise time

Unlike in the underdamped case, the step response of an overdamped oscillator never reaches its final value x(), it approaches it asymptotically. When using the same definition of rise time as in the underdamped case, you could only conclude that tr= 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% of x() to 90% of x().

To compute the rise time, we need to find t1 and t2 such that: (27)x(t1)=0.1 x()x(t2)=0.9 x()tr=t2t1 This results in equations of the form (28)1λ1/λ21eλ1t+1λ2/λ11eλ2t=γ1, where γ is either 0.1 or 0.9.
Note that the coefficients 1λ1/λ21 and 1λ2/λ11 do not depend on ωn, only on the damping factor ζ.

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 ζ increases, λ1 decreases in magnitude, while λ2 increases. The inverses of λ1,2 are the time constants τ1,2 of the exponential terms in x(t). As λ2 increases, the corresponding time constant τ2=1/λ2 becomes very small. This means that the second term c2eλ2t 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 λ1 decreasing and λ2 increasing, c1 approaches 1, and c2 goes to zero, as shown in the following figure.

images/lambdas-overdamped.svgImage source code

This means that if ζ is large enough, the second term of x(t) can be ignored: x(t)c1eλ1t+c3 An approximation of t1 and t2 would therefore be the solution of 1λ1/λ21eλ1t~=γ1. Moving the denominator to the right-hand side and taking the logarithm from both sides: (29)t~1,2=log((γ1)(λ1λ21))λ1=log((γ1)(λ1λ21))ζωn+ωnζ21=1ωnlog((γ1)(λ1λ21))ζ+ζ21. The corresponding rise time is then (30)t~r=t~2t~1=1ωnlog((0.91)(λ1λ21))ζ+ζ211ωnlog((0.11)(λ1λ21))ζ+ζ21=1ωnlog(10.910.1)ζ+ζ21=1ωnlog(9)ζζ21.

We can simplify the denominator of equation (30): 1ζζ21=ζ+ζ21(ζζ21)(ζ+ζ21)=ζ+ζ21ζ2ζζ21+ζζ21ζ2+1=ζ+ζ212ζ(ζ1) In conclusion, the approximation of rise time by ignoring the second exponential is (31)t~r=1ωnlog(9)(ζ+ζ21). When ζ is large, this simplifies to (32)t~r2log(9)ωnζ(ζ1)4.39ωnζ.

To solve equation (28), you could use Newton's or Halley's method, since you can easily write out x˙(t) and x¨(t). However, it's faster to use fixed-point iteration to solve: t=1λ1log(γ11λ2/λ11eλ2tλ1/λ21)

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 ζ.

images/risetime-overdamped.svgImage source code

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

To get a better look at the differences between the true solution and the approximation, we subtract both of them from their asymptote. That results in the following figure.

images/risetime-overdamped-error.svgImage source code

Both curves have a similar shape. A possible approach to refine our approximation for the rise time could be to start from t~r (equation (31)), and vary some of the parameters, minimizing the error between this approximation and the numerical solution. As an additional constraint, their behavior for ζ should be the same, so they both have an asymptote 2log(9)ζ. (33)minimizeδ,c,p,hf(δ,c,p,h)=t^r(ζ;δ,c,p,h)tr(ζ)2=log(9)(ζ+(ζ+δ)pcp)+htr(ζ)2s.t.limζt^r(ζ)2log(9)ζ=0. The previous approximation t~r is a special case of t^r(ζ;δ,c,p,h), with δ=0,c=1,p=2,h=0.
The constraint completely determines the value of h, so only three parameters remain.

For example, using t~r as an initial guess, optimizing over the interval ζ(1,10] yields an optimizer δ=0.764,c=1.732,p=2.178,h=1.678. The result is shown in the following figure. (34)trlog(9)(ζ+(ζ+0.764)2.1781.7322.178)1.678

images/risetime-overdamped-optimize.svgImage source code

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 t1,2 (equation 28). To get the 1% settling time, γ=0.99.

The same kind of approximation is used as well, by ignoring the second exponential term. (35)t~s=1ωnlog(1001λ1/λ2)(ζ+ζ21)2log(100)ωnζ(ζ1)9.21ωnζ

images/settling-time-overdamped.svgImage source code images/settling-time-overdamped-error.svgImage 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.

images/critically-damped.svgImage source code

Rise time

Like the overdamped regime, the 100% 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%90% rise time is defined in the same way as for the overdamped case (equation (27)), and leads to similar transcendental equations: (36)eλt+λteλt=γ1 Rearranging this equation for fixed-point iteration: (37)t=1λlog(γ1λt1) Solving this equation numerically yields (38)tr3.35790856147781ωn.

Settling time

The formulation of the settling time is of same the form of equation (36). The numerical result for the 1% settling time (γ=0.99) is (39)ts6.638352067993811ωn.