# ODE $x''(t)+\eta x'(t)+x(t)=0$ with the $\eta$ extremely small

I want to numerically solve the damped oscillator equation

$$\frac{d^2x}{dt^2}+\eta\frac{dx}{dt}+x=0 .$$

Usually it's very easy to solve this kind of ODE analytically and numerically. But now the damping coefficient is extremely small, $\eta=10^{-10}$. To get the obvious damping pattern, the value of oscillating time $t$ should be bigger than $10^{10}$. I use Mathematica to solve it. It can solve from $t=0$ to $t=10^6$ easily, but it's still far from what I want to get. Calculating for a large domain may be just impossible in principle. (I think this difficulty is universal for other languages and softwares like Fortran, C, etc.)

I try to rescale the domain, like set $t=e^\tau$ and then modify the ODE correspondingly, so the domain becomes $\tau=[0,10]$ rather than $t=[0,10^{10}]$. But this method fails. It turns out that it will cost the same amount of time.

Does anyone have some good ideas or similar experience to handle this problem?

PS: Actually the ODE I want to calculate is more complicated than this damped oscillator, and it cannot be solved analytically, so I have to numerically solve it. But these two ODEs have the same structure.

• I added the tag finite-difference, assuming that Mathematica use this method for the numerical solution, is that right? – nicoguaro Feb 24 '17 at 4:52
• I believe this is the sort of problem for which multiscale methods were invented. Otherwise, the timestep can only be as large as the oscillation period at most, to resolve the oscillation. – Kirill Feb 24 '17 at 4:57
• @ I seldom use numerical methods. I never heard "multiscale method" before. I googled it and find some PDFs that may be helpful to my problem. Thank you. My problem is exactly a multiscale problem. – youyou Feb 24 '17 at 5:44
• Look into exponentially fitted methods. – David Ketcheson Feb 25 '17 at 6:26
• (1) For this particular equation, you know the period and the exponential decay of the amplitude (Wikipedia, damped oscillator), so even without solving the system we can sketch how the solution looks like. (2) For another system, if you want to solve every single oscillation, I see no other way than to increase the numerical precision and use something with more significant digits than 64bit floats. – Sampo Smolander Feb 28 '17 at 19:21

The problem you have is what is called "stiff": it has two time scales, namely one for the oscillation (which is ${\cal O}(1)$) and one for the damping (which is ${\cal O}(\eta^{-1})$). If $\eta$ is very small, then these time scales are very different.

Stiff ODEs are difficult to solve because in order to solve them accurately, you have to resolve the short time scale with the time step, which yields very large number of time steps to resolve the long time scale. You can, however, use an implicit solver with long time steps (sufficient to resolve the long time scale) and get a qualitatively correct answer (although it will no longer show the oscillatory behavior).

There is a large literature on stiff ODEs and solution approaches. You will find much information if you look for it.

• I hadn't thought of it this way before, so correct me if I'm wrong, but I thought even an implicit solver for stiff problems would still need to resolve individual oscillations. For the typical stiff system, the r.h.s. eigenvalues are large, negative, and real, as opposed to large and almost pure-imaginary. Indeed that's why I mentioned multiscale methods up in my comment, because I kind of thought that's why they were necessary, as opposed to just using the usual methods for stiff problems. – Kirill Feb 24 '17 at 21:42
• Yes, implicit methods are going to remain stable. But just because the numerical solution you get is stable doesn't mean that it's accurate. – Wolfgang Bangerth Feb 26 '17 at 16:43
• I would look at exponential time differencing (Cox & Matthews 2001 for starters). They treat linear terms with exact integration which makes it very well suited for stiff systems. – KyleW Feb 28 '17 at 5:06
• Stiff methods solve the long term behavior, but they don't solve the fast, short term transients correctly. It's just that for the long term solution, getting the exact form and shape of the fast transients correctly does not matter, for the class of problems stiff methods are used. If OP wants to get every single oscillation solved correctly, then stiff methods won't produce a correct solution. – Sampo Smolander Feb 28 '17 at 19:25

Let us write the equation in state-space form

$$\frac{d}{dt}\begin{bmatrix}x_{1}(t)\\ x_{2}(t) \end{bmatrix}+\begin{bmatrix}\eta/2 & -\omega\\ \omega & \eta/2 \end{bmatrix}\begin{bmatrix}x_{1}(t)\\ x_{2}(t) \end{bmatrix}=0 \tag1$$

where the natural frequency is $\omega \triangleq \sqrt{1-\eta^2/4}$. It is easy to verify that the characteristic polynomial is indeed $p(\lambda) = \lambda^2 + \eta \lambda + 1$ by construction. When $\eta \ll 1$, we have $\omega \approx 1$, so the solution can always be approximated to great accuracy as

