MATH2033 Introduction to Scientific Computation
— Coursework 4 —
Submission deadline: 15:00 Friday 3 May 2024
Note: This is currently version 4 of the Coursework 4 pdf.
No more questions will be added.
This coursework contributes 10% towards the overall mark for the module.
Rules:
• Each student is to submit their own coursework.
• You are allowed to work together and discuss in small groups (2 to 3 people), but you must write your
own coursework and program all code by yourself.
• You are allowed and certainly encouraged to use the forum on the MATH2033 Moodle page to ask
questions to obtain clarification of the coursework questions, as well as general Python queries.
• However, please ensure that you are not revealing any answers to others. Also, please ensure that
you are not revealing any code that you wrote. Doing so is considered academic misconduct.
• When in doubt, ask Richard Rankin during MATH2033 classes, during his office hours or by email
(richard.rankin@nottingham.edu.cn).
• Please be informed of the Policy on academic misconduct (including plagiarism, false authorship,
and collusion).
Coursework Aim and Coding Environment:
• In this coursework you will develop Python code related to algorithms that solve numerical integration
and numerical ODE problems, and you will study some of the algorithms’ behaviour.
• You should write (and submit) plain Python code (.py), and you are strongly encouraged to use the
Spyder IDE (integrated development environment). Hence you should not write ipynb files, and you
should not use Jupyter.
How and Where to Run Spyder:
• Spyder comes as part of the Anaconda package (recall that Jupyter is also part of Anaconda). Here
are two options on how to run Spyder:
(1) You can choose to install Anaconda on you personal device (if not done already).
(2) Anaconda is installed on the computers in some UNNC computer labs.
Helpful Resources:
• Python 3.12.0 documentation – Spyder IDE website – NumPy user guide – Matplotlib Quick start
guide – Analytical and Computational Foundations (MATH1028 UNNC) (FCH1 22-23)
• You are expected to have basic familiarity with Python, in particular: logic and loops, functions, NumPy
and matplotlib.pyplot. Note that it will always be assumed that the package numpy is imported as np,
and matplotlib.pyplot as plt.
Some Further Advice:
• Write your code as neatly and readable as possible so that it is easy to follow. Add some comments
to your code that indicate what a piece of code does. Frequently run your code to check for (and
immediately resolve) mistakes and bugs.
Some Further Information:
• Your work will be mostly marked by running your functions and comparing their output against the
correct results.
Department of Mathematical Sciences Page 1 of 9
Getting Started:
• Download the contents of the “Coursework 4 Pack” folder from the MATH2033 Moodle page into a
single folder.
(More files may be added in later weeks.)
• Open Spyder.
I Numerical Integration
Suppose that f ∈ C2[a, b]. Let h = b− a
n
and let xi = a + ih for i = 0, 1, . . . , n.
The composite trapezium rule with n subintervals is∫ b
a
f(x) dx ≈ h
2
(
f(a) + 2
n−1∑
i=1
f(xi) + f(b)
)
which uses the values of f at n+ 1 evenly spaced points.
0 This is a warm up exercise. No marks will be awarded for it. Its solution will
appear on Moodle. The solution to this problem will be required in questions that
follow.
The file warmup.py contains an unfinished function with the following first line:
def composite_trapezium(a,b,m,f):
The function should return as output:
integral approx CTR - a floating point number of type numpy.float64.
• Complete the composite_trapezium function so that it returns the approx-
imation to
∫ b
a
f(x) dx obtained by the composite trapezium rule with f being
evaluated at m evenly spaced points.
• Once you have completed composite_trapezium, test it by running
main_warmup.py. You should obtain:
integral_approx_CTR =
1.9835235375094549
Coursework 4 Page 2 of 9
1(a) Suppose that f ∈ C4[a, b]. Let h = b− a
2n
and let xi = a+2ih for i = 0, 1, . . . , n.
The composite Simpson’s rule with n subintervals is∫ b
a
f(x) dx ≈ h
3
(
f(a) + 4
n−1∑
i=0
f(xi + h) + 2
n−1∑
i=1
f(xi) + f(b)
)
which uses the values of f at 2n+ 1 evenly spaced points.
The file numerical_integration.py contains an unfinished function with the
following first line:
def composite_simpson(a,b,m,f):
The function should return as output:
integral approx CSR - a floating point number of type numpy.float64.
• Complete the composite_simpson function so that it returns the approxima-
tion to
∫ b
a
f(x) dx obtained by the composite Simpson’s rule with f being evalu-
ated at m evenly spaced points.
• Once you have completed composite_simpson, test it by running main.py.
In this case, you should obtain:
integral_approx_CSR =
2.0001095173150043
(b) The file numerical_integration.py contains an unfinished function with the
following first line:
compute_errors_integral(M,f,a,b,true_val):
The function should return as output:
– errors_trapezium - a numpy.ndarray with its length determined by the
length of the input list M;
– errors_simpson - a numpy.ndarray with its length determined by the length
of the input list M;
– fig - a matplotlib.figure.Figure.
• Complete the compute_errors_integral function to compute the absolute
errors in the approximations obtained with the composite trapezium rule and the
composite Simpson’s rule for different numbers (stored in M) of evaluations of f
where the true value of the integral is known and provided in true_val. Hence,
your function should call composite_trapezium and composite_simpson.
• errors_trapezium and errors_simpson store the approximation errors
Coursework 4 Page 3 of 9
computed for different numbers of evaluations of f.
• The function should also plot the absolute errors against the number of evalua-
tions m for the composite trapezium rule and the composite Simpson’s rule using
a log-log scale and return this figure in fig.
• Once you have written compute_errors_integral, test it by running
main.py. In this case, you should obtain:
errors_trapezium =
[4.29203673e-01 1.64764625e-02 1.64496113e-04 1.64493434e-06]
errors_simpson =
[9.43951024e-02 1.09517315e-04 1.08245040e-08 1.08224540e-12]
[15 / 40]Marks can be obtained for your composite_Simpson and
compute_errors_integral functions for generating the required output, for
certain set(s) of inputs. The correctness of the following may be checked:
• The type of the outputs;
• The shape of the outputs;
• The values of the outputs;
• The plot produced in fig.
Note that in marking your work, different sets of inputs may be used.
Coursework 4 Page 4 of 9
I Numerical Integration - double integration
Let’s explore the double integral of f(x, y) over a rectangular region [a, b]× [c, d]:∫ d
c
∫ b
a
f(x, y) dx dy.
2
(a) To illustrate, we’ll start by using the composite trapezium rule in a simple scenario
with just 3 ∗ 3 = 9 evaluation points. This straightforward example is designed
for easy comprehension. It’s crucial to highlight that the principles we’ll cover
can be easily adapted to handle any number of subintervals. Let’s break down
the x-interval [a, b] into two subintervals with endpoints x0 = a, x1 = a + h and
x2 = a+2h. While keeping y constant, we can then easily approximate the inner
integral as: ∫ b
a
f(x, y) dx ≈ h
2
[
f(x0, y) + 2f(x1, y) + f(x2, y)
]
. (1)
Further, we divide the y-interval [c, d] into two subintervals with endpoints y0 = c,
y1 = c + k and y2 = c + 2k. By employing the composite trapezium rule to
approximate the integral of each term on the right hand side of (1), we obtain:
h
2
[
k
2
(
f(x0, y0) + 2f(x0, y1) + f(x0, y2)
)
+2
k
2
(
f(x1, y0) + 2f(x1, y1) + f(x1, y2)
)
+
k
2
(
f(x2, y0) + 2f(x2, y1) + f(x2, y2)
)]
.
Now, let’s adapt this approach for mx*my evaluation points.
In numerical_integration.py you will find an unfinished function with the fol-
lowing first line:
def double_integral_trapezium(f,rangex ,rangey ,mx,my,true_val):
The function should return as output:
– double_CTR - a floating point number of type numpy.float64.
– error_CTR - a floating point number of type numpy.float64.
• Complete the function so that it approximates the double integral of a two-
variable function over a rectangular region by utilizing the composite trapezium
rule. The output double_CTR should have the value of this approximation. The
output error_CTR should have the value of the absolute error in this approx-
imation, where the true value of the double integral is known and provided in
true_val.
• rangex: This is a list with two elements defining the lower and upper limits of
Coursework 4 Page 5 of 9
integration for the variable x.
• rangey: Similar to rangex, this is a list with two elements that sets the lower
and upper integration bounds for the variable y.
• The input mx is an integer that determines the number of evenly spaced points
in the interval [a, b] that will be the x-coordinates of the points (x, y) at which f
will be evaluated.
• The input my is an integer that determines the number of evenly spaced points
in the interval [c, d] that will be the y-coordinates of the points (x, y) at which f
will be evaluated.
• In total, f will be evaluated at mx*my points.
• Hint : within double_integral_trapezium you could define an inner inte-
gral function, which computes the composite trapezium rule along the x-axis for
a fixed value of y. It then uses this inner integral as part of another composite
trapezium rule computation for the y-axis, effectively nesting two applications of
the rule to approximate the double integral. To complete the calculation, your
function should call composite_trapezium.
• Once you have written double_integral_trapezium, test it by running
main.py. In this case, you should obtain:
Double Integral Result for CTR: 3.4375
(b) The double integral
∫ d
c
∫ b
a
f(x, y) dx dy can be approximated using different
methods. In 2(a), we utilized the composite trapezium rule. Now, in 2(b), let’s
explore an alternative method: the composite Simpson’s rule. In this case, with
3 ∗ 3 = 9 evaluation points, (1) is redefined as∫ b
a
f(x, y) dx ≈ h
3
[
f(x0, y) + 4f(x1, y) + f(x2, y)
]
. (2)
Extending this method for double integration involves approximately integrating
each term on the right hand side of (2) using the composite Simpson’s rule,
resulting in the following expression:
hk
9
[(
f(x0, y0) + 4f(x0, y1) + f(x0, y2)
)
+ 4
(
f(x1, y0) + 4f(x1, y1) + f(x1, y2)
)
+
(
f(x2, y0) + 4f(x2, y1) + f(x2, y2)
)]
.
Next, generalize this approach to accommodate mx*my evaluation points.
In numerical_integration.py you will find an unfinished function with the fol-
lowing first line:
Coursework 4 Page 6 of 9
def double_integral_simpson(f,rangex ,rangey ,mx,my,true_val):
The function should return as output:
– double_CSR - a floating point number of type numpy.float64.
– error_CSR - a floating point number of type numpy.float64.
• Complete the function so that it approximates the double integral of a two-
variable function over a rectangular region by utilizing the composite Simpson’s
rule. The output double_CSR should have the value of this approximation. The
output error_CSR should have the value of the absolute error in this approx-
imation, where the true value of the double integral is known and provided in
true_val.
• Other inputs are as in 2(a).
• Once you have written double_integral_simpson, test it by running
main.py with your own tests.
[12 / 40]Marks can be obtained for your double_integral_trapezium and
double_integral_simpson functions for generating the required output, for
certain set(s) of inputs. The correctness of the following may be checked:
• The type of the outputs;
• The values of the outputs.
Note that in marking your work, different sets of inputs may be used.
Coursework 4 Page 7 of 9
I Runge-Kutta Schemes
Suppose we wish to solve the system of m first-order initial value problems
dy
dt
= f(t,y), a ≤ t ≤ b, (3)
y(a) = y0. (4)
3 In numerical_integration.py you will find an unfinished function with the fol-
lowing first line:
def runge_kutta(a,b,f,N,y0 ,m,method):
The function should return as output:
– t - a numpy.ndarray of shape (N+1,);
– y - a numpy.ndarray of shape (m,N+1).
• Complete the function so that for a given m-vector function f(t,y) (stored in
f), it will find approximate solutions to (3)-(4) at time points
{
ti = a + ih
}N
i=0
(output in t), where N is the number of time steps to take and h = (b − a)/N .
The initial value y0 is supplied in y0.
• On return, the jth row of y should contain the approximations to the jth com-
ponent of y at all the time points.
• Finally, method selects the Runge-Kutta method to be used and can take the
following values:
– method = 1 - Forward Euler;
– method = 2 - Modified Euler;
– method = 4 - RK4.
• Also add a brief description at the top of your definition, the so-called
docstring. This description should become visible whenever one types:
help(ni.runge_kutta) or ni.runge_kutta?.
Coursework 4 Page 8 of 9
• Test your runge_kutta function by running main.py, which should give the
following output:
In all cases you should obtain:
t = [0. 0.2 0.4 0.6 0.8 1. ]
Method 1:
y = [[ 1. 1.2 1.48 1.875328 2.4469248 3.29496678]
[ 0. 0. -0.048 -0.18528 -0.4653696 -0.960192 ]
[ 1. 1. 1.0016 1.02016 1.0781824 1.19919718]]
Method 2:
y = [[ 1. 1.24 1.5997568 2.15615444 3.04208295 4.47877757]
[ 0. -0.024 -0.1319232 -0.38732212 -0.8841379 -1.76914977]
[ 1. 1.0008 1.0149984 1.06320417 1.16660459 1.33548735]]
Method 4:
y = [[ 1. 1.24630727 1.62092974 2.20871137 3.15704679 4.71318699]
[ 0. -0.02720673 -0.14334098 -0.41704869 -0.95227519 -1.9148448 ]
[ 1. 1.00181147 1.01795655 1.06935282 1.17695518 1.34906871]]
[13 / 40]Marks can be obtained for your runge_kutta function for generating the required
output, for certain set(s) of inputs. The correctness of the followingmay be checked:
• The type of the outputs;
• The shape of the outputs;
• The values of the outputs;
• The output of help(ni.runge_kutta).
Note that in marking your work, different sets of inputs may be used.
I Finished!
Submission Procedure:
To submit, upload your numerical integration.py file through the Coursework 4 as-
signment activity in the Coursework 4 section of the MATH2033 Moodle page.
Coursework 4 Page 9 of 9