400 likes | 408 Views
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 1 st -order approximations are one-sided approximations, full lines in figure:
E N D
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
Basics • Accuracy • Stability • Dispersion • Fourier Finite Difference Outline
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 2x0, when x0 • 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 kx are well represented by sin(kx) • 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
Basics • Accuracy • Stability • Dispersion • Fourier Finite Difference Outline
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 1j 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
Basics • Accuracy • Stability • Dispersion • Fourier Finite Difference Outline
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
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. FD 5 -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. FD 3 -point 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. FD 5 -point 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