Numerical Modeling in Applied Geophysics by benbenzhou


									Numerical Modeling in Applied Geophysics

                Hansruedi Maurer
    Institute of Geophysics, ETH Zürich, Switzerland
 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:
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 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
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 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
 Root finding algorithms
 Numerical differentiation
 Numerical integration

                       IDEA League              10
Getting started: Number representation in a

 Numerical modelling is concerned with finite
 precision representations of numbers on a
 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
  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 for further information

                                              IDEA League                           13
Cancellation           From

   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)
                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 )
        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
                                                 f '( x )
                             a −1
     f '( x ) = ax                               f ( x)

 Relativ error of f(xe) is
                                    f ( x ) − f ( xe ) x a − xe
                                           f ( x)          xa

                     Numerically moderately stable!

                                            IDEA League           18
 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:


                          IDEA League           23
 Spline interpolation

Expand linear interpolation to

(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

                             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
                                  Construct secant through points 3
                                  and 4, etc.
                                  Always retains the most recent prior

                    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

                        Approximate f(x) by f’(x) at initial
                        guess (x1)

                        Find zero-crossing of resulting
                        tangent (x2)

                        Compute f’(x2), find zero-crossing,

                        Method may fail in the presence of
                        several local minima

                        Comparable with least-squares

                    IDEA League                                 33
 Numerical differentiation

 Taylor series approximation of f(x)
            f ( x + h) = f ( x) +        h + h.o.t. ≈ f ( x ) + f '( x )h

Rearranging terms yields
                         f ( x + h) − f ( x)
            f '( x ) =
                                                         One-sided finite difference
Alternatively one may write                              approximations of f’(x)
                         f ( x) − f ( x − h)
           f '( x ) =

                                               IDEA League                             34
Numerical differentiation

Central finite difference approximation is more accurate
                             f ( x + h) − f ( x − h)
                f '( x ) =

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

 Gaussian quadrature addresses this problem by
 incorporating uneven sampling and weigthing

                         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

 Actual integration includes evaluation of the sum on the
 right-hand side

                                 IDEA League                       41

To top