Download Numerical Analysis Homework 2: Iterative Techniques Solutions - Prof. Erin K. Mcnelis and more Assignments Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!
MATH 640 – Numerical Analysis Homework # 2: Iterative Techniques Solutions
- For each of the iterative methods we’ve discussed: Jacobi, Gauss-Seidel, SOR, Conjugate Gradient Method. Turn in a copy of these m-files. Jacobi Iteration Method
function [final_approx, num_it] = MyJacobi(A,b,x0,tol,max_it);
[num_row,num_col] = size(x0); if num_col > 1 x0 = x0’; end;
D=diag(diag(A)); L = -1tril(A,-1); U = -1triu(A,1);
T = inv(D)(L+U); c = inv(D)b;
converged = 0; i=1;
x(:,1) = x0; old_x = x0;
while (i <= max_it & ~converged) new_x = T*old_x + c; x(:,i+1) = new_x; if norm(new_x - old_x) < tol final_approx = new_x; num_it = i; converged = 1; fprintf(’Jacobis Method Converged! Yippee!\n’); break; else old_x = new_x; i=i+1; end; end;
if (i > max_it & ~converged) fprintf(’Maximum number of iterations were met without converging.’); fprintf(’Final approximation was:\n’); final_approx = x(:,i); num_it = i-1; disp(final_approx); end;
Gauss-Seidel Iteration Method (omitting repeated material of Jacobi)
function [final_approx, num_it] = MyGaussSeidel(A,b,x0,tol,max_it);
D=diag(diag(A)); L = -1tril(A,-1); U = -1triu(A,1);
T = inv(D-L)(U); c = inv(D-L)b;
SOR Iteration Method (omitting repeated material of Jacobi)
function [final_approx, num_it] = MySOR(A,b,x0,omega,tol,max_it);
D=diag(diag(A)); L = -1tril(A,-1); U = -1triu(A,1);
T = inv(D-omegaL)((1-omega)D + omegaU); c = omegainv(D-omegaL)*b;
Conjugate Gradient Iteration Method (omitting repeated material of Jacobi)
function [final_approx,num_it] = MyConjGrad(A,b,x0,tol,max_it);
old_x = x0; old_r = b - Aold_x; old_v = old_r; t = dot(old_r,old_r)/dot(old_v,Aold_v);
while (i <= max_it & ~converged) new_x = old_x + told_v; if norm(new_x - old_x,inf) < tol final_approx = new_x; : break; else new_r = b - Anew_x; s = dot(new_r,new_r)/dot(old_r,old_r); new_v = new_r + sold_v; t = dot(new_r,new_r)/dot(new_v,Anew_v);
old_x = new_x; old_r = new_r; old_v = new_v;
i=i+1; end; end;
- Solve
Amat ∗ u = bvec from Homework #1 again, using the SOR method, but with each of the following sets of velocity values: (a) Case 1: vx = 3, vy = 3
Room Temperature with m = n = 8 Room Temperature with m = n = 16
0 5
10 15
20 25
0 5
10 15 20 2530
35
40
45
50
55
60
65
70
Boundary 4
Case 1: (SOR) Room Temperature with m=n=8 and ω =1.
Boundary 1
Temperature
0 5
10 15
20 25
0 5
10 15 20 2530
35
40
45
50
55
60
65
70
Boundary 4
Case 1: (SOR) Room Temperature with m=n=16 and ω =1.
Boundary 1
Temperature
(b) Case 2: vx =
ln((x − 0)^2 + (y − 7 .5)^2 + e)
ln((x − 25)^2 + e)
, vy =
ln((y − 12)^2 + e)
Room Temperature with m = n = 8 Room Temperature with m = n = 16
0 5
10 15
20 25 0 5 10
15 20 2530
35
40
45
50
55
60
65
70
Boundary 4
Case 2: (SOR) Room Temperature with m=n=8 and ω =1.
Boundary 1
Temperature
0 5
10 15
20 25 0 5 10
15 20 2530
35
40
45
50
55
60
65
70
Boundary 4
Case 2: (SOR) Room Temperature with m=n=16 and ω =1.
Boundary 1
Temperature
(c) Case 3: vx = y − 12. 5 , vy = −(x − 12)
Room Temperature with m = n = 8 Room Temperature with m = n = 16
0 5
10 15
20 25 0 5 10 15 20 2530
35
40
45
50
55
60
65
70
Boundary 4
Case 3: (SOR) Room Temperature with m=n=8 and ω =1.
Boundary 1
Temperature
0 5
10 15
20 25 0 5 10 15 20 2530
35
40
45
50
55
60
65
70
Boundary 4
Case 3: (SOR) Room Temperature with m=n=16 and ω =1.
Boundary 1
Temperature
(d) Case 4: vx = 3 cos(xy), vy = 3 sin(xy)
Room Temperature with m = n = 8 Room Temperature with m = n = 16
0 5
10 15
20 25 0 5 10
15 20 2530
35
40
45
50
55
60
65
70
Boundary 4
Case 4: (SOR) Room Temperature with m=n=8 and ω =1.
Boundary 1
Temperature
0 5
10 15
20 25 0 5 10
15 20 2530
35
40
45
50
55
60
65
70
Boundary 4
Case 4: (SOR) Room Temperature with m=n=16 and ω =1.
Boundary 1
Temperature
Case Number of Iterations m = n = 8 m = n = 16 One 8 10 Two (original) 42 72 Two (new numerator) 9 12 Three) 92 242 Four 211 441
- Let g(x) = xT^ Ax − 2 xT^ b =< x, Ax > − 2 < x, b > We have shown in class that if A is a symmetric positive definite matrix then the vector x that solves Ax = b also minimizes g(x). Suppose that A is not symmetric positive definite, and the conjugate gradient method applied to the linear system with coefficient matrix A converges. Show that it converges to x that solves
(Ax + AT^ x) = 2b
Proof:
Since the conjugate gradient method converges for the n × n matrix A, we know it is the so- lution, x∗, to
min x g(x) = min x (xT^ Ax − 2 xT^ b) = g(x∗)
thus, the gradient of g(x) is equal to the zero vector at x∗. By the product rule we see that this means
dg dxi
= eTi Ax + xT^ Aei − 2 eTi b = 0 for i = 1, 2 , · · · , n.
Since xT^ Aei is a scalar, then
xT^ Aei = (xT^ Aei)T^ = eTi AT^ x.
Thus, for all i
eTi Ax + eTi AT^ x − 2 eTi b = 0 eTi (Ax + AT^ x) − 2 eTi b = 0 eTi ((Ax + AT^ x) − 2 b) = 0
i.e. the ith^ entry of (Ax + AT^ x) − 2 b is zero for all i = 1, 2 , · · · , n.
This implies that
(Ax + AT^ x) − 2 b = 0
and hence,
(Ax + AT^ x) = 2 b
Thus, the solution x to the conjugate gradient method is also the solution to the linear equation (Ax + AT^ x) = b.
