Computational Mathematics II

winter semester 2011-12

Syllabus

1. Iterative methods for solving
   • linear systems of equations: sparse formats of matrices, stationary iterative methods (incl. Jacobi, Gauss-Seidel, SOR, Richardson), projection methods (incl. steepest descent method, minimal residual method), Krylov subspace methods (incl. CG and GMRES), preconditioning (ILU, spectral equivalence of matrices, PCG, preconditioned GMRES, termination conditions), structural methods (reduction to Schur matrix), multigrid methods (only 2-grid was discussed) Basic references: posiiton vi. (Saad - Iterative... ) and 4. (Kelley - Iterative....)
   • nonlinear systems of equations - Banach iterations, multidimensional Newton methods, the chord method, Newton method with Jacobian approximated by finite differences, inexact Newton methods. Basic reference: position 4. (Kelley - Iterative....)
2. Eigenproblems - methods for finding (approximations) of eigenvalues and eigenvectors of a matrix - similar matrices, transformation of a matrix into a similar matrix in the Hessenberg form by using Householder reflections, the power and inverse iterations, the Rayleigh quotient iteration, pure QR and shifted QR methods, "divide and conquer" method, the Hyman method. The basic reference: position 8. (Trefethen, Bau, NLA) and 7. (Stoer at al, Intro. to Numer. Anal.- in particular Hyman met.)
3. Multidimension integration - nothing was done - the dimensionality curse, the Monte Carlo (MC) and quasi MC methods should have been discussed

Lecture notes

References

Text books

1. James W. Demmel, Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics (SIAM), Philadelphia 1997.
2. Peter Deuflhard, Andreas Hohmann. Numerical analysis in modern scientific computing, vol. 43 in Texts in Applied Mathematics. Springer-Verlag, New York, 2nd edition, 2003. An introduction.
3. J.M. Jankowscy, M. Dryja. Przegląd metod i algorytmów numerycznych, tom I i II. Biblioteka inżynierii oprogramowania. Wydawnictwo Naukowo-Techniczne, Warszawa, 1995.
4. C. T. Kelley. Iterative methods for linear and nonlinear equations, vol. 16 of Frontiers in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1995 - covers iterative methods for linear equations (CG, PCG and GMRES), and solving nonlinear equations
5. A. Kiełbasiński, H. Schwetlick. Numeryczna algebra liniowa. Wydawnictwa Naukowo-Techniczne, 1992.
6. D.Kincaid, W.Cheney, Numerical Analysis, 2nd editione, Brooks/Cole, 1996. A general textbook for numerical methods.
7. J. Stoer, R. Bulirsch. Introduction to numerical analysis. Translated from the German by R. Bartels, W. Gautschi and C. Witzgall. Third edition. Texts in Applied Mathematics, 12. Springer-Verlag, New York, 2002
8. Lloyd N. Trefethen, David Bau, III, Numerical linear algebra. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997.

Monographs or advanced text books

1. John E. Dennis Jr., Robert B. Schnabel. Numerical methods for unconstrained optimization and nonlinear equations. Prentice-Hall Series in Computational Mathematics. Prentice-Hall Inc., En- glewood Cliffs, N.J., 1983.
2. Peter Deuflhard. Newton methods for Nonlinear Problems. Affine Invariance and Adaptive Algori- thms. Springer International, 2002.
3. Eugene G. Dyakonov. Optimization in solving elliptic problems. CRC Press, Boca Raton, FL, 1996. Translated from the 1989 Russian original, Translation edited and with a preface by Steve McCormick.
4. Gene H. Golub, Charles F. Van Loan. Matrix computations. Johns Hopkins Studies in the Mathe- matical Sciences. Johns Hopkins University Press, Baltimore, MD, 3rd ed., 1996.
5. J. M. Ortega, W. C. Rheinboldt. Iterative solutions of nonlinear equations in several variables. Academic Press, New York, 1970.
6. Yousef Saad. Iterative methods for sparse linear systems. Society for Industrial and Applied Ma- thematics, Philadelphia, PA, 2nd ed., 2003. On-line
7. Yousef Saad, Numerical Methods for Large Eigenvalue Problems. Society for Industrial and Applied Ma- thematics, Philadelphia, PA, 2nd ed., 2011. On-line
8. A. A. Samarski, J. S. Nikołajew. Metody rozwiązywania równań siatkowych. PWN 1988.
9. Barry F. Smith, Petter E. Bjorstad, William D. Gropp. Domain decomposition. Cambridge Univer- sity Press, Cambridge, 1996. Parallel multilevel methods for elliptic partial differential equations.
10. Andrea Toselli, Olof Widlund. Domain decomposition methods - algorithms and theory, wolumen 34 serii Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 2005.
11. J. F. Traub. Iterative Methods for the Solution of Equations. Englewood Cliffs, New York, 1964.

