i
n - 1:
| -[K(ui+1 - ui )/h - K(ui - ui-1 )/h] + h C ui = h F(ih) with u0 given, |
UCES
Course: Methods and Analysis in UCES
Applied Area: cooling fin's continuum model
Method: Jacobi and SOR algorithms
Assessment: strict diagonal dominant matrices, existence, convergence
Prerequisites:
Objectives:
General Information: Although iterative methods are not practical (use the tridiagonal algorithm) for tridiagonal problems, they are very useful for more complicated problems that arise from 2D and 3D physical problems.
Contact Person: R. E. White, NCSU, white@math.ncsu.edu
Revision Date: 8-3-94
Copyright ©1994, R. E. White. All rights reserved.
Models of heat flow in more than one direction will generate large and nontridiagonal matrices. Alternatives to the full version of Gaussian elimination, which requires large storage and number of operations, are the iterative methods. These usually require less storage, but the number of iterations needed to approximate the solution can vary with the tolerance in the approximation. In this lecture we present the most elementary iterative methods: Jacobi, Gauss-Seidel and successive over relaxation (SOR). These methods are useful for sparse (many zero components) matrices where the nonzero patterns are very systematic. Another iterative method called preconditioned conjugate gradient (PCG) is particularly useful, but it is beyond the scope of these notes.
We will consider the cooling fin problem from the previous lecture, but here we will use the iterative methods to solve the algebraic system. Also we will study the effects of varying the parameters of the fin such as thickness, T, and width,W.
Let C = (2W + 2T)/(TW) c and
Let ui be an approximation of u(ih) where
i
n - 1:
| -[K(ui+1 - ui )/h - K(ui - ui-1 )/h] + h C ui = h F(ih) with u0 given, |
| -[c(usur - un)/h - K(un - un-1 )/h] + (h/2) C un = (h/2) F(nh). |
The matrix form of this is
Here with n = 4 we have
In order to motivate the definition of these algorithms, consider the
following
Since the diagonal component is the largest, an approximation can be gotten by letting ui-1 and ui+1 be either previous guesses or calculations, and then computing the new ui from the above equation.
Jacobi Method: Let u0 = [0, 0, 0] be the initial guess.
One repeats this until there is little change for all the nodes i.
Gauss-Seidel Method: Let u0 = [0, 0, 0] be the initial guess.
This method will vary from the Jacobi method because we will always use the most recently computed values.
| u2 | = [ (1 + 12/27)/3, (1 + 39/81 + 13/27)/3, (1 + 53/81)/3] |
| = [117/243, 159/243, 134/243]. |
The formula for the next iteration m+1 for node i is
Note the m+1 on the right side.
One repeats this until there is little change for all the nodes i.
The Gauss-Seidel algorithm usually converges much faster than the Jacobi
method. Even though we can define both methods for any matrix, the method may
or may not converge. Even if it does converge, it may do so very slowly and
prove to be of no practical use. Fortunately, for many problems similar to
heat conduction, these methods and their variation do effectively converge.
One variation of the Gauss-Seidel method is the successive over relaxation
(SOR) method, which has an acceleration parameter
. Here the choice
of the parameter should be between 1 and 2 so that convergence is as
rapid as possible. For very special matrices there are formulae for such
optimal , but generally the optimal are approximated by
virtue of experience. Also the initial guess should be as close as possible to
the solution, and in this case one may rely on the nature of the solution as
dictated by the physical problem that is being modeled.
Jacobi Algorithm. Consider Ax = d, m be the Jacobi iteration and i be the node.
for m = 0, maxit for i = 1,nendloop test for convergence endloop.
SOR Algorithm. Consider Ax = d, m be the SOR
iteration, i be the node and
< 2.
for m = 0, maxit
for i = 1,n
endloop
test for convergence
endloop.
There are a number of tests for convergence. One common test is determine if at each node the absolute value of two successive iterates is less the some given small number. This does not characterize convergence! Consider the following sequence of partial sums of the harmonic series which will blowup:
Note that xm+1 - xm = 1/(m+1) goes to zero and yet xm goes to infinity. Thus, the above convergence test could be deceptive.
In most applications of these interative methods the matrix is sparse.
Consequently, one tries to use the zero pattern to reduce the computations in
the summations. It is very important to not do the parts of the summation
where the components of the matrix are zero. Also, it is not usually necessary
to store all the computations. In Jacobi's algorithm you only need two
This section contains Maple code for three types of calculations. First, the cooling fin problem of the previous section is reconsidered with the tridiagonal algorithm replaced by the Jacobi iterative method. The output in graphical form has the space along the rows and the Jacobi iterations along the columns. Second, we use the SOR iterative method and note more rapid convergence than the Jacobi method. Third, we do some numerical experiments with variable thickness, T, of the fin, and note the temperature rises as the thickness increases.
Maple Code for the Cooling Fin.
Input Data:
> with(linalg):
> n:=10:
> T:=.1:
> W:=10.:
> L:=1.:
> csur:=.1:
> usur:= 70.:
> uleft:= 160:
> K:=.001:
Procedure BVPF is Defined with the Jacobi Method:
> BVPF:=proc(n,K,T,W,L,csur,usur,uleft)
> a:=vector(n,0):
> b:=vector(n,0):
> c:=vector(n,0):
> d:=vector(n,0):
> h:=L/n:
> C:=csur*2*(W + T)/(T*W):
> for i from 1 to n-1 do
> a[i]:=2*K + h*h*C:
> c[i]:=-K:
> b[i]:=-K:
> d[i]:=h*h*C*usur:
> od:
> d[1]:= d[1]+K*uleft:
> a[n]:= 2*K + h*h*C + 2*h*csur:
> b[n]:=-2*K:
> d[n]:= h*h*C*usur + 2*csur*usur*h:
> JAC(n,a,b,c,d,uleft):
> end:
Procedure JAC is Defined for the Jacobi Method:
> JAC:=proc(n,a,b,c,d,usur)
> u:=matrix(n,100,usur);
> eps:=.1:
> for m from 1 to 99 do
> numi:= 0:
> u[1,m+1]:=(d[1] - c[1]*u[2,m])/a[1]:
> error := abs(u[1,m+1] - u[1,m]):
> if error < eps then
> numi:=numi +1:
> fi:
> for i from 2 to n-1 do
> u[i,m+1]:=(d[i] - b[i]*u[i-1,m] - c[i]*u[i+1,m])/a[i]:
> error := abs(u[i,m+1] - u[i,m]):
> if error < eps then
> numi:=numi +1:
> fi:
> od:
> u[n,m+1]:=(d[n] - b[n]*u[n-1,m] )/a[n]:
> error := abs(u[n,m+1] - u[n,m]):
> if error < eps then
> numi:=numi +1:
> fi:
> if numi = n then break fi:
> od:
> end:
BVPF is Executed and the Output Graphed:
> with(plots):
> BVPF(10,.001,.1,10.,1.,.001,70.,160):
> uc:=(vector(100,uleft)):
> ut:=transpose(u):
> uu:=transpose(augment(uc,ut));
> matrixplot(uu,style=wireframe,axes=framed);
 | Figure: Jacobi with m = 41 for Convergence |
