# Finite difference for 1D wave equation: why the spike initial data results in a noisy output?

I am using a second-order finite difference in space and time approximation for the 1D wave equation.

No source but initial data: $$I(x)=\mathrm{e}^{-400 (x-0.5)^2}$$. Velocity $$c=1$$, $$nx=501$$, $$nt=1001$$, $$x \in[0,1]$$ , $$t \in[0,1]$$, $$dt=0.001$$.

After 0.2 seconds I get

However, if I change the initial condition to a spike $$I(x)=1$$ for $$x=0.5$$ and $$0$$ for $$x \ne 0.5$$, I get a noisy, low amplitude response that I include next.

Why is my output (response to a spike) so noisy and where did my amplitude go? I expect to see a traveling impulse (perhaps with lots of small wiggles) but distinctive spikes in the middle with amplitude 0.5 each (as in the smooth case).

I believe it is not a bug in my code. I tried other programs from the web and got the same results. I tested the accuracy of my program against an analytic solution and got an error smaller than $$10^{-4}$$. There is not an issue with the CFL condition since my $$dt$$ is small enough for the velocity $$c=1$$ and spatial sampling $$dx$$. I tried sampling rates 10 times smaller. Nothing in particular changed. I believe it is something inherent to the finite difference method but I still do not get it. I have the program in Python and can include it here if necessary.

The issue has nothing to do with boundary conditions because the impulse has not even reached the boundary yet.

Any hint about what could be going on here would be much appreciated.

I think you're probably seeing artifacts that are due to numerical dispersion. In brief, in the discrete case different (spatial) frequencies of a wave function will propagate at different phase/group velocities. This stands in contrast to the continuous case, where the all frequency components travel with exactly the same speed (c). In the first simulation, the initial condition is apparently sufficiently smooth that this effect isn't noticeable (it can be decomposed using only a narrow band of spatial frequencies, that have similar speeds). However, with impulsive initial conditions like the second simulation, every spatial frequency is (equally) present, leading to a "spreading" or dispersion of the signal due to their different propagation speeds.

If memory serves, in 1D you can pick a "magic" timestep, exactly at the courant criterion, and this effect disappears (all spatial frequencies have same velocity c, no spreading/dispersion). Probably worth trying that, even if for no other reason than verifying the claim that this is a numerical dispersion artifact.

Unfortunately there is no such "magic" timestep in 2D or 3D. In fact, this problem is even more exaggerated in higher ND, because there is not only this chromatic dispersion, but also dispersion anisotropy (waves going in different directions exhibit different velocities). High order discretizations (in space and time) are the best way to combat these dispersion problems.

• OMG you are a life saver. Thanks so muuuuchhh! I have been struggling with this for two weeks. I made CFL=1 and Yes! the spike splits into two clean spikes with half amplitude traveling in opposite directions. It is weird that when they meet again, the amplitude does not add up to 1 as in the Gaussian example. This is a minor thing that I still need to fix (any ideas?). Again, you hit the point exactly right. Now I need to study "dispersion" to understand better this. – Herman Jaramillo Oct 4 '18 at 2:21
• You can derive the dispersion relationship using von-Neumann-like techniques: impose a planewave solution (with an unknown wavenumber) upon the samples, then impose periodic boundary conditions in space and time, and solve the resulting eigenproblem to determine what wavenumber has to be, in order satisfy the FD relations. See csc.kth.se/utbildning/kth/kurser/DN2255/ndiff12/Lecture5.pdf, section 3.2 (starts on page 9) for ideas. A mentor can also walk you through this, or feel free ask another detailed question about the reference if you get stuck. – rchilton1980 Oct 5 '18 at 15:19
• I was thinking about taking the Fourier transform of the finite difference equation. This will transform my $u(x,t)$ into $U(k,\omega)$ (which is the dispersion relation) then I can find the group $\partial \omega/\partial k$ and phase velocities $\omega/k$ and see why your "magic" number is CFL=1. Does this make sense? Thanks. – Herman Jaramillo Oct 5 '18 at 16:16
• Yep it does, it's basically the same procedure. – rchilton1980 Oct 5 '18 at 17:16
• "High order discretizations (in space and time) are the best way to combat these dispersion problems." There is some truth to that, but the bottom line is that you need to resolve your initial data better. – David Ketcheson Oct 15 '18 at 7:03

The credits of this should go to @rchilton1980. I only want to validate the hypothesis of dispersion as the result of the response of a spike to the finite difference operator.

$$\bf{ Basic \; Theory: }$$
Physical dispersion accounts for the phenomena of waves moving in a way that different frequencies or wavenumbers travel at different velocities. We can make an analogy with a group of cliclist traveling along a flat road at constant speed. Once they hit the mountain it could happen that the heaviest cliclists will ride slower and the group will spread (or disperse) along the road. If we think of a wavelet as the supperposition (synthesis) of many frequencies (and wavenumbers) it could happen that the velocity of each frequency is different. However we need to be more precise on the meaning of velocity here. We define two velocities:

$$\bf{Phase \; velocity}$$

$$$$v_f= \frac{\omega}{k}$$$$

$$\bf { Group \; Velocity: }$$

$$$$v_g= \frac{\partial \omega}{\partial k}$$$$

The previous definitions indicate, in a way, that the angular frequency $$\omega$$ should be a function of the wavenumber $$k$$. That is, $$\omega=\omega(k)$$. To better understand this let us consider the 1D wave equation with constant wavespeed $$c$$.

