Transcript High Performance Computing 811

Computational Methods
in Physics
PHYS 3437
Dr Rob Thacker
Dept of Astronomy & Physics (MM-301C)
[email protected]
Today’s Lecture

Interpolation & Approximation I

Overview of ideas


Interpolation vs fitting
Polynomial interpolation
Lagrange interpolation
 Hermite interpolation
 Cubic splines (introduction)

General idea of interpolation

Suppose we are given data at discrete points
only, and we wish to estimate values between
these known points
f(x)
x1
x2
x3 X
x4
x5
Difference between interpolation &
fitting


Interpolation passes through each datum
Fitting seeks to produce a curve that approximates the
data
f(x)
x1
x2
x3
x4
x5
Global versus piecewise interpolation


We can fit a single polynomial of degree n to n+1 data points
Alternatively we can use a series of polynomials to link points together –
frequently called splines

Splines will usually be cubic polynomials – I’ve used linear here to emphasize the piecewise
nature
f(x)
x1
x2
x3
x4
x5
Lagrange interpolation



To fit a polynomial of degree n we need n+1 points
Consider first interpolating between two points,
f(x2),f(x3) using a function F(x)
The Taylor expansion about a point x gives



f(x2)=F(x)+(x2-x)F’(x)+…
f(x3)=F(x)+(x3-x)F’(x)+… (I’ve highlighted what we know)
If we use a polynomial of degree 1, p(x), then the
Taylor series truncates after p’(x)


f(x2)=p(x)+(x2-x)p’(x)
f(x3)=p(x)+(x3-x)p’(x)
--- (1)
--- (2)
Eliminate p’(x)

After a short derivation (ex) we find that eqns (1) and (2) give
( x  x3 )
( x  x2 )
p( x) 
f ( x2 ) 
f ( x3 )
( x2  x3 )
( x3  x2 )

On varying the expansion point x, but keeping fixed end points
x2,x3 this corresponds to a straightline – as expected (check the
following for yourself):



If x=x3 => p(x)=f(x3)
If x=x2 => p(x)=f(x2)
Linear fits aren’t that helpful though, and bringing in more
points allows us to up the order of the interpolation
Quadratic fit using three points

If we again Taylor expand put assume the function is a
polymial, p(x) of degree 2, then




f(x1)=p(x)+(x1-x)p’(x)+(x1-x)2p’’(x)/2 --- (3)
f(x2)=p(x)+(x2-x)p’(x)+(x2-x)2p’’(x)/2 --- (4)
f(x3)=p(x)+(x3-x)p’(x)+(x3-x)2p’’(x)/2 --- (5)
With a bit more algebra (fairly lengthy) we can eliminate
both p’(x) and p’’(x) using (3),(4),(5) to get
p( x) 
( x  x2 )( x  x3 )
( x  x1 )( x  x3 )
( x  x1 )( x  x2 )
f ( x1 ) 
f ( x2 ) 
f ( x3 )
( x1  x2 )( x1  x3 )
( x2  x1 )( x2  x3 )
( x3  x1 )( x3  x2 )
You should see a pattern emerging here…
General (Lagrange) Polynomial fit

If we have n points, a polynomial of degree n-1 can be
constructed using the formula
n
p( x)   l j ,n ( x) f ( x j )
j 1
where
l j ,n ( x) 

( x  x1 )( x  x2 )...( x  x j 1 )( x  x j 1 )...( x  xn )
( x j  x1 )( x j  x2 )...( x j  x j 1 )( x j  x j 1 )...( x j  xn )
This is called the interpolation polynomial in Lagrange
form, but it is frequently referred to as the Lagrange
polynomial
Properties of the lj,n(x)


The lj,n(x) are clearly polynomials, and a called basis polynomials
since they sum to find the interpolating polynomial
If we let x=xi in the formula (so that we are picking one of the
interpolation points) then
l j ,n ( xi ) 
( xi  x1 )( xi  x1 )...( xi  x j 1 )( xi  x j 1 )...( xi  xn )
( x j  x1 )( x j  x1 )...( x j  x j 1 )( x j  x j 1 )...( x j  xn )
1
0
if i  j (since the numerator & denominato r are equal)
if i  j (since one of numerator terms must be zero)
Thus
l j ,n ( xi )   ij where  ij is the Kronecker delta, and hence
n
p ( xi )    ij f ( x j )  f ( xi ) which is what we would expect.
j 1
Problems with polynomial
interpolation

As the number of points increases so does the degree
of the polynomial



Cubic fits tend to be about optimal


