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 **D**_{x} 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
**D**_{x}^{2}, which is just applying
**D**_{x} 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
**D**_{x}^{2} 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
2^{nd} derivative has a value of 2.47.

Now, let's apply the formula for the 2^{nd}-difference approximation
at the same point, *x*=3. The estimate of the 2^{nd} derivative
with**D**_{x}^{2}
= [*P(x+dx)-2P(x)+P(x-dx)*]/*dx*^{2}
is **centered** at a value of *x*, which will be 3.

On the left, the 3^{rd} sample on the green line is at a minimum in the sine
function. The 2^{nd} derivative is supposed to be about 2.47, and
**D**_{x}^{2} 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 2^{nd} derivative
close to 2.47. However, this coarser sampling does not produce nearly as accurate
a **D**_{x}^{2} estimate, coming back with a
value of 1. The 2^{nd} difference estimate is highly inaccurate, despite
the sampling being sufficient for the Nyquist Criterion.

Note that computing **D**_{x}^{2} at *x*=2,
at the sine wave's zero crossing, will always yield the correct value, zero, for the
2^{nd} derivative there. This is true for both the green and red samplings.
In fact, D_{x}^{2}
is near-perfect where the 2^{nd} derivative is small, and underestimates
the 2^{nd} 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 10*dh* < *lambda*_{min} .
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.

**Assignment:**

- Given the criterion for avoiding grid dispersion,
10
*dh*<*lambda*_{min}, derive an expression that gives you the*dh*needed to meet this criterion for a highest desired frequency f_{max}. This expression will involve a velocity-- is this the minimum or the maximum velocity found in the model? - If one had a code that was computing only P-wave propagation in a crustal
model, with a
*dh*that avoided grid dispersion, and switched to computing only shear-wave propagation across the same realistic model, would*dh*have to increase or decrease? Could you guess by how much? - Briefly describe a few examples of lithospheric geologic features that affect the
Courant Condition on
*dt*, and then describe the types of features that affect the need to set*dh*small enough to avoid grid dispersion. - To compute a seismogram representing a given duration, the cost of computation
goes up as
*dt*gets smaller, since more time steps are needed. Each time step in general requires computation effort proportional to the number of grid nodes in the model, which is*nx*times*ny*times*nz*. Thus the computational effort required for a model of a certain size goes up as*dh*decreases, as well. Can you think of a set or sets of lithospheric features, which, if included in a model, cause the computational cost to skyrocket by demanding both small*dt*and small*dh*? - Download the test3.in E3D input file. Examine it
in a text editor. Note that we are not using ModelAssembler for this lab, so all
velocities are set in test3.in . Given the frequency specified on the
*source*line of the input file, and the other settings in the file, what is the minimum number of grid samples per wavelength for the computation set up here?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. - Edit the
*source*line of the test3.in file to increase the source frequency by a factor of five. How many points per wavelength will this computation have as a minimum? Run this model and observe the grid dispersion artifacts, which become more obvious as the computation proceeds. Again rename the LSM and UNLV output files to save them. - Plot the four seismogram files for the two runs on the LSM and UNLV stations
in SAC, label them
completely, including which are from which station and which run, and point out
some possible features of grid dispersion on the plot.
- The original source line in test3.in is not a efficient as it might be. How high can you raise the frequency and still avoid grid dispersion? Set this value and make another run. Print out a comparison of the seismograms from these two runs (both without grid dispersion).

Resources:

- E3D modeling code manual, defining inputs and outputs
- test3.in E3D input file, defining how E3D will run.