sci.math.numanalysis
SUNA44, 1D problem condition numbers
=====================================
A poster called "Condition Number" was submitted to "sci.math.numanalysis".
The response to this query is compiled as the file "CNN.sum" (Condition Number
News.summary), at the same anonymous ftp site. Now we seek for applications of
this theory. Let's start with the onedimensional example below.
Suppose one has a mesh of 1D linear and identical elements:
o  o  o  o  o  o  o  ....
1 2 3 4 5 6 7
The matrix of a single (1D diffusion) element is:
 +1 1 
 1 +1 
Consequently, the grand matrix becomes, after assembly:
 +X 1  Here the numbers "X" are equal to "1".
 1 +2 1  If we add Dirichlet boundary conditions
 1 +2 1  at both endpoints of the mesh, while
 1 +2 1  accounting for them with the "matrix":
 1 +2 1 
 .........   1  , then X := 2 too.
 1 +X  
In order to determine the condition number of such matrices, the determinant of
the following one must be determined:
 +2L 1 
 1 +2L 1 
 1 +2L 1 
 1 +2L 1 
 1 +2L 1 
 ........... 
 1 +2L 
The first case can be done by hand:
 +2L 1  gives (2  L).(2  L)  1 = 4  4.L + L^2  1 =
 1 +2L  = (L  1).(L  3) = 0 Hence L = 1 , 3 .
Condition number = 3/1 = 3 .
After this, it already becomes somewhat complicated. Use has been made of MAPLE
while performing the next exercises. Comments are indicated with "#".
\^/ MAPLE V
._\ /_. Copyright (c) 19811990 by the University of Waterloo.
\ MAPLE / All rights reserved. MAPLE is a registered trademark of
<____ ____> Waterloo Maple Software.
 Type ? for help.
> with(linalg);
[ stuff deleted ]
> A := array([[2L,1,0],[1,2L,1],[0,1,2L]]);
#
[ 2  L 1 0 ]
[ ]
A := [ 1 2  L 1 ]
[ ]
[ 0 1 2  L ]
> solve(det(A),L);
1/2 1/2
2, 2 + 2 , 2  2
#
# Condition number = [2 + sqrt(2)] / [2  sqrt(2)] = 5.8284263 ...
#
> A := array([[2L,1,0,0],[1,2L,1,0],[0,1,2L,1],[0,0,1,2L]]);
#
[ 2  L 1 0 0 ]
[ ]
[ 1 2  L 1 0 ]
A := [ ]
[ 0 1 2  L 1 ]
[ ]
[ 0 0 1 2  L ]
> solve(det(A),L);
1/2 1/2 1/2 1/2
3/2 + 1/2 5 , 3/2  1/2 5 , 5/2 + 1/2 5 , 5/2  1/2 5
#
# Condition number = [2.5 + 0.5*sqrt(5)] / [1.5  0.5*sqrt(5)] = 9.4721335 ...
#
> A := array([[ 2L, 1, 0, 0, 0 ],
> [ 1, 2L, 1, 0, 0 ],
> [ 0, 1, 2L, 1, 0 ],
> [ 0, 0, 1, 2L, 1 ],
> [ 0, 0, 0, 1, 2L ]]);
#
[ 2  L 1 0 0 0 ]
[ ]
[ 1 2  L 1 0 0 ]
[ ]
A := [ 0 1 2  L 1 0 ]
[ ]
[ 0 0 1 2  L 1 ]
[ ]
[ 0 0 0 1 2  L ]
> solve(det(A),L);
1/2 1/2
1, 2, 3, 2 + 3 , 2  3
#
# Condition number = [2 + sqrt(3)] / [2  sqrt(3)] = 13.9282202 ...
#
# After this, the fun is over: _too_ complicated expressions occur.
>
> quit;
As a function of the matrix rank, condition numbers grow bigger & bigger:
Rank = 2 3 4 5 6
C.N. = 3 5.828 9.472 13.928 ??.??? .....
As a matter of fact, more can be done by hand, while evaluating the determinant
of the above matrices. Let us define:
Det (n) = the determinant of a matrix like above, with rank n .
Then it is a matter of routine to prove the following recurrence relationship:
Det(n) = Det(n1).(2  L)  Det(n2)
The difficult part remains, however, to solve the algebraic equations for L .
Time to invoke numerical means. The following (Fortran) program was written:
program eigen
*
* fc eigen.f lnag.fl lveclib
* a.out > eigen.out
*
parameter (NN=100)
double precision d(NN),e(NN),eps
*
do 100 N=2,NN
do 10 i=1,N
d(i)=2. ; e(i)=1
10 continue
e(1)=0.
*
data eps / 1.d10 /
call f02avf(N,eps,d,e,ifail) ! NAG
*
write(*,*) N,d(N)/d(1)
*
100 continue
end
And the outcomes are:
Matrix order Condition Number
 
2 3.00000000000000 
3 5.82842712474619  in accordance with
4 9.47213595499957  the MAPLE results.
5 13.9282032302755 
6 19.1956693580892
7 25.2741423690882
8 32.1634374775260
9 39.8634581890617
10 48.3741500787084
20 178.064274610860
30 388.812134493367
40 680.617070022187
50 1053.47899119953
60 1507.39787487881
70 2042.37371293282
80 2658.40650192184
90 3355.49624017367
100 4133.64292681876
200 16373.2418993205
300 36718.5355745369
400 65169.5239518397
500 101726.207032200
600 146388.584830023
700 199156.657373040
800 260030.424586180
900 329009.886771312
1000 406095.043797044
^^^^^^
Five or six significant digits lost for order 1000. Is that realistic?

* Han de Bruijn; Applications&Graphics  "A little bit of Physics * No
* TUD Computing Centre; P.O. Box 354  would be NO idleness in * Oil
* 2600 AJ Delft; The Netherlands.  Mathematics" (HdB). * for
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood