Is it possible to write explicitly the exact solution for forced damped harmonic oscillator?

Question

Preamble

Consider a damped harmonic oscillator, with his well know differential equation \begin{equation*} m \ddot{x} + c \dot{x} + kx=0 \end{equation*} and let's find the solution that satisfies $x(0)=x_0$ and $\dot{x} (0) = v_0$. The problem admit a lovely exact solution (exploiting the considerations I found in books, I derived the complete explicit solution I write below). We have three cases. In the particular case in which $c=2\sqrt{mk}$, we have \begin{equation*} x(t) = \left[ x_0 + \left( v_0 + x_0 \sqrt{\frac{k}{m}} \right) t \right] e^{- \sqrt{\frac{k}{m}} t} \end{equation*} Otherwise we have to consider separately the two cases $c \gtrless 2\sqrt{mk}$, and to don't make heavier equations we have to introduce this parameters (we can observe that $\alpha$ and $\beta$ are length and depend to initial conditions too, while $\gamma$ e $\xi$ are adimensional and depend only by physical characteristics of the system) \begin{equation*} \alpha \equiv \frac{x_0}{2} \qquad \beta \equiv \frac{m v_0}{c} \qquad \gamma \equiv \frac{4mk}{c^2} \qquad \xi \equiv \frac{ct}{2m} \end{equation*} In the case $c>2\sqrt{mk}$ (i.e. $\gamma<1$) we have \begin{equation*} \begin{split} \displaystyle x(t) = & \frac{ \alpha \left(\sqrt{1 - \gamma } + 1 \right) + \beta}{\sqrt{1 - \gamma }} \cdot \exp \left[ {\left({-1 + \sqrt{1 - \gamma}} \right) \cdot \xi} \right] + \\ \displaystyle &\frac{ \alpha \left(\sqrt{1 - \gamma } - 1 \right) - \beta}{\sqrt{1 - \gamma }} \cdot \exp \left[ {\left({-1 - \sqrt{1 - \gamma}} \right) \cdot \xi} \right] \end{split} \end{equation*} While if $c<2\sqrt{mk}$ (i.e. $\gamma>1$) the solution is \begin{equation*} x(t) = 2 \alpha \cdot \frac{{\sin \left[ \tan^{-1} \left( \frac{\sqrt{\gamma - 1}}{1+\frac{\beta}{\alpha}} \right) + \sqrt{\gamma - 1} \cdot \xi \right]}} {\sin \left[ \tan^{-1} \left( \frac{\sqrt{\gamma - 1}}{1+\frac{\beta}{\alpha}} \right) \right]} \cdot {\exp(-\xi)} \end{equation*} Note that in this expressions, time dependence is included in $\xi$.

Question

Is it possible to do the same with the sinusoidal forced case? We have \begin{equation*} m \ddot{x} + c \dot{x} + kx= F \cos(\omega t + \phi) \end{equation*} but what is the solution of this differential equation that satisfies the condition $x(0)=x_0$ and $\dot{x} (0) = v_0$? I can't write it, even in the simpler case with $\phi = 0$.

Answer 1

You just need to find a particular solution of that equation, then general solution will be a linear combination of the particular solution and the general solution of the homogeneous differential equation. Well, let $\phi = 0$, to simplify.

Let a complex equation: $m\ddot z + c\dot z + kz = Fe^{i\omega t}$. Your equation is just the real part of this. If you make the substitution $z(t) = Ae^{i\omega t}$ in the differential equation, you will find the value of $A$.

Now you can write $A$ this way: $A = |A|e^{i\theta}$. Now your particular solution is just the real part of $z(t)$: $$ x_p(t) = |A|\cos(\omega t + \theta). $$

Therefore, you can have a particular solution explicity using the values of $|A|$ and $\theta$ you have found. Now, do a linear combination of the general solution, and plug the initial condition values for the initial position and initial velocity. Then you have the general solution you are looking for.

Answer 2