LAB

Program of labs and octave scripts with solution of some problems

• Lab 1 (10 X 2011)- introduction to octave; direct methods for solving system of linear equations - LU factorization
  mo2basic.m - a simple octave script with basic operations like matrices multiplications, basic structures like loops etc
  
• Lab 2 (Oct 24th, 2011) - direct methods for solving linear systems (LU and L'L factorizations), sparse matrices, Jacobi method
  • Check how functions lu() i chol() work. Using chol() check if tridiagonal matrix with -1 on sub- and super diagonals and zero on the main diagonal is positive definite. Compare times of solving system with T_N i.e. tridiagonal matrix with 2 on the diagonal and -1 on sub- and super diagonals, using \ and lu() or chol(). Create this matrix in sparse format using sparse() octave function NOT creating any full matrix. Compare time of solving the system with T_N in the sparse or full T_N format. Find the max N for which octave solver works (both for T_N sparse/full)
  Examples and exercises in section 1.3.1 of lect notes.
  • Test octave script to example 1.8 in lect note i.e we want to test the convergence of Jacobi method for A=B+pI with p a parameter. We choose p in order to get a diagonally dominant matrix
  • Solve exercise 1.9, i.e check how the value of p influence the convergence rate of the Jacobi method; Test the time of solving the system with respect to N; test if for N=100 the Jacobi method is it reasonable to use it (comparing to the direct method)
  • exercise 1.10: add or modify the lines of code to the script of example 1.8 in which the rate of convergence in max norm is tested; (the rate of the norm reduction) and compare with the theoretical results. Test if the symmetry does influence the convergence rate.
  • Test octave script to example 1.9 in lect note we want to test the Jacobi method applied to the system with T_N+p*I - (T_N tridiagonal with 2 on the main diagonal, and -1 on the remaining diagonals) for p>=0 - we will test the time for solving the system by the Jacobi method and compare it with the time of the direct solver.
  • exercise 1.11: explain why for large p - the convergence is faster; increase N and repeat the the tests: which method wins for large N? Is it dependent on p?; Decrease tol to the level fl error e.g. take TOL=1e-16 - which method is now better?
• Lab 3 (Nov 14th, 2011) - Jacobi method, cont., Gauss-Seidel, SOR methods and Richardson methods. Projection methods: steepest descent method.
  Test examples 1.10, 1.12; solve exercises 1.13 and 1.15 in section 1.3.1 of lect notes:
  octave script to example 1.10: (tests of Jacobi and Gauss-Seidel methods for B+pI - B random matrix; p parameter);
  octave script to example 1.12: (tests of Jacobi, Gauss-Seidel and SOR methods for T_N+pI - T_N tridiagonal [-1,2,-1]; p parameter)
  • Implement Richardson method e.g. modifying script of example 1.10. Test Richardson method with parameter [10,5,1,1/5,1/0,1/100] for B=A^TA - B random matrix as in example 1.8.
  • Test Richardson in 2nd norm for T_N N=1000 with known solution for different values of parameter, and then with optimal parameter (why do we know it?) for T_N for N=10,100,1000, 10000 etc Test convergence of Richardson for T_N and optimal parameter but in energy norm. Does the results confirm the theory? Compare times with the octave direct solver (backslash). Test Richardson method for T_N+p*I for varying values of parameters and random initial starting vectors.
  • Test example 1.13 in section 1.5.1:
    octave script to example 1.13 (tests of Jacobi, Gauss-Seidel and steepest descent methods for B'B+pI - B random matrix; p parameter)
  • Test steepest descent method x_{k+1}=x_k+ p_k r_k for p_k=(r_k^Tr_k)/(r_k^TAr_k) for the same matrices as in the tests of Richardson method i.e. T_N for N=10,100,1000 etc and A'A with a random A. Compare the results. (SDM is implemented in example 1.13)
  • Implement and test minimal residuum method x_{k+1}=x_k+ p_k r_k for p_k=(r_k^TA^Tr_k)/(r_k^TA^TAr_k) for the same matrices as in the tests of Richardson method and B random matrix and T_N+ (1/N)*T_f where T_f - tridiagonal with zeros on the main diagonal, ones on super-diagonal and minus ones on sub-diagonal. Test also the method for the matrix for A=[0,1;1,0] and x_0=0, b=[1,0]^T or in general for N-cycle permutation matrix in R^N : Av_k=v_{k+1} - v_k kth unit vector with v_{N+1}=v_1.
• Lab 4 (Nov 28th,2011) Tests of CG method.
  • Test octave function pcg(), i.e. read help and run it for matrix [2,1;1,2] with a random right hand side.
  • Run script from example 2.1 from Section 2.1.1 from lecture notes:
    script with the code of example 2.1 , Continuation of example 1.13. We compare 4 methods: Jacobi, SOR, steepest descent and cg for a matrix A = B' B + pI, where p > 0 , B is sparse and random. We increase p which improves also condition of A . SOR parameter was set as \omega= 1.3 .
  • Exercise 2.7. Check in example 2.1 if the convergence speed of cg and the steepest descent method depends on cond(A). To check cond(A) use cond() octave function.
  • Exercise 2.8. Modify code of example 2.1, to check if A is nonsymmetric or non-definite, then the convergence rate of cg will change (and if so, how it is changed). take A = B + pI for p > 0 (no symmetry) selected in such a way, that: A_{sym} > 0 and A = B' B + pI for p < 0 such that A ahs both positive and negative eigencvalues.
  • Run script from example 2.2:
    script with the code of example 2.2 , We want to compare 4 methods: Jacobi, SOR, steepest descent and cg for a matrix T_N (1D laplacian matrix). As a parameter for SOR we take the optimal value cf. example 1.11.
  • Modify the code of example 2.2 in order to test CG for T_N+pI for p=0,5,10,20,50 for N=25 and 50.
  • Modify the code of example 2.2 in order to test cg for P_N - a matrix of 2D laplacian discretized on a uniform mesh script with code which creates this matrix
• Lab 5 (Dec 12th,2011) Tests of multigrid
  • Multigrid for T_N - test smoother - take random s a solution, compute b=T_Ns and use 3-4 iterations of Richardson with parameter a=0.25 with random x0 and plot the errors: x0-s and x-s (x obtained by 2-3 R. iterations). Compute now plots of eigenvalues of I- aT_N for a =1,0.5,0.25,0.125 for N=100. Find using eigs() eigenvalues and eigenvalues of T_N - plot eigenvectors for the largest and the smallest eigenvalues i.e. Q(:,1) and Q(:,N).
  • Download to your computer octave m-file with a function: ExtensionCreate.m which creates extension matrix E:R^{N_0}-->R^N (N_0>N). Create E for N=100 and N_0 = 25 and plot graphs of columns of E (1st, 20th etc)
  • Download m-file twogrid.m and test for T_Nx=b with N_0=N/2. for random b - plot the graph of norms of residual vectors for N=24,100,1000. Is the convergence rate dependent on N?
• Lab 6 (January 2nd, 2012) Tests of GMRES method.
  • Download to your computer octave m-file with a function: gmres() . Test this octave function gmres(), i.e. read help and run it for matrix [2,1;3,2] with a random right hand side.
  • GMRES for spd matrices. Change script from example 2.1 adding tests for gmres(). (Section 2.1.1 from lecture notes:
    script with the code of example 2.1 , Continuation of example 1.13.) Thus we compare 5 methods: Jacobi, SOR, steepest descent, cg and gmres for a matrix A = B' B + pI, where p > 0 , B is sparse and random. We increase p which improves also condition of A . SOR parameter was set as \omega= 1.3 .
  • Modify code of example 2.1, to check the convergence fo GMRES for random nonsymmetric function. Take B - random 25X25 and then test GMRES for the systems Ax=b with:
    1. A = B + pI for p > 0 (no symmetry but positive definite for larger p)
    2. A = B'B-pI for p>0 (so A is symmetric but not necessarily positive definite)
    3. A - orthogonal e.g. find QR factors of random B and take A=Q
  • Modify script of example 2.2:
    script with the code of example 2.2 , We want to compare 5 methods: Jacobi, SOR, steepest descent, cg and GMRES for a matrix T_N +pI (T_N - 1D laplacian matrix). As a parameter for SOR we take the optimal value cf. example 1.11. Test for p=0,5,10,20,50 for N=25 and 50.
  • Modify the code of example 2.2: in order to test gmres for T_N + (1/N)T_f + p I for p=0,5,10,20,50 for N=25 and 50, where T_f - tridiagonal with zeros on the main diagonal, ones on super-diagonal and minus ones on sub-diagonal.
  • Modify the code of example 2.2: in order to test gmres for P_N +pI - a matrix of 2D laplacian discretized on a uniform mesh script with code which creates this matrix
  • Test fsolve.m for $x+y=1;x^2+y^2=1$ and for Allen Cahn equation $-u''+\delta u(1-u^2)=f$ with $u=0$ on boundary. After FDM discretization we get: $h^{-2}T_N U + \delta U.*(1-U^3)=F$ - $delta=1,1e-2,1e2 etc$
  • Implement Newton with FD approximation of Jacobian.
  • Test Newton for $x+y=1;x^2+y^2=1$ and for Allen Cahn equation $-u''+\delta u(1-u^2)=f$ with $u=0$ on boundary. After FDM discretization we get: $h^{-2}T_N U + \delta U.*(1-U^3)=F$
  • Test Newton with Jacobian approximated by FD on the problem $x+y=1;x^2+y^2=1$ for different values of FD parameter.
  Some scripts (in gzipped tar file): modified script of example 2.1 with tests of GMRES, the newton method (classic and wit hJacobian approximated by FD) with test of order of convergence etc: lab6.tgz
• Lab 7 (January 16th, 2012) Tests of Newton, power iteration, inverse iteration, Rayleigh quotient iteration, QR method etc
  • Test if fsolve() finds the solution of the system (x,a)|-->[Ax-ax;0.5*x'*x-1]=0 for A=A^T e.g. A=[2,1;1;2]
  • Download to your computer octave m-file with a function: power() . Test this octave function power(), i.e. read help and run it for matrix A=[2,1;1;2] and for A^{-1}. Modify the code in order to print Rayleigh quotient r_k and x_k in each iteration and the errors |r_k-\lambda_1| and \|x_k-q_1\| (for known eigenpair (\lambda_1,q_1)).
  • Repeat tests of power.m but for NONSYMMETRIC matrix with known eigenvalues e.g. A=CDC^{-1} with D=diag(3,2,1) and C a random nonsingular matrix.
  • Modifying power.m write a m-file: inverse.m with inverse iteration (add an extra parameter - for the inverse method) Test for A=[2,1;1;2] and a matrix with known eigenvalues: M,2,1 for M=3,5,10,100. Print out the errors of eigenvector and eigenvalue.
  • Modifying power.m write a m-file: Rayleigh.m with Rayleigh quotient iteration (inverse iteration with changing parameter - in each iteration we take r_k Rayleigh quotient as a shift) Test for A=[2,1;1;2] and a matrix with known eigenvalues: M,2,1 for M=3,5,10,100. Print out the errors of eigenvector and eigenvalue.
  • Modifying power.m that it works for a matrix A such that it has symmetric eigenvalues of the largest magnitude e.g. A symmetric and it has as eignvalues M,-M,2,1 for M=3,5,10,100.
  • Implement QR method with or without shift (use qr() function of octave) as a shift take Rayleigh shift i.e. (A_k)_{N N}.