A First Course in Numerical Methods

Written for undergraduate and beginning graduate students, this book provides comprehensive coverage of modern techniques in scientific computing. The book provides an in-depth treatment of fundamental issues and methods, as well as the reasons behind the success and failure of numerical software. Topics include numerical algorithms, round-off errors, nonlinear equations in one variable, and linear algebra.

MATLAB is used to solve numerous examples in the book. In addition, a supplemental set of MATLAB code files is available for download.

Product Details

• Series: Computational Science and Engineering (Book 7)
• Paperback: 574 pages
• Publisher: Society for Industrial & Applied Mathematics (June 22, 2011)
• Language: English
• ISBN-10: 0898719976
• ISBN-13: 978-0898719970
• Product Dimensions: 6.8 x 1.2 x 9.7 inches
• Shipping Weight: 2.2 pounds
  

1. function [p,n] = bisect(func,a,b,fa,fb,atol)
   
2. %
   
3. % function [p,n] = bisect(func,a,b,fa,fb,atol)
   
4. %
   
5. % Assuming fa = func(a), fb = func(b), and fa*fb < 0,
   
6. % there is a value root in (a,b) such that func(root) = 0.
   
7. % This function returns in p a value such that
   
8. %      | p - root | < atol
   
9. % and in n the number of iterations required.
   
10. %
    
11. 
    
12. % check input
    
13. if (a >= b) || (fa*fb >= 0) || (atol <= 0)
    
14.   disp(' something wrong with the input: quitting');
    
15.   p = NaN; n=NaN;
    
16.   return
    
17. end
    
18. 
    
19. n = ceil ( log2 (b-a) - log2 (2*atol));
    
20. for k=1:n
    
21.   p = (a+b)/2;
    
22.   fp = feval(func,p);
    
23.   if abs(fp) < eps, n = k; return, end
    
24.   if fa * fp < 0
    
25.     b = p;
    
26.   else
    
27.     a = p;
    
28.     fa = fp;
    
29.   end
    
30. end
    
31. p = (a+b)/2;

复制代码

1. function [p] = bisect_recursive (func,a,b,fa,fb,atol)
   
2. %
   
3. % function [p] = bisect_recursive (func,a,b,fa,fb,atol)
   
4. %
   
5. % Assuming fa = func(a), fb = func(b), and fa*fb < 0,
   
6. % there is a value root in (a,b) such that func(root) = 0.
   
7. % This function returns in p a value such that
   
8. %      | p - root | < atol
   
9. 
   
10. p = (a+b)/2;
    
11. if b-a < atol
    
12.   return
    
13. else  
    
14.   fp = feval(func,p);
    
15.   if fa * fp < 0
    
16.     b = p;
    
17.     fb = fp;
    
18.   else
    
19.     a = p;
    
20.     fa = fp;
    
21.   end
    
22.   p = bisect_recursive (func,a,b,fa,fb,atol);
    
23. end

复制代码

1. function x = ainvb(A,b)
   
2. %
   
3. % function x = ainvb(A,b)
   
4. %
   
5. % solve Ax = b
   
6. 
   
7. [p,LU] = plu (A);
   
8. y = forsub (LU,b,p);
   
9. x = backsub (LU,y);

复制代码

1. function x = backsub (A,b)
   
2. %
   
3. % function x = backsub (A,b)
   
4. %
   
5. % Given an upper triangular, nonsingular n by n matrix A and
   
6. % an n-vector b, return vector x which solves Ax = b
   
7. 
   
8. n = length(b); x = b;
   
9. x(n) = b(n) / A(n,n);
   
10. for k = n-1:-1:1
    
11.   x(k) = ( b(k) - A(k,k+1:n)*x(k+1:n) ) / A(k,k);
    
12. end

复制代码

1. function y = forsub (A,b,p)
   
2. %
   
3. % function y = forsub (A,b,p)
   
4. %
   
5. % Given a unit lower triangular, nonsingular n by n matrix A,
   
6. % an n-vector b, and a permutation p,
   
7. % return vector y which solves Ay = Pb
   
8. 
   
9. n = length(b);
   
10. 
    
11. % permute b according to p
    
12. b = b(p);
    
13. 
    
14. %forward substitution
    
15. y = b;
    
16. for k = 2:n
    
17.   y(k) = b(k) - A(k,1:k-1) * y(1:k-1);
    
18. end

复制代码

1. function [p,A] = plu (A)
   
2. %
   
