The second equation defines the Laplacian, which is just a second spatial derivative of the scalar field P(x,y,z) (partial derivatives over the spatial axes only, ignoring the time t axis). The Laplacian appears in our differential equation for heat flow, and it also appears in differential equations for wave propagation, called wave equations.
Now, one way to solve such differential equations is to find a closed-form analytic expression that solves a simplified case of the differential equation. We looked at an analytic solution for the heat flow equation in one dimension. The way to solve these equations without simplifying them is to finite difference them. This works well for fields that are discretized and exist as discrete samples in computer memory. We start with the definition of the analytic derivative, the third equation at right.
The fourth equation simply removes the limit and keeps dx finite, or greater than zero. The bold Dx is an operator, one that approximates the analytic partial-derivative operator. It says ``apply the single finite difference in x to the field P(x)'' that may vary in the x direction.
The fifth equation shows how to derive the second finite-difference operator Dx2, which is just applying Dx twice to take the difference of the difference and approximate the second spatial derivative.
The sixth line is the result, defining the second-difference operator
Dx2 that we can just plug into differential
equations containing the Laplacian, along with the first difference, to solve them
explicitly. Finite-difference computer programs just loop along time t, in
steps dt, to build up a computed wavefield (or temperature field) for the
time range we desire. As you saw in lab 2, if we make dt too large in an
attempt to save computation time, the faster velocities make the discretized
wavefield aliased in time, and the program blows up.
At right is a blue sinusoid of amplitude one and wavelength lambda equal to 4 in x. On the left it is sampled with the green line at dx=1 for 4 samples per cycle (or wavelength). On the right it is sampled with the red line at dx=2 for 2 samples per cycle. The coarser sampling and the red line are not aliased, since the Nyquist Criterion says 2 samples per cycle perfectly represent the wave.
The equations on the lower right of the figure here show that the second derivative of this sinusoid is a scaled sinusoid itself. At a value of x equal to three-quarters of the wavelength lambda, which is x=3, the analytic 2nd derivative has a value of 2.47.
Now, let's apply the formula for the 2nd-difference approximation at the same point, x=3. The estimate of the 2nd derivative withDx2 = [P(x+dx)-2P(x)+P(x-dx)]/dx2 is centered at a value of x, which will be 3.
On the left, the 3rd sample on the green line is at a minimum in the sine function. The 2nd derivative is supposed to be about 2.47, and Dx2 produces an estimate of 2.0. For geophysical work, this is sometimes accurate enough.
The second point on the red line, sampling the sine at 2 points per cycle, is at an identical minimum in the sine wave, and ought to still give a 2nd derivative close to 2.47. However, this coarser sampling does not produce nearly as accurate a Dx2 estimate, coming back with a value of 1. The 2nd difference estimate is highly inaccurate, despite the sampling being sufficient for the Nyquist Criterion.
Note that computing Dx2 at x=2, at the sine wave's zero crossing, will always yield the correct value, zero, for the 2nd derivative there. This is true for both the green and red samplings. In fact, Dx2 is near-perfect where the 2nd derivative is small, and underestimates the 2nd derivative more and more at it gets larger.
This factor of two error at the maxima will have drastic effects on the wave propagation
in E3D. Waves with high frequencies, having the shortest wavelengths, will seem to slow
down and even stop propagating once they get below six spatial samples per wavelength.
This slowing of high-frequency parts of the wavefield is reminiscent of the
dispersion (spread) of the velocities of surface waves according to their frequency.
The computational artifact is thus known as grid dispersion, since it
depends on the spatial sampling dh of the model grid.
For reliable computations, Shawn Larsen and other authors of finite-difference
wave modeling codes suggest that the shortest wavelengths in the calculation cover
at least 10 grid points-- in other words,
that 10dh < lambdamin .
At left is an example E3D computation visual (a map with m=0)
that has at least 10 spatial samples per
wavelength. At right is a visual from the same model, where the source was set to
pump a higher frequency into the grid that results in fewer than 6 points per cycle.
The grid dispersion generates false long durations of shaking. You can also notice
a ``moire'' pattern of waves that seem to be propagating parallel to the axes rather
than radiating circularly from the source. As the computation
continues, you will see waves propagating so slowly they appear to be stationary.
Run this model with the command:
/quake/rg/e3d/bin/e3d test3.in &Watch the visual to verify that no obvious grid dispersion occurs. As you did in Lab 2, rename the LSM and UNLV x-component output files to save them.