Computer assignment 3/7/2007
Penta-diagonal system solver
Math 501, CSUF winter 2007
Team: Nasser M. Abbasi
CSUF, Math Dept.
Part
1. Derive special Gauss-elimination
Part
2. Derive special backward-substitution algorithm
Part 1. Derive special Gauss-elimination
Derive special Gauss-elimination strategy to transfer the resulting penta-diagonal system into an upper semi-penta system
The following algorithm was designed and developed to handle a general banded A matrix to solve the problem of solving Ax=b.
This algorithm works on Matrices which contain only one band of specific width such as those found in tri-diagonal and penta-diagonal matrices.
The main idea of the algorithm is to locate submatrices within the main matrix A, so as to process those by applying the standard Gaussian elimination algorithm on them.
The algorithm locates these submatrices which are bounded below and to the right by the first zero entry. Starting at first pivot in A(1,1), looking down and locating the first zero entry to determine the lower bound, and then looking right from that location to locate the first zero entry. This determines the boundaries of the submatrix.
This process is repeated by shifting one row down and one column to the right, and each time a new submatrix boundaries are located as described above, and Gaussian elimination is called to process this new lower submatrix.
Hence we travel down the main matrix from the top left corner to the bottom right corner, processing small submatrices along the way. The b vector is updated all the time. Hence in each step, we create a new separate Ax=b with its own A and b variables extracted from the original A and original b variables. Notice that no data copying is involved, and the data is processed in-place.
The advantage of this algorithm is that it will work on any central banded Matrix A, tri-diagonal, penta-diagonal and larger bands.
The algorithm is illustrated in the following diagram
Part 2. Derive special backward-substitution algorithm
Derive special backward-substitution algorithm to solve the resulting upper semi-penta diagonal matrix.
The resulting U matrix from part(1) is a banded matrix. Hence a special backward substitution algorithm was devised to take advantage of the sparseness of this U matrix.
Only the non-zero entries in each row are used to solve for x during the process of back substitution.
This is in comparison with the standard back substitution routine written for solving general Ax=b, which processed all entries in the upper triangular matrix regardless if the entries contain zero or not.
The following diagram illustrates the algorithm
Part 3. Long operation counts
Elimination process: each submatrix is of size 3x5. There are n-1 such submatrices
Each submatrix requires 2 divisions (for the multipliers), and 6 multiplications. (3 per row, we have 2 rows). Hence each submatrix requires 8 ops. Hence the total is 10(n-1)
For the special backsub: each x requires 2 multiplications and one division, and there are (n-1) rows to process. Also there is the first division for the x(n). Hence the total is 2(n-1)+1.
Adding the elimination process with the backsub, we obtain
8(n-1)+2(n-1)+1
= 10(n-1) + 1
Part 4. Implement in Matlab
Files used
|
File name |
Purpose |
|
nma_pentaSolve.m |
Main
driver. Called to solve Ax=b when A is penta-banded. Locates matrices within
the main matrix A and calls nma_gaussian_elimination
on each |
|
nma_backSub.m |
special
back substitution for special banded U matrices |
|
nma_gaussian_elemination.m |
Gaussian elimination routine |
|
nma_penta_test.m |
script that calls nma_pentaSolve
repeatedly with different A,b and compares the
results with Matlab "\" solver to check
for correctness. |
See appendix for source code listing.
Part 5 Test code
A Matlab script was written which tested the above implementation using different A,b input. Each test was verified against Matlab "\" solver. The following is the output of the test.
==========>test 1
A =
15
-2 -6 0
-2
12 -4 -4
-6
-4 19 -9
0
-1 -9 21
b =
300
0 0 0
=======>our result
x =
27.16548702392990
11.42568250758342
14.10515672396360
6.58914728682170
=======>Matlab result
ans =
27.16548702392989
11.42568250758342
14.10515672396360
6.58914728682170
==========>test 2
A =
15
8 -6 0
0 0
-2
12 -4 -4
0 0
-6
-4 19 -9
4 0
0
-1 -9 21
6 7
0
0 9 10
11 8
0
0 0 10
-2 2
b =
300
0 0 0
1 2
=======>our result
x =
1.0e+002 *
0.14783336170627
0.22913057224154
0.17509083392106
0.43838420195044
0.60556360806440
-1.57635740168778
=======>Matlab result
ans =
1.0e+002 *
0.14783336170627
0.22913057224154
0.17509083392106
0.43838420195044
0.60556360806440
-1.57635740168779
==========>test 3
A =
15
8 -6 0
0 0 0
-2
12 -4 -4
0 0 0
-6
-4 19 -9
4 0 0
0
-1 -9 21
6 7 0
0
0 9 10
11 8 3
0
0 0 10
-2 2 4
0
0 0 0
-2 2 4
b =
300
0 0 0
1 2 6
=======>our result
x =
22.60711151041847
10.45280071875325
20.45484640105050
-0.40000000000000
-53.69705242060892
75.01839040740875
-62.85772141400877
=======>Matlab result
ans =
22.60711151041846
10.45280071875324
20.45484640105048
-0.40000000000000
-53.69705242060884
75.01839040740865
-62.85772141400874
==========>test 4
A =
15
8
0
12
b =
300
0
=======>our result
x =
20
0
=======>Matlab result
ans =
20
0
==========>test 5
A =
15
8 -6 0
0 0 0
0
-2
12 -4 -4
0 0 0
0
-6
-4 19 -9
4 0 0
0
0
-1 -9 21
6 7
0 0
0
0 9 10
11 8 3
0
0
0 0 10
-2 2 4
3
0
0 0 0
-2 2 4
7
0
0 0 0
0 4 8
9
b =
300
0 0 0
1 2 6 10
=======>our result
x =
1.0e+002 *
0.13293913687916
0.13247435914403
0.00898032105660
0.32197318793590
1.01366621229967
-1.80430528168548
-0.00214695022696
0.81493296983974
=======>Matlab result
ans =
1.0e+002 *
0.13293913687916
0.13247435914403
0.00898032105660
0.32197318793589
1.01366621229967
-1.80430528168548
-0.00214695022696
0.81493296983974