week 2

1
week 2

• last week introduction
• Numerical errors from representation doubles
• machine dependant
• Machine accuracy, about 17 digits for double, 7
  for float
• Round-off errors
• Truncation errors
• Minimize with optimal routines
• Stability exercise 1 (due next week)
• Binnen van Bouwe, Rik
• Test for stability when applicable
• Routines www.nikhef.nl/henkjan/CODE

2
chapter 3 Interpolation

• Interpolation of functions
• fixed number of points known
• from measurements
• from calculations. e.g continuous wave Faddeev
  equation, takes too much time to calculate per
  variable value
• calculate points in between
• error?
• interpolation in between
• extrapolation -gt outside range dangerous
• However Pade formalism - analyticity
• interpolation
• polynomials
• trigonometric -gt Fourier analysis (later in
  course)
• rational functions (pade)
• other types of polynomials (e.g. Legendre.....)

3
interpolation

• very weak pole
• mocks all interpolation schemes close to pole

4
interpolation

• very weak pole
• mocks all interpolation schemes close to pole
• function approximation
• approximate a function by an easier calculable
  function
• calculated points of your own choosing!
• e.g. Chebyshev polynomials, gaussian quadrature
• fitting functional form is known.
• not part of the course
• We may add it later

5
interpolation

• interpolation
• use tabulated points around x-value of interest
• interpolated function shifts at tabulated point!
• discontinuous derivatives
• spline- smooth derivatives, stiffer
• cubic spline often used (spline)
• higher-order better?

6
Interpolation
7
interpolation

• higher order interpolation fine for smooth
  functions
• worse for rapidly-changing derivatives
  (especially polynomial)
• most simple case polynomial interpolation
• Lagrange
• goes by construction through n points
• solved with Nevilles algorithm
• gives error estimate

8
Nevilles algorithm

• zeroth order
• P1 y1
• P2 y2
• ... Pn yn
• higher orders P12, P123, P1234 .....
• higher orders recursively obtained
• e.g. P12 (x-x2)y1 (x1-x)y2/(x1-x2)
• keep track of difference between Pi..im and the
  lower order. This gives error estimate.

9
Nevilles algorithm

• Di,m1 (xim1-x)(Ci1,m - Di,m)/(xi- xim1)
• Ci,m1 (xi -x)(Ci1,m-Di,m)/(xi - xim1)
• error estimate last difference added
• differences with previous order remembered
  recursive formula.
• higher-order interpolation can be easily
  recursively updated
• Routine polint
• input arrays of x and y values
• give it appropriate offset e.g. x_15ltxltx_16
  x14, n4

10
Rational interpolation

• -poles
• (complex plane ruins convergence)
• rational functions will stay good, as long as
  there are enough powers in denominator
• Pade approximation, chapter 5.12
• R(x) P(x)/Q(x)
• by definition Q0 1
• n-th order goes through n1 points munu1n
• Bulirsch and Stoer Neville-type algorithm
• diagonal mu (1) nu, as many powers in P(x) as
  in Q(x)
• also error estimate by comparison with previous
  order

11
Pade extrapolation

• A Pade approximant is the rational function that
  has a power series expansion that agrees to a
  given power series to the highest known order
  (MN or M1N)
• very useful if you can calculate a function and
  several orders of the derivative.

12
Pade Approximation, example

• Deep-inelastic scattering structure functions
  are determined and interpreted in terms of
  quark/gluon contents in the nucleon. (Spin and
  momentum of the nucleon) hermes.eps
• more sea-quarks resolved at higher momentum
  transfer.
• structure functions can be evolved from one value
  of x,Q2 to another Q2 via DGLAP, third order
  QCD calculations. Pade extrapolation is used.
• DGLAP evolutions, AltarelliParisi

13
Pade approximation, ch. 5.12
14
Spline

• linear interpolation
• zero second derivative, discontinous first on
  grid points
• Cubic spline continuous on second derivative,
  smooth on first
• How to proceed add 3rd-order polynomial, such
  that y linearly varies from yi to yi 1 .
• Polynomial constructed to be zero at grid points

15
spline

• This yields for the 2nd derivative
• second derivative linearly changes between xj and
  xj1
• However Y not known.
• Require that first derivative y is continuous
  across boundary
• use this to get y
• equate yj from j-1,j and j,j1
• this yields for j1..N-2

16
Spline

• N-2 equations, N unknowns
• set y0 and yN-1 zero
• this is called natural cubic spline
• or set either y(0) and y(N-1) to values to
  give a specified value of y(0) and/or y(N-1)
• cubic spline tri-diagonal. Each y is coupled
  to its two neighbours.
• O(N) calculations to solve the N equations.
• for general schemes one would expect O(N3)
  calculations
• spline call only once. Store parameters. Splint
  routine gives interpolated values.
• spline typically available in operating system
  (e.g. linux)

17
Exercise 2008
18
Exercise 2008
19
Interpolation in multiple dimensions

• much tougher to obtain good interpolation
• linear interpolation in two directions often
  done.
• bicubic splines is an alternative. Smooth
  derivatives
• take out factors that you know! E.g. if a
  function f(a,b) has a factor exp(a), divide it
  out.
• If derivatives can be calculated, higher accuracy
  may be obtained (also for bicubic spline).

20
interpolation in 2 dimensions

• bilinear interpolation

21
interpolation in 2 dimensions
22
Quadratic roots

• We have discussed
• Analytically, you can get the roots in 2 ways
• Right procedure