MTHS2008 Advanced Mathematical Modelling for Process Engineers - COURSEWORK 2: DOUBLE NUMERICAL INTEGRATION AND FINITE
1 MTHS2008
SPRING SEMESTER 2023-2024
ASSESSED COURSEWORK 2 – DOUBLE NUMERICAL INTEGRATION AND FINITE
DIFFERENCE METHODS
Submission deadline: 5pm, Friday 22 March 2023
This coursework contributes 20% towards the overall grade 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 4 people), but you must
write all code by yourself.
You must submit only the .py files requested in this question paper. Details of the required filenames are given within the questions. You are strongly encouraged to use the Spyder IDE (integrated development environment). Hence you should not write IPython Notebooks (.ipynb), and you should not use Jupyter.
You may adapt any code that we developed in class (i.e. you do not have to start from scratch).
In each template file, the packages and libraries required to complete the question are listed. You are not permitted to use any libraries or packages beyond those given for each question.
MTHS2008 Turn Over
Marks breakdown:
2 MTHS2008
This coursework is out of 100 marks:
• Outputs in Q1a – 20 marks;
• Outputs in Q1b – 30 marks;
• Outputs in Q2a – 20 marks;
• Outputs in Q2b – 20 marks;
• Commenting and structure – 10 marks.
To obtain the maximum marks for commenting and structure:
-
Your comments should be helpful to the reader; they should help make your program easier to navigate. Note that having too many comments is potentially as bad as having too few comments! The style of commenting we used in the example programs is what you should aim for.
-
Your program structure should be: imports at the top (if any), followed by function definitions, followed by anything else. You ought to remove or comment out any lines of code used to call and test your functions prior to submission; however, you will not be penalised if you leave them in, provided that they are in the correct place (e.g. not sandwiched in-between two function definitions).
Guidelines for submitting the coursework correctly:
For full marks, your functions need to return the correct outputs for particular input arguments. Suggestions for tests you can use to check your code with are given but you are encouraged to test your code using your own examples also.
Please be aware that we will test your functions to check that they work and produce the desired output(s), both with the test given in the question and with different undisclosed data.
If your modules have filenames that differ to what we have explicitly asked for, then you risk losing marks.
• Therefore, please do not add your username or student ID number to your filename. Your functions must have the same inputs and outputs as those specified in the question(s),
and in the same order as specified.
MTHS2008 Turn Over
3 MTHS2008 1. Question 1a: Plotting region of integration for a double integral
Using template_question_1a.py (on Moodle) as a template, write a module containing a function plot_region(a,b,y1,y2,n) to create a figure containing the plot of the region of integration for the double integral
𝑏 𝑦2(𝑥)
𝑓(𝑥,𝑦)d𝑦d𝑥. 𝑎 𝑦1(𝑥)
You function should return the figure as an output. To do so, give the figure a ‘handle’ when you create it, i.e. assign it to a variable. For example
fig = plt.figure()
You can then return the figure to the user using the command return fig.
Your module file must be named q1a.py and your function should return one output - the figure containing the plot of the region of integration.
Test: To test that your code is working as expected, plot the region of integration for the double integral
2𝑥
∫ ∫ 𝑥𝑦2d𝑦d𝑥
1 −𝑥
with 𝑛 = 12. You should be able to verify the accuracy of your plot by sketching the region
of integration by hand.
∫∫
[20 marks] Using template_question_1b.py (on Moodle) as a template, write a module containing a
function double_int_approx(a,b,y1,y2,f,n,m) to approximate the double integral
𝑏 𝑦2(𝑥)
Question 1b: Approximating a double integral
𝑓(𝑥,𝑦)d𝑦d𝑥, 𝑎 𝑦1(𝑥)
using the quadrature rule where Simpson’s rule has been applied in both coordinate directions. The quadrature rule was given in the first Python session of Spring semester.
The function should take as its parameters
-
𝑎, 𝑏, 𝑦1(𝑥), 𝑦2(𝑥) as the integral limits that define the region of integration;
-
𝑓 as the integrand;
-
𝑛 and 𝑚 are the number of strips to be used in the Simpson’s approximation in the 𝑥-direction and 𝑦-direction respectively.
Your function should assert that any user gives values of 𝑛, 𝑚 that are even as required by Simpson’s rule.
Your module file must be named q1b.py and your function should return one output - an approximation of the double integral.
Test: To test that your code is working, compare your approximation of ∫2 ∫𝑥 𝑥𝑦2 d𝑦 d𝑥
∫∫
to the true value of the double integral: 62. 15
1 −𝑥
MTHS2008
Turn Over
[30 marks]
4 MTHS2008 2. Question 2a: Poisson’s equation with a Neumann boundary conditions
Using template_question_2a.py (on Moodle) as a template, write a module containing a function poisson_neumann_approx(N,a,b,alpha,beta,f) that uses central finite differences to find numerical solutions to the 1D Poisson equation on [𝑎, 𝑏] subject to boundary conditions:
−𝑑2𝑢 = 𝑓(𝑥), 𝑑2𝑥
Your module file must be named q2a.py and your function should return one output - an array containing approximations 𝑢0, 𝑢1, 𝑢2, ⋯, 𝑢𝑁−2, 𝑢𝑁−1, 𝑢𝑁.
Test: To test that your code is working as expected, compute approximations for − 𝑑2𝑢 = 4𝜋2 cos(2𝜋𝑥),
[25 marks]
𝑢′(𝑏) = 𝛽.
𝑢(𝑎) = 𝛼,
The parameter 𝑁 acts as the number of sub-intervals considered in the mesh of [𝑎, 𝑏].
𝑑2𝑥
The true solution to this PDE is 𝑢(𝑥) = cos(2𝜋𝑥).
𝑢(0) = 1,
𝑢′(1) = 0.
Question 2b: Time independent advection-diffusion-reaction equation
Using template_question_2b.py (on Moodle) as a template, write a module containing a function advection_1d_solver(N,a,b,f) that uses central differencing approximations to find numerical solutions to the 1D time-independent advection-diffusion-reaction equation on [0, 1] with boundary conditions:
−𝑎𝑑2𝑢 + 𝑏𝑑𝑢 = 𝑓(𝑥), 𝑑𝑥2 𝑑𝑥
𝑢(0) = 0,
where 𝑎 > 0 and 𝑏 are scalar constants. The parameter 𝑁 acts as the number of sub-intervals
𝑢(1) = 0, considered in the mesh of [0, 1].
Your module file must be named q2b.py and your function should return one output - a plot of your solution approximation.
Test: To test that your code is working as expected, plot your approximate solution for −𝑑2𝑢 + 2𝑑𝑢 = 1,
𝑑2𝑥 𝑑𝑥
𝑢(0) = 0, 𝑢(1) = 0.
ThetruesolutiontothisPDEis𝑢(𝑥)=𝑥+ 2
1 2(𝑒2−1) (
1−𝑒2𝑥 . )
MTHS2008
End
[25 marks]
