I  Lagrange Polynomial Interpolation

Recall, given p + 1 distinct points {i }i(p)=0 , the Lagrange interpolating polynomials

are:

p

pp (x) =       f (i )Li (x).                                         (1)

i=0

0 Introductory question, not for credit

•  Let f (x) = cos(πx) + x and let 0  = - , 1  = 0 and 2  = .

•  Find the (quadratic) Lagrange polynomials L0 (x), L1 (x) and L2 (x) based on the points 0 , 1  and 2 .  (You can do this either by hand, or by using Python). Use Python to plot these 3 Lagrange polynomials on one figure and make sure they look correct.

•  Use (1) to construct the polynomial interpolant p2 (x) of f (x). (Again, you can do this either by hand, or using Python).  Once you have this interpolant, use Python to plot it and compare against f (x) over the interval [-3, 3].

•      Go   to   the   Geogebra   demonstration https://www.geogebra.org/m/ bwmpekfaand compare your answers with the one obtained there when a = - and b = .  Investigate using higher polynomials degrees and try different func-

tions f (x) over different intervals.

1 (This is a warm up exercise. No marks will be awarded for it.  Its solution will be made available on the 7th February, however, you are thoroughly encouraged to attempt this question in advance of the questions for credit being released. The solution to this problem will be required in questions that follow.)

The file lagrange_polynomials .py contains an unfinished function with the following definition:

def   lagrange _ poly (p , xhat ,n ,x , tol )

This function returns:

–  lagrange_matrix - a numpy .ndarray of shape (p+1,n)

–  error_flag - an integer.

•  Complete the function so that given a set of p + 1 distinct nodal points {i }i(p)=0 (stored in xhat) and a set of n evaluation points {xj } (stored in x), it will return a (p + 1) x n matrix stored in lagrange_matrix, where the ijth entry of the matrix is Li (xj ).

•   In addition, your function should also perform a check to make sure that the nodal points {i }i(p)=0  are distinct.  If the points are distinct, then the output vari- able error_flag should be set to 0, otherwise error_flag should be set to 1. Floating point numbers x and y should be considered equal if lx - yl < tol.

•  No other checks on the inputs are required.

•        Once    you    have    written    lagrange_poly,     test    it    by    running lagrange_polynomials.py. You should obtain:

lagrange_matrix  =

[[    6 .5625     2 .1875     0 .3125    -0 .0625     0 .0625    -0 .3125    -2 .1875]  [-11 .8125    -2 .1875     0 .9375     0 .5625    -0 .3125      1 .3125     8 .4375]  [    8 .4375      1 .3125    -0 .3125     0 .5625     0 .9375    -2 .1875  -11 .8125]  [  -2 .1875    -0 .3125     0 .0625    -0 .0625     0 .3125     2 .1875     6 .5625]]

error_flag  =  0

2 In approximations .py you will find the following function definition:                          [6 / 40]

def   p o l y _ i n t e r p o l a t i o n (a ,b ,p ,n ,x ,f , produce _ fig )

The function should return as output:

–  interpolant - a numpy .ndarray of shape (n,);

–  fig - a matplotlib .figure .Figure when the Boolean input produce_fig is True and None if produce_fig is False.

•  Complete the function so that it evaluates the pth order polynomial interpolant pp (x) of a function f at a set of points {xj } The nodal interpolation points should be uniformly spaced over the interval [a, b], including the endpoints.

•  The function must call lagrange_poly from Q1 with tol  =  1 .0e-10.

•   If produce_fig is True then the function should also plot (on the same set of axes) the function f evaluated at the points {xj } and the interpolant pp (x)

evaluated at the same points.  (Note, when testing this, it is assumed that the points in x form an increasing sequence.)

•       Once   you   have   written   poly_interpolation   test   it   by   running

approximations.py. Your output should be:

interpolant  =

[2 .64872127  2 .79438712  2 .83153269  2 .80097589  2 .74353466  2 .7000269 2 .71127054  2 .81808351  3 .06128371  3 .48168907]

def   i n t e r p o l a t i o n _ e r r o r s (a ,b ,n ,P , f )

The function should return as output:

–  error_matrix - a numpy .ndarray of shape (n, );

–  fig - a matplotlib .figure .Figure.

•       Complete   the   function   to   compute   (approximations   to)   the   error maxxe[a,b] lppj (x) - f (x)l for  a  range  of  polynomial degrees P =  {pj } The  numpy .ndarray  error_matrix  should  contain  the  errors  when  a  uni- form  set  of  interpolating  nodes  is  used.    Hence,  your  function  should  call poly_interpolation.

Hint: In order to find an approximation to the error, evaluate the error lpp (x) - f (x)l for 2000 equally spaced points over [a, b] and take the maximum of those.

•  The function should also plot the errors against {pj }, using plt .semilogy and return this figure in fig.

•   Add a description at the top of your interpolation_errors function as a “docstring” that when help(interpolation_errors) is run will comment and explain the results when P = {1, 2, 3, . . . , 10} and

(a)  f (x) = e2x , [a, b] = [0, 1];

(b)  f (x) = , [a, b] = [-5, 5];

•  There is no need to include the usual function description in docstring.

•      Once   you   have   written   interpolation_errors  test   it   by   running compute_errors.py. You should obtain:

error_matrix  =

[0 .17863278  0 .04634435  0 .01368204  0 .00431683  0 .0014182  ]

I  Interpolation of Multivariate functions

Lagrange interpolation can also be used for interpolation of multivariate functions.