Procedure SOR is Defined for Successive Over Relaxation Method:
> SOR:=proc(n,a,b,c,d,usur,w)
> u:=matrix(n,100,usur):
> eps:=.1:
> for m from 1 to 99 do
> numi:= 0:
> u[1,m+1]:=(d[1] - c[1]*u[2,m])/a[1]:
> u[1,m+1]:=(1 - w)*u[1,m] + w*u[1,m+1]:
> error := abs(u[1,m+1] - u[1,m]):
> if error < eps then
> numi:=numi +1:
> fi:
> for i from 2 to n-1 do
> u[i,m+1]:=(d[i] - b[i]*u[i-1,m+1] - c[i]*u[i+1,m])/a[i]:
> u[i,m+1]:=(1 - w)*u[i,m] + w*u[i,m+1]:
> error := abs(u[i,m+1] - u[i,m]):
> if error < eps then
> numi:=numi +1:
> fi:
> od:
> u[n,m+1]:=(d[n] - b[n]*u[n-1,m+1] )/a[n]:
> u[n,m+1]:=(1 - w)*u[n,m] + w*u[n,m+1]:
> error := abs(u[n,m+1] - u[n,m]):
> if error < eps then
> numi:=numi +1:
> fi:
> if numi = n then break fi:
> od:
> end:
BVPFS is Executed with SOR (w = 1.0) in Place of JAC:
> BVPFS(10,.001,.1,10.,1.,.001,70.,160,1.00):
> uc:=(vector(100,uleft)):
> ut:=transpose(u):
> uu:=transpose(augment(uc,ut)):
> matrixplot(uu,style=wireframe,axes=framed);
 |