This implies many turns and a lot of “structure” between
interpolation points in high order polynomials
High order polynomials will have a lot of “wiggles” and if
points do not change smoothly they can “overshoot” in
between datums
4 data points can surround x symmetrically (2 on each side)
Increasing the number of points adds more data from
further away from a given interpolation point

This doesn’t necessarily improve the accuracy of the fit
locally
Example of “wiggles”
Piecewise linear fit
6th order polynomial fit
From http://dmpeli.mcmaster.ca/Matlab/Math1J03/LectureNotes/Lecture3_3.htm
Segue - numerically implementing
the general formula


Implementing the general formula is not exactly trivial
It helps to notice that the lj,n(x) can be written as a
product
( x  x1 )( x  x2 )...( x  x j 1 )( x  x j 1 )...( x  xn )
x  xi
l j ,n ( x)  

( x j  x1 )( x j  x2 )...( x j  x j 1 )( x j  x j 1 )...( x j  xn )
i 1,i  j x j  xi
n


One can then think about using a loop to do the
multiplication of the terms in the product
However, the best way to calculate large interpolation
polynomials is to use Neville’s Algorithm
Segue - Neville’s Algorithm


We won’t go into the details, however you can take a
look (if you are interested) at Numerical Recipes section 3.1
Neville’s algorithm is better because it builds up the
interpolation polynomial in a recursive fashion


It uses a similar table idea to building up coefficients in a
binomial expansion
The employed recursion cuts down the total number of
operations required and helps reduce concerns about
rounding error
Hermite Interpolation

If we have both the function value, and the
derivative we can fully constrain a higher order
polynomial through two points
f(x)
f(x2) & f ’(x2)
f(x1) & f ’(x1)
x1
x
x2
(x1<x<x2)
Interpolating a cubic to two points

Given the general form for a cubic
p( x)  ax 3  bx 2  cx  d
the derivative is
p' ( x)  3ax 2  2bx  c
fitting the two points gives
f ( x1 )  ax13  bx12  cx1  d
f ( x2 )  ax23  bx22  cx2  d
fitting the derivative s gives
f ' ( x1 )  3ax  2bx1  c
2
1
f ' ( x2 )  3ax22  2bx2  c
So we have 4 equations
with 4 unknowns to solve.
Solve for p(x)

While again somewhat lengthy, we can solve for p(x) in terms of
the known values
(subscript denotes for x1,x2)
(3x1  2 x  x2 )( x  x2 ) 2
(3x2  2 x  x1 )( x  x1 )2
p1, 2 ( x) 
f
(
x
)

f ( x2 )
1
( x1  x2 )3
( x2  x1 )3
( x  x1 )( x  x2 ) 2
( x  x2 )( x  x1 ) 2

f ' ( x1 ) 
f ' ( x2 )
2
2
( x1  x2 )
( x2  x1 )

One important property of Hermite interpolation is that if we fit
a different polynomial p2,3(x) between x2 & x3 then it will have
the same derivative at x2 as p1,2(x2) i.e. p’1,2(x2)=p’2,3(x2)



Interpolation between successive pairs of points is continuous
The derivatives are also continuous
Hermite interpolation is thus considered to be smooth
Diagram of two Hermite
interpolation polynomials
Second polynomial fit
p2,3(x)
First polynomial fit
p1,2(x)
f(x)
f(x1) & f ’(x1)
f(x3) & f ’(x3)
f(x2) & f ’(x2)
Derivatives match at
x2
x1
x2
x3
Advantages of Hermite interpolation
over Lagrange interpolation