Suppose we have a function f : R2  → R that we would like to approximate over the rectangular region [a, b] x[c, d]. To interpolate with a polynomial of degree p, we can define a set of p + 1 distinct nodal points in both the x and y directions: {i }i(p)=0 and {yˆi }i(p)=0 , respectively. f can then be interpolated by the function

p      p

pp (x, y) =             f (i , yˆj )L北,i (x)Ly,j (y),

i=0  j=0

and {Ly,j (y)}j(p)=0 are the Lagrange polynomials associated with the y nodal points.

4 In approximations .py you will find a function with the following definition:              [8 / 40]

def   p o l y _ i n t e r p o l a t i o n _ 2 d (p ,a ,b ,c ,d ,X ,Y ,n ,m ,f , produce _ fig )

The function should return as output:

–  interpolant - a numpy .ndarray of shape (m,n);

–  fig - a matplotlib .figure .Figure when the Boolean input produce_fig is True and None if produce_fig is False.

•  Complete the function so that it evaluates the pth order polynomial pp (x, y) of a function at a set of grid points stored in the m x n arrays X and Y. Here, X and Y are assumed to have been created using the numpy .meshgrid function. Hence, the ijth component of interpolant should contain the value of pp (Xij , Yij ). The nodal interpolating points in x and y should be uniformly spaced over the intervals [a, b] and [c, d] respectively, including the endpoints.

•  The function must call lagrange_poly from Q1 with tol  =  1 .0e-10.

•  If produce_fig is True then the function should produce a contour plot of the interpolant pp (x), which is returned in fig.

•     Once  you   have  written  poly_interpolation_2d  test   it   by   running approximations.py. You should obtain:

interpolant  =

[[2 .71828183  3 .03773178  4 .23949621  7 .3890561  ]

[0 .99999993  1 .11751899  1 .55962339  2 .71828164]

[2 .71828183  3 .03773178  4 .23949621  7 .3890561  ]]

I  Numerical Differentiation

Recall that, given a point x and a set of distinct points {j }i(p)=0 , such that i  = x for some 0 < i < p we can find an approximation to f\ (x) as follows:

p

f\ (x) ≈ p(x) =       f (xi )L(x),

i=0

where pp (x) is the interpolant and {Li }i(p)=0  are the usual Lagrange polynomials as- sociated with the points {i }i(p)=0 .

tion:                                                                                                                                    [8 / 40]

def   d e r i v _ l a g r a n g e _ p o l y (p , xhat ,n ,x , tol )

The function returns:

–  deriv lagrange matrix - a numpy.ndarray of shape (p+1,n)

–  error flag - an integer

•  Complete the function so that given a set of p + 1 distinct nodal points {i }i(p)=0 (stored in xhat) and a set of n evaluation points {xj } (stored in x), it will return a (p + 1) x n matrix stored in deriv_lagrange_matrix, where the ijth entry of the matrix is L(xj ).

• Hint: Derive a formula for L(x) by hand using the product rule then implement this. Do not use SymPy.

•   In addition, your function should also perform a check to make sure that the nodal points {i }i(p)=0  are distinct.  If the points are distinct, then the output vari- able error_flag should be set to 0, otherwise error_flag should be set to 1. Floating point numbers x and y should be considered equal if lx - yl < tol.

•  No other checks on the inputs are required.

•      Once   you   have   written   deriv_lagrange_poly   test   it   by   running lagrange_polynomials.py. You should obtain:

deriv_lagrange_matrix  =

[[-17 .875    -7 .435    -1 .315     0 .485    -2 .035    -8 .875]

[  41 .625    13 .905    -0 .855    -2 .655     8 .505    32 .625]

[-32 .625    -8 .505     2 .655     0 .855  -13 .905  -41 .625]

[    8 .875     2 .035    -0 .485      1 .315     7 .435    17 .875]]

error_flag  =  0

6 In approximations .py you will find a function with the following definition:               [6 / 40]

def   a p p r o x im a t e _ d e r i v a t i v e (x ,p ,h ,k , f )

The function returns:

–  deriv_approx - a floating point number of type numpy .float64

•   Complete the function so that it will evaluate the derivative of the pth order polynomial interpolant of a function f at a point x. This value is an approximation to f\ (x).  The interpolating points should be equally separated so that {j }j(p)=0 are such that

j  = 0 + jh,    j = 0, . . . , p,

for h > 0 and should be positioned so that k  = x, where 0 < k < p. That is, x coincides with one of the nodal points.

•    The  function  should  make  use  of  deriv_lagrange_poly  from Q5 with tol  =  1 .0e-10.

• Once you have written deriv_approx, test it by running approximations .py. The outputs should be

k  =  0,  deriv_approx  =  -2 .150454361444488

k  =  1,  deriv_approx  =  -2 .1405845045118443

k  =  2,  deriv_approx  =  -2 .1405845045118417

k  =  3,  deriv_approx  =  -2 .1504543614444867

7 In compute_errors .py you will find a function with the following definition:               [6 / 40]

def   d e r i v a t i v e _ e r r o r s (x ,P ,m ,H ,n ,f , fdiff )

The function returns:

– E - a numpy .ndarray of shape (m,n)

–  fig - a matplotlib .figure .Figure

•  Complete the function so that for a set of m even polynomial degrees P = {pi }1  and a set of n widths H = {hj }, it will return an m x n matrix E , where

Eij  = lf\ (x) - pi,hj (x)l,    0 < i < m - 1, 0 < j < n - 1

and ppi,hj (x) is the polynomial interpolant of order pi with interval width hj  such that x = pi/2 . i.e. pi,hj (x) is the centred pi + 1 point approximation to f\ (x).

•  Input fdiff is a function that computes the exact derivative f\ (x).

•  The function must call derivative_approximation from Q6.