3. % function [p,A] = plu (A)
   
4. %
   
5. % Perform LU decomposition with partial pivoting.
   
6. % Upon return the coefficients of L and U replace those
   
7. % of the input n-by-n nonsingular matrix A. The row interchanges
   
8. % performed are recorded in the 1D array p.
   
9. 
   
10. n = size(A,1);
    
11. 
    
12. % initialize permutation vector p
    
13. p = 1:n;
    
14. 
    
15. % LU decomposition with partial pivoting
    
16. for k = 1:n-1
    
17. 
    
18.   % find row index of relative maximum in column k
    
19.   [val,q] = max(abs(A(k:n,k)));
    
20.   q = q + k-1;
    
21. 
    
22.   % interchange rows k and q and record this in p
    
23.   A([k,q],:)=A([q,k],:);
    
24.   p([k,q])=p([q,k]);
    
25. 
    
26.   % compute the corresponding column of L
    
27.   J=k+1:n;
    
28.   A(J,k) = A(J,k) / A(k,k);
    
29.   
    
30.   % update submatrix by outer product
    
31.   A(J,J) =  A(J,J) - A(J,k) * A(k,J);
    
32.   
    
33. end

复制代码

1. function [p,A] = plu_scaled (A)
   
2. %
   
3. % function [p,A] = plu_scaled (A)
   
4. %
   
5. % Perform LU decomposition with scaled partial pivoting.
   
6. % Upon return the coefficients of L and U replace those
   
7. % of the input n-by-n nonsingular matrix A. The row interchanges
   
8. % performed are recorded in the 1D array p.
   
9. 
   
10. n = size(A,1);
    
11. 
    
12. % find scales, initialize permutation vector p
    
13. s = max(abs(A'))';
    
14. p = 1:n;
    
15. 
    
16. % LU decomposition with partial pivoting
    
17. for k = 1:n-1
    
18. 
    
19.   % find row index of relative maximum in column k
    
20.   [val,q] = max ( abs(A(k:n,k)) ./ s(k:n) );
    
21.   q = q + k-1;
    
22. 
    
23.   % interchange rows k and q and record this in p
    
24.   A([k,q],:)=A([q,k],:);
    
25.   p([k,q])=p([q,k]);
    
26. 
    
27.   % compute the corresponding column of L
    
28.   J=k+1:n;
    
29.   A(J,k) = A(J,k) / A(k,k);
    
30.   
    
31.   % update submatrix by outer product
    
32.   A(J,J) =  A(J,J) - A(J,k) * A(k,J);
    
33.   
    
34. end

复制代码

1. % Example 5.11 : comparing scaled and unscaled GEPP for a special example
   
2. 
   
3. n = 100; h = 1/(n-1); K = 100;
   
4. A = zeros(n,n);
   
5. for i = 2:n-1
   
6.   A(i,i) = -2/h^2 - K;
   
7.   A(i,i-1) = -1/h^2; A(i,i+1) = -1/h^2;
   
8. end
   
9. A(1,1) = 1; A(n,n) = 1;  % end definition of A
   
10. 
    
11. xe = ones(n,1);          % exact solution of 1's  
    
12. b = A*xe;                % corresponding right hand side
    
13. 
    
14. %solve using ainvb
    
15. xu = ainvb(A,b);
    
16. err_ainvb = norm(xu-xe)
    
17. 
    
18. %solve using scaled GEPP version ainvb_scaled
    
19. xs = ainvb_scaled(A,b);
    
20. err_ainvb_scaled = norm(xs-xe)

复制代码

1. % Example 5.13 : Cholesky for random matrices
   
2. 
   
3. n = 500;
   
4. C = randn(n,n); A = C'*C;
   
5. xe = randn(n,1); % the exact solution
   
6. b = A*xe;        % generate right hand side data
   
7. 
   
8. R = chol(A);     % Cholesky factor
   
9. % the following line is for compatibility with forsub
   
10. D = diag(diag(R)); L = D \ R'; bb = D \ b; p = 1:n;
    
11. y = forsub(L,bb,p);            % forward substitution R'y = b
    
12. x = backsub(R,y);              % backward substitution Rx = y
    
13. rerx = norm(x-xe)/norm(xe)     % error by Cholesky
    
14. 
    
15. xd = ainvb(A,b);               % ignore spd and use partial pivoting
    
16. rerxd = norm(xd-xe)/norm(xe)   % error by general routine

复制代码