Part I: how to become an expert at solving differential equations in five minutes…
What is Finite Difference?
In computational mathematics, finite-difference (FD) methods are numerical methods for solving differential equations by approximating them with difference equations, in which finite differences approximate the derivatives. Many outstanding texts have stimulated the development of the calculus of finite differences. One of the first, presented in a form suited to the needs of students and teachers, is that by the English mathematician George Boole (1815–1864), who published the Treatise on Differential Equations in 1859, followed in 1860 by the Treatise on the Calculus of Finite Differences.
The history of calculus does not begin with the findings of the English mathematician and physicist Sir Isaac Newton (1642–1726) and the German philosopher and mathematician Gottfried Wilhelm Leibniz (1646–1716). Their calculus was the culmination of centuries of work by many mathematicians. Newton and Leibniz argued over who had first invented calculus: Newton made his discoveries in 1664–1666, but his findings were not published until 1693. Leibniz worked between 1672 and 1676, but published the results in 1684–1686, before Newton. In 1711, the controversy was taken to court. The Royal Society appointed a commission to investigate the charges, which found Leibniz guilty of plagiarism; perhaps not that surprising since Newton was the president of the Royal Society. Years after Leibniz’s death, the mathematical community came to realize that Newton and Leibniz had made their discoveries independently.
History of Solving Finite Difference
In 1928, the German-born American mathematician Richard Courant (1888–1972) published the theoretical fundament for the solution of problems of mathematical physics by means of finite differences. Among other things, Courant and coworkers defined the FD approximation for the wave equation, and the famous Courant–Friedrichs–Lewy (CFL) condition – a necessary condition for convergence while solving partial differential equations numerically by the method of finite differences.
Considerable progress in FD methods was made during and after WWII, when practical applications became possible with the use of computers. John von Neumann (1903–1957), the Hungarian-American mathematician, physicist, and computer scientist, developed the von Neumann stability analysis, known as Fourier stability analysis, used to check the stability of finite difference schemes as applied to linear partial differential equations. In 1949 he worked with meteorologists Jule Charney (1917–1981) and Ragnar Fjørtoft (1913–1998) on numerical weather prediction, and von Neumann’s Fourier method was given a rigorous treatment in their joint publication in the periodical Tellus (Charney et al., 1950).
For hyperbolic equations the FD method has played and is still playing a dominant role, starting with the work of, e.g., Friedrichs, Lax, and Wendroff in the 1950s and 1960s (Lax, 1954; Lax, 1961; Lax and Wendroff, 1964). Standard references on FD methods are the books of Forsythe and Wasow (1960) and Richtmyer and Morton (1967).
Finite Difference Methods in Seismology
Today, FD methods are the dominant approach to numerical solutions of partial differential equations (Grossmann et al., 2007). For wave equations solved in seismics, the FD method is accurate and robust (Fornberg, 1988; Etgen et al., 2009; Robertsson and Blanch, 2011; Ikelle and Amundsen, 2018). We note that the FD method was introduced in 1966 by Yee to discretize the differential form of Maxwell’s equations. Alterman and Karal (1968), Boore (1970, 1972), and Kelly et al. (1976) then applied the FD method applied to seismology and seismics. They used a displacement formulation with conventional grids, which yielded instability problems in models with high-velocity contrasts. Virieux (1984, 1986) introduced the stable staggered-grid velocity-stress FD schemes, following Madariaga (1976) who introduced the staggered-grid formulation for dynamic modeling of earthquake ruptures.
High School Calculus
Equation Sheet 2.
Calculus (from Latin calculus, literally ‘small pebble’; the small stone the ancient Romans used in counting and gambling) is the mathematical study of continuous change. One way to think of calculus is as the study of functions of time or of space. There are two different types of calculus:
Differential calculus
Differential calculus divides functions into different pieces and tells us how they change from one moment to the next.
Integral calculus
Integral calculus joins or integrates the small pieces together and tells us how much of something is made by a series of changes.
In high school, you learn that differential calculus is concerned with finding the instantaneous rate of change (i.e., derivative) of a function’s value, with respect to changes within the function’s arguments. The derivative of a function with respect to a variable is denoted, after Leibniz, by: (see Equation Sheet 1).
The Wave Equation
Equation Sheet 2.
The laws of physics are generally written down as differential equations, or relations involving rates at which things happen; that is, derivatives. A differential equation expresses a relationship between a function and its derivatives. If we know the function and its derivatives at a particular point or time, then this information, together with the differential equation, can be used to determine the function over its entire domain.
One example is the wave equation. It is a second-order linear partial differential equation which describes how waves, such as sound, light and water waves, travel. The French mathematician Jean le Rond d’Alembert (1717–1783) was the first to find and solve the 1D wave equation – now known as d’Alembert’s solution in his honor. Ten years later, Leonhard Euler (1707– 1783) solved the three-dimensional wave equation.
See Equation Sheet 2.
Solving the wave equation means finding p in terms of t and x. It turns out that solving the wave equation, as most other differential equations, can be quite hard, and there is no general method that solves every differential equation. We will generally focus on finding a numerical solution to the wave equation by using computers.
Solving the Wave Equation by Computer
Figure 1: To compute the wave field u at time-space location (t+Δ t , x) (black dot) we use the wavefield at two previous times t and t– Δ t at three space locations x– Δ x, x, x+Δ x (red and blue dots).
To solve the wave equation by computer, we discretize the time and space variables (see Figure 1), replace both second derivatives in the wave equation with finite differences, and obtain:
This solution is a time marching scheme, where the next value of the wavefield at the discrete time t+Δ t is computed from current values known at time t and the previous time t– Δ t . For every next time step, the source is activated with its current value.
The choice of time step Δ t is dependent on the grid size Δ x. It turns out that as we increase the time step while keeping the grid size fixed, the FD method eventually becomes unstable. Likewise, if we decrease the grid size while keeping the time step fixed, we run into stability problems. To know how to change the time step with changes in grid size in order to maintain stability, we revert to the famous CFL condition for stability, which for our wave equation has the form:
Cmax Δt / Δx < 1, where Cmax is the maximum velocity in the model.
There is little need to know much more at this point. You can now leap on the train to a job as a geophysicist in the oil industry with whistles and fanfare as you chug out of the station.
Solution for Constant Velocity
When the velocity is constant, and the source is a point monopole source, located say at x=0, the wave equation has the analytical solution of an outgoing wave from the source:
u (t,x) = S(|x| – ct),
where S is the time integral of the source signature s in the wave equation. The pressure wave is not altered as distance x increases. (In 3D, however, the pressure wave decreases in amplitude with distance.)
Let’s check if the FD solution to the wave equation reproduces the analytical solution. The source is located at x=0. It is given a time signature that is a time derivative of a Ricker signature, sometimes called the Mexican hat wavelet:
where fM is the peak frequency of the signal, and t0 is time where the zero-phase Mexican hat has its maximum amplitude. Parameters are c0=3,000 m/s, Δt=0.00058 s, Δx=2.5m, and fM=10 Hz. The output is displayed in Figure 2 as time traces at four different recording locations from the source. Indeed, the waveform is the Mexican hat wavelet travelling outwards from the source without change.
Figure 2: In a homogeneous 1D medium, the FD solution of the wave equation produces a wave that moves outwards from the source without change. The four Mexican hat signals are recorded at four different locations. This result agrees with d’Alembert’s solution, discovered in 1746.
Solution for Non-Constant Velocity
We let the velocity take the form of the long-wavelength sine-model
and run the FD scheme with velocities c0=2,000 m/s, Δ c=500 m/s and λ=250m. We set the pressure to zero at the end of the computational domain and display the pressure wavefield in Figure 3 as a 2D color plot over distance and time. Observe that the wave changes sign after hitting the pressure-free ends of the domain.
Figure 3: FD solution of the 1D wave equation in the long-wavelength sine model. The source is located at the center distance of the model (*), and the wavefield is recorded at all distances and times. At the endpoints of the model, the wavefield is fixed to zero, so that the wave is mirrored into the model domain with opposite sign in the amplitude (colors change from dominantly red to blue).
Coming soon in the next issue…
In Part 2 of this series, we will discuss the numerical errors that one can make in FD modeling and bring you the latest research on how to eliminate errors caused by the time-stepping of the wavefield.
