Consider the following ODE over $\mathbb{R}^+$ with $x(t)\in\mathbb{R}^n$ for large $n$:

\begin{equation} \ddot x + K x = 0, \quad K=\begin{bmatrix} 2 & -1 &\dots \\ -1 & 2 &-1&\dots \\ \vdots & \ddots & \ddots & \ddots \\ 0 & \dots & -1 & 2 & -1 \\ 0 &\dots &\dots & -1 & 1\end{bmatrix} \end{equation}

Physically, $x$ can be understood as the displacement of $n$ nodes of a spring-mass chain, fixed at one end (note that $K_{n,n}=1$, that's not a typo). This is what happens for an initial condition $x(0)=[0,\dots,0]$ and $\dot x(0)=[0,\dots,0,1]$, i.e. for a initial velocity of the mass at the end.

enter image description here

It appears that there is something looking like a propagating wave (going to the left), with a constant velocity.

Question How to estimate this apparent velocity as a function of $n$?

Numerically, it is easy to estimate it, for example by plotting $(i,t,x_i(t))$ where $i\in \{ 1,\dots,n\}$, here for $n=100$:

enter image description here

We can clearly see a straight line in the plot $(i,t)$. On the left of this line, $x$ is (almost) zero: there is no displacement (information has not reached these points yet).

Note that the eigenvalues of $K$ are $$\omega_k=2\Big(1-\cos\Big(\dfrac{(2k-1)\pi}{2n+1}\Big)\Big)$$ and its eigenvectors are given by $$\phi_k^{(i)}=\sin\Big(\dfrac{i(2k-1)\pi}{2n+1}\Big)$$ where $i$ is the component index and $k$ the index of the eigenvector.

The ODE can be solved in the first order form, posing $y^\top = [x,\ \dot x]$:

\begin{equation} \dot y = Ay, \quad A = \begin{bmatrix} 0 & I \\ -K & 0 \end{bmatrix} \end{equation} so the solutions are $y(t)=\exp(tA)y(0)$. The question reduces to "extracting" the apparent velocity from $\exp(At)e_{2n}$ ($e_{2n}$ is the last vector of the canonical basis) but I can't see how. Digging in that direction, we can show that with $$q_i(t)=\dfrac{1}{\omega_i} \phi_i^{(n)} \sin(\omega_i t)$$ and $q^\top =[q_1,\dots,q_n]$, the solution for the initial velocity $y=e_{2n}$ is $$ x(t)=\phi^{-1} q(t)$$ so the information in inside $\phi^{-1}q$.

  • $\begingroup$ If the propagation speed really is finite, then the boundary conditions don't matter, and you can replace $x$ with a doubly-infinite sequence $\mathbb Z\to\mathbb R$ satisfying $\ddot x_i + (-x_{i-1} + 2x_i - x_{i+1}) = 0$. Then it's easier to take the Fourier transform and analyze the behavior of modes of different wavenumbers. $\endgroup$ – Rahul Sep 2 '17 at 20:02
  • $\begingroup$ @Rahul Edited to include the eigenvalues and eigenvectors. What do you mean by "analyze the behavior"? $\endgroup$ – anderstood Sep 2 '17 at 23:29
  • $\begingroup$ Seeing that your $K$ is basically the discrete Laplacian, I would've expected the speed to be "1"; so I am somewhat confused by the result of your numerical simulation. $\endgroup$ – Willie Wong Sep 3 '17 at 1:50
  • $\begingroup$ Also, you will not be able to see the finite speed information from the eigenvectors and eigenvalues, the same way that you cannot see the finite speed information from the Fourier transform of the linear wave equation. $\endgroup$ – Willie Wong Sep 3 '17 at 1:53
  • $\begingroup$ @WillieWong As I understand, you guessed the velocity by recognizing the discretized wave equation $\ddot x = \Delta x$. What if you had never heard about it? I mean, the "information" should be somewhere in the ODE?! $\endgroup$ – anderstood Sep 3 '17 at 2:03

In the $i$th row of the differential system, the term $(x_{i+1} - 2 x_{i} + x_{i-1})/\Delta\xi^2$ is an order-2 central finite-difference approximation of the space derivative $\partial^2 x/\partial \xi^2$, where $\xi$ is a space coordinate such that $x_i(t) \simeq x(i\, \Delta \xi ,t)$. Thus, the differential system $\ddot{x_i} + K x_i = 0$ may be viewed as a finite-difference spatial discretization of the wave equation $$ \frac{\partial^2 x}{\partial t^2} - c^2 \frac{\partial^2 x}{\partial \xi^2} = 0 \, , $$ which speed of sound in $\xi$-$t$ coordinates is $c = \Delta \xi$ /s.

Now, let us assume that $x = \exp\left({\text{i}(\omega t - k \xi)}\right)$ is a monochromatic wave. Injecting this Ansatz in the wave equation, we obtain the expression of the physical wave number $k = \omega/c$, i.e. the dispersion relation. Injecting the Ansatz $x = \exp\left({\text{i}(\omega t - \kappa \xi)}\right)$ in the corresponding discrete equation $$ \ddot{x}_i - c^2 \frac{x_{i+1} - 2 x_{i} + x_{i-1}}{{\Delta\xi}^2} = 0 \, , $$ we obtain the relation $$ \kappa \Delta\xi = \arccos\left( 1 - \frac{(k\Delta\xi)^2}{2} \right) $$ satisfied by the numerical wave number $\kappa$. A series expansion as $k\Delta\xi \to 0$ gives the numerical dispersion relation $$ \kappa \simeq k + \frac{k^3 \Delta\xi^2}{24} + O\left(k^4 \Delta\xi^3\right) . $$


We may decompose $x(t)$ into a linear combination of eigenvectors, $x(t)=\sum_{k=1}^n\alpha_k(t)\phi_k$.

Each coefficient satisfies $\ddot\alpha_k(t)+\omega_k\alpha_k(t)=0$, so $\alpha_k(t)=a_k\exp(j\sqrt{\omega_k}t)+b_k\exp(-j\sqrt{\omega_k}t)$ for some complex $a_k,b_k$. I'm using $j=\sqrt{-1}$ because we are already using $i$ for the spatial index.

Similarly, each eigenvector can be written as $\phi_k^{(i)}=\frac1j\bigl(\exp(j\kappa_ki)-j\exp(-j\kappa_ki)\bigr)$, where $\kappa_k=\pi(2k-1)/(2n+1)$.

Multiplying them together, each eigenmode $\alpha_k(t)\phi_k^{(i)}$ can be expressed as a sum of terms of the form $c\exp\bigl(j(\pm\sqrt{\omega_k}t\pm\kappa_ki)\bigr)$, or equivalently, $c\exp\bigl(j\kappa_k\bigl(\pm i\pm(\sqrt{\omega_k}/\kappa_k)t\bigr)\bigr)$. Thus each eigenmode is a superposition of travelling waves of wavenumber $\kappa_k$ and velocity $\pm\sqrt{\omega_k}/\kappa_k$.

In particular, for $k\ll n$, we have $\sqrt{\omega_k}/\kappa_k\approx1$, consistent with the analysis from the wave equation. However, there is some dispersion at higher wavenumbers, which explains the trailing high-frequency oscillations.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.