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

Question

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.

Answer 1

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.

Answer 2

Let us write the equation in state-space form

\begin{equation} \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 \end{equation}

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{equation} \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 \end{equation}

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{equation} \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 \end{equation}

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

\begin{equation} \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 \end{equation}

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.

Answer 3

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.