VIEWS: 9 PAGES: 41 POSTED ON: 11/16/2011 Public Domain
Numerical Modeling in Applied Geophysics Hansruedi Maurer Institute of Geophysics, ETH Zürich, Switzerland maurer@aug.ig.erdw.ethz.ch Your course instructor Name: Hansruedi Maurer Research Group: Applied and Environmental Geophysics Research interests: -Inversion of various geophysical data sets -Statistical experimental design -Application of geophysical tomography in multidisciplinary projects Contact information: Tel: 044 633 6838; Fax: 044 633 1065; email: maurer@aug.ig.erdw.ethz.ch Address: Institute of Geophysics, Schafmattstrasse 30, 8093 Zürich Office: HPP O7 IDEA League 2 Course schedule Date Content Module Lecturer Mon 25.2. Introduction, basics Mod M1 H. Maurer Tue 26.2. Potential field modeling Mod M2 H. Maurer Wed 27.2. Potential field modeling (exercise) Mod M3 H. Maurer Thur 28.2. Layered Earth modeling Mod M3 H. Maurer Fr 29.2. Layered Earth modeling (exercise) Mod M3 H. Maurer Mon 3.3. Explicit finite-difference schemes Mod M4 H. Maurer Tue 4.3. Implicit finite differences Mod M5 H. Maurer Wed 5.3. Implicit finite differences (exercises) Mod M5 H. Maurer Thur 6.3. Finite elements Mod M6 H. Maurer Fri 7.3. Finite elements Mod M6 H. Maurer Mon 10.3. Sensitivities Mod M7 H. Maurer Tue 11. – Student projects Mod M8 H. Maurer Wed 13.3. Thu 14.3. Presentation student projects Mod M8 H. Maurer Doctoral students Fri 14.3. – Finalizing exercises, exam preparation H. Maurer Wed 19. 3. Doctoral students Thu 20.3. Exam IDEA League 3 Lectures Lectures start at 08.15 in HCI E2 Typically, we will teach for about 90 minutes followed by a 30 minute break. Afterwards, there will be some more teaching or exercises. Modules finish either at 9.45 or 11.45. IDEA League 4 Exercises Most modules include exercises. Most (all) of them require MATLAB. 75% of all the exercises must be solved. You can team up for the exercises, but maximum group size is two! Exercises should be handed in no later than April 1. Completed exercises should be sent by email to the lecturer responsible. IDEA League 5 Student projects Student projects are either programming exercises or literature studies. They will be conducted by teams of 2 students. Results of the projects will be presented to the class in form of ~15 min talks followed by discussions. Grading of the student project will be based primarily on the presentation. Themes of student projects will be presented later in the course. IDEA League 6 Exam Exam will last from 8.30 to 11.30. Questions are set up such that they can be answered within 2 hours. No extra supporting material is allowed. IDEA League 7 Earning the credits To receive the 3 credit points you are requested to complete 75% of the exercises satisfactorily, make an acceptable presentation of your project, achieve a sufficient mark (> 4.0) in the exam. Overall grading of the course: 25% exercises 25% student project 50% final exam IDEA League 8 Course notes You can download the course notes from the EVA platform. Access information will be provided by email. It is highly recommended to print the handouts of the slides (3 per page) prior to the lectures (PDF files with handouts are available on EVA). You are responsible for printing your course notes. On EVA you will also find a literature list including books and key papers on which this course is based. IDEA League 9 Module 1: Some basic issues of computational geophysics Number representation in a computer Cancellation and error propagation Interpolation Root finding algorithms Numerical differentiation Numerical integration IDEA League 10 Getting started: Number representation in a computer Numerical modelling is concerned with finite precision representations of numbers on a computer. Inappropriate choice of a data representation may result in Inaccurate results Excessive memory requirements Suboptimal computation speed IDEA League 11 Integer numbers Data type Size Range Char 1 byte (8 bit) -127 to 127 -27-1 to 27-1 Short integer 2 bytes (16 bit) -32‘768 to 32‘768 -215-1 to 215-1 Long integer 4 bytes (32 bits) -2‘147‘483‘648 to 2‘147‘483‘648 -231-1 to 231-1 Operations within the valid range of a data type will be perfectly accurate Integer operations are generally much faster than floating point operations Range of individual data types can be extended (shifted), when positivity can be ensured Some programming languages also offer 8 byte integers Range of long integers limits addressable memory of 32 bit operating systems to about 2 Gigabytes. IDEA League 12 Floating point numbers Data type Size Range Float (Real4) 4 bytes (32 bit) ~ +/- 10-38 to 1038 Double float (Real8) 8 bytes (64 bits) ~ +/- 10-308 to 10308 Operations involving floating point numbers are not exact 4 byte numbers have 6 significant digits and 8 byte floating point numbers have approx. 15 significant digits (Def. MachEps) Numerical modeling in Geophysics considers quantities that extend over several magnitudes 6 significant digits are often not enough! There exist various representations of floating point numbers See http://en.wikipedia.org/wiki/Floating_point for further information IDEA League 13 Cancellation From http://www.cs.mcgill.ca/~chang/teaching/cs350/content.php Relative error “worst case scenario“ IDEA League 14 Error propagation Given a function f(x). How does an error in the argument x (xe instead of x) influences the function f? Taylor approximation of f(xe) arround f(x) f ( xe ) = f ( x ) + f '( x )( x − xe ) + f ''( x )( x − xe ) 2 + … f ( xe ) ≈ f ( x ) + f '( x )( x − xe ) f ( x ) − f ( xe ) f ' ( x ) x − xe x − xe Relativ error of f(xe) is ≈ x =K f ( x) f ( x) x x f ( x ) − f ( xe ) K is the so called condition number f ( x) f '( x ) K= = x x − xe f ( x) x If K is large, the problem is said to be ill-posed! IDEA League 15 Error propagation: subtraction f ( x) = a − x K= f '( x ) x= −x a−x f '( x ) = −1 f ( x) f ( x ) − f ( xe ) xe − x Relativ error of f(xe) is = f ( x) a−x If a is close to x numerical cancellation may occur! IDEA League 16 Error propagation: multiplication f ( x ) = ax K= f '( x ) x =1 f '( x ) = a f ( x) Relativ error of f(xe) is f ( x ) − f ( xe ) x − xe = f ( x) x Numerically stable! IDEA League 17 Error propagation: exponentiation f ( x) = x a K= f '( x ) x=a a −1 f '( x ) = ax f ( x) Relativ error of f(xe) is f ( x ) − f ( xe ) x a − xe a = f ( x) xa Numerically moderately stable! IDEA League 18 Interpolation Many options are available in literature Most popular variants include Polynomial interpolation and Spline interpolation Here we consider only 1D examples. Extensions to higher dimensions are generally straightforward IDEA League 19 Polynomial interpolation (Vandermonde approach) IDEA League 20 Polynomial interpolation (Lagrange approach) IDEA League 21 Polynomial interpolation Use Vandermonde approach, when many interpolations for a given set of (x,y) pairs need to be performed Use Lagrange approach (resp. more clever implementations of it), when only a few interpolations need to be performed IDEA League 22 Spline interpolation Find interpolation scheme that is smooth in the first derivative f’(x) and continuous in the second derivative f’’(x). Linear interpolation: with IDEA League 23 Spline interpolation Expand linear interpolation to with (it can be shown that y’’ = Ay’’j + B y’’j+1) Problem: y’’(x) is not known Solution: By requesting that y’(x) is continuous everywhere, a system of equations can be set up that yields y’’(x). IDEA League 24 Root finding algorithms Find x such that f(x) = 0 Many applications in numerical modeling Related to optimization problems (f’(x) = 0 minimum of f(x)) Algorithms should be efficient and/or robust We consider Bracketing and bisection Secant and regula falsi method Brent’s method Newton method Script follows “Numercial Recipes” from Press et al. (see also http://www.nrbook.com/a/) IDEA League 25 Root finding algorithms - Bracketing Find an interval [a,b], in which the function f(x) has at least one root In some cases such an interval may not exist! IDEA League 26 Root finding algorithms - Bisection method 1. Find interval [a,b] in which f(x) contains a root (f(a)*f(b) < 0) 2. Determine midpoint of [a,b] and check, in which part of the interval the sign change occurs 3. Find midpoint of chosen part of the interval 4. Repeat steps 2 and 3 until sub-intervals become sufficiently small IDEA League 27 Root finding algorithms - Bisection method Advantage: ALWAYS finds a root (if it exists) Disadvantage: Convergence may be slow IDEA League 28 Root finding algorithms - Secant method Construct secant through initial points 1 and 2 Evaluate f(x) at zero crossing of secant point 3 Construct secant through points 2 and 3 Evaluate f(x) at zero crossing point 4 Construct secant through points 3 and 4, etc. Always retains the most recent prior estimate IDEA League 29 Root finding algorithms - Method of false position Similar to secant method, but retains previous point with opposite sign IDEA League 30 Root finding algorithms - Pathological case Here, secant and regula falsi will provide a poor convergence! IDEA League 31 Root finding algorithms - Brent’s method Combines advantages of bisection and regula falsi/secant method Utilizes quadratic interpolation Replaces midpoint of bisection method by more clever estimate Special care is required, when Q is close to zero (From “Numerical Recipes”) IDEA League 32 Root finding algorithms - Newton method Requires first derivative f’(x) to be known. Approximate f(x) by f’(x) at initial guess (x1) Find zero-crossing of resulting tangent (x2) Compute f’(x2), find zero-crossing, etc. Method may fail in the presence of several local minima Comparable with least-squares inversions IDEA League 33 Numerical differentiation Taylor series approximation of f(x) ∂f f ( x + h) = f ( x) + h + h.o.t. ≈ f ( x ) + f '( x )h ∂x Rearranging terms yields f ( x + h) − f ( x) f '( x ) = h One-sided finite difference Alternatively one may write approximations of f’(x) f ( x) − f ( x − h) f '( x ) = h IDEA League 34 Numerical differentiation Central finite difference approximation is more accurate f ( x + h) − f ( x − h) f '( x ) = 2h but computational effort is usually higher, since f(x) is mostly required as well! IDEA League 35 Numerical integration Numerical solution of More complicated than numerical differentiation Many variants exists Here we discuss Rectangle, trapzoidal and Simpson rule Gaussian quadrature and modifications thereof IDEA League 36 Numerical integration – Rectangle rule IDEA League 37 Numerical integration – Trapezoidal rule IDEA League 38 Numerical integration – Simpson’s rule IDEA League 39 Numerical integration – Gaussian quadrature Rectangle, trapezoidal and Simpon’s rule consider equally spaced sampling points. This may lead to inefficiencies or numerical inaccuracies when the functions are not well behaved (e.g. the occurrence of integrable singularities) Gaussian quadrature addresses this problem by incorporating uneven sampling and weigthing schemes. IDEA League 40 Numerical integration – Gaussian quadrature Rewrite integrand as a product with a weighting function W(x) Weighting coefficients need to be determined at suitable sampling points xj Sampling points xj are determined via the roots of a polynomial function Actual integration includes evaluation of the sum on the right-hand side IDEA League 41