3.4.3.   Model.

The above partial differential equation is usually associated with boundary conditions which may or may not have derivatives of the solution. For the present we will assume the temperature is given as a constant on the boundary of (0,L)×(0,W), L = W = 1, K = 1, f(x,y) = 0 and T = 2. So, we are considering

-uxx - uyy + cu = cusur .

(2)

Let u_ij be the approximation of u(ih, jh) where h = dx = dy = 1.0/n. Approximate the second order derivatives by centered finite differences.

Finite Difference Model of (2)

-[(u_i+1,j - u_ij )/h - (u_ij - u_i-1,j )/h]/h - [(u_i,j+1 - u_ij )/h - (u_ij - u_i,j-1 )/h]/h + cu_ij = cu_sur

There are (n - 1)² equations for (n - 1)² unknowns u_ij . This can be put into matrix form by multiplying the equations by h² and listing the unknowns starting with the smallest y values and moving from the smallest to the the largest x values. The first grid row of unknowns is denoted by U₁ = (u₁₁ , u₂₁ , u₃₁ , u₄₁ )^T for n = 5.

Figure: Grid for 2D Diffusion

Hence, the block form of the above system is

The above block tridiagonal system is a special case of the following block tridiagonal system where all entries are N×N matrices (N = 4). The block system has N² blocks, and so there are N² unknowns. If the full version of the Gaussian elimination algorithm were used, it would require (N² )³ /3 = N⁶/3 operations.

.

3.4.4.   Method.

We present three alternatives to the Gaussian elimination algorithm. The first method is a block version of the tridiagonal algorithm. It requires only order 3N⁴ operations and less storage than the full version of the Gaussian elimination algorithm. The other two algorithms are variations on the SOR method and can require, for proper choice of SOR parameter, fewer operations and less storage.

In the block tridiagonal algorithm all entries are matrices. The "divisions" for the point form must be replaced by matrix solves, and one must be careful to preserve the proper order of matrix multiplication. The derivation of the following is similar to the derivation of the point form.

Block Tridiagonal Algorithm.

    (1) = A(1)
    solve (1)*(1) = C(1)
    solve (1)*Y(1) = D(1)
    for i = 2, n
        (i) = A(i)- B(i)*(i-1)
        solve (i)*(i) = C(i)
        solve (i)*Y(i) = D(i) - B(i)*Y(i-1)
    endloop
    X(n) = Y(n)
    for i = n - 1,1
        X(i) = Y(i) - (i)*X(i+1)
    endloop.

The block or line version of the SOR algorithm also requires a matrix solve step in place of a "division." Note the matrix solve has a point tridiagonal matrix for the problem in (2).

Block SOR Algorithm.

    for m = 1,maxit
        for i = 1,n
            solve A(i) Xtemp = D(i) - B(i)*X(i-1) - C(i)*X(i+1)
            X(i) = (1-w)*Xtemp + w*X(i)
        endloop
        test for convergence
    endloop.

The point version of SOR for the problem in (2) can be very effectively implemented because we know exactly where and what the nonzero components are in each row of the coefficient matrix. In the following note how the test for convergence is done and how the storage is kept to a minimum.

Point SOR Algorithm for -(u_xx + u_yy ) + cu = f(x,y) on = (0,1)×(0,1) and u given on Boundary of .

    define w, n, eps
    h = 1/n
    tol = eps*h*h
    for m = 1,maxit
        numi = 0
        for j = 1,n-1
            for i = 1,n-1
                utemp = (h*h*f(i*h,i*h) + u(i,j-1) + u(i-1,j) +
                               u(i+1,j) + u(i,j+1))/(4 + h*h*c)
                utemp = (1 - w)*u(i,j) + w*utemp
                error = abs(utemp - u(i,j))
                u(i,j) = utemp
                if error < tol then
                   numi = numi + 1
                endif
             endloop
        endloop
        if numi = (n-1)*(n-1) go to end
    endloop
    end.

3.4.5.   Implementation.

Maple Code for a Cooled Plate and Point SOR.

   Input Data:
        > with(linalg):
	> with(plots):
        > c:=.4:
        > usur:= 70:
        > eps:=.1:
        > n:= 8:
        > w:=1.4:
        > f:=(x,y)->c*usur:

   Procedure BVP2D and SOR is Defined:
        > BVP2D:=proc(f,c,n,w,eps):
        > u:=matrix(n+1,n+1,200):
        > h:=1/n:
        > maxit:=400:
        > tol:=eps*h*h:
        > for m from 1 to maxit do
        >     numi:=0:
        >     for j from 2 to n do
        >         for i from 2 to n do
        >             utemp:= evalf((h*h*f(i*h,j*h) + u[i,j-1] + u[i-1,j]
                                         + u[i+1,j] + u[i,j+1])/(4 + h*h*c)):
        >             utemp:= (1-w)*u[i,j] + w*utemp:
        >             error:= abs(utemp - u[i,j]):
        >             u[i,j]:= utemp:
        >             if error < tol then 
                            numi:=numi + 1
                      fi:
        >         od:
        >     od:
        >     if numi = (n-1)*(n-1) then break fi:
        > od:
        > end:

   Output Data:
        > BVP2D(f,1.,8,1.4,.1):
        > matrixplot(u,style=wireframe,axes=boxed);

