Transcript High Performance Computing 811

Computational Methods
in Physics
PHYS 3437
Dr Rob Thacker
Dept of Astronomy & Physics (MM-301C)
[email protected]
Project titles/descriptions

Please email me your project title if you have not
done so already

I would also like a one paragraph description of
what you intend to do, to be emailed to me next
week (Tuesday 3rd)
Why?
This will help you to plan things out now!


Today’s Lecture

Interpolation & Approximation II

More on cubic splines
Clamped versus natural
 Matrix techniques
 Tridiagonal solvers

Approximating derivatives
 Richardson Extrapolation

Recap: Origin of clamped versus
natural

Setting p´j,j+1(xj)=p´j-1,j(xj) in the last lecture we found
that
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 j 1 
 h j
(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”
Case 1: Clamped Spline

Specify first derivatives at end points only


p´1 & p´n are known
Go back to equation for p´ (8) at x1 to get
p '1, 2 ( x1 ) 
p ' '2  p ' '1
p  p1 h1 p ' '2 h1 p ' '1
( x1  x1 ) 2  p' '1 ( x1  x1 )  2


2h j
h1
6
3
p2  p1 h1 p' '2 h1 p ' '1



h1
6
3
 2h1 p ' '1  h1 p ' '2  6
p2  p1
 6 p'1
h1
(11)
Similarly , for j  n : hn 1 p ' 'n 1 2hn 1 p' 'n  6
pn  pn 1
 6 p 'n
hn 1
(12)
Collect equations

We now have a complete set of n equations in
the n unknowns p´´j:
 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 j 1 
 h j
2h1 p' '1  h1 p' '2  6
p2  p1
 6 p'1
h1
(11)
pn  pn 1
hn 1 p' 'n 1 2hn 1 p' 'n  6 p'n 6
hn 1

We now write this in matrix form
(12)
(10)
Matrix form for clamped splines
This matrix is said to have tridiagonal form: the only
non zero elements are between the two dotted lines.
2h1

 h1
 0








h1
0
0
0
2(h1  h2 )
h2
0
h2
2(h2  h3 ) h3

0


  p ' '1 


  p ' '2 
  p' '3 


   
  








2(hn  2  hn 1 ) hn 1 p' 'n 1 


p
'
'
hn 1
2hn 1   n 





hn  2
0
Note typo in DeVries, equation 3.38
1


(
p

p
)

p
'
2
1
1


h1


1
1


( p3  p2 )  ( p2  p1 )


h2
h1


1
1
( p4  p3 )  ( p3  p2 )


h3
h2













 1

1
( pn  pn 1 ) 
( pn 1  pn  2 )

hn  2
 hn 1

1


p 'n 
( pn  pn 1 )


hn 1


6
Eqn. (13)
Case 2: Natural Spline

p´´1 & p´´n are both set to zero


Provides two additional equations to the set (10)
Matrix form is then
0
0
0
1

h2
0
0 2(h1  h2 )
0
h2
2(h2  h3 ) h3











0
  p' '1 
  p' ' 
 2 
  p' '3 


   
  



  
2(hn  2  hn 1 ) 0  p' 'n 1 


p
'
'
0
1  n 
0





hn  2
0
0


1
1


(
p

p
)

(
p

p
)
3
2
2
1


h2
h1


1
1


( p4  p3 )  ( p3  p2 )


h3
h2











 1

1

( pn  pn 1 ) 
( pn 1  pn  2 )
hn  2
 hn 1



0
6
Eqn. (14)
How to calculate p(x)


Given either (13) or (14) we have all the pj necessary to
ensure continuous derivatives at the datums
Algorithm to find f(x) given x and set of pj:




Find j such that xj<x<xj+1
Warn user if x<x1 orx>xn (cannot extrapolate outside region)
Solve tridiagonal system to find all p´´j from the pj values
Once all p´´j are found then construct p(x) using equation (7)
and calculate value of p(x) (=f(x))
Tridiagonal Linear Systems

Consider the following matrix equation
 b1

 a2
0








c1
0
0
b2
c2
0
a3
b3

c3


0
0






an 1 bn 1
0
an
  x1   r1 

 r 
x
2

  2 
  x3   r3 

  


  
     

  


   
cn 1   xn 1  rn 1 

  
bn   xn   rn 
bj j=1,…,n occupy the Main Diagonal
aj j=2,…,n occupy the Sub-diagonal
cj j=1,…,n-1 occupy the Super-diagonal
Consider Gaussian Elimination

Let’s perform Gaussian elimination on the second row:
 b1

 a2
0








c1
0
0
b2
a3
c2
b3

0
c3


0
0






an 1 bn 1 cn 1
0
an
bn
r1 

r2 
r3 

 
 

 
rn 1 

rn 
Step 1: Need to subtract
a2/b1 of row 1 from row 2
to eliminate a2 from the
second row
“(2)=(2)- (a2/b1)×(1)”
After performing first elimination
c1
0
b1

a
 0 b2  2 c1 c2
b1

a3
b3
0









0
0
0
c3





0
 1

0
0








c1
0
 2 c2
0
c3


0
b3


an 1 bn 1 cn 1
0
an
bn
0
a3


0







r1 









rn 1 

rn 
r1
a
r2  2
b1
r3


an 1 bn 1 cn 1
0
an
bn
To make notation easier, let
1b1, 1=r1
2=b2-(a2/1)c1, 2=r2-(a2/1)1
1 

2 
r3 

 
 

 
rn 1 

rn 
Next perform elimination on
the third line
“(3)=(3)-(a3/2)×(2)”
Further substitution
 1

0
0








c1
0
0
2
c2
0
3
0
c3



0
0






an 1 bn 1 cn 1
0
an
bn
1 

2 
3 





rn 1 

rn 



Where we have substituted
3=b3-(a3/2)c2, 3=r3-(a3/2)2
It should be clear that we can do the same pattern of elimination
and then relabelling throughout the matrix.
After completing elimination steps
 1

0
0








c1
0
0
2
c2
0
3
0
c3



0
1
2
3
0






0  n 1 cn 1
0 0
n





 
 

 
 n 1 

 n 
We have reduced the tridiagonal
system to a very simple form.
The very last row can immediately
be solved.
Writing out the last row explicitly:
 n xn   n  xn  n /  n
We can now use back substitution to
solve for all the remaining values.
Applying back substitution

Writing out the penultimate line:
 n 1 xn 1  cn 1 xn   n 1
 xn 1 

1
 n 1
 n1  cn1 xn 
It should be clear that in general
xn j 
1
 n j

n j
 cn j xn j 1 ,
j  1,..., n  1
Summary of substitution process


Once we have the form of the initial triadiagonal
matrix:
Substitute 1=b1 1=r1 then
 j  bj 

aj
 j 1
c j 1 and
 j  rj 
aj
 j 1
 j 1 j  2,..., n
Once this is done, solve xn=n/n and then use back
substitution to get remaining xj
Basic algorithm for Cubic Spline
interpolation


(1) Use data (xj,f(xj)) to evaluate p´´(xj) where
j=1,…,n depending on whether user wants
clamped or natural spline
Calculate the interpolation polynomial around
the two points that bracket the desired x value
using the polynomial form (7)
Approximating derivatives

If we take the algebraic definition of a derivative
df ( x)
f ( x  h)  f ( x )
 lim
h 0
dx
h

We would naturally generalize this to numerical approximation
using
f ( x  h)  f ( x )
f 'h ( x ) 
(A)
h

Where h would be a separation between points
Note this is calculating a series of derivatives at the points f(x)
Just how good is this estimate?


Estimating the accuracy

If we Taylor expand about the point x, then
h2
h3
f ( x  h)  f ( x)  hf ' ( x) 
f ' ' ( x) 
f ' ' ' ( x)  ...
2!
3!
Rearrange to define
f ( x  h)  f ( x ) h
h2
f f ' ( x) 
 f ' ' ( x) 
f ' ' ' ( x)  ... (B)
h
2
6
Take this to be the error term which we write
as E(h). In this case the lowest power that
appears is h, so E(h)~h, or in “Big-O” notation
E(h) is O(h).
Remember h is small, so that h2<h, hence O(h2) is more accurate than O(h)
Forwards vs Backwards difference

The difference defined in (A)
is called a forward difference

The derivative is clearly
defined in the direction of
increasing x


We can equivalently define a
backwards difference

x-h
x
x+h
Using x & x+h
The derivative is defined in the
direction of decreasing x

Using x & x-h
h2
h3
f ( x  h)  f ( x)  hf ' ( x) 
f ' ' ( x) 
f ' ' ' ( x)  ...
2!
3!
f ( x )  f ( x  h) h
h2
 f b ' ( x) 
 f ' ' ( x) 
f ' ' ' ( x)  ... (C)
h
2
6
Improving the error


It should be clear that both the forward & backward
differences are O(h)
However, by averaging them we can eliminate the term
that is O(h) there by defining the centered difference:
1
 f ' f ( x)  f 'b ( x)
2
f ( x  h)  f ( x )  f ( x )  f ( x  h) h
h

 f ' ' ( x)  f ' ' ( x)
2h
2
2
h2
Note that all odd ordered

f ' ' ' ( x)  ....
6
terms cancel too. Further the
central difference is symmetric
f ( x  h)  f ( x  h)
f 'c ( x ) 
 O(h 2 )
about x, forward & backward
2h
f 'c ( x ) 
differences are clearly not.
Higher order differences

We can also derive approximations to the second
derivative by subtracting the backward difference from
the forward (i.e. setting f´f(x)-f´b(x)=0):
f ( x  h)  f ( x )  f ( x )  f ( x  h) h
h
 f ' ' ( x)  f ' ' ( x)  O(h 3 )
h
2
2
f ( x  h)  2 f ( x )  f ( x  h)
2
 f ' ' ( x) 

O
(
h
)
2
h
0  f ' f ( x)  f 'b ( x) 

This formula is also accurate to O(h2)
Richardson Extrapolation

This is a very widely applicable idea due to L.F.
Richardson (1881-1953)



We will see the idea again in the section on numerical
integration of functions
The algorithm may seem somewhat unusual at first, but it is
deceptively simple and elegant
The idea is that if we know the order of the underlying
error we can systematically improve the accuracy of a
numerical estimate

Makes sequences converge faster too
Richardson made seminal contributions to weather forecasting and turbulence
Outline

Let’s consider the “three point”* central difference
formula
f ( x  h)  f ( x  h) h 2
f 'h ( x ) 

f ' ' ' ( x)  O(h 4 ) (1)
2h
6
The same formula applies if we use a step size 2h :
f ( x  2h)  f ( x  2h) 4h 2
f '2 h ( x ) 

f ' ' ' ( x)  O(h 4 ) (2)
4h
6

Notice that we can subtract ¼ of (2) from (1) to
eliminate the leading term of the error
*Three point since it requires and produces information at three points
Eliminating the leading order error

We can now define a new estimate of the derivative using these
two formulae:
3 ( 4)
1
f ' ( x)  f 'h ( x)  f '2 h ( x) (Prefactor of 3/4 since f 'h ( x), f '2 h ( x) measure the same thing )
4
4
f ( x  2h)  8 f ( x  h)  8 f ( x  h)  f ( x  2h)
 f ( 4) ' ( x) 
 O(h 4 ) (3)
12h
Similarly for the second derivative
 f ( x  2h)  16 f ( x  h)  30 f ( x)  16 f ( x  h)  f ( x  2h)
4
f ( 4)''(x) 

O
(
h
) (4)
2
12h


These are “five point” central difference formulae, with the
furthest points receiving the least weight
I have added a (4) superscript to denote order h4 accurate
There is no reason to stop!

Given our new formula for f´(x) we could apply exactly the same
approach again


Since the factor of 2 in the change of step size will contribute 24
we need to subtract 1/24 of f´2h(x) from f´h(x)


The leading error term is now O(h4)
Notice in the first step we only had to subtract 1/22=1/4 because the
error term was O(h2)
Define
15 ( 6)
1 ( 4)
f ' ( x )  f ( 4 ) 'h ( x ) 
f '2 h ( x )
16
16

 f ( 6 ) ' ( x) 
16  ( 4)
1 ( 4)

f
'
(
x
)

f
'
(
x
)
h
2
h

15 
16

And so on… but we can’t continue indefinitely


Expressions rapidly get very complex
Better to work numerically…
A note on step size

The key point is that from one step to the next we let h double
in size (we could allow h to halve in size as well!)


This is independent of the number of points in the difference formula, or
the positions at which data are taken from
The key point is the relative change in the value of h, and not the value of
h itself


You can even think of resetting the value of h each time you want to
extrapolate to the next highest level of approximation
e.g. in the definition of f(6)´(x) we would use
f ( x  2h)  8 f ( x  h)  8 f ( x  h)  f ( x  2h)
 O(h 4 ) (3)
12h
f ( x  4h)  8 f ( x  2h)  8 f ( x  2h)  f ( x  4h)
( 4)
f 2 h ' ( x) 
 O(( 2h) 4 ) (3)
24h
f h ' ( x) 
( 4)

However the error in f2h(4)´(x) is (2h)4 so 24 times larger
Generalized notation


Let D1(h) be the three point central difference operator
we considered earlier (eqn (1))
Since E(D1(h))~h2

=> D1(2h) has 4 times the error of D1(h)
D1 (2h)  D1 (h) D1 (2h)  D1 (h)
 error 

3
22  1

Thus we can subtract this calculation of the error from
D1(h) to create a more accurate difference operator
(D2(h))
 D (2h)  D1 (h)  4 D1 (h)  D1 (2h)
D2 (h)  D1 (h)   1 2

(5)

2 1
3


i=2 step


Consider now the error in D2(h)
Since E(D2(h))~h4

=> D2(2h) has 16 times the error of D2(h)
D2 (2h)  D2 (h) D2 (2h)  D2 (h)
 error 

15
24  1

Thus we can subtract this calculation of the error from
D2(h) to create a more accurate difference operator
(D3(h))
 D2 (2h)  D2 (h)  16 D2 (h)  D2 (2h)
D3 (h)  D2 (h)  

(6)
4

2 1
15


ith step

Clearly this process generalizes
E ( Di (h)) ~ h 2i  Di (2h) contains 2 2i times error of Di(h)
Di (2h)  Di (h)
 error 
2 2i  1
2i
 Di (2h)  Di (h)  2 Di (h)  Di (2h)
 Di 1 (h)  Di (h)  

(7)
2i
2i

2 1
2 1


Numerical implementation

In practice, we implement the scheme by
actually halving the step size rather than doubling
it


That allows us to use old values of Di
For example, consider f(x)=xex, and we want to
find f’(x) at x=2
By hand: f’(x)=ex+xex=ex(1+x)
 Thus f’(2)=3e2=22.16716830…
 This is solely used as a check of our convergence

Answer=22.16716830
Start with h=0.4


Step 1: D1=23.16346429, used f(1.6) and f(2.4)
Step 2: Set h=0.2



Step 3: Set h=0.1




D1=22.41416066, used f(1.8) and f(2.2)
D2=22.16439278, where we have used D1 from step 1 in eqn (5)
D1=22.22878688, used f(1.9) and f(2.1)
D2=22.16699562, where we have used D1 from step 2 in eqn (5)
D3=22.16716914, where we have used D2 from step 2 in eqn (6)
Step 4: Set h=0.05




D1=22.18256486, used f(1.95) and f(2.05)
D2=22.16715752, where we have used D1 from step 3 in eqn (5)
D3=22.16716831, where we have used D2 from step 3 in eqn (6)
D4=22.16716830, where we have used D3 from step 3 in eqn (7)
Notice that D2 on step 2 (h=0.2) is a better estimate of the derivative than
D1 from step 3 (h=0.1). This shows the power of error removal process
Answer=22.16716830
Extrapolation table

These results can be neatly expressed as follows
Step
h
D1
D2
D3
1
0.4
23.16346429
2
0.2
22.41416066 22.16439278
3
0.1
22.22878688 22.16699562 22.16716914
4
0.05
D4
22.18256486 22.16715752 22.16716831 22.16716830
I have added arrows to indicate how values used in one calculation contribute to the next.
Deciding when to stop


When we use R.E. we are using it because we don’t
know the answer
We must therefore decide on the number of digits of
accuracy we require



Let nacc = number of digits of accuracy required
Need to consider both absolute and relative errors
On step i, define the difference between current value
and previous to be the absolute error


abserr=|Di-Di-1|
If (-log10(abserr).ge.real(nacc)) then done (A)
Also consider relative error

Define the relative/fractional error by
fracerr=|(Di-Di-1)/Di|
 If (-log10(fracerr).ge.real(nacc)) then done (B)


Need to consider both of these conditions (A &
B) at the same time
The relative error ensures convergence when Di is
large
 The absolute error is tripped when Di are small and
the relative error can potentially blow-up

Summary

Fitting cubic splines involves solving a
triadiagonal system


Tridiagonal solvers are efficient, and use back
substitution to calculate solutions
The accuracy of finite difference approximations
to derivatives can be greatly improved using
Richardson Extrapolation

R.E. relies upon cancelling the leading order error in
any expansion
Next Lecture


LU decompositions for matrix algebra
Least squares polynomial fitting