| Figure: SOR ( = 1.0)
with m = 21 for Convergence |
BVPFS is Executed with SOR (w = 1.4) in Place of JAC:
> BVPFS(10,.001,.1,10.,1.,.001,70.,160,1.40):
> uc:=(vector(100,uleft)):
> ut:=transpose(u):
> uu:=transpose(augment(uc,ut)):
> matrixplot(uu,style=wireframe,axes=framed);
 |
| Figure: SOR ( = 1.4)
with m = 11 for Convergence |
BVPFS is Executed for Variable T:
> BVPFS(10,.001,.05,10.,1.,.001,70.,160,1.40):
> L:=[[0,160],seq([i*h,u[i,m]],i=1..10)]:
> BVPFS(10,.001,.1,10.,1.,.001,70.,160,1.40):
> LL:=[[0,160],seq([i*h,u[i,m]],i=1..10)]:
> BVPFS(10,.001,.15,10.,1.,.001,70.,160,1.40):
> LLL:=[[0,160],seq([i*h,u[i,m]],i=1..10)]:
> plot({L,LL,LLL},x=0..1);
 |
| Figure: Temperatures versus Space: T = .05, .10, .15 |
The following code uses Matlab, and it illustrates how the number of
iterations needed for convergence of the SOR algorithm depends on the choice of
the SOR parameter w. This is done for the special problem
(i.e., that for which the number of iterations needed
for convergence is smallest) will depend on the problem and the number of
unknowns.
 |
| Figure: Iterations for Convergence versus
|
In the previous section we noted some short falls of the cooling fin model with diffusion in only one direction. We shall see that the new models for such problems will have more complicated matrices. They will not be tridiagonal, but they will still have large diagonal components relative to the other components in each row of the matrix.
Definition. Let A = [aij ]. A is called strictly diagonally dominant if and only if for all i
The 3×3 example in the beginning of this section and the matrix from the cooling fin are strictly diagonally dominate matrices. The next section will contain some more examples which are not tridiagonal. Also, in the next section we shall prove that algebraic systems which have strictly diagonally dominant matrices do have solutions, and they are unique.
Existence Theorem. Consider Ax = d and assume A is strictly diagonally dominant. If there is a solution, it is unique. Moreover, there is a solution.
Jacobi and Gauss-Seidel (SOR with
= 1.0)
Consider Ax = d. If A is strictly diagonally dominant, then for all x0 both the Jacobi and the Gauss-Seidel algorithms will converge to the solution.
Proof. Let xm+1 be from the Jacobi iteration and Ax = d. The component forms of these statements are
Let the error at node i be
. Record the number of iterations
required for convergence. For a fixed n and and eps find the
such that the number of iterations required
for convergence is a minimum.