$$\begin{eqnarray*} \frac{\partial^2 u}{\partial x^2} -\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} = 0. \end{eqnarray*}$$ If we take a double Fourier transform (from time $$t$$ to frequency $$\omega$$, and from space $$x$$ to wavenumber $$k$$) we find

$$\begin{eqnarray*} \left ( -k^2 + \frac{\omega^2}{c^2} \right ) U(k,\omega) = 0, \end{eqnarray*}$$ where $$U(k,\omega)$$ is the Fourier transform of $$u(x,t)$$. We get, from here

$$\begin{eqnarray*} \omega = \pm c k \end{eqnarray*}$$ The "$$+$$" sign indicates a wave traveling to the right, while The $$-$$'' sign indicates a wave traveling to the left. From now on, we will ignore the sign on this discussion. Then

$$\begin{eqnarray*} \frac{\omega}{k} = \frac{\partial \omega}{\partial k} = c \end{eqnarray*}$$ So we find that both phase and group velocity coincide with the wavespeed $$c$$. So there is no physical dispersion.

However, when we find a finite difference approximation we could introduce dispersion (numerical dispersion). That is, both phase and group velocities could change with frequency (or wave number).

Let us consider the following finite difference scheme for the wave equation without source.

$$\begin{eqnarray} w_{i j+1} = \gamma^2 w_{i-1 j} + 2(1-\gamma^2) w_{ij} + \gamma^2 w_{i+1 j} - w_{i j-1}, \end{eqnarray}$$ where $$\gamma$$ is the CFL number $$c \Delta t/\Delta x$$.

We take the Fourier transform of this expression, and using the shift property of the Fourier transform to find

$$\begin{eqnarray*} W_{ij} \left ( \mathrm{e}^{\mathrm{i} \omega \Delta t} + \mathrm{e}^{-\mathrm{i} \omega \Delta t} \right ) = \left ( \gamma^2 \mathrm{e}^{-\mathrm{i} k \Delta x} + 2(1 - \gamma^2) W_{ij} + \gamma^2 \mathrm{e}^{\mathrm{i} k \Delta x} W_{ij} \right ) \end{eqnarray*}$$ where $$W_{ij}$$ is the Fourier transform of $$w_{ij}$$. That is, we find

$$\begin{eqnarray*} \left [ 2 \cos \omega \Delta t - 2 \gamma^2 \cos k \Delta x + 2(1 - \gamma^2) \right ] W_{ij} = 0. \end{eqnarray*}$$ We isolate $$\omega$$ from here so that

$$\begin{eqnarray*} \omega \Delta t = \arccos \left ( \gamma^2 \cos k \Delta x - (1 - \gamma^2) \right ), \end{eqnarray*}$$ That is, $$\begin{eqnarray*} v_f = \frac{\omega}{k} = \frac{1}{k \Delta t} \arccos \left ( \gamma^2 \cos k \Delta x - (1 - \gamma^2) \right ), \end{eqnarray*}$$ Clearly $$v_f$$ changes with wavenumber $$k$$. That is there is numerical dispersion. The group velocity is given by

$$\begin{eqnarray*} \frac{\partial \omega}{\partial k} = -\frac{1}{\Delta t} \frac{ -\gamma^2 \sin k \Delta x} { \sqrt{1 - \gamma^2 \cos k \Delta x - (1 - \gamma^2) }}. \end{eqnarray*}$$ Again, here group velocity $$v_g$$ changes with $$k$$

There is a special case given by $$\gamma=1$$. Here

$$\begin{eqnarray*} w \Delta t = \arccos(\cos k \Delta x)= k \Delta x \end{eqnarray*}$$ or

$$\begin{eqnarray*} \frac{\omega}{k} = \frac{\Delta x}{\Delta t} \\ \frac{\partial \omega}{\partial k} = \frac{\Delta x}{\Delta t} \\ \end{eqnarray*}$$ and since $$\gamma=c \Delta t/\Delta x =1$$, then $$c=\Delta x/\Delta t$$ and the phase and group velocities coincide with $$c$$. So there is no dispersion.

If $$\gamma < 1$$ then we can expect dispersion as shown in the original question.

Here I include the plots after a few fractions of second of propagation.

You can only expect to get a sensible answer if your data can be resolved by the grid. Shannons Sampling Theorem is central to all digital signal processing and states that no signal with wavelength less than 2h can be represented on a grid with spacing h (and for most methods the minimum wavelength is much more than this). Your spike is severely underresolved.

The magic CFL number, unfortunately does not work for nonlinear problems with variable wavespeed or as rchilton observes, in higher dimensions, so it is of limited practical use use

• @Phillip Roe : Thanks for your input. What is the relation between the sampling theorem and dispersion? In other words, can we say that aliasing and dispersion are the same thing? The concepts are totally different because one is related to sampling while the other is related to group and phase velocities. However they seem to be connected. How? Thanks. – Herman Jaramillo Oct 5 '18 at 13:17
• Im not sure what connections you are seeing. And I dont think you ever said exactly what method you are using. Second-order finite-difference in space and time sounds like "leapfrog" $(u_j^{n+1}-u_j^{n-1})+\nu(u_{j+1}+^n-u_{j-1}^n)=0$ Yes? This is not the most typical FD scheme. In general, a scheme has both dissipation (wrong amplitude) and dispersion (wrong phase) The leapfrog scheme is symmetrical and has no dissipation. All of the different frequencies that are the Fourier components of your spike propagate without damping at different speeds to create the mess that you generated. – Philip Roe Oct 6 '18 at 15:56