[SOLVED] Homework 4. amath 301 beginning scientific computing

Written Exercises
(10 points total)
Exercise 1. (Component Skill 4.1)
Consider the following ridiculously simple linear system:
(
1 0
0 2) (x1
x2
) = (
1
1
) .
a. What is the solution of this linear system? You don’t have to show your work.
b. Set up a Neumann iteration with initial guess x0 = (0, 0)
T . Calculate x1, x2, x3, and x4. Show
your work.
c. Based on your calculations in (b), what will xk equal for k > 0 when (i) k is odd and (ii) k is
even? You don’t have to show your work.
d. Based on your answer to (c), does the Neumann iteration converge?
e. Let M represent the matrix used in your Neumann iteration. Does Mn → 0 as n → ∞? What
does this tell you about the convergence of the Neumann iteration method? Is this consistent
with your answer to (d)?
Exercise 2. (Component Skill 4.2)
For each linear system, decide whether Jacobi iteration is guaranteed or not guaranteed to converge to the solution. Justify each of your decisions in a sentence or two.
a.
⎛
⎜
⎝
3 1 0
1 3 1
0 1 3
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
b.
⎛
⎜
⎝
−3 1 0
1 −3 1
0 1 −3
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
c.
⎛
⎜
⎝
1 3 0
3 1 3
0 3 1
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
d.
⎛
⎜
⎝
1 −3 0
−3 1 −3
0 −3 1
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
e.
⎛
⎜
⎝
2 1 0
1 2 1
0 1 2
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
2
Coding Exercises
(10 points total)
Exercise 1. (Component Skill 4.2)
The following MATLAB function file executes Jacobi iteration:
1 function [ x1 ] = JacobiHW4 (A ,b , x0 ,N )
2 % INPUTS : A is a nxn coefficient matrix
3 % b is a nx1 column vector of knowns
4 % x0 is the nx1 initial guess of the solution
5 % N is the number of iterations
6 % OUTPUTS : x1 is the nx1 solution
7
8 D = diag ( diag ( A ) ) ;
9 L = tril (A , -1) ;
10 U = triu (A ,1) ;
11
12 for j = 1: N
13 x1 = D ( -( L + U ) * x0 + b );
14 x0 = x1 ;
15 end
16
17 end
a. Use the function above to solve
⎛
⎜
⎝
3 1 0
1 3 1
0 1 3
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
(1)
by Jacobi iteration. Be sure to put the code for the function at the very bottom of your hw4
script file. Take the initial guess to be x0 = (0, 0, 0)
T and iterate N = 10 times. Assign the
value of x10 to the variable A1. Use format long to show as many decimal places as possible.
b. What is the exact solution of (1) according to the MATLAB backslash command? Assign
this solution to the variable A2.
c. Let x∗ be the exact solution of (1). Define the error of the kth Jacobi iterate xk as the
maximum difference in absolute value between corresponding components of xk and x∗:
E = max
1≤i≤3
{∣xk,i − x∗,i∣} .
Assign the variable A3 to the value of E for your x10 computed in (a). Use format long to
show as many decimal places as possible.
3
Answers:
Exercise 2. (Component Skill 4.3)
Modify the Jacobi function file given in Coding Exercise 1 to create a new function file in MATLAB
that executes Gauss-Seidel iteration. Call your new function file GaussSeidelHW4. Retain the same
inputs and output as before.
a. Use your new function to solve
⎛
⎜
⎝
3 1 0
1 3 1
0 1 3
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
(2)
by Gauss-Seidel iteration. Be sure to put the code for your function at the very bottom of
your hw4 script file. Take the initial guess to be x0 = (0, 0, 0)
T and iterate N = 10 times.
Assign the value of x10 to the variable A4. Use format long to show as many decimal places
as possible.
b. Let x∗ be the exact solution of (2). Define the error of the kth Gauss-Seidel iterate xk as the
maximum difference in absolute value between corresponding components of xk and x∗:
E = max
1≤i≤3
{∣xk,i − x∗,i∣} .
Assign the variable A5 to the value of E for your x10 computed in (a). How does this error
compare to your error in Coding Exercise 1(c)? Use format long to show as many decimal
places as possible.
4
Answers:
Exercise 3. (Component Skill 4.4)
Modify the Jacobi function file given in Written Exercise 3 to create a new function file in MATLAB
that executes SOR iteration. Call your new function file SORHW4. Retain the same inputs as before,
but add one more input ω for the relaxation parameter. Keep the output the same.
a. Use your new function to solve
⎛
⎜
⎝
3 1 0
1 3 1
0 1 3
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
(3)
by SOR iteration with ω = 1.09. Be sure to put the code for your function at the very bottom
of your hw4 script file. Take the initial guess to be x0 = (0, 0, 0)
T and iterate N = 10 times.
Assign the value of x10 to the variable A6. Use format long to show as many decimal places
as possible.
Remark: For this linear system, ω = 1.09 is very close to the optimal choice.
b. Let x∗ be the exact solution of (3). Define the error of the kth SOR iterate xk as the maximum
difference in absolute value between corresponding components of xk and x∗:
E = max
1≤i≤3
{∣xk,i − x∗,i∣} .
Assign the variable A7 to the value of E for your x10 computed in (a). How does this error
compare to your errors in Coding Exercises 1(c) and 2(b)? Use format long to show as many
decimal places as possible.
5
Answers:
Exercise 4. (Component Skill 4.2)
Rather than specify the number of Jacobi iterations, we could specify a desired error threshold
to achieve instead. In particular, if the kth Jacobi iterate is within this specified threshold of error
of the true solution, we terminate our Jacobi iterations. It’s a great idea, but the problem is that
we usually don’t know the exact solution of our linear system and, hence, can’t exactly measure
the error of our kth iteration.
What we do instead is we keep track of what’s called the Cauchy error. This is the maximum
difference in absolute value between two successive Jacobi iterations:
C = max
i
∣xk,i − xk−1,i∣.
Here, the index i keeps track of the components of our iterates. If the Cauchy error is decreasing,
our iterates must be converging closer and closer to the exact solution of the linear system.
The following represents a MATLAB function file that executes Jacobi iteration until the Cauchy
error falls below a certain threshold:
1 function [ x1 , N ] = JacobiCauchyHW4 (A ,b , x0 , threshold )
2 % INPUTS : A is a nxn coefficient matrix
3 % b is a nx1 column vector of knowns
4 % x0 is the nx1 initial guess of the solution
5 % threshold is the Cauchy error threshold
6 % OUTPUTS : x1 is the nx1 solution
7 % N is the number of iterations
8
9 D = diag ( diag ( A ) ) ;
10 L = tril (A , -1) ;
11 U = triu (A ,1) ;
12
13 x1 = D ( -( L + U ) * x0 + b );
14 err = max ( abs ( x1 – x0 ) ) ;
15 N = 1;
6
16
17 while err > threshold
18 xtemp = x1 ;
19 x1 = D ( -( L + U ) * x1 + b );
20 x0 = xtemp ;
21 err = max ( abs ( x1 – x0 )) ;
22 N = N +1;
23 end
24
25 end
a. Use the function above to solve
⎛
⎜
⎝
3 1 0
1 3 1
0 1 3
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
by Jacobi iteration. Be sure to put the code for the function above at the very bottom of your
hw4 script file. Take the initial guess to be x0 = (0, 0, 0)
T and the error threshold to be 10−10
.
Assign the number of iterations necessary to achieve this error threshold to the variable A8.
b. Modify the function above to solve
⎛
⎜
⎝
3 1 0
1 3 1
0 1 3
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
by Gauss-Seidel iteration. Call your new function GaussSeidelCauchyHW4. Be sure to put the
code for your function at the very bottom of your hw4 script file. Take the initial guess to be
x0 = (0, 0, 0)
T and the error threshold to be 10−10. Assign the number of iterations necessary
to achieve this error threshold to the variable A9. Compare your answer with (a).
c. Modify the function above to solve
⎛
⎜
⎝
3 1 0
1 3 1
0 1 3
⎞
⎟
⎠
⎛
⎜
⎝
x1
x2
x3
⎞
⎟
⎠
=
⎛
⎜
⎝
1
1
1
⎞
⎟
⎠
by SOR iteration. Call your new function SORCauchyHW4. Make sure you add an extra input
for the relaxation parameter ω. Be sure to put the code for your function at the very bottom
of your hw4 script file. Take the initial guess to be x0 = (0, 0, 0)
T , the relaxation parameter to
be ω = 1.09, and the error threshold to be 10−10. Assign the number of iterations necessary to
achieve this error threshold to the variable A10. Compare your answer with (a) and (b).
7
Answers:
8