I found expressions for the solution, and I write them below. With the help of a software I checked that they really are the solutions, I'm 100% sure: they solve the differential equation \begin{equation*} m \ddot{x} + c \dot{x} + kx= F \cos(\omega t ) \end{equation*} and satisfy the initial conditions (initial position $x_0$ and initial speed $v_0$). But obviously the job isn't done, for the simple reason that the solution I write are really ugly: it must be possible to write this formulas in a better way. We should do it for aesthetic reasons, obouvisly, but also to make them more easy to handle, and so more useful. The beautiful work I did with the not forced case (see the question) allow users to find easily the equation of motion of the particle (exploring limiting cases too is simple). I tried to do the same with the forced case but I wave the white flag, I can't do it (I won't work more about this problem). I hope someone will do. Anyway I write out the complete solution (althoung in this ugly form).

Case 1: $c=2\sqrt{m k}$

\begin{equation*} x(t) = \frac{e^{-\frac{c t}{2 m}} \left( A+B+C \right) }{B_2} \end{equation*} where \begin{equation*} \begin{split} &A= \left( A_1 + A_2 \right) F \\ &B= \left( B_1 + B_2 \right) x_0 \\ &C= 2 m (m \omega^2 + k)^2 t v_0 \\ &A_1= 2 m e^{\frac{c t}{2m}} \left( c \omega \sin \left( \omega t\right) -m \omega^2 \cos\left( \omega t\right) +k \cos\left( \omega t\right) \right) \\ &A_2 = -c \left( m \omega^2 + k \right) t+2 m^2 \omega^2-2 k m \\ &B_1 = c {\left( m \omega^2+k\right) }^2 t \\ &B_2 = 2 m ( m \omega^2 + k )^2 \end{split} \end{equation*}

Case 2: $c>2\sqrt{m k}$

\begin{equation*} x(t) = \frac{e^{-\frac{U}{2 m}-\frac{c t}{2 m}} \left( x_0 W \left( e^{\frac{U}{m}} +1\right) +\sqrt{c^2-4 k m} \left( F L +v_0 E+x_0 D\right) +F G\right) }{2W} \end{equation*} where \begin{equation*} \begin{split} & L = \left( c m \omega^2+c k\right) e^{\frac{U}{m}}-c m \omega^2-c k \\ & D = cV \left( 1-e^{\frac{U}{m}} \right) \\ & E = 2mV \left( 1-e^{\frac{U}{m}} \right) \\ & G = K e^{\frac{U}{m}}+K+G_2+G_1 \\ & K = O \left( m \omega^2-k\right) \\ & W = O V \\ & G_1 = 2 c O \omega \sin\left( \omega t\right) Z \\ & G_2 = -2 O \left( m \omega^2-k\right) \cos \left( \omega t\right) Z \\ & O = 4 k m-c^2 \\ & Z = e^{\frac{U}{2 m}+\frac{c t}{2 m}} \\ & U = \sqrt{c^2-4 k m} t \\ & V = {m}^2 \omega^{4}-2 k m \omega^2+c^2 \omega^2+{k}^2 \end{split} \end{equation*}

Case 3: $c<2\sqrt{m k}$

\begin{equation*} x(t) = \frac{e^{-\frac{c t}{2 m}} \left( S+R+\sqrt{4 k m-c^2} \left( Q+P \right) +N \right) }{\sqrt{4 k m-c^2} \left( m^2 \omega^4+\left( c^2-2 k m \right) \omega^2+k^2 \right) } \end{equation*} where \begin{equation*} \begin{split} P=&\left( c \omega e^{\frac{c t}{2 m}} \sin \left( \omega t \right) +\left( k-m \omega^2 \right) e^{\frac{c t}{2 m}} \cos \left( \omega t \right) +\left( m \omega^2-k \right) \cos \left( \frac{\sqrt{4 k m-c^2} t}{2 m} \right) \right) F \\ Q=&\left( m^2 \omega^4+\left( c^2-2 k m \right) \omega^2+k^2 \right) \cos \left( \frac{\sqrt{4 k m-c^2} t}{2 m} \right) x_0 \\ R=&-c\left( m \omega^2 + k \right) \sin \left( \frac{\sqrt{4 k m-c^2} t}{2 m} \right) F \\ S=&c \left( m^2 \omega^4+\left( c^2-2 k m \right) \omega^2+ k^2 \right) \sin \left( \frac{\sqrt{4 k m-c^2} t}{2 m} \right) x_0 \\ N=& m \left( 2 m^2 \omega^4+\left( 2 c^2 -4 k m \right) \omega^2+2 k^2 \right) \sin \left( \frac{\sqrt{4 k m-c^2} t}{2 m} \right) v_0 \end{split} \end{equation*}

