Course assignment ENG2079/ENG5322: Introduction to computer programming (Additional Hints Provided) February 21, 2021 Instructions You should submit the following through the ENG2079/ENG5322 Moodle submission system by 4pm, 26th March 2021: 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. Newton-Raphson method 2. Conjugate gradient 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: NEWTON RAPHSON METHOD 1. Testing: You should write appropriate code that performs tests on your program. 1 Problem 1: Newton Raphson method A common problem in engineering analysis is this: given a function f(x) determine the values of x for which f(x) = 0. The solutions (values of x) are known as the roots of the equation f(x) = 0, or the zeroes of the function f(x). The equation y = f(x) (1) contains three elements: an input value x, an output value y, and the rule f for computing y. The Newton-Raphson algorithm is the best-known method of finding roots for a good reason: it is simple and fast. The only drawback of the method is that it uses the derivative f ′(x) of the function as well as the function f(x) itself. Therefore, the Newton-Raphson method is usable only in problems where f ′(x) can be readily computed. The Newton-Raphson formula can be derived from the Taylor series expansion of f(x) about x: f(xi+1) = f(xi) + f ′(xi)(xi+1 xi) +O(xi+1 xi)2 (2) where O(z) is to be read as ”of the order of z”, i.e., the higher order terms. If xi+1 is a root of f(x) = 0, the above equation becomes: 0 = f(xi) + f ′(xi)(xi+1 xi) +O(xi+1 xi)2 (3) Assuming that xi is a close to xi+1, we can drop the last term and solve for xi+1. The result is the Newton-Raphson formula: xi+1 = xi f(xi) f ′(xi) (4) If x denotes the true value of the root, the error in xi is Ei = x xi. It can be shown that if xi+1 is computed from Eq. 4, the corresponding error is: Ei+1 = f ′′(xi) 2f ′(xi) E2i (5) indicating that Newton-Raphson method converges quadratically (the error is the square of the error in the previous step). As a consequence, the number of significant figures is roughly doubled in every iteration, provided that xi is close to the root. The graphical representation is shown in Fig. 1. The algorithm for the Newton-Raphson is simple and can be easily implemented in Python. It repeatedly applies Eq. 4, starting with an initial value x0, until the convergence criteria is reached: |xi+1 xi| < (6) where being the error tolerance. The algorithm for Newton-Raphson is then simple: 2 1 PROBLEM 1: NEWTON RAPHSON METHOD Figure 1: Graphical interpretation of the Newton-Raphson formula which approproximates f(x) by the straight line that is tangent to the curve at xi, and in doing so, xi+1 is aat the intersection of the x-axis and the tangent line. Figure 2: The pseudo-algorithm for Newton-Raphson method. You are going to compute the root of f(x) = x3 10x2 + 5 = 0 that lies close to x = 0.7. The derivative of the function is f ′(x) = 3x2 20(x), so that the Newton-Raphson formula in Eq. 4 is x← x f(x) f ′(x) = x x 3 10x2 + 5 3x2 20x = 2x3 10x2 5 x(3x 20) (7) It takes only two iterations to reach five decimal place accuracy: x← 2(0.7)3 10(0.7)2 5 0.7[3(0.7) 20] = 0.73536 x← 2(0.73536)3 10(0.73536)2 5 0.73536[3(0.73536) 20] = 0.73460 (8) In addition, you are required to find the roots of the following equations, f(x) = x3 2x 5 = 0 (9) that has a root near x = 2. And, f(x) = x 2sinx = 0 (10) that has a root near x = 1.5. In the second case, show a scenario where the Newton-Raphson will not work by choosing a different starting point.You can improve solving Eq. 10 by rewriting 2/x 1/sinx = 0, and if x0 is any number in (0, pi), Newton quickly takes us to the root. Reformulation as (sinx)/x 1/2 = 0 also works nicely. 3 1.1 Bonus question for extra credit 2 PROBLEM 2: CONJUGATE GRADIENT METHOD 1.1 Bonus question for extra credit A loan of A pounds is repaid by making n equal monthly payments of M pounds, starting a month after the loan is made. It can be shown that if the monthly interest rate is r, then Ar = M ( 1 1 (1 + r)n ) (11) A car loan of ￡10,000 was repaid in 60 monthly payments of ￡250. Use the Newton Method to find the monthly interest rate correct to 4 significant figures.You can seek help from https: //www.mathway.com/popular-problems/Calculus/526998 to calculate the derivative. 1.2 Hint To program this in python, you can go about defining two function: def f(x): return def df(x): return And can then make a function, def newtonRaphson(x,tol=1.0e-9): Additionally, you can plot f(x) and show that the function is converging. 2 Problem 2: Conjugate gradient method Overview of task: Write a program that solves equations of the form Ax = b using Conjugate Gradient Method. Use this to solve the equation: 4 1 1 1 4 2 1 2 4 x1x2 x
3 = 12 1 5 (12) 2.1 Derivation The problem involves finding the vector x that minimizes the scalar function f(x) = 1 2 xTAx bTx (13) 4 2.1 Derivation 2 PROBLEM 2: CONJUGATE GRADIENT METHOD where the matrix A is symmetric and positive definite. Because f(x) is minimized when its gradient f = Ax b, we see that the minimization is equivalent to solving Ax = b (14) Gradient methods solve this problem of minimization through iteration starting with an initial estimate of vector x0 and then iterate for k computations using xk+1 = xk + αksk (15) such that the step length αk is chosen so that xk+1 minimizes f(xk+1) in the search direction sk such that xk+1 satisfies A(xk + αksk) = b (16) If we now introduce a residual rk = b Axk (17) and pre-multiplying both sides by sTk and solving for αk gives αk = sTk rk sTkAsk (18) While intuitively we can choose sk = f = rk as this is the direction of the largest negative change in f(x) (called method of steepest descent), it converges very slow. A better choice is to use conjugate gradient method that uses the search direction sk+1 = rk+1 + βksk (19) The constant βk is then chosen so that the two successive search directions are conjugate to each other, i.e., sTk+1Ask = 0 (20) Subsituiting Eq. 20 into Eq. 19 gives: βk = rTk+1Ask sTkAsk (21) The algorithm is then as follows: It can be shown that the residual vectors r1, r2, r3, … are mutually orthogonal i.e. ri.rj = 0, i 6= j Please note that while writing the algorithm, if you want to multiply a transpose of a vector with a regular vector, say xTy, you can use in numpy, np.dot(x,y). While implementing the algorithm, set the tolerance to very low tol=1.0e-9. 5 2.2 Hint 2 PROBLEM 2: CONJUGATE GRADIENT METHOD 91 2.7 Iterative Methods Here is the outline of the conjugate gradient algorithm: Choose x0 (any vector will do). Let r0 ← b Ax0. Let s0 ← r0 (lacking a previous search direction, choose the direction of steepest descent). do with k = 0, 1, 2, . . .: αk ← s T k rk sTkAsk . xk+1 ← xk + αksk . rk+1 ← b Axk+1. if |rk+1| ≤ ε exit loop (ε is the error tolerance). βk ← rTk+1Ask sTkAsk . sk+1 ← rk+1 + βksk . It can be shown that the residual vectors r1, r2, r3, . . . produced by the algorithm are mutually orthogonal; that is, ri · r j = 0, i = j . Now suppose that we have carried out enough iterations to have computed thewhole set ofn residual vectors. The resid- ual resulting from the next iteration must be a null vector (rn+1 = 0), indicating that the solution has been obtained. It thus appears that the conjugate gradient algorithm is not an iterative method at all, because it reaches the exact solution after n compu- tational cycles. In practice, however, convergence is usually achieved in less than n iterations. The conjugate gradient method is not competitive with direct methods in the solution of small sets of equations. Its strength lies in the handling of large, sparse systems (where most elements of A are zero). It is important to note that A enters the algorithm only through its multiplication by a vector; that is, in the form Av, where v is a vector (either xk+1 or sk). If A is sparse, it is possible to write an efficient subrou- tine for the multiplication and pass it, rather than A itself, to the conjugate gradient algorithm. ! conjGrad The function conjGrad shown next implements the conjugate gradient algorithm. The maximum allowable number of iterations is set to n (the number of unknowns). Note that conjGrad calls the function Av that returns the product Av. This func- tion must be supplied by the user (see Example 2.18). The user must also supply the starting vector x0 and the constant (right-hand-side) vector b. The function returns the solution vector x and the number of iterations. ## module conjGrad ’’’ x, numIter = conjGrad(Av,x,b,tol=1.0e-9) Conjugate gradient method for solving [A]{x} = {b}. Figure 3: Pseudo algorithm for conjugate gradient method. 2.2 Hint The algorithm requires calculating r, α, and β (if you can work back from the worked example) along with modulus of r, then you can break down this whole algorithm to very simple ma- trix/vector multiplications/addition/subtraction that you have to do in a loop whether it is for or a while loop: a) Multiplying a matrix A with x and subtracting from a vector b to get r. This can all be done using numpy arrays, A = np.array([[ 4, -1 ,1], [ -1, 4 ,-2], [ 1, -2 ,4]]) b = np.array([12, -1, 5]) x = np.ones(3) then r = Ax b is nothing but r = np.dot(A,x) – b b) You also have to solve for α and β. If I take α, then the numerator bit is nothing but multiplying the transpose of s with r (can be done using np.dot). The denominator is transpose of s with As which you can again solve using np.dot provided if you can precalculate A times s using np.dot. numerator = np.dot(s,r) As = np.dot(A,s) denominator = np.dot(s,As) alpha = numerator / denominator You can similarly take hints from above to solve for β c) xk+1 = xk + αksk, i.e., calculating next value of x is nothing but x = x + alpha * s or in short form 6 2.2 Hint 2 PROBLEM 2: CONJUGATE GRADIENT METHOD x+= alpha * s We have also covered the += operator in the lectures. You can similarly solve for s. d) To calculate |r|, you can simply use rnorm = np.dot(r,r) and then use an if statement to see if rnorm is less than a tolerance say 1e-5 to break from the loop. If I have said not to use numpy library function then it means not to use an equivalent function such as np.linalg.solve() to cheat which, we already covered in the lectures and act as a one line solution def conjugate_gradient_cheating(A,b) x=np.linalg.solve(A,b) return x instead I want you to implement your own algorithm: def conjugate_gradient_snazzy(A,b) return x 7