from numpy import * from numpy.linalg import * from myarray import * ## plot, show, axes, clf from pylab import * def condition_number( A ): ''' Returns the condition number of a matrix. The condition number is the ratio of absolute largest to absolute smallest eigenvalues (for symmetric matrices). ''' eigenvalues = eig( A )[0] abseigenvalues = abs( eigenvalues ) return max( abseigenvalues ) / min( abseigenvalues ) def lesson1(): ''' Prints the condition numbers of an identity matrix and an identity matrix with scaled rows. ''' A = identityf(2) print 'condition_number( [1 0; 0 1] ):', condition_number( A ) A[0,:] *= 1e20 A[1,:] *= 1e-20 print 'condition_number( [1e20 0; 0 1e-20] ):', condition_number( A ) def lesson2a(): ''' Plots all eigenvalues of our curve editing system. Note the 'S' shape. Keep in mind that the condition number is the largest divided by the smallest. ''' import direct_linear steps = 101 undeformed = zerosf(steps+1) h = 1 A,b = direct_linear.curve_editing_1d_system( undeformed, h ) #A,b = direct_linear.curve_editing_1d_system_fixed( {}, undeformed, h ) ## Plot the n eigenvalues eigenvalues = eig(A)[0] eigenvalues.sort() plot( eigenvalues ) plot( eigenvalues, 'go' ) ylabel('eigenvalue') title('eigenvalues of the curve editing system') def lesson2b(): ''' Plots the condition number of our curve editing system as a function of 'h', the resolution of our discrete curve approximation. ''' import direct_linear steps = 21 undeformed = zerosf(steps+1) fixed = {} ## Calculate the condition number as a function of h hs = arange( .01, .05, .001 ) cs = [ condition_number( direct_linear.curve_editing_1d_system_fixed( fixed, undeformed, h )[0] ) for h in hs ] plot( hs, cs ) plot( hs, cs, 'go' ) xlabel('h'), ylabel('condition number') title('h versus condition number, curve editing system') def improve_condition_number( A, b ): ''' Improves the condition number of a system by scaling all coefficients in an equation so that the diagonal entry is 1. Returns the improved system A,b. NOTE: Systems where the diagonal entry is not the largest coefficient in a row (equation) should be pivoted to make this true. ''' A = A.copy() b = b.copy() assert( A.shape[0] == A.shape[1] ) assert( A.shape[0] == b.shape[0] ) for row in range( A.shape[0] ): diag = A[ row, row ] oodiag = 1./diag A[ row, : ] *= oodiag b[ row ] *= oodiag return A,b def residual( A, x, b ): ''' Returns the residual of a system of equations Ax = b. The residual of a system of equations Ax = b is defined as |Ax-b|. If there were no roundoff error, this would be 0. ''' r = dot( A, x ) - b return sqrt( sum( r ** 2 ) ) def lesson3(): ''' Print the original and improved condition numbers of our curve editing system. Print the residual of the solutions obtained from the original system and the improved system. ''' import direct_linear steps = 101 undeformed = zeros(steps+1) fixed = {steps//2: 1.0} h = .01 A,b = direct_linear.curve_editing_1d_system_fixed( fixed, undeformed, h ) x = solve( A, b ) #x = dot( inv(A), b ) c = condition_number( A ) print 'condition number:', c print 'residual:', residual( A, x, b ) Ap,bp = improve_condition_number( A, b ) xp = solve( Ap, bp ) #xp = dot( inv(Ap), bp ) cp = condition_number( Ap ) print 'condition number, improved:', cp print 'residual, improved:', residual( Ap, xp, bp ) def main(): import matplotlib ion() lesson1() figure(1) lesson2a() figure(2) lesson2b() lesson3() show() if __name__ == '__main__': main()