$$\begin{bmatrix}x_{1}(t)\\ x_{2}(t) \end{bmatrix}=\exp\left(-t\begin{bmatrix}0 & -1\\ 1 & 0 \end{bmatrix}\right)\begin{bmatrix}y_{1}\\ y_{2} \end{bmatrix}\tag2$$

with $y_{1}=x_{1}(0), y_{2}=x_{2}(0)$. Indeed this is the exact solution when $\eta=0$. Now, with a finite $\eta$, this solution is not completely correct, so we let $y_{1},y_{2}$ vary with time to yield

$$\begin{bmatrix}x_{1}(t)\\ x_{2}(t) \end{bmatrix}=\exp\left(-t\begin{bmatrix}0 & -1\\ 1 & 0 \end{bmatrix}\right)\begin{bmatrix}y_{1}(t)\\ y_{2}(t) \end{bmatrix},\tag3$$

with $y_{1}(0)=x_{1}(0)$, $y_{2}(0)=x_{2}(0)$. These variables can be intuitively understood as “corrections” to the initial guess above. Then using the laws of matrix exponentials, we derive an ODE for the correction terms

$$\frac{d}{dt}\begin{bmatrix}y_{1}(t)\\ y_{2}(t) \end{bmatrix}+\begin{bmatrix}\eta/2 & -(\omega-1)\\ (\omega-1) & \eta/2 \end{bmatrix}\begin{bmatrix}y_{1}(t)\\ y_{2}(t) \end{bmatrix}=0.\tag4$$

This latter ODE is non-stiff, and can be solved using any simple integration method. With $y_{1}(t)$ and $y_{2}(t)$ known, $x(t)$ is easily recovered.

To summarize, the principal idea is to start with an excellent initial guess, and to write the true solution as a time-dependent correction multiplied by this initial guess. Solving for the correction is then a far easier problem (read: less stiff) than solving for the solution directly.

The particularly strategy above is known as an “exponential integrator” (more specifically, the Lawson type), and assumes that the nonlinear dynamics can be well approximated by a linear one. For this very specific problem of a damped oscillation, it is also known as a phasor transform.

• Are there nonlinear dynamics involved? All the equations seem to be linear – DanielRch Feb 28 '17 at 15:36
• OP posted a linear ODE but alluded to the "real" problem being nonlinear. – Richard Zhang Feb 28 '17 at 15:40
• Thank you for the answer. I have a question about how you go from (3) to (4). I simply differentiate (3) by $t$, but I cannot get (4) because the matrix {{0,-1}{1,0}} does not commute the matrix {{0,0}{0,$\eta$}} – youyou Mar 3 '17 at 0:27
• Yikes, you're completely right. I will fix this (at the risk of making the solution look a bit artificial), but for real life problems, I would suggest the following reference: journals.cambridge.org/article_S0962492910000048 – Richard Zhang Mar 3 '17 at 20:00
• While the previous answer was not technically correct, you will find that since the commutator is so small (in any matrix norm), it still works really well. In practice, physicists and engineers would freely assume that the commutator is exactly zero, though that does make mathematicians a little uneasy :-) – Richard Zhang Mar 3 '17 at 20:13

I would give spectral methods a trial. For your one-dimensional equation

$$\frac{d^2x}{dt^2}+ \eta \frac{dx}{dt} + x =0\,,$$

this basically means the following:

1. Choose some basis, e.g. some Lagrange polynomials over a grid.
2. Calculate the matrix elements of the derivatives (and possibly other terms).
3. Diagonalize the matrix (which is possible for the function above, as it's symmetric), and use the eigenvalues and eigenvectors to construct the matrix exponential, with which one can propagate the solution to very long times.

Points 1. and 2. are identical as if you would use an explicit or implicit method. However, point 3. gives you the exact result for the chosen discretization, by which you can take very large steps. This is in sharp contrast to the mentioned explicit and implicit methods (which so to say only approximate the matrix exponential).

Using spectral methods, the stiffness then appears through the fact that the eigenvalues differ by some orders of magnitude.

• Thank you for the suggestion. The ODE I want to solve is in fact more complicated:$$\frac{d^2x}{dt^2}+ \eta f(x) \frac{dx}{dt} + g(x) =0\,,$$ with $f(x)\sim\frac{1}{\sqrt{x}}$ and $g(x)\sim\sqrt{x}$. Is the spectral method general enough to solve this kind of equation? – youyou Mar 4 '17 at 2:09
• As long as it's one-dimensional, no problem. Just take the usual description you already have -- if it's not all too large, it should work quite well. What is your current discretization, by the way? – davidhigh Mar 4 '17 at 10:24