Basics of Finite Differences: Stability, Accuracy, Dispersion, Fourier (100 characters)

Simple FD Gerard T. Schuster

Basics • Accuracy • Stability • Dispersion • Fourier Finite Difference Outline

Basics of Finite Differences • The discretization of differential operators is not unique • 1st-order approximations are one-sided approximations, full lines in figure: • 2nd-order difference has higher accuracy but does not involve values at grid point xj dotted line in figure • Demands on numerical stability as well as numerical accuracy can not be satisfied in all cases

Explicit Methods • Advance numerical solution in time from the old time level t(n) to t(n+1) with spatial difference operatorH • Explicit methods: • Spatial derivatives are calculated using the values at the old time level U(n) • Extrapolating the solution into the future • Time step t is restricted, e.g. CFL-condition • Easy to implement, parallelization possible • Most numerical methods, operator splitting

Accuracy I • The accuracy of difference approximations can be tested by simple Taylor-series expansion • For the right-sided approximation we obtain • Adding the two expressions we see the more accurate central differences

Accuracy II • right-sided, 1st - order • left-sided, 1st - order • central difference, 2nd - order central difference on a non-equidistant grid

2nd Derivatives • A simple approximation of the 2nd-derivative on an equidistant grid can be obtained by summing up the left- an right-sided finite difference expressions, subtract 2fj • Higher order approximations involve more grid points, e.g. 4-order expression • Higher order: Be careful with the corresponding boundary conditions

Accuracy of 2nd Derivatives • Simplest approximation of the 2nd-derivative on an equidistant grid • Again, simple test function: f = C eikx and same procedure as before • See again 2x0, when x0 • Approximation is 2nd-order, waves are dispersive with ~k2 • Such difference scheme for diffusion problems, Laplace and Poisson equation

Accuracy of waves • Compare finite difference expression with differential operator, e.g. central difference and 1st-derivative • Simple test function: f = C eikx • Only small values of kx are well represented by sin(kx) • k has to be small or the wavelength 2/k has to large • The smaller k the less diffusion, transport of waves is dispersive • Errors are always seen on the smallest scales, zigzag-pattern

Stability • The stability of any numerical method is measured by the amplification of numerical errors in time • nis the error at time step tn and g is the so-called amplification factor (or amplification matrix for systems of equations): • Numerical stability requires that the amplification of errors remains limited, i.e. • Only for simple, linear discretizations allow a rigorous analysis of the amplification matrix, non-linear systems yield to non-linear numerical schemes

Calculation of Stability • Taking a 1-dimensional interval X = [X1,X2], dividing it into (J-1)-subintervals where a function f(x) is defined on this interval X at distinct grid points xj with 1j J and f(xj)=fj • Fourier decomposition of f(x) on X, assume an equidistant grid with x • Fourier modes are orthogonal, therefore single modes can be analyzed to calculate stability of the numericalscheme

Example: Von Neumann Calculation of Instability ∂^2P/∂t^2=-c∂(∂P/∂t)/∂x = c^2∂^2P/∂x^2 (2) ∂P/∂t=-c∂P/∂x (1) t+1 t t t ∂P/∂t = (P - P )/dt ;∂P/∂x=(P - P )/(2dx) (3) j j j+1 j-1 Substitute (3) into (1) we get t+1 t t t P = P -cdt/(2dx)(P - P ) (4) j j j+1 j-1 ikjdx -iwt Substitute P = e e into (3) we get ik(j+1)dx ik(j-1)dx -iwt ikjdx -iwt ikjdx -iw(t+dt) e e = e -cdt/(2dx)(e - e )e (5) ikjdx – iwt Divide by e we get – iwdt g= e = 1 - icdt/(dx)sin{kdx} ; |g|=sqrt{1+(cdt/dxsin{kdx})^2} Amplification factor because Pjn=gnPj0 Stability demands |g|<1, which is impossible. Above is unstable because errors grow with time

Example: Von Neumann Calculation of Stability ∂^2P/∂t^2=-c∂(∂P/∂t)/∂x = c^2∂^2P/∂x^2 (2) ∂P/∂t=-c∂P/∂x (1) t+1 t t t ∂P/∂t = (P - P )/dt ;∂P/∂x=(P - P )/(2dx) (3) j j j+1 j-1 Substitute (3) into (1) we get t t t P = 0.5(P + P ) j j +1 j-1 t+1 t t t P = P -cdt/(2dx)(P - P ) ; (4) j j j+1 j-1 ikjdx -iwt Substitute P = e e into (4) we get eikjdx-iw(t+dt)= {cos(kdx) - icdt/(dx)sin[kdx] }eikjdx-iwt; or dividing by eikjdx-iwt Amplification factor because Pjn=gnPj0 – iwdt g= e = cos(kdx) - icdt/(dx)sin[kdx] ; |g|=sqrt{cos(kdx)^2 +(cdt/dxsin{kdx})^2} Stability demands |g|<1, which it is if cdt/dx<1. Above is conditionally stable because errors damps with time if CFL holds (a necessary but not suficient condition for stability). As along as accuracy of FD does not exceed highest order derivative, CFL for symmetric stencils is both sufficient and necessary.

