ENG5322M程序讲解、Python编程语言辅导、programming程序设计讲解

日期：2020-11-12 08:36

Course assignment
ENG5322M: Introduction to computer programming
(Additional Hints Provided)
October 26, 2020
Instructions
You should submit the following through the ENG5322M Moodle submission system by 4pm,
27th November 2020:
1. Python notebook (.ipynb format) used for to solve each of the assignment problems.
This can be a separate notebook for each problem or a single file.
You will be marked according to the following criteria:
? Program correctness (60%): Does your program solve the problem given to you?
? Program readability (20%): Have you written your program in such a way that it is
easy to follow? Have you included appropriate comments?
? Testing (20%): Have you done appropriate testing to show your program works as
intended?
You are given two problems to solve for your project assignment which will be outlined in more
detail shortly. The two problems are focussed on:
1. Numerical integration
2. Multivariate least square method
Notebook format
I expect to see you explain your code throughout your notebook both using text and images,
much like I have down in the lectures during this course. As a quick guide, the default type of a
cell is code, but if you change it to ’markdown’ (look in the toolbar) and then hit shift and enter
you will see it become text in the notebook.
1
1 PROBLEM 1: NUMERICAL INTEGRATION
1. Testing: You should write appropriate code that performs tests on your program.
1 Problem 1: Numerical integration
Overview of task: Write a program that calculates the integral of a function using both the
trapezoidal rule and the rectangular rule. That is, we want to calculate
I =
Z b
a
f(x) dx (1)
which is the area under the curve as shown graphically in Fig. 1.
Figure 1: Graphical representation of Eq. (1)
Specifics: You are required to write a program that first defines variables of the the limits of
the integral (i.e. the values a and b in Eq. (1) ) and the number of integral divisions (explained
shortly). Please note, it is not acceptable to simply use a library function to compute these
integrals - I want to see you implement the algorithms for both the rectangular and trapezoidal
rule.
Your program should then calculate the integral for f(x) = x
2
, f(x) = sin x and f(x) = e
?x
sin x
using both rules. That is, calculate
To calculate these values, we can approximate them using numerical integration. A key point
to understand is that we almost always approximate the integral using these methods - there is
usually an error.
Finally, I want you to output a file for each integral and each rule which gives the result for N=2
(i.e. 2 integral divisions) up to N=100.
2
1.1 Rectangular rule 1 PROBLEM 1: NUMERICAL INTEGRATION
Figure 2: Rectangular rule for number of divisions N = 2
1.1 Rectangular rule
This works by dividing the area under the curve into rectangles as shown in Fig. 2
Let’s take an example of how to apply this to calculate the integral R 1
0
x
3dx for N = 2. We first
calculate the step size h = (b ? a)/N = (1 ? 0)/2 = 0.5. We now apply the rectangular rule to
get
approximate integral = f(x1)h + f(x2)h
= f(0) × 0.5 + f(0.5) × 0.5
= 0 × 0.5 + 0.125 × 0.5
= 0.0625 (2)
We can also do the same for N = 4 (See Fig. 3). h = (1 ? 0)/4 = 0.25.
approximate integral = f(x1)h + f(x2)h + f(x3)h + f(x4)h
= f(0) × 0.25 + f(0.25) × 0.25 + f(0.5) × 0.25 + f(0.75) × 0.25
= 0 × 0.25 + 0.015625 × 0.25 + 0.125 × 0.25 + 0.421875 × 0.25
= 0.140625 (3)
And for N = 8 we get 0.19140625, N = 50 we get 0.2401, N = 100 we get 0.245025 and
N = 1000 we get 0.24950025.
If you were to work this integral out by hand, you would see that the exact answer is 0.25. Notice
how we converge to this answer as we increase N.
The general formula for the rectangular rule is as follows:
approximate integral = h [f(x1) + f(x2) + · · · + f(xN )] (4)
where h = (b ? a)/N.
3
1.2 Trapezoidal rule 1 PROBLEM 1: NUMERICAL INTEGRATION
Figure 3: Rectangular rule for number of divisions N = 4
1.2 Trapezoidal rule
The trapezoidal rule approximates the integral by constructing trapezoids under the function.
This is shown in Fig. 4 for N = 2 and Fig. 4 for N = 4. It is clear that this provides a more
accurate approximation the integral.
The general rule for trapezoidal rule is as follows. Letting h = (b ? a)/N, the integral can be
approximated as:
approximate integral = h
2
[f(x1) + 2f(x2) + 2f(x3) + 2f(xN ) + f(xN+1)] (5)
Let’s apply this to the previous example of f(x) = x
3
. For N = 2 with h = (1 ? 0)/2 = 0.5 we
have
approximate integral = h
2
[f(x1) + 2f(x2) + f(x3)]
=
0.5
2
[f(0) + 2f(0.5) + f(1.0)]
= 0.25[0 + 0.25 + 1]
= 0.3125 (6)
and for N = 4 with h = (1 ? 0)/4 = 0.25 we have
approximate integral = h
2
[f(x1) + 2f(x2) + 2f(x3) + 2f(x4) + f(x5)]
=
0.25
2
[f(0) + 2f(0.25) + 2f(0.5) + 2f(0.75) + f(1.0)]
= 0.125 [0 + 2 × 0.015625 + 2 × 0.125 + 2 × 0.421875 + 1]
= 0.265625 (7)
and for N = 8 we have 0.25390625, N = 50 we have 0.2501, N = 100 we have 0.250025 and
N = 100 we have 0.25000025.
1.3 Testing
The obvious way to test your program is to compare your results against integrals you have
performed by hand. For example, in the results shown above we know the exact result to be 0.25
4
2 PROBLEM 2: MULTIVARIATE LEAST SQUARES METHOD
Figure 4: Trapezoidal rule for number of divisions N = 2
Figure 5: Trapezoidal rule for number of divisions N = 4
and we can see that they are converging to the exact result as N is increased.
Another test might be to integrate a function that is a constant. E.g. what is the area under the
function f(x) = 1?
2 Problem 2: Multivariate least squares method
We want to find the regression coefficients a0, a1,...,am for the following equation:
yi = a0 + a1x1i + a2x2i + a3x3i + ... + amxmi (8)
To solve this, let us suppose, we take a case for m = 2 and add an error term i to our equation:
yi + i = a0 + a1x1i + a2x2i (9)
5
2 PROBLEM 2: MULTIVARIATE LEAST SQUARES METHOD
We can then obtain an objective function S to minimize as:
(12)
Minimizing S with respect to a0, a1, and a2, we get:
?S
?a0
= 2X(a0 + a1x1i + a2x2i ? yi) (13)
?S
?a1
= 2X(a0 + a1x1i + a2x2i ? yi)x1i (14)
?S
?a2
= 2X(a0 + a1x1i + a2x2i ? yi)x2i (15)
Equating these to zero, we obtain three equations in three unknowns, a0, a1, a2:
a0 + a1 < x1 > +a2 < x2 >=< y > (16)
a0 < x1 > +a1 < x1x1 > +a2 < x1x2 >=< x1y > (17)
a0 < x2 > +a1 < x1x2 > +a2 < x2x2 >=< x2y > (18)
where < . > represents average quantity defined as:
We can represent in matrix form to use matrix methods as:
or XA = Y . The matrix of unknown coefficients, A, is obtained by finding the inverse of X,
so that A = X?1Y .
6
2.1 Hint 2 PROBLEM 2: MULTIVARIATE LEAST SQUARES METHOD
Extending to m variables, we get:
(21)
Overview of task: The assignment is to use the above logic to find the coefficients a0,a1,a2
for the data given in Table 1 (also provided as a tab-delimited text file). The exact solution is
y = ?2.145 ? 3.117x1 + 1.486x2.
2.1 Hint
To program this efficiently in python, we can make use of a trick. Imagine we get the values of x
and y from file, and since we are reading one line at a time, rather than calculating the averages,
we can simply maintain the sums by taking 1/n out and cancelling it on both sides:
? (22)
where i is an individual record (row). For example, when i = 1, x1i = 0.408, x2i = 0.58, and
yi = ?2.39. Similarly, when i = 2, x1i = 0.429, x2i = 0.51, and yi = ?2.53, and so on.
The summations are as follows:
Pn
i=1 x1i = 0.408 + 0.429 + 0.492 + 0.529 + ... + 1.655 + 1.684 + 1.897
Pn
i=1 x1ix2i = (0.408 × 0.58) + (0.429 × 0.51) + (0.492 × 0.53) + ... + (1.655 × 1.2) + (1.684 ×
0.96) + (1.897 × 1.24)
If you can import all the function from numpy library from numpy import *, and if you have a
variable m as the number of independent variables, and if you can define x=zeros((m+1,m+1)),
y=zeros((m+1,1)), and xi=[0.0]*(m+1), then you can populate both matrices x and y in a
nested-loop with two equations, x[j,k] += xi[j]*xi[k] and y[j,0] += yi*xi[j], provided if
you have xi[0]=1, xi[1] the value of first predictor (first column), xi[2] the value of second
predictor (second column), and so on. yi is dependent variable (last column) in Table 1. Finally,
you can simply convert numpy arrays x and y to numpy matrix so that you can get the inverse
automatically: X=mat(x.copy()), Y=mat(y.copy()), and A=X.I*Y.
Also, the data file provided is tab-delimited. Look up how to split a string using .split()
function.
7
2.1 Hint 2 PROBLEM 2: MULTIVARIATE LEAST SQUARES METHOD
Table 1: Data to be fitted
x1 x2 y
0.408 0.58 -2.39
0.429 0.51 -2.53
0.492 0.53 -3.38
0.529 0.6 -2.72
0.569 0.58 -2.95
0.677 0.64 -3.32
0.703 0.61 -3.45
0.841 0.66 -3.81
0.911 0.73 -3.9
0.936 0.72 -4.11
0.96 0.57 -4.24
0.978 0.84 -4.81
0.993 0.69 -4.1
1.038 0.75 -4.28
1.051 0.76 -4.05
1.119 0.85 -4.23
1.134 0.75 -4.58
1.16 0.85 -4.4
1.209 0.72 -5.05
1.272 0.82 -4.79
1.303 0.86 -4.76
1.353 0.98 -4.65
1.367 0.9 -5.18
1.388 0.71 -5.5
1.425 1.04 -5.01
1.453 0.9 -5.15
1.484 0.77 -5.8
1.503 0.93 -5.33
1.536 0.81 -5.23
1.563 0.83 -6.11
1.655 1.2 -5.34
1.684 0.96 -6.13
1.897 1.24 -6.45
8
2.1 Hint 2 PROBLEM 2: MULTIVARIATE LEAST SQUARES METHOD
Figure 6: Pictorial representation of the overall process
9