I have attached the file for the assignment.
Department of Physics: 1CL Assignment 1 PHY1038 Mathematical and Computational Physics Assignment 3 Finite Difference Methods SUBMISSION DEADLINE: 16:00 Tuesday 24th March 2020 1. INTRODUCTION In physics problems, we often need to calculate the first derivatives of smooth continuous functions, that is, their rate of change. It is however not always possible or perhaps difficult to write down a simple analytic form for the function f(x) although it may be perfectly possible to work out the value of the function for any required value of the variable x. For instance, f(x) may be obtained from a numerical integration or by solving numerically a differential equation. Similarly, we may have a function that we would like to integrate but which either has no analytic closed-form integral, or we may only know numerical values of the function at a series of points - perhaps some measured data values. In these circumstances, it is convenient to have reliable techniques, or algorithms, with which to estimate derivatives and definite integrals of the function f(x) from our knowledge of the function’s values at a finite number of points x. The aim of this assignment is to become familiar with the use of such numerical techniques for calculating numerical derivatives and integrals. 2. NUMERICAL APPROXIMATION OF DERIVATIVES 2.1 Forward Difference Method One estimate of the derivative (i.e. the slope of the tangent) of a function f at the point x is obtained by taking a small value of h and calculating Dfor (x,h) = [ f (x + h) − f (x)] h . Provided the step h is small enough, then for reasonably smooth functions this should provide a reasonable estimate for f ′(x). 2.2 Centred Difference Method Another reasonable estimate of the derivative of the function f at the point x comes from choosing a small value of the step length h and then evaluating Dcen (x,h) = [ f (x + h /2) − f (x − h /2)] h . Since this uses knowledge of the function symmetrically about the point x at which one is trying to find the derivative, one might expect that Figure 1: graphical representation of derivative Figure 2: f(x) = 1/(x2+3x+2) this approximation would be more accurate than the forward difference method. This is the case, and you should be able to show this by the use of Taylor series expansions of the formulas above for small values of h. We will test the accuracy of these two methods by obtaining the numerical first derivative of f (x) = 1 x 2 + 3x + 2 , a case for which we can also calculate the derivative exactly. The function y=f(x) is plotted in Figure 2 with a negative slope for all positive x. Step 1. Write a Python function functionf(x) which returns values of the function f(x) that is to be differentiated. Step 2. Write a function dfor(x,h) and a function dcen(x,h) which evaluate the forward and centred difference approximations given above to compute the first derivative at x using a step-length parameter h. Step 3. Write a function that returns the exact derivative of the function f(x). Compute the exact derivative on paper first and then code it. Step 4. In your program ask the user to input values of x and h, evaluate Dfor(x,h) and Dcen(x,h), and the errors in these approximations compared to the exact derivative. Run your program for the values x=0.0, x=1.0 and x=2.0, all with h=0.1. Step 5. Repeat the calculations of Step 4, but with h=0.01 and note the changes in the errors due to this factor of 10 reduction in the step size h. How does the error depend on h? Put all the functions in a module which should be imported in your program, as described in Unit 5 of last semester’s Python Workbook. 2.3 NUMERICAL SECOND DERIVATIVES Now we have an approximate way of calculating first derivatives numerically with reasonable (and controllable) accuracy. Using repeated application of the centred difference first derivative we can now obtain an expression for the second derivative; Dcen 2 (x,h) = f (x + h /2) − f (x − h /2) h Dcen (x + h /2) − Dcen (x − h /2) h Using the expressions we already have for Dcen, these give, upon substitution, Dcen 2 = f (x + h)− f (x)( )− f (x)− f (x − h)( ) h2 = f (x + h)− 2 f (x)+ f (x − h) h2 Either of these two forms (in terms of Dcen or directly in terms of the function) can be used, but, since you have already coded Dcen it makes sense to use this function. Step 6. Add to the module a function d2cen(x,h) which uses your existing function dcen(x,h) to estimate the second derivative as Dcen 2 (x,h) = Dcen (x + h /2,h) − Dcen (x − h /2,h) h . As well, add to the module a function which calculates the exact second derivative of the function f(x). Modify your program to evaluate the approximation to the second derivative, and compare it to the exact, analytic second derivative. Run your program with x=0.0, x=π/4 and x=π/2, and with steps h=0.1 and h=0.01, as for the first derivative. Note the difference in error between the different step sizes. 3. NUMERICAL INTEGRATION The numerical evaluation of one– or multi–dimensional integrals is very common in physical problems. A wide variety of techniques (algorithms) are available to evaluate such integrals to any required precision. It is far easier to obtain accurate integrals than derivatives as we shall see. We consider the evaluation of definite integrals in one dimension on the range a ≤ x ≤ b, that is, the area under the graph of f(x) between x = a and x = b, = b a dxxfI )( . Rather than deal with the integral from a to b in one bite, we will divide the required interval (b-a) into N equal sub-intervals, called bins, each of width h=(b-a)/N, and consider the area under the graph from each bin. By controlling the number of intervals N, or equivalently the step size h, we aim to obtain an accurate estimate of the integral without having to use too large a value of N. You will evaluate an integral numerically using two different methods – the Trapezium rule and Simpson’s rule. Each will be investigated for its accuracy considering a function whose integral can be computed exactly so that we can compare our results with the analytic solution. 3.1 THE TRAPEZIUM RULE (linear approximation) Consider the area under the graph of the function f in the range of values from x to x+h, that is the values x+z where 0 ≤ z ≤ h. The trapezium rule arises when we consider one of the simplest approximate ways to describe the function f(x) within this small range of x values, namely approximating f(x+z) by a straight line which passes through the points (x,f(x)) and (x+h,f(x+h)) (see the figure). The equation of this straight line is f lin (x + z) = zf (x + h) h − (z − h) f (x) h which is clearly linear in z and passes through the end points, as can be seen by substituting z=0 and z=h. The area under the straight line is then the area of the trapezium, which by simple geometry we find to be h[f(x)+f(x+h)]/2. If we then add up the area of all the trapeziums between our limits we get f (x)dx = h 2a b ([ f (a) + f (a + h)]+ [ f (a + h) + f (a + 2h)]+ +[ f (a + [N − 2]h) + f (a + [N −1)h)]+ [ f (a + (N −1]h) + f (b)]) = h 2 [ f (a) + 2 f (a + h) + 2 f (a + 2h) + + 2 f (a + [N −1]h) + f (b)]. This yields the trapezium rule approximation to the integral Itrap = h 2 [ f (a) + f (b)]+ h f (a + jh) j=1 N−1 Step 7 Write a Python program which evaluates the integral 1 x 2 + 3x + 2 dx 0 1 using the trapezium rule expression Itrap. Your program should contain a function, f(x), which defines the expression to be integrated, variables a and b for the limits of integration and variable for the number N of bins. Compare the result from the trapezium rule with the exact analytic result: I = 1 x 2 + 3x + 2 dx = ln x +1 x + 2 a b a b = ln b +1 b + 2 − ln a +1 a + 2 . Compare the answers for N = 10, N = 20 and N = 40. Use a = 0 and b = 1. In each case, evaluate the error E = I - Itrap. 3.2. SIMPSON’S