AMATH242/CS371 Spring 2024 Assignment III

Release date: Friday, June 28th

Due date: Monday, July 15th, 11:59 pm

Due date with your one time extension: Thursday, July 18th, 11:59 pm

•  Questions below are either theoretical or computational. For the computational questions you may use any pro- gramming language you prefer except symbolic ones like Maple or Mathematica.

• This assignment should be submitted to Crowdmark. Besides your written/typed solutions you must also submit your code. Please make sure your upload has the correct orientation and ordering.

•  If you want to use your one time 3-day extension, please email the instructor before the original due date.

• You are not allowed to post this assignment on sites like stackexchange.com, chegg.com, etc. This will be checked.

Offenders will face penalties from the Math Faculty (e.g. suspension).

Total points: 20

1. (0points)Please sign the Academic Integrity Checklist on the lastpage ofthis pdf. Ifyou do not sign the Academic Integrity Checklist you will receive a 0 for this assignment.

2. (Theoretical, 3 points) When applying Gauss-Seidel iteration to solve a linear system Ax→ = b(→) and when A is one

of the two matrices below, we expect convergence of the iteration. You’re asked to argue why the convergence. Use two different arguments for convergence of Gauss-Seidel iteration applied to A1  and A2 . (You’re allowed to use a software/calculator to compute, for example, an inverse of a matrix, if needed.)

3. (Computational, 7 points) We will be implementing Jacobi and Gauss-Seidel methods in this question, along with exploring the vectorization*technique to accelerate both iterations.

(1)  First, implement Jacobi, Gauss-Seidel methods based on their pseudocodes Algorithms 3.5, 3.6 in [YW] Lec

12. This also means you’re required to NOT use the matrix-based iteration. For example, for Jacobi’smethod, you should NOT create matrices D, L, U from A and write an one-line code for the following:

