ECE 2240: Numerical Methods for Engineers Winter 2017 Lab 5: Numerical Differentiation and Integration Apr 3-7, 2017 The purpose of this lab is to: 1) Implement numerical methods for approximating derivatives. 2) Implement numerical integration techniques. Important: 10% of your lab mark will be based on evaluation during your specified lab period. Please be sure to attend your lab section. Report: Hand-in a report of your work that clearly indicates your name and student number as well as the name of your partner if applicable. What should be included is clearly laid out at the end of each question. Numerical Differentiation and Integration Objectives/Introduction During this lab you are to implement and evaluate numerical differentiation and integration methods. Question 1 You are to construct a discrete approximation to the function f (x) = sin(x) over the interval x [0, 2] at n equally spaced points. The easiest way to do this is to use the Matlab function linspace. Type 'help linspace' for details. Given this discrete approximation to f (x) = sin(x) you now seek to compute the derivative of f (x) at each of the n sample points. To accomplish this you will write one function and one main program: 1. A function numericalDerivative: function [df ] = numericalDerivative(x,f,type) The inputs are the x locations in x, the function values f corresponding to x and a character array that denotes the type of derivative you seek. Supported types should 5-1 Lab 5: Numerical Differentiation and Integration Apr 3-7, 2017 be 'F1', 'B1', 'C1' for the 2-point forward, backward and centered difference respectively as well as 'F2', 'B2' and 'C2' for the three point forward, backward and centered differences. You can check for string equality easily. For example strcmp(type,'F1'). The output is the derivative df of f (x) with respect to x evaluated at the x locations. Remark: You will need to handle the end-points carefully: for example at the first x-location you will not be able to apply a backward difference even if 'B1' is requested. For these cases use a different rule (just at those points where the requested rule fails) that maintains the accuracy of the approximation that has been requested. 2. A main program Lab5Q1: This main program should be used to tabulate the error in the approximations for each of the 6 derivatives for the following cases n = 10, n = 100, n = 1000 and n = 10000. We will compute the error as: error = ||x(df f 0 )||2 where f 0 is an array of the same size as df that contains the analytic value of the derivative at the sample points. Note that we are using the two-norm for vectors in this case. We include the step-size x to account for the fact that as n changes the number of entries used to compute the norm will change (this makes it analogous to a pulse integration). The main program should also produce the following plots: on one figure plot the derivative function C1 versus the analytic derivative for n = 10 and 1000; on a separate figure plot the derivative function F 2 versus the analytic derivative for n = 10 and 100. Hand-In: The code for the the function and the main program, a table of the requested simulation results and the requested plots. Make sure that the code is well documented and formatted. Question 2 You are to construct a discrete approximation to the function f (x) = sin(x) over the interval x [0, ] at n equally spaced points. Given this discrete approximation, you now seek to compute the integral of f (x) using composite integration rules. To accomplish this you will write one function and one main program: 1. A function numericalIntegral: function I = numericalIntegral(x,f,type) The inputs are the x locations in x, the function values f corresponding to x and a character array that denotes the type of derivative you seek. Supported types should 5-2 Lab 5: Numerical Differentiation and Integration Apr 3-7, 2017 be 'P', 'T', 'S3', for the pulse, trapezoidal and Simpson's 1/3 rule respectively. The output is the integral I of f (x) over the range of x provided. Remark: The function computes the integral over the entire range of x it receives as input. You may (and should) apply composite rules directly. However, in the case of Simpson's 1/3 rule the composite formula requires an even number of intervals. Your code should handle this case by applying a trapezoidal rule to the last interval when there is an odd number of intervals (or if you're really keen you can use Simpson's 3/8 rule). 2. A main program Lab5Q2: This main program should compute the integral of f (x) over the entire range for n = 20 and n = 100 using all three methods. Compute and tabulate the error in the integral over the entire range for each of the integration methods for both values of n. Also, generate the integral up to different values of x so that you can plot it. This can be accomplished by calling the numericalIntegral function using different ranges of x and f in a loop. For example, using x(1 : p) and f (1 : p) will give the integral value for the limits x(1) to x(p). That it, it will approximate: Zx(p) I(x(p)) = f (x)dx. a=x(1) Plot the integral for n = 20 using 'P' and 'S3' integration schemes compared to the analytic solution (you may choose to use the scatter function if the number of integral evaluations is small and it makes the plots clearer). Hand-In: The code for the function and the main program, a table of the requested simulation results and the requested plot. Make sure that the code is well documented and formatted. Question 3 Using the function you developed in Question 2 we will implement a main program that shows how Romberg integration works. To accomplish this you will write one main program: 1. A main program Lab5Q3: This main program should again compute the integral of f (x) = sin(x) over the interval x [0, ]. You will first compute this integral using a composite trapezoidal rule with n = 2, 4, 8 and 16 intervals (previously n was the number of points!!) to provide four initial integral estimates Ij,k at level k = 1, namely: I1,1 , I2,1 , I3,1 and I4,1 . You do not need to minimize the number of function evaluations, i.e., it is not necessary to pull the function evaluations for j = 1 from the j = 4 list etc. Given these 4 initial estimates, use Romberg integration to produce I1,4 by first constructing 3 level 2 estimates followed by 2 level 3 estimates and then the 5-3 Lab 5: Numerical Differentiation and Integration Apr 3-7, 2017 final level 4 estimate. Provide a table showing the value of the results for all estimates. Use the Matlab command format long to get enough significant figures to clearly show the improvement in accuracy. Hand-In: The code for the main program, and a table of the requested integration results. Make sure that the code is well documented and formatted. Question 4 From the previous question, based on the number of significant digits you get correct in I4,1 (compared to the analytic solution), how many intervals are required to achieve the same level of accuracy (in terms of significant digits) for a composite trapezoidal rule? A composite Simpson's 1/3 rule? (You may find the answers by trial and error.) Assuming that you only compute a given function value once (meaning you just need to count the number of function evaluations required to compute the largest j at k = 1), how much speedup (in terms of function evaluations) does Romberg integration provide in both cases? Note: When comparing values assume that the next digit is rounded. Therefore if your value from Romberg integration is 2.00003 (which has 5 correct significant digits as the analytic solution is 2.0), a composite trapezoidal rule that provides an answer of 1.99999 would have the correct number of significant digits because if this is rounded we get 2.0000. 5-4