Figure: Temperature of Cooled Plate, c = 1.0 and x and y Trace Curves

        > numi;
        49
        > m;
        19
        > BVP2D(f,10.,8,1.4,.1):
        > matrixplot(u,style=wireframe,axes=boxed);

Figure: Temperature of Cooled Plate, c = 10.0 and Contour Curves

The following Matlab code for the numerical solution of

-u_xx -u_yy = sin(x) sin(y) on (0,1)×(0,1) with u = 0 on the boundary

illustrates the SOR algorithm and how it depends on the choice of the SOR parameter, .

Matlab Code and SOR.

   function [w,soriter] = sor2d(n,w)
   h = 1./n;
   tol =.1*h*h;
   maxit = 500;
   u = zeros(n+1);
   f = zeros(n-1);
   for i=2:n
       for j=2:n
             f(i,j) = h*h*pi*pi*sin(pi*(i-1)*h)*sin(pi*(j-1)*h);
         end
     end
     for iter =1:maxit
         numi = 0;
         for j = 2:n
             for i = 2:n
                 utemp = (f(i,j) + u(i,j-1) + u(i-1,j) + u(i+1,j) + u(i,j+1))*.25;
                 utemp = (1. - w)*u(i,j) + w*utemp;
                 error = abs(utemp- u(i,j));
                 u(i,j) = utemp;
                 if error<tol
                    numi = numi + 1;
                 end
             end
         end
         if numi==(n-1)*(n-1) 
            break;
         end
     end
     w;
     soriter = iter;

Output in Graphical Form:

     >>for i = 1:20
           [x(i),y(i)] = sor(20,1. + (i-1)*.05);
          end
     >>plot(x,y)

Figure: Iterations for Convergence versus

3.4.6.   Assessment.

In the above first 2D heat diffusion model we kept the boundary conditions simple. However, in the 2D model of the cooling fin one should consider the heat that passes through the edge portion of the fin. This is similar to what was done for the cooling fin with diffusion in just the x direction. There the heat flux at the tip was given by the boundary condition Ku_x (L) = c(u_sur - u(L)). For the 2D steady state cooling fin model we have similar boundary conditions on the edge of the fin that is away from the mass to be cooled.

0 = f(x,y) + (2c/T)(u_sur - u) + (Ku_x(x,y))_x + (Ku_y (x,y))_y ,

    Ku_x (L,y) = c(u_sur - u(L,y)),

    -Ku_y (x,0) = c(u_sur - u(x,0)),

    Ku_y(x,W) = c(u_sur - u(x,W)) and

    u(0,y) = given temperature of the mass to be cooled.

Figure: Cooling Fin with Boundary Conditions

The finite difference model must have additional equations for the cells near the edge of the fin. So, in the above figure there are 10 equations for the (dx dy) cells, 5 equations for the (dx/2 dy) cells, 4 equations for the (dx dy/2) cells and 2 equations for the corner (dx/2 dy/2) cells. For example, the cells at the right most side with (dx/2 dy) the finited difference equations are for i = n and 1 < j < n

0 = f(L,j dy) (dx/2 dy T) dt

      + (2 dx/2 dy) dt c(u_sur - u_nj )

      + (dx/2 T) dt [K(u_n,j+1 - u_nj )/dy - K(u_n,j - u_n,j-1 )/dy]

      + (dy T) dt [c(u_sur - u_nj ) - K(u_nj - u_n-1,j )/dx].

The other portions of the boundary are similar.

Another question is related to the existence of solutions. Note the diagonal components of the coefficient matrix are much larger than the off-diagonal components. For each row the diagonal component is strictly larger than the sum of the absolute value of the off-diagonal components. In the finite difference model for (2) the diagonal component is 4/h² + c and the four off-diagonal components are -1/h².

Definition. Let A = [a_ij ]. A is called strictly diagonally dominant if and only if for all i

Examples. The one and two dimensional models of heat diffusion give strictly diagonally dominant matrices.

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.

Proof. Let x and y be two solutions. Then A(x - y) = 0. If x - y is not zero, then we can choose i so that |x_i - y_i | = max |x_j - y_j| > 0.

Divide by (x_i - y_i ) and use the triangle inequality to contradict the strict diagonal dominance.

The existence follows from the uniqueness. Since the matrix is square, we can find elementary matrices such that Ax = d is equivalent to EAx = Ux = Ed where U is upper triangular and E is a product of elementary matrices. If U has a zero diagonal component, then there will be more than one solution, which cannot happen. So, U must have all nonzero diagonal components, and we can construct via back substitution a solution of Ux = Ed, that is, a solution of Ax = d.

3.4.7.   Homework.