A series of H.I. polynomials will be continuous
at the datums, as will the derivatives (i.e. smooth)
While a series of L.I. polynomials is continuous
at the datums the derivatives will not be
continuous
H.I. can fit a cubic to `local’ data (i.e. x1<x<x2)
L.I. requires four points to fit a cubic
(x1<x2<x<x3<x4)
However, even if we don’t have derivatives it would be nice to fit a series of
polynomials with continuous slopes at each datum
Cubic Splines



Let’s examine how we can create a series of cubic
polynomials with equal derivatives at the datums
Choose x, such that xj<x<xj+1 where j=1,…,n
Around x, we may write p(x) as
p( x)  a j ( x  x j )3  b j ( x  x j ) 2  c j ( x  x j )  d j
(1)
as in the discussion on Lagrange interpolat ion set x  x j then
p(x j )  p j  d j  f(x j )
(2)
If we set x  x j 1 then
p ( x j 1 )  p j 1  a j ( x j 1  x j ) 3  b j ( x j 1  x j ) 2  c j ( x j 1  x j )  p j
where we have substitute d for d j  p j using (2)
If we let x j 1  x j  h j then
p ( x j 1 )  p j 1  a j h 3j  b j h 2j  c j h j  p j
(3)
Next consider derivatives


We just derived two equations (3) & (2) for the value of
the cubic fit at xj and xj+1
Taking the first & second derivative of p(x) we find
p' ( x)  3a j ( x  x j ) 2  2b j ( x  x j )  c j

and
p' ' ( x)  6a j ( x  x j )  2b j
If we again equate at xj and xj+1 but using 2nd deriv:
p' ' ( x j )  p' ' j  2b j so b j  p'' j / 2 (4) and p' ' ( x j 1 )  p' ' j 1  6a j h j  2b j
 p' ' j 1  6a j h j  p' ' j
1 p' ' j 1  p' ' j
a j 
6
hj
solve for a j to get
(5)
Next find the cj

We’ve just shown
1 p' ' j 1  p' ' j
aj 
6
hj
(5), b j  p'' j / 2
(4)
we substitute into
p( x j 1 )  p j 1  a j h 3j  b j h 2j  c j h j  p j
(3)
to get
1
p j 1  ( p' ' j 1  p' ' j )h 2j  p' ' j h 2j / 2  c j h j  p j
6
p j 1  p j 1
 cj 
 h j p' ' j 1 2h j p' ' j 
(6)
hj
6
Substitute back to get p(x)

We now substitute (2),(4),(5) and (6) back into
our first equation for p(x) to get
p( x) 
p' ' j 1  p' ' j
6h j
(x  x j ) 
3
p' ' j
2
(x  x j )
2
 p j 1  p j h j p' ' j 1 h j p' ' j 



( x  x j )  p j
6
3 
 h j
(7)
To get any further we now need some information about p’’j+1 and p’’j
Setting p´´j+1 and p´´j

We have no information about the derivative of f(x) though



We just have the series of f(x) values
What requirements can we impose?
Recall we want the derivatives of neighbouring polynomials to match at the
datums (xj)
f(x)
f(x3)
f(x1)
f(x2)
x1
x2
Derivatives match at
x2
x3
Apply continuity in the derivative


We now have to consider adjacent polynomials
Label:



pj,j+1 to be the fit between xj and xj+1
pj-1,j to be the fit between xj-1 and xj
Continuity of the derivative implies that p´j,j+1 and p´j-1,j
should be equal at xj, the general p’ forms are:
p' j , j 1 ( x) 
p' j 1, j ( x) 
p' ' j 1  p' ' j
2h j
p' ' j  p' ' j 1
2h j 1
( x  x j ) 2  p' ' j ( x  x j ) 
p j 1  p j
hj
( x  x j 1 ) 2  p' ' j 1 ( x  x j 1 ) 

h j p' ' j 1
p j  p j 1
h j 1
6


h j p' ' j
h j 1 p' ' j
6
3

(8)
h j 1 p' ' j 1
3
(9)
Equate at xj

Setting p´j,j+1(xj)=p´j-1,j(xj) the first few terms of (8)
vanish, while terms in (x-xj-1) in (9) cancel with hj-1:
p ' j , j 1 ( x j ) 
p ' j 1, j ( x j ) 
p j 1  p j
hj

p ' ' j  p ' ' j 1
2
h j p ' ' j 1
6

h j p' ' j
3
h j 1  p ' ' j 1 h j 1 
p j  p j 1
h j 1

h j 1 p ' ' j
6

h j 1 p ' ' j 1
3
Equate and rearrange to get
 p j 1  p j p j  p j 1 
h j 1 p ' ' j 1 2(h j  h j 1 ) p ' ' j  h j p ' ' j 1  6 


h
h


j
j 1
(10)
Where j=2,3,…,n-1
Finding the p´´j values

We now have a system of n-2 equations


There are n unknown values though


Recall limits of j=2,3,…,n-1
The p’’j values range through j=1,…,n
To have enough information to solve for all the values
we have to provide information about the end points –
two cases

Specify derivatives at end points


“Clamped spline”
Set p’’1=0 and p’’n=0

“Natural spline”
Summary




Interpolation shouldn’t be confused with fitting,
although the terms are frequently used interchangeably
Lagrange interpolation will fit an order n polynomial to
n+1 data points
Hermite interpolation uses derivatives to solve for a
cubic fit with only 2 datums
Cubic spline interpolation gives a smooth piecewise
interpolation scheme but requires conditions to be set
for the end points
Next Lecture

Matrix forms for spline fits
Clamped spline
 Natural spline


Tridiagonal solvers