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
equired 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
easonably 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 e
ors 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 e
ors due to this factor of 10 reduction in the step size h. How does the e
or
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
easonable (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 e
or 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)
etween x = a and x = 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
([ 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 XXXXXXXXXXf (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
a
= ln
+1
+ 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 e
or E = I - Itrap.
3.2. SIMPSON’S