sci.math.num-analysis SUNA44, 1-D problem condition numbers ===================================== A poster called "Condition Number" was submitted to "sci.math.num-analysis". 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 one-dimensional example below. Suppose one has a mesh of 1-D linear and identical elements: o ----- o ----- o ----- o ----- o ----- o ----- o ---- .... 1 2 3 4 5 6 7 The matrix of a single (1-D 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 end-points 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: | +2-L -1 | | -1 +2-L -1 | | -1 +2-L -1 | | -1 +2-L -1 | | -1 +2-L -1 | | ........... | | -1 +2-L | The first case can be done by hand: | +2-L -1 | gives (2 - L).(2 - L) - 1 = 4 - 4.L + L^2 - 1 = | -1 +2-L | = (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) 1981-1990 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([[2-L,-1,0],[-1,2-L,-1],[0,-1,2-L]]); # [ 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([[2-L,-1,0,0],[-1,2-L,-1,0],[0,-1,2-L,-1],[0,0,-1,2-L]]); # [ 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([[ 2-L, -1, 0, 0, 0 ], > [ -1, 2-L, -1, 0, 0 ], > [ 0, -1, 2-L, -1, 0 ], > [ 0, 0, -1, 2-L, -1 ], > [ 0, 0, 0, -1, 2-L ]]); # [ 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(n-1).(2 - L) - Det(n-2) 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.d-10 / 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 * E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood