1 / 40

Simple FD

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:

sguzman
Download Presentation

Simple FD

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Simple FD Gerard T. Schuster

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

  3. 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

  4. 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

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

  6. 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

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

  8. 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

  9. 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

  10. 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

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

  12. 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

  13. 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

  14. 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

  15. 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.

  16. 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

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

  18. 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

  19. Snapshot Example

  20. Seismogram Dispersion

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

  22. 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

  23. 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

  24. 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

  25. 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

  26. 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

  27. 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

  28. 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

  29. 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

  30. 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

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

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

  33. 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

  34. 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

  35. 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

  36. 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

  37. 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

  38. 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

  39. 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

  40. 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

More Related