Example: Von Neumann Calculation of Stability t+1 t t t P = P -cdt/(2dx)(P - P ); (4) j j j+1 j-1 t+1 t t t t t t P = P -cdt/(2dx)(P - P ) + 0.5dx^2(P -2P + P )/dx^2 t t t j j j+1 j-1 j+1 j j-1 P = 0.5(P + P ) j j +1 j-1 This is a Laplacian, and emulates a diffusion term that damps errors; this is why Lax worked. Less accuracy but better stability

Problems: Dispersion Dispersion: The numerical approximation has artificial dispersion, in other words, the wave speed becomes frequency dependent. You have to find a frequency bandwidth where this effect is small. The solution is to use a sufficient number of grid points per wavelength. True velocity

Snapshot Example

Seismogram Dispersion

Basics • Accuracy • Stability • Dispersion • Pseudo-Spectral or Fourier Finite-Difference Modeling Outline

What is a Pseudo-spectral Method? Spectral solutions to time-dependent PDEs are formulated in the frequency-wavenumber domain and solutions are obtained in terms of spectra (e.g. seismograms). This technique is particularly interesting for geometries where partial solutions in the -k domain can be obtained analytically (e.g. for layered models). In the pseudo-spectral approach - in a finite-difference like manner - the PDEs are solved pointwise in physical space (x-t). However, the space derivatives are calculated using orthogonal functions (e.g. Fourier Integrals, Chebyshev polynomials). They are either evaluated using matrix-matrix multiplications or the fast Fourier transform (FFT). Numerical Methods in Geophysics The Fourier Method

MATLAB Code FD Approx. to 1-D Acoustic Wave Equation a = (c*dt/dx)^2; dp2 fort=1:nt f = force(:,:,t); dfuture = dpres - 2* dpast + a*del(dpres); dfuture = dfuture + f; dpast=dpresent; dpresent=dfuture; Dp2=fft(p);Dp2=k^2*Dp2;dp2=ifft(Dp2); end

Fourier Derivatives .. let us recall the definition of the derivative using Fourier integrals ... ... we could either ... 1) perform this calculation in the space domain by convolution 2) actually transform the function f(x) in the k-domain and back Numerical Methods in Geophysics The Fourier Method

The Fast Fourier Transform ... the latter approach became interesting with the introduction of the Fast Fourier Transform (FFT). What’s so fast about it ? The FFT originates from a paper by Cooley and Tukey (1965, Math. Comp. vol 19 297-301) which revolutionised all fields where Fourier transforms where essential to progress. The discrete Fourier Transform can be written as Numerical Methods in Geophysics The Fourier Method

The Fast Fourier Transform ... this can be written as matrix-vector products ... for example the inverse transform yields ... .. where ... Numerical Methods in Geophysics The Fourier Method

The Fast Fourier Transform ... the FAST bit is recognising that the full matrix - vector multiplication can be written as a few sparse matrix - vector multiplications (for details see for example Bracewell, the Fourier Transform and its applications, MacGraw-Hill) with the effect that: Number of multiplications full matrix FFT N2 2Nlog2N this has enormous implications for large scale problems. Note: the factorisation becomes particularly simple and effective when N is a highly composite number (power of 2). Numerical Methods in Geophysics The Fourier Method

The Fast Fourier Transform Number of multiplications Problem full matrix FFT Ratio full/FFT 1D (nx=512) 2.6x105 9.2x103 28.4 1D (nx=2096) 94.98 1D (nx=8384) 312.6 .. the right column can be regarded as the speedup of an algorithm when the FFT is used instead of the full system. Numerical Methods in Geophysics The Fourier Method

Acoustic Wave Equation - Fourier Method let us take the acoustic wave equation with variable density the left hand side will be expressed with our standard centered finite-difference approach ... leading to the extrapolation scheme ... Numerical Methods in Geophysics The Fourier Method

Acoustic Wave Equation - Fourier Method where the space derivatives will be calculated using the Fourier Method. The highlighted term will be calculated as follows: multiply by 1/ ... then extrapolate ... Numerical Methods in Geophysics The Fourier Method

Acoustic Wave Equation - 3D .. where the following algorithm applies to each space dimension ... Numerical Methods in Geophysics The Fourier Method

Fourier Method - Dispersion and Stability ... with the usual Ansatz we obtain ... leading to Numerical Methods in Geophysics The Fourier Method

Fourier Method - Dispersion and Stability What are the consequences? a) when dt << 1, sin-1 (kcdt/2) kcdt/2 and w/k=c -> practically no dispersion b) the argument of asin must be smaller than one. Numerical Methods in Geophysics The Fourier Method

Fourier Method - Comparison with FD - 10Hz Example of acoustic 1D wave simulation. FD 3 -point operator red-analytic; blue-numerical; green-difference Numerical Methods in Geophysics The Fourier Method

Fourier Method - Comparison with FD - 10Hz Example of acoustic 1D wave simulation. Fourier operator red-analytic; blue-numerical; green-difference Numerical Methods in Geophysics The Fourier Method

Fourier Method - Comparison with FD - 20Hz Example of acoustic 1D wave simulation. Fourier operator red-analytic; blue-numerical; green-difference Numerical Methods in Geophysics The Fourier Method

Fourier Method - Comparison with FD - Table Difference (%) between numerical and analytical solution as a function of propagating frequency Simulation time 5.4s 7.8s 33.0s Numerical Methods in Geophysics The Fourier Method

Basics of Finite Differences: Stability, Accuracy, Dispersion, Fourier (100 characters)