An example

This work look a mess but it isn't: it allow to get elegant solution for physical problems (provided you have a computer to help you, if expressions won't be simplified, as I hope someone will do, calculus by hands will be prohibitive!). But let's consider a concrete example. Suppose we have exactly (I don't transcribe units, simply suppose they are consistent) $c=1$, $k=2$, $\omega = 3$, $F=4$, $x_0 = 5$, $v_0 = 6$, $m=2{.}5$ we get this equation of motion (we are in the third case) \begin{equation*} \begin{split} x(t)=& \frac{59703 e^{-\frac{t}{5}} \sin \left( \frac{\sqrt{19} t}{5} \right) }{1717 \sqrt{19}}+\frac{8913 e^{-\frac{t}{5}} \cos \left( \frac{\sqrt{19} t}{5} \right) }{1717}+\\ & \frac{48 \sin \left( 3 t \right) }{1717}-\frac{328 \cos \left( 3 t \right) }{1717} \end{split} \end{equation*} We see that when $t$ is big (after transient) the motion is approximately given by \begin{equation*} x(t) = -\frac{8 }{\sqrt{1717}} \cos \left( 3 t+\tan^{-1} \left( \frac{6}{41}\right) \right) \end{equation*} I obtained this from $x(t) = \frac{48 \sin \left( 3 t \right) }{1717}-\frac{328 \cos \left( 3 t \right) }{1717}$, but looking the formula $x(t)$ of case 3 we can write a general formula for this motion after transient: I find strange that it is completely independent by initial condition: \begin{equation*} x_{at}(t) = \frac{\left( c \omega \sin \left( \omega t\right) +\left( k-m {\omega}^2\right) \cos \left( \omega t\right) \right) F}{m^2 {\omega}^{4}+\left( {c}^2 -2 k m\right) {\omega}^2 +{k}^2 } \end{equation*} But we know exactly all the motion, not only this limiting state: here I plot $x=x(t)$

If we wait in doing the final substitution $m=2{.}5$, we can get the two variable function that tell us how varies the motion if we vary the mass $m$ (the previous graph is the section of this below at $m=2{.}5$ and with $t$ until 25):

With derivative of $x(t)$ we can find speed and acceleration as function of $t$:

In short, knowing analytic solution allow us to find exactly any kind of information: for example we can find that after 8 seconds (or unit of time you are using) the position is about $1{.}727745920373746$ and after 15 seconds acceleration is $0{.}3102092667738571$, or more precisely \begin{equation*} \frac{72 \left( 41 \cos \left( 45\right) - 6 \sin \left( 45\right) \right) }{1717} - \frac{24 \left( 6133 \sin \left( 3 \sqrt{19}\right) +2332 \sqrt{19} \cos \left( 3 \sqrt{19}\right) \right) }{8585 \sqrt{19} e^3 } \end{equation*} We can also find exactly other kind of information, for example we can estimate (necessarily numerically in this case, but only because we have transcendental equations: errors are very small and we can estimate them) that the particle meet the $x=0$ position, for the first time, when $t \approx 2{.}9818321933898$ (and for the first time it is at rest when $t=0{.}90687408051175$). I defy anyone to find better data using a numerical ODE solver: with my formulas I can find exact position at any time, while solving ODE numerically give more approximate data than find the value of a formula (and it is not easy estimate the error). Of course all this digit have no sense in a real physical problem (physical constant and conditions are not exactly know, linear law for friction is only a rough approximation) but I feel that exact analytic solution of a model (a physical model with all its limits, like every other one) is very gratifying (what about the artificial but beautiful problem we can find in electrostatic for example?). I hope someone will write my exact but ugly solution, in the right way (it MUST exists a RIGHT way to write it).