ACMS 40390: Fall 2013 - Numerical Analysis

Course Information

Instructor: Zhiliang Xu, zxu2@nd.edu Hayes-Healy 226

Office phone number: (574)631-3423

Class time: M, W, F, 8:20am - 9:10am.

Classroom: Edward J. DeBartolo Hall 129

Teaching assistant: Tian Jiang, tjiang@nd.edu, Hayes-Healy 215

Office hours:

Xu: W 4:00pm - 6:00pm or by appointment and drop-in (when available), HH226
Jiang: R 4:00pm - 5:00pm, HH215

Textbook: Richard L. Burden and J. Douglas Faires, Numerical Analysis, 9th Edition.

Prerequisites:

MATH 20750 or MATH 20860 or MATH 30650 or ACMS 20750 or PHYS 20452

The course requires a moderate amount of programming. C, C++ or Fortran programming languages are preferred. However, you may also use software programs including Matlab, Mathematica.

Syllabus is here.

Tips about programming:

Important Dates
First test review Monday, 09/23 in class
First in class test Wednesday, 09/25 in class
Second test review Monday, 10/28 in class
Office Hours 4:30pm - 6:30pm, Tuesday, 10/29 Hurley 154
Second in class test Wednesday, 10/30 in class
Final test review 5:00pm - 7:00pm, Tuesday, 12/17 Hayes-Healy, 231
Office Hrs for Final Test 4:00pm - 6:00pm, Wednesday, 12/18 226 Hayes-Healy
Final test 8:00am - 10:00am, 12/19 129 DeBartolo

Tutorial Session
Computer Lab 1 Tueday 09/03, 4:00pm - 5:30pm DeBartolo 140
Computer Lab 2
(Bisection method; Fixed point iteration) 		Tueday 09/17, 4:30pm - 5:30pm Hayes-Healy 117
Computer Lab 3
(Newton divided difference; Hermite interpolation) 		Monday 10/07, 4:30pm - 5:30pm Hayes-Healy 127
Computer Lab 4
(Runge-Kutta methods) 		Monday 11/11, 4:30pm - 5:30pm 127 Hayes-Healy
Computer Lab 5
(Multistep methods) 		Tuesday 11/19, 4:00pm - 5:00pm 127 Hayes-Healy
Computer Lab 6
(Gaussian elimination) 		Friday 11/22, 3:00pm - 4:00pm 127 Hayes-Healy