(new)  = −D−1 (L + (old)  + D−1 b(→) .

Instead, you should update components of the vector with loops:

and this is exactly what Algorithm 3.5 is doing.

The following requirements apply to your implemention throughout this question:

• Use vector 2-norm to evaluate the residuals;

•  Use relative tolerance rel  =  10−3 , i.e., the stopping criterion should be: ‖r(⃗)(k)‖ / ‖r(⃗)(0) ‖ ≤ 10−3 ;

• Use the zero vector (a vector with all its components being zeros) as the initial guess;

•  Output the final iteration count when stopping criterion is met.

Apply your Jacobi and Gauss-Seidel codes to the two linear systems A ⃗x = b(⃗) and Abigger ⃗(x)  = b(⃗)bigger , with their data provided on LEARN (A .txt,b.txt,A bigger.txt,b bigger.txt‡ ). Here is how you load matrix/vector data into MATLAB:

A = load ( ' A .txt ' ) ; b = load ( ' b .txt ' ) ;

And here is how you load them into Python:

import num py a s np A = np . loadtxt ( ' A .txt ' ) b = np . loadtxt ( ' b . txt ' , nd min =2)

You are also asked to record the execution time of your code when solving each system. Here is how you do it if, for example, your Jacobi iteration is written as a function jacobi itr(A,b,guess,tol) in MATLAB§ :

% The execution of codes between tic and to c is timed by MATLAB %  and   the elapsed time is printed out ( see the documentation for examples ) tic number _ it r  =  jacobi _ it r (A ,b , guess , to l ) ; % the final   iteration  count   is  the output toc

And here is how you do it in PythonT :

import time # ...... # other codes here # ...... start = time . time () nun mber_it r  = jacobi_it r (A ,b , guess , to l ) T = time . time () -  start print ( ' Elapsed time is ' ,  T , ' seconds . ' )

Finally, tabulate your results as below and make a brief comment comparing your results between Jacobi and Gauss-Seidel: which iteration is faster? And how much faster?

	A ⃗x = b(⃗)		Abigger ⃗(x) = b(⃗)bigger
	Elapsed time (s)	Iteration count	Elapsed time (s)	Iteration count
Jacobi method
Gauss-Seidel method

(2)  Show that we can rewrite the Jacobi iteration formula as follows:

xi(new)  =  ( bi  − j=Σ1 =i aij xj(old)) ⟹ xi(new)  = xi(old)  +  ( bi  − Ai ⃗(x)(old)) ,

where Ai denotes the ith-row of A. (This is all you need to do for this subquestion. The remaining information is for you to read and is useful for the next subquestion.)

For the Gauss-Seidel iteration, we have (see [YW] Lec 12, above Algorithm 3.6)

xi(new)  =   ( bi − aij xj(new) − aij xj(old)) .

Since Gauss-Seidel successively updates the vector component and hence works on only one vector, the formula above becomes (by dropping superscripts “(new)” and “(old)”)

This formula is precisely what is being implemented in Algorithm 3.6.  Now, with the identical steps to rewrite Jacobi’s formula, we can show:

(3) In the previous subquestion, we obtained a vectorized form of the j-loop in Algorithms 3.5 and 3.6.  In MATLAB (and NumPy), converting loops (that operate on scalars) to their vectorized form. (that operate on matrices/vectors) often lead to much faster codes. This is thanks to that much effort in developing MATLAB (and NumPy) focuses on optimizing operations involving matrices and vectors. With j-loop vectorized, the pseudocode in Algorithm 3.5 becomes:

Initial guess: ⃗(x)(0) k = 0; r(⃗)(0) = b(⃗) − A ⃗x(0) while ‖r(⃗)(k)‖ / ‖r(⃗)(0) ‖ > rel do

Similarly, the pseudocode in Algorithm 3.6 becomes:

Initial guess: ⃗(x) r(⃗)(0) = b(⃗) − A ⃗x r(⃗) = r(⃗)(0) while ‖r(⃗)‖ / ‖r(⃗)(0) ‖ > rel do for i = 1 ∶ n do xi = xi + (bi − Ai ⃗(x))/aii end r(⃗) = b(⃗) − A ⃗x end

Now, modify your codes from the first subquestion and run them for the same tests. Tabulate your results as below:

	A ⃗x = b(⃗)		Abigger ⃗(x) = b(⃗)bigger
	Elapsed time (s)	Iteration count	Elapsed time (s)	Iteration count
Jacobi method
Gauss-Seidel method

The iteration counts should be the same as in the previous table (we didn’t change the iterations; we only changed the way they are calculated). Make a brief observation on how much of an acceleration you gained by vectorizing the j-loop.

4.  (Theoretical, 4 points) In this question, we interpolate and approximate the natural logarithm function ln(x). When evaluating ln(x), you may use a calculator. When presenting your final results, keep three significant digits.

(1) Find the quadratic interpolating polynomial of ln(x) in the Lagrange form.  The data points are (1, ln(1)), (2, ln(2)), (3, ln(3)).  Compute an approximation of ln(2.9)  by evaluating the interpolating polynomial at x = 2.9 and compute the relative error of the approximation.

(2)  Find the Hermite interpolating polynomial of the function ln(x).  The x-value of the three data points are x   =   1, 2, 3  (same as the previous subquestion).   Compute an approximation of ln(2.9)  by evaluating the Hermite interpolating polynomial at x  =  2.9 and compute the relative error of the approximation.

5.  (Computational, 6 points) In this question, you’re asked to interpolate the Runge function f(x) = on

[−1, 1] with equidistant nodes and Chebyshev nodes, respectively:

• The n + 1 equidistant nodes on [−1, 1]: xk   = − 1, with k  = 0, 1, 2, … , n.

• The n + 1 Chebyshev nodes on (−1, 1): xk ), with k  = 0, 1, 2, … , n.

(1)  Given n + 1 nodes {xk }k(n)=0  (be it equidistant or Chebyshev), write down the Lagrange form of the ℙn   inter-

polating polynomial for the Runge function. You need to explicitly write out all Lagrange polynomials, i.e., lk (x)’s by their definitions.

(2)  Implement the ℙn  interpolating polynomial in subquestion 1 with the following requirements:

•  Your program should take as an input a set of n +  1 nodes;

•  For a given set of n + 1 nodes, your program should evaluate the resulting interpolating polynomial at

400 equidistant nodesⅡ in [−1, 1];

•  Your program should also evaluate the Runge function at these 400 nodes;

•  Your program should compute the errors between the interpolating polynomial and the Runge function at these 400 nodes and find the maximum error (in terms of absolute value);

•  Your program should generate a figure with two plots: one for the Runge function, the other the inter- polating polynomial. You may use the data you obtained on the 400 nodes for plotting as well. Here is an example of how the figure should look:

Feed your program with equidistant and Chebyshev nodes for n = 4, 8, 12. Submit the six plots and tabulate the maximum errors as below (with three significant digits):

	equidistant	Chebyshev
n = 4
n = 8
n = 12

Make the following observations on the tabulated results:

•  Between the two choices of nodes, which one is better when it comes to the maximum error of the interpolating polynomial?

• As n grows, how does the maximum error behave for each set of nodes? Briefly explain why we would expect the behavior. for each case.

ⅡThis vector of nodes can be created using, for example, linspace(-1,1,400)