Computational Methods in Physics PHYS 3437

1. Computational Methods in Physics PHYS 3437 Dr Rob Thacker Dept of Astronomy & Physics (MM-301C) thacker@ap.smu.ca

2. Today’s Lecture • Interpolation & Approximation I • Overview of ideas • Interpolation vs fitting • Polynomial interpolation • Lagrange interpolation • Hermite interpolation • Cubic splines (introduction)

3. f(x) x4 x1 x3 x5 x2 X General idea of interpolation • Suppose we are given data at discrete points only, and we wish to estimate values between these known points

4. f(x) x4 x1 x3 x5 x2 Difference between interpolation & fitting • Interpolation passes through each datum • Fitting seeks to produce a curve that approximates the data

5. f(x) x4 x1 x3 x5 x2 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

6. 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) --- (1) • f(x3)=p(x)+(x3-x)p’(x) --- (2)

7. Eliminate p’(x) • After a short derivation (ex) we find that eqns (1) and (2) give • 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

8. 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 You should see a pattern emerging here…

9. General (Lagrange) Polynomial fit • If we have n points, a polynomial of degree n-1 can be constructed using the formula • This is called the interpolation polynomial in Lagrange form, but it is frequently referred to as the Lagrange polynomial

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

11. Problems with polynomial interpolation • As the number of points increases so does the degree of the polynomial • 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 • Cubic fits tend to be about optimal • 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

12. Example of “wiggles” Piecewise linear fit 6th order polynomial fit From http://dmpeli.mcmaster.ca/Matlab/Math1J03/LectureNotes/Lecture3_3.htm

13. 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 • 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

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

15. 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) x2 x1 x (x1<x<x2)

16. Interpolating a cubic to two points • Given the general form for a cubic So we have 4 equations with 4 unknowns to solve.

17. 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) • 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

18. Diagram of two Hermite interpolation polynomials Second polynomial fit p2,3(x) First polynomial fit p1,2(x) f(x) f(x3) & f’(x3) f(x1) & f’(x1) f(x2) & f’(x2) Derivatives match at x2 x1 x3 x2

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

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

21. 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 • If we again equate at xj and xj+1 but using 2nd deriv:

22. Next find the cj • We’ve just shown

23. Substitute back to get p(x) • We now substitute (2),(4),(5) and (6) back into our first equation for p(x) to get To get any further we now need some information about p’’j+1 and p’’j

24. 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) Derivatives match at x2 x1 x3 x2

25. 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:

26. 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: Where j=2,3,…,n-1

27. Finding the p´´j values • We now have a system of n-2 equations • Recall limits of j=2,3,…,n-1 • There are n unknown values though • 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”

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

29. Next Lecture • Matrix forms for spline fits • Clamped spline • Natural spline • Tridiagonal solvers