Week 1:

  • Wednesday (08/28): Sections 1.1 HW:
    Computer project: Setup your account for computer projects. Compile and run the C++ code on the machine "darrow" and submit your work following the instruction at Instruction and help on Unix systems . The code for practising is ex0-1.cpp .

    Section 1.1 Review of Calculus

    Section 1.2 Arithematic

    Matlab code for finding extreme values of functions .
    Computer Lab (Bring your own laptop): When: 4:00pm - 5:00pm, 09/03 (Tu). Where: 140 DeBartolo.
    Increase AFS space quota: email oithelp@nd.edu to ask for 200 MB AFS disk space if you are NOT engineering students.
  • Friday (08/30): Sections 1.1, 1.2 HW: Section 1.1: 1(b,c), 2(b), 4(c, you can also use Matlab to solve for the equation f'(x) = 0. To do this, take the above Matlab code and modify.), 9, 16 (Use Matlab or Maple), 19 ( "within 10^-6" means that the error bound should be less than 10^-6), 21.

    Short review of binary numbers.

    A Short Tutorial on The Binary System

    In the decimal system, we use the digits 0-9 to represent numbers, and things are organized into columns. Take the decimal number 285 for example: 

        H | T | O
    2 | 8 | 5
    such that "H" is the hundreds column, "T" is the tens column, and "O" is the ones column. So the decimal number "285" is 2-hundreds plus 8-tens plus 5-ones.

    The ones column stands 10^0, the tens column stands 10^1, the hundreds column stands 10^2 and so on, so 

          10^2|10^1|10^0
        2 |  8 |  5
    the number 285 is really {(2*10^2)+(8*10^1)+(5*10^0)}.

    The binary system works under the exact same way as the decimal system, except it operates in the base 2 rather than the base 10. In other words, instead of columns being

           10^2|10^1|10^0
    they are 
            2^2|2^1|2^0

    Instead of using the digits 0-9, we only use 0-1 digits.

    Examples: What will the binary number 1011 be in the decimal system?

    Answer: By using the columns, we have 

            2^3|2^2|2^1|2^0
         1 | 0 | 1 | 1
    Therefore 
      1011=(1*2^3)+(0*2^2)+(1*2^1)+(1*2^0)
      = (1*8) + (0*4) + (1*2) + (1*1)
  = 11 (in decimal system)

    Week 2:

  • Monday (09/02): Section 1.2 HW: Section 1.2: 1(a,b,c), 3(b,d), 4(b,c,d), 12, 15(a,b), 16(use the number given in 15(a,b)).
  • Wednesday (09/04): Sections 1.2 HW: Section 1.2: 13, 17, 22, 24

    Section 1.3 Algorithms and Convergence

    Summary of Pseudo-code Language Constructions:

    An algorithm is an ordered sequence of unambiguous and well-defined instructions that performs some tasks.

    Three Categories of Algorithmic Operations

    1. sequential operations(Sequence) - instructions are executed in order
    2. conditional  ("question asking") operations - a control structure that asks a true/false question and then selects the next instruction based on the answer
    3. iterative operations (loops) - a control structure that repeats the execution of a block of instructions

    Computation/Assignment 

            set the value of "variable" to :"arithmetic expression" or
            "variable" equals "expression" or
            "variable" = "expression"

    Input/Output

            get "variable", "variable", ...
            display "variable", "variable", ...

    Conditional

                if  "condition" then
                            (subordinate) statement 1 
                             etc ...
                else
                            (subordinate) statement 2
                            etc ...

    Iterative

                while "condition" 
                            (subordinate) statement 1
                            (subordinate) statement 2 ...

                for "iteration bounds" 
                            (subordinate) statement 1
                            (subordinate) statement 2 ...

  • Friday (09/06): Sections 1.2, 1.3 HW: Section 1.2: 19(a), 21, 28.
    Matlab code for bisection method .
    The sample code for the bisection method is ALG021_bisection.cpp .

    Week 3:

  • Monday (09/09): Sections 1.3 HW: Section 1.3: 1(b), 6(a,c), 7(a, d), 9, 10, 14, 15.
  • Wednesday (09/11): Sections 1.3, 2.1 HW: Section 2.1: 2(a), 9, 11, 14, 19(Computer project).
    Use Alg. 2.1 (the bisection method) to solve Exercise 2.1.19.
    Submit the code in acms40390/hw1 with ex1-1.cpp.

    Section 2.1 The Bisection Method

    Matlab code for bisection method .
    The sample C code for the bisection method is ALG021_bisection.cpp .
  • Friday (09/13): Sections 2.1, 2.2 HW: Section 2.2: 1(a,c), 2(do 1(a,c)), 16(a), 18(b)

    Section 2.2 Fixed Point Iteration

    Main Matlab code for fixed-point iteration , Matlab code of Routine for fixed-point .
    The sample C code for the fixed-point iteration method: ALG022_fixed_pt_iter.cpp

    Week 4:

  • Monday (09/16): Sections 2.2, 2.3.
    HW: Section 2.2: 8(Use Thm 2.3 to show g(x) has a unique fixed point. Use Corollary 2.5 estimate the number of iterations required to achieve 10^-4 accuracy in abs. error.), 20(a).
    Section 2.3: 2 (write down each step).
    Computer project:
    Use Alg. 2.2 (fixed point iteration) to solve Exercise 2.2.6.
    Submit the code in acms40390/hw2 with ex2-1.cpp.

  • Section 2.3 Newton's Method

  • Algorithm of the Newton's method.

    Matlab demo code for Newton's method

    The sample code for the Newton's method is ALG023_Newton.cpp
  • Algorithm of the Secant method.

    Matlab demo code for Secant method

    The sample code for the Secant method is ALG024_Secant.cpp
  • Wednesday (09/18): Sections 2.3, 2.4
    HW: Section 2.3: 3(a), 4(a). Section 2.4: 6, 7, 8(a).
    Computer project:
    Use Newton's method and Secant method respectively to solve Exercise 2.3.26. Which one converges faster?
    Submit the code in acms40390/hw2 with ex2-2.cpp (for Newton's method) and ex2-3.cpp (for Secant method) respectively.
  • Friday (09/20): Sections 2.4, 2.5
    HW: Section 2.4: 10.
    Section 2.5: 4, 5, 13(b).
    Computer project:
    Use the modified Newton's method to solve Exercise 2.4.5.
    Submit the code in acms40390/hw2 with ex2-4.cpp.

  • Section 2.4 Error Analysis for Iterative Methods

    modified_Newton_241.cpp

  • Section 2.5 Accelerating Convergence

    ALG026_ALG026_Steffensen.cpp

    Week 5:

  • Monday (09/23): In class review
  • Wednesday (09/25): Exam
  • Friday (09/26): Sections 2.6, 3.1
    HW: Section 2.6: 2(a), Let x0=1, compute f(x0) by Horner's method and find the expression of Horner's decomposition, which is the last eqn of Theorem 2.19. 2(h), Let x0=2, find the expression of Horner's decomposition, which is the last eqn of Theorem 2.19.
    Section 3.1: 2(a,c), 4(analyze exercises 2(a,c)).

  • Section 2.6 Zeros of polynomials and Horner's method

  • Algorithm for the Horner's method is here.
    The sample code for the Horner's method is ALG027_Horner_2.cpp

    Week 6:

  • Monday (09/30): Sections 3.1, 3.3 HW: Section 3.1: 6(a, justify how you choose points), 8(a), 14, 19 (a. Use the Matlab code).
    Section 3.3: 2.

  • Section 3.1 Interpolation and the Lagrange Polynomial


    The Matlab code for the Lagrange Interpolating Polynomial is Lagrange Interpolating Polynomial and Lagrange basis Polynomial .
  • To use this sample code to solve Example 3.1.1, do the following:
    Download these two files and in Matlab add paths where we save these two files;
    In the command window, excute the following statement to see the curve of the 2nd degree Lagrange polynomial for interpolating points taken from 1/x:
    xi = [2.0, 2.5, 4.0];
    fi = 1./xi;
    lp = lagrange_2(xi, fi);
    xx = 0.5:0.01:5;
    plot(xx, polyval(lp, xx), xi, fi, 'or', xx, 1./xx, '--g');
    To evaluate the value of the interpolating polynomial at a number x, use the command: polyval(lp, x).
  • Remark: the return values from calling " lagrange_2(xi, fi)" are coefficents of the interpolating polynomial from the highest degree to zeroth degree in the form of P_n(x) = a_n*x^n + a_{n-1} *x^{n-1} + ... + a_0. (a_n, ..., a_0) are return values.

  • Section 3.3 Divided Differences

    Matlab function for computing Newton divided differences table
    C code of Newton divided differences
    The C code needs the input data.
  • Wednesday (10/02): Sections 3.3, 3.4 HW: Section 3.3: 7, 16, 17, 19.
    Section 3.4: 1(a,b).
    Computer project:
    Use Alg. 3.2 (Newton's Divided-Difference) to solve Exercise 8 of Section 3.3.
    Submit the code in acms40390/hw3 with ex3-1.cpp.
    Hint: Use the above C code of Newton divided differences to start.

  • Section 3.4 Hermite Polynomial; Section 3.5 Cubic Splines

    The Matlab code for the Hermite Polynomial to interpolate function 1/x
    The code needs the input data.
    To run this code, open it in Matlab and click on "run" button. Then you will see options to input data set. If you want to input data from a file. The format of the data set in the file should follow the one in the sample data file.

    C code of Hermite Interpolation
    The Hermite C code needs the input data.

  • Friday (10/04): Sections 3.4, 3.5 HW: Section 3.4: 2(a,c, use Alg. 3.3 which is divided difference approach), 4(a,c), 6(a).

    Computer project:
    Use Alg. 3.3 (Divided-Difference) to solve Exercise 10 of Section 3.4.
    Submit the code in acms40390/hw3 with ex3-2.cpp.
    Hint: Use the above C code of Hermite Interpolation to start.

    The cubic spline Matlab code
    The input data for Ex3.5.3 is here .
    To run this cubic spline code, open it in Matlab and click on "run" button. Then you will see options to input data set. If you want to input data from a file. The format of the data set in the file
    Oscillation at the edges of an interval when using high degree polynomial interpolation ( known as Runge's phenomenon ).
    Consider the function: f(x) = 1/(1+25*x*x) on [-1, 1].
    Run the following Matlab code:
    xi = -1.0:0.1:1.0;
    fi = 1./(1.0+25.0*xi.*xi);
    lp = lagrange_2(xi, fi);
    xx = -1.0:0.05:1.0;
    plot(xx, polyval(lp, xx), xi, fi, 'or', xx, 1./(1.0+25*xx.*xx), '--g');
  • Interpolation oscillates toward the end of the interval.

    Week 7:

  • Monday (10/07): Sections 3.5, 4.1 HW: Section 3.5: 2, 11.
  • Wednesday (10/09): Sections 4.1, 4.3 HW: Section 4.1: 1(a), 6(a,c), 8(a,c), 13, 15(do exercise 1(a)), 26(use the most accurate formula for each of time points).

  • Section 4.1 Numerical Differentiations

  • Friday (10/11): Section 4.3 HW: Section 4.3: 2(a,b,c), 4, 6, 8, 10, 12 (for exercises 4, 6, 8, 10, 12, do problems 2(a,b,c)), 14, 16, 17.

  • Section 4.3 Elements of Numerical Integration

    Week 8:

  • Monday (10/14): Section 4.4 HW: Section 4.4: 1(a,b,c), 3, 5 (for exercises 3, 5, do problems 1(a,b,c)), 9, 11, 22 (use composite Simpson's rule).

  • Section 4.4 Composite Numerical Integration

  • Wednesday (10/16): Sections 4.6, 4.7 HW: Section 4.6: (a,c), 2(do 1(a, c)), 9. Section 4.7: 1(a, b, c).

  • Section 4.6 Adaptive Quadrature Methods


    The sample code for the adaptive quadrature is ALG043.cpp .
  • Friday (10/18): Sections 4.7, 4.8 HW: Section 4.7: 2, 3 (do problems 1(a,c) for exercises 2 and 3), 5, 7 (do n = 3 case).
    Section 4.8: 5(Only use Gaussian quadrature to solve problems 1(a, b)).

  • Section 4.7 Gaussian Quadrature

  • Section 4.8 Multiple Integration and MC Integration


    Matlab code for MC simulation of Buffon's needle problem

    Matlab code for MC integration

    Week 9:

  • Fall break

    Week 10:

  • Monday (10/28): Exam review
  • Wednesday (10/30): Exam
  • Friday (11/01): MC Integration, Sec. 5.1

    Week 11:

  • Monday (11/04): Sec. 5.1, 5.2
    HW: Section 5.1: 1(a,b,c), 2(a,b), 4. Section 5.2: 2(a,b), 4 (a,b), 9 (note: part a could be done by using the Matlab code), 11(a) (note: turn in the matlab plots for different h values), 12.

  • Section 5.1 Theory of IVP; Section 5.2 Euler's Method

  • Euler's Method Matlab code .
  • Wednesday (11/06): Sec. 5.3
    HW: Section 5.3: 2(a,b), 4(a,b), 9(a, b, c, d. Note: modify the sample Matlab code to solve a and c; turn in the Matlab plots of the solutions and the hard copy of the Matlab code.).

  • Section 5.3 Higher-Order Taylor Method

  • Taylor method of order two to solve Example 1(section 5.3) (Matlab code)

  • Friday (11/08): Sec. 5.4
    HW: Section 5.4: 2(a,c), 4(a, implement a Matlab code to solve. Turn in the solution plot which consists of both approximate and exact solutions. Hint: take the sample 4th order Runge-Kutta Matlab code and modify it), 6(do problems 2(c)), 14 (do problems 2(c)), 31.

    Computer project:
    1. Solve problem 5.4.28(a) by the Runge-Kutta method of order 4. Submit the code in acms40390/hw5 with ex5-1.cpp (Hint: Use the modified Euler's method C code as the base and see how different stages are implemented in the 4th order Matlab code)
    2. Solve problem 5.4.28(a) by the Heun's Runge-Kutta method. Submit the code in acms40390/hw5 with ex5-2.cpp (Hint: Modify the sample modified Euler method C code)

  • Section 5.4 Runge-Kutta Methods

    4th order Runge-Kutta method to solve Example 3(section 5.4) (Matlab code)

    Modified Euler's method (C++ code)

    Week 12:

  • Monday (11/11): Sec. 5.6
    HW: Section 5.6: 2(a,b) (For 2(a,b), use two-step and four-step Adams-Bashforth explicit methods respectively. You do not need to compute solutions over the entire time interval. The two-step explicit method computes w2, w3, and w4; The four-step explicit method computes w4, w5 and w6. You can use the Matlab code for the 4th order RK method to compute the starting values for each method respectively),
    4(Do exercise 1(a) and only computes w2 and w3), 13, 14.

  • Section 5.6 Multistep Methods


    The Adams-Bashforth four step explicit method Matlab code .
    The Adams fourth-order predictor-corrector method Matlab code .
    The Adams-Bashforth four step explicit method C++ code .
    The 4th order predictor-Correction method C++ code .
  • Wednesday (11/13): Sec. 5.6, 5.10
    HW: Section 5.6: 6(solve problem 2(b) by the 4th order predictor-corrector method, only compute approximations w4, w5 and w6. You can use the Matlab code for RK4 to get the initial approximations). Section 5.10: 1.
    Computer project:
    Use Alg. 5.4 (predictor-corrector method) to compute approximate solution to Extercise 5.6.3(b). Compare with exact solution. Submit the code in acms40390/hw6 with code name ex6-1.cpp.

  • 5.10 Stability


    C++ code of unstable 2-step method for solving y' =0, y(0) = 9.4 .
    C++ code of AB 4-step method for solving y' =0, y(0) = 9.4 .
  • Friday (11/15): Sec. 5.10, 5.11
    HW: Section 5.10: 4(d. Hint: for analyzing consisteny by local truncation error, do 3rd order Taylor expansion for y_{i+2} and y_{i+1} about y_i respectively. And see what you have left with for the difference equation whose approximations are replaced by exact solutions after some cancellation.), 7, 8.
    Section 5.11: 2(a,d, use the Matlab code to do the computation), 10, 12, 15.

  • 5.11 Stiff Equations

  • RK4 to solve a single stiff ODE (Matlab code)

  • RK4 to solve system of stiff ODEs (Matlab code)

    Week 13:

  • Monday (11/18): Sec. 5.11, 6.1, 6.2
    HW: Section 6.1: 3(a), 6(a,c), 9.

  • Section 6.1 & Section 6.2


    Algorithm for Gaussian elimination with backward substitution method.
    Gaussian elimination with backward substitution for solving linear systems of eqns code .
    Gaussian elimination with partial pivoting C++ code .

  • Section 6.3 - Section 6.5

  • Wednesday (11/20): Sec. 6.2, 6.3
    HW: Section 6.2: 4(do problems 2(a,c)), 14(do problem 10(a)), 15(do problem 9(b)), 20 (do problems 10(a,b)), 31. Section 6.3: 6(a,c).

    Computer project:
    Use Alg. 6.3 (Gaussian elimination with scaled partial pivoting) to solve exercise 10.c of Section 6.2 (Hint: Modify the code for Gaussian elimination with partial pivoting). Submit the code in acms40390/hw7 with ex7-1.cpp.

    Algorithms for LU, LL^t, LDL^t and Crout method for tridiagonal matrices


    LU factorization C++ code
    LDL^t factorization C++ code
    LL^t Cholesky factorization C++ code
    Crout factorization for tridiagonal linear system C++ code
  • Friday (11/22): Sec. 6.4, 6.5, 6.6
    HW: Section 6.4: 6. Section 6.5: 4(b,c), 9(b,c). Section 6.6: 2(a,b,c).

  • Section 6.6 Special Types of Matrices

    Week 14:

  • Monday (11/25): Sec. 6.6, 7.3.
    HW: Section 7.1: 1, 3(a,b), 5(b). Section 7.3: 2(b,c).
    Computer project:
    Use Alg. 7.1 (Jacobi iterative method) to solve Exercise 7.3.2(d) with TOL = 10^{-6}. Use the stopping criterion ||x^k-x^{k-1}||/||x^k|| < TOL in L_2 norm to stop iterations. Submit the code in acms40390/hw8 with ex8-1.cpp.

  • Sections 7.1 & 7.2

  • Section 7.3 Jacobi and Gauss-Siedel Methods



    Jacobi iteration method C++ code .
    Gauss Seidel iteration C++ code .

    Week 15:

  • Monday (12/02): Sec. 7.3, 7.4
    HW: Section 7.3: 4(do exercise 2(b,c)), 9(a,c).
    Computer project:
    Use Alg. 7.2 (Gauss-Seidel iterative method) to solve Exercise 7.3.2(d) with TOL = 10^{-3} in l_2 norm. Submit the code in acms40390/hw8 with ex8-2.cpp.
  • Wednesday (12/04): Sec. 7.4, 7.5
    HW: Section 7.4: 2(b,c), 10(b, Modify the sample code to solve. Use the stopping criterion given in Step 4 of SOR Alg. 7.3). Section 7.5: 4(a,c), 7.

  • Section 7.4 SOR Methods & 7.5 Error Bounds


    SOR C++ code
  • Friday (12/06): Sec. 8.1, 8.2
    HW: Section 8.1 (Write matlab codes to solve these problems. Also turn in the codes with homework): 5(b,c), 9.
    Section 8.2: 4(Use Legendre polynomials to construct least square approximation of degree 2 for function in Exercise2(c)), 14.

  • Section 8.1 Discrete least squares approximation

    Matlab code of polynomial least squares fitting .

  • Section 8.2 Orthogonal Polynomials and least squares approximation

  • Section 8.3 Chebyshev Polynomials

    Week 16:

  • Monday (12/09): Sec. 11.1, 11.2.
    HW: Section 11.1: 2, 6. Section 11.2: 3(a,b), 5(a).

  • Section 11.1-11.4 Methods for solving boundary value problem (BVP)


    Matlab code for solving linear BVP
    Matlab code for solving nonlinear BVP