Project of Numerical Analysis
February 18, 2014
Consider the Dirichlet boundary-value problem:
−∆u ≡ −uxx− uyy = 2π2sin πx sin πy, for (x, y) ∈ Ω, (1) u(x, y) = 0 (x, y) ∈ ∂Ω,
for Ω := {x, y|0 < x, y < 1} ⊆ R2 with boundary ∂Ω, which has the exact solution
u(x, y) = sin πx sin πy, and is shown in Figure 1.
Figure 1: Exact solution.
1 Center difference discretization
To solve (1) by means of a difference methods, one replaces the differential operator by a difference operator. Let
Ωh:= {(xi, yi)|i, j = 1, . . . , n},
∂Ωh:= {(xi, 0), (xi, 1), (0, yj), (1, yj)|i, j = 0, 1, . . . , n + 1},
where xi = ih, yj = jh, i, j = 0, 1, . . . , n + 1, h := n+11 , n ≥ 1, is an integer.
From the Taylor’s theorem, we have u(xi+ h) = u(xi) + u0(xi)h +h2
2 u00(xi) +h3
6 u000(xi) +h4
24u(4)(ξ1) u(xi− h) = u(xi) − u0(xi)h +h2
2 u00(xi) −h3
6 u000(xi) +h4
24u(4)(ξ2), where ξ1 is between xi and xi+ h and ξ2 is between xi and xi− h. Hence
u00(xi) =u(xi+ h) − 2u(xi) + u(xi− h)
h2 −h2
12u(4)(ξ)
=u(xi+1) − 2u(xi) + u(xi−1)
h2 −h2
12u(4)(ξ), where ξ is between xi− h and xi+ h. Similarly,
∂2u
∂x2(xi, yj) = u(xi+1, yj) − 2u(xi, yj) + u(xi−1, yj)
h2 −h2
12
∂4u
∂x4(ξi, yj),
∂2u
∂y2(xi, yj) = u(xi, yj+1) − 2u(xi, yj) + u(xi, yj−1)
h2 −h2
12
∂4u
∂x4(xi, ηj), where ξi ∈ (xi−1, xi+1) and ηj∈ (yj−1, yj+1). It implies that
∂2u
∂x2(xi, yj) +∂2u
∂y2(xi, yj)
=u(xi, yj−1) + u(xi−1, yj) − 4u(xi, yj) + u(xi+1, yj) + u(xi, yj+1) h2
−h2 12
∂4u
∂x4(ξi, yj) +∂4u
∂x4(xi, ηj)
.
Let uij denote an approximated value of function u at the grid point (xi, yj) for i, j = 1, . . . , n + 1. Then
− uxx(xi, yj) − uyy(xi, yj) ≈ −ui,j−1− ui−1,j+ 4ui,j− ui+1,j− ui,j+1
h2 with an error O(h2) and the equation
−uxx(xi, yj) − uyy(xi, yj) = 2π2sin πxisin πyj≡ fij
can be replaced by the following equation
−ui,j−1− ui−1,j+ 4ui,j− ui+1,j− ui,j+1
h2 = fij (2)
for i, j = 1, . . . , n.
For j = 1, we have
−u1,0− u0,1+ 4u1,1− u2,1− u1,2=h2f1,1, (3a)
−u2,0− u1,1+ 4u2,1− u3,1− u2,2=h2f2,1, (3b) ...
−un−1,0− un−2,1+ 4un−1,1− un,1− un−1,2=h2fn−1,1, (3c)
−un,0− un−1,1+ 4un,1− un+1,1− un,2=h2fn,1. (3d) By the boundary condition, it holds that
u1,0=u2,0 = · · · = un,0= 0, (4a)
u0,1=un+1,1= 0. (4b)
Substituting (4) into (3), we get
4u1,1− u2,1−u1,2 = h2f1,1, (5a)
−u1,1+ 4u2,1− u3,1−u2,2 = h2f2,1, (5b) ...
−un−2,1+ 4un−1,1− un,1−un−1,2 = h2fn−1,1, (5c)
−un−1,1+ 4un,1−un,2 = h2fn,1. (5d) Let, for j = 1, . . . , n,
u:,j=
u1,j
u2,j
... un,j
, f:,j=
f1,j
f2,j
... fn,j
, A1=
4 −1
−1 . .. . .. . .. . .. −1
−1 4
∈ Rn×n.
Then (5) can be rewritten as following matrix form:
A1 −In
u:,1 u:,2
= h2f:,1. For j = 2, . . . , n − 1, using u0,j = un+1,j= 0, we have
−u1,j−1+ 4u1,j− u2,j−u1,j+1 = h2f1,j,
−u2,j−1− u1,j+ 4u2,j− u3,j−u2,j+1 = h2f2,j, ...
−un−1,j−1− un−2,j+ 4un−1,j− un,j−un−1,j+1 = h2fn−1,j,
−un,j−1− un−1,j+ 4un,j−un,j+1 = h2fn,j.
Above equations can be represented as following matrix form:
−In A1 −In
u:,j−1
u:,j
u:,j+1
= h2f:,j.
For j = n, using u1,n+1= u2,n+1= un,n+1= 0, we have
−u1,n−1+ 4u1,n− u2,n = h2f1,n,
−u2,n−1− u1,n+ 4u2,n− u3,n = h2f2,n, ...
−un−1,n−1− un−2,n+ 4un−1,n− un,n = h2fn−1,n,
−un,n−1− un−1,n+ 4un,n = h2fn,n. Above equations can be represented as following matrix form:
−In A1
u:,n−1
u:,n
= h2f:,n.
Therefore, (2) with boundary conditions is equivalent to a linear system
Au = h2f (6)
with
A =
A1 −In
−In A1 . .. . .. . .. −In
−In A1
∈ Rn2×n2, (7)
and
A1=
4 −1
−1 . .. . .. . .. . .. −1
−1 4
, u =
u:,1
u:,2
... u:,n
, f =
f:,1
f:,2
... f:,n
.
2 Project for direct method
(a) Use Algorithms 1, 2 and 3 (Gaussian elimination) to reduce A in (7) to an upper triangular matrix and modify the entries of b accordingly. Compare and plot the CPU times for reducing A to upper triangular with various n by using these three algorithms. (Use “tic” and “toc” functions in MATLAB to estimate the CPU times.)
Require: Nonsingular matrix A and right hand side vector b.
Ensure: This algorithm implements the Gaussian elimination procedure to re- duceAtoupper triangularand modify the entries ofbaccordingly.
1: for k = 1, . . . , n − 1 do
2: Let p be the smallest integer with k ≤ p ≤ n and apk6= 0.
3: If @ p, then stop.
4: If p 6= k, then perform (Ep) ↔ (Ek).
5: for i = k + 1, . . . , n do
6: Compute t = A(i, k)/A(k, k);
7: Set A(i, k) = 0;
8: Update b(i) = b(i) − t × b(k);
9: for j = k + 1, . . . , n do
10: Update A(i, j) = A(i, j) − t × A(k, j);
11: end for
12: end for
13: end for
Algorithm 1: Gaussian elimination
Require: Nonsingular matrix A and right hand side vector b.
Ensure: This algorithm implements the Gaussian elimination procedure to re- duce A to upper triangular and modify the entries of b accordingly.
1: for k = 1, . . . , n − 1 do
2: Let p be the smallest integer with k ≤ p ≤ n and apk6= 0.
3: If @ p, then stop.
4: If p 6= k, then perform (Ep) ↔ (Ek).
5: for i = k + 1, . . . , n do
6: Compute t = A(i, k)/A(k, k);
7: Set A(i, k) = 0;
8: Update b(i) = b(i) − t × b(k);
9: UpdateA(i, k + 1 : n) = A(i, k + 1 : n) − t × A(k, k + 1 : n);
10: end for
11: end for
Algorithm 2: Vector version of Gaussian elimination
Require: Nonsingular matrix A and right hand side vector b.
Ensure: This algorithm implements the Gaussian elimination procedure to re- duce A to upper triangular and modify the entries of b accordingly.
1: for k = 1, . . . , n − 1 do
2: Let p be the smallest integer with k ≤ p ≤ n and apk6= 0.
3: If @ p, then stop.
4: If p 6= k, then perform (Ep) ↔ (Ek).
5: Computet = A(k + 1 : n, k)/A(k, k);
6: SetA(k + 1 : n, k) = 0;
7: UpdateA(k + 1 : n, k + 1 : n) = A(k + 1 : n, k + 1 : n) − t × A(k, k + 1 : n);
8: Updateb(k + 1 : n) = b(k + 1 : n) − b(k) × t.
9: end for
Algorithm 3: Matrix version of Gaussian elimination
(b) Use backward substitution to solve the upper triangular linear system in (a). Plot the CPU times for solving such linear system with various n.
(c) Compare the CPU times for using left matrix divide “A \ b” in MATLAB with that in (a) and (b).
(d) Store the matrix A with sparse format. Plot the CPU times for generating matrix A and solving the associated linear systems by left matrix divide
“A \ b” with various n.
3 Project for iterative method
(e) Use Jacobi method to solve linear system (6).
Given an initial vector x(0), rewrite the linear system as:
a11x(k)1 + a12x(k−1)2 + a13x(k−1)3 + · · · + a1nx(k−1)n = b1
a21x(k−1)1 +a22x(k)2 + a23x(k−1)3 + · · · + a2nx(k−1)n = b2
... an1x(k−1)1 + an2x(k−1)2 + an3x(k−1)3 + · · · +annx(k)n = bn. If we decompose the coefficient matrix A as
A = L + D + U,
where D is the diagonal part, L is the strictly lower triangular part, and U is the strictly upper triangular part, of A, then we derive the iterative formulation for Jacobi method:
x(k)= −D−1(L + U )x(k−1)+ D−1b.
• Use Algorithm 4 with initial vector x(0) = [1, · · · , 1]> to solve linear system (6). Plot the CPU times and iteration numbers k for solving such linear system with various n.
Require: Given x(0), tolerance T OL, maximum number of iteration M . Ensure: The solution x.
1: Set k = 1.
2: Compute x = −D−1(L + U )x(0)+ D−1b.
3: while k ≤ M and kx − x(0)k2≥ T OL do
4: Set k = k + 1, x(0)= x;
5: Compute x = −D−1(L + U )x(0)+ D−1b;
6: end while
Algorithm 4: Jacobi method
(f) Use Gauss-Seidel method to solve linear system (6).
Given an initial vector x(0), rewrite the linear system as:
a11x(k)1 + a12x(k−1)2 + a13x(k−1)3 + · · · + a1nx(k−1)n = b1
a21x(k)1 +a22x(k)2 + a23x(k−1)3 + · · · + a2nx(k−1)n = b2
a31x(k)1 +a32x(k)2 +a33x(k)3 + · · · + a3nx(k−1)n = b3
... an1x(k)1 +an2x(k)2 +an3x(k)3 + · · · +annx(k)n = bn.
This improvement induce the Gauss-Seidel method. The iteration of the Gauss-Seidel method is defined as follows:
x(k)= −(D + L)−1U x(k−1)+ (D + L)−1b.
Require: Given x(0), tolerance T OL, maximum number of iteration M . Ensure: The solution x.
1: Set k = 1.
2: Compute x = −(D + L)−1U x(0)+ (D + L)−1b.
3: while k ≤ M and kx − x(0)k2≥ T OL do
4: Set k = k + 1, x(0)= x;
5: Compute x = −(D + L)−1U x(0)+ (D + L)−1b;
6: end while
Algorithm 5: Gauss-Seidel method
1. Use MATLAB functions “triu(A,1)” and “tril(A,-1)” to extract the strictly upper and lower triangular parts of A, respectively.
2. Use Algorithm 5 with initial vector x(0)= [1, · · · , 1]> to solve linear system (6). Plot the CPU times and iteration numbers k for solving such linear system with various n.
3. Compare the results produced by Jacobi and Gauss-Seidel methods.
(g) Use SSOR method to solve linear system (6).
Given an initial vector x(0), rewrite the linear system as:
a11x(k)1 + a12x(k−1)2 + a13x(k−1)3 + · · · + a1nx(k−1)n = b1 a21x(k)1 +a22x(k)2 + a23x(k−1)3 + · · · + a2nx(k−1)n = b2 a31x(k)1 +a32x(k)2 +a33x(k)3 + · · · + a3nx(k−1)n = b3
... an1x(k)1 +an2x(k)2 +an3x(k)3 + · · · +annx(k)n = bn.
Let the approximate solution x(k,i) produced by Gauss-Seidel method be defined by
x(k,i)=h
x(k)1 , . . . , x(k)i−1, x(k−1)i , . . . , x(k−1)n iT and
r(k)i =h
r(k)1i , r2i(k), . . . , rni(k)iT
= b − Ax(k,i)
be the corresponding residual vector. Then the ith component of r(k)i is
rii(k)= bi−
i−1
X
j=1
aijx(k)j −
n
X
j=i+1
aijx(k−1)j − aiix(k−1)i ,
so
aiix(k−1)i + rii(k)= bi−
i−1
X
j=1
aijx(k)j −
n
X
j=i+1
aijx(k−1)j = aiix(k)i .
Consequently, the Gauss-Seidel method can be characterized as choosing x(k)i to satisfy
x(k)i = x(k−1)i +r(k)ii aii
.
Relaxation method is modified the Gauss-Seidel procedure to
x(k)i = x(k−1)i + ωr(k)ii aii
= x(k−1)i + ω aii
bi−
i−1
X
j=1
aijx(k)j −
n
X
j=i+1
aijx(k−1)j − aiix(k−1)i
= (1 − ω)x(k−1)i + ω aii
bi−
i−1
X
j=1
aijx(k)j −
n
X
j=i+1
aijx(k−1)j
(8)
for certain choices of positive ω. These methods are called for ω < 1: under relaxation,
ω = 1: Gauss-Seidel method, ω > 1: over relaxation.
Over-relaxation methods are called SOR (Successive over-relaxation). To determine the matrix of the SOR method, we rewrite (8) as
aiix(k)i + ω
i−1
X
j=1
aijx(k)j = (1 − ω)aiix(k−1)i − ω
n
X
j=i+1
aijx(k−1)j + ωbi,
so that if A = L + D + U , then we have
(D + ωL)x(k)= [(1 − ω)D − ωU ] x(k−1)+ ωb.
Theorem 1 (Ostrowski-Reich) If A is positive definite and the relax- ation parameter ω satisfying 0 < ω < 2, then the SOR iteration converges for any initial vector x(0).
Let A be symmetric and A = D + L + LT. The idea is in fact to imple- ment the SOR formulation twice, one forward and one backward, at each iteration. That is, SSOR method defines
(D + ωL)x(k−12) = (1 − ω)D − ωLT x(k−1)+ ωb, (9) (D + ωLT)x(k) = [(1 − ω)D − ωL] x(k−12)+ ωb. (10) Define
Mω: = D + ωL,
Nω: = (1 − ω)D − ωLT. Then from the iterations (9) and (10), it follows that
x(k)= Mω−TNωTMω−1Nω x(k−1)+ ω Mω−TNωTMω−1+ Mω−T b
≡ T (ω)x(k−1)+ M (ω)−1b, where
M (ω) = 1
ω(2 − ω)(D + ωL)D−1 D + ωLT . 1. Take x(0)= [1, · · · , 1]> as an initial vector.
2. Use MATLAB functions “triu(A,1)” and “tril(A,-1)” to extract the strictly upper and lower triangular parts of A, respectively.
3. Fixed n = 100 and uniformly took 40 values for the parameter ω in the interval (0, 2), show the iteration numbers and CPU times of SSOR iterative method for each ω. Find the optimal value ω∗ of the parameter ω.
4. Compare the iteration numbers and CPU times for Jacobi, Gauss- Seidel and SSOR(ω∗) iterative methods with various n.
(h) Use conjugate gradients method to solve linear system (6).
1. Use MATLAB function pcg without any preconditioner:
[x, flag, relres, iter] = pcg(A, b, tol, maxit) 2. Use MATLAB function pcg with a given preconditioner:
[x, flag, relres, iter] = pcg(A, b, tol, maxit, M), [x, flag, relres, iter] = pcg(A, b, tol, maxit, M1, M2), [x, flag, relres, iter] = pcg(A, b, tol, maxit, [], M2), [x, flag, relres, iter] = pcg(A, b, tol, maxit, MFUN).
(i) Jacobi method: A = D + (L + U ), M = D xk+1= −D−1(L + U )xk+ D−1b (ii) Gauss-Seidel: A = (D + L) + U , M = D + L
xk+1= −(D + L)−1U xk+ (D + L)−1b.
(iii) SSOR: A = D + L + LT, M = M (ω)
x(k)= Mω−TNωTMω−1Nω x(k−1)+ M (ω)−1b, where
M (ω) = 1
ω(2 − ω)(D + ωL)D−1 D + ωLT . (iv) M may be a function handle MFUN returning M−1x
[x, flag, relres, iter] = pcg(A, b, tol, maxit, ...
@(x)precSSOR(x,omega,mtxLower,mtxdiag) 3. Compare the iteration numbers and CPU times for pcg by using
different preconditioner with various n.
