ACMS 40390: Fall 2014 - Numerical Analysis

Course Information

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

Office phone number: (574)631-3423

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

Classroom: 229 Hayes Healy Center

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

Office hours:

Xu: M 11:30am - 1:30pm or by appointment and drop-in (when available), 242HH
Jiang: T 3:00pm - 5:00pm, 154 Hurley Hall

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/22 in class
First in class test Wednesday, 09/24 in class
Second test review Monday, 10/27 in class
Second in class test Wednesday, 10/29 in class
Final test review 5:00pm - 6:30pm, Wednesday, 12/17/2014 TBA
Final exam 8:00am - 10:00am, Friday, 12/19/2014 229 Hayes-Healy

Tutorial Session
Computer Lab 1 Thursday 08/28, 6:00pm - 7:00pm 129 Hayes Healy

Week 1:

  • Wednesday (08/27): 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

    Computer Lab (Bring your own laptop): When: 6:00pm - 7:00pm Thursday (08/28). Where: 129 Hayes Healy.
    Increase AFS space quota: email oithelp@nd.edu to ask for 200 MB AFS disk space if you are NOT engineering students.
  • Friday (08/29): Section 1.1 HW: Section 1.1: 1(b,c), 2(b), 4(a, b), 9, 19, 22(b).

    Section 1.2 Arithmetic

    Matlab code for finding extreme values of functions .

    Week 2:

  • Monday (09/01): Section 1.2 HW: Section 1.2: 15(a,c), 16(use the number given in 15(a,c)).
  • Wednesday (09/03): Section 1.2 HW: Section 1.2: 1(g), 3(b), 4(c), 12 (a,b,d (only find relative error of (b))), 19(a), 22 .
  • Friday (09/05): Section 1.2 HW: Section 1.2: 26(a,b(i)), 28.

    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 3:

  • Monday (09/08): Section 1.3 HW: Section 1.3: 6(b,d), 7(a, d), 14, 15.

    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 ...

  • Wednesday (09/10): Section 2.1 HW: Section 2.1: 2(a), 10, 14 (Do not need to find the approximation to the root), 16.

    Section 2.1 The Bisection Method

    Matlab code for bisection method .
    The sample code for the bisection method is ALG021_bisection.cpp .
  • Friday (09/12): Sections 2.1, 2.2 HW: Section 2.2: 1.
    Computer project (due on 09/19):
    Use Alg. 2.1 (the bisection method) to solve Exercise 2.1.20.
    Submit the code under folder acms40390/hw1 with file name ex1-1.cpp.

    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/15): Section 2.2 HW: Section 2.2: 2(do 1(a)), 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.), 16.
  • Wednesday (09/17): Sections 2.2, 2.3 HW: Section 2.2: 6(construct the fixed-point function and show that it generates a convergent sequence of approximations, namely, verify against Thm 2.4.). Section 2.3: 6(a)(Let P0= 1.5. Use Newton's method to find P3), 13 (Derive the equation to be solved. Then Use P0 = 1 and find P3).
    Computer project (due on 09/26):
    Use Alg. 2.2 (fixed point iteration) to solve Exercise 2.2.23.
    Submit the code in acms40390hw/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 The code is set to solve cos(x)-x=0

  • Algorithm of the Secant method.

    Matlab demo code for Secant method solving exercise2.3.6.a

    The sample code for the Secant method is

    ALG024_Secant.cpp The code is set to solve cos(x)-x=0


  • Friday (09/19): Sections 2.3, 2.4
    HW: Section 2.3: 4(a). Section 2.4: 6, 8(a)

    Computer project (Due: Friday Sep 26, 2014):
    Use Newton's method and Secant method respectively to solve Exercise 2.3.24. Use |P_{n}-P_{n-1}|/|p_n|<10^-4 as the stopping criterion.
    Submit the code in acms40390hw/hw2 with ex2-2.cpp (for Newton's method) and ex2-3.cpp (for Secant method) respectively.

  • Section 2.4 Error Analysis for Iterative Methods

    modified_Newton_241.cpp The code is set to solve e^x-x-1.0=0


    Week 5:

  • Monday (09/22): Review
  • Wednesday (09/24): Exam
  • Friday (09/26): Sections 2.4, 2.5
    HW: Section 2.4: 4(use modified Newton's method to solve exercise 2.a to obtain p_2)
    Section 2.5: 13(a) (generate elements till p_3 hat).

  • Section 2.5 Accelerating Convergence

    ALG026_ALG026_Steffensen.cpp

    Week 6:

  • Monday (09/29): Sections 2.5, 2.6
    HW: Section 2.5: 4, 6.
    Section 2.6: 2(a), Let x0=1, compute f(x0) and f ' (x0) (derivative of f at x0) by Horner's method.

  • 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
  • Wednesday (10/01): Sections 2.6, 3.1 HW: Section 3.1: 2(c), 14.

  • 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.
  • Friday (10/03): Sections 3.1 HW: Section 3.1: 6(a, justify how you choose nodes), 8(a), 10.

    Week 7:

  • Monday (10/06): Sections 3.3 HW: Section 3.1: 17. Section 3.3: 2, 16.

    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/08): Sections 3.3, 3.4 HW: Section 3.3: 10, 14, 17.
    Computer project (Due: Friday Oct 17, 2014):
    Use Alg. 3.2 (Newton's Divided-Difference) to solve Exercise (8) of Section 3.3.
    Submit the code in acms40390hw/hw3 with ex3-1.cpp.
    Hint: Use the above C code of Newton divided differences to start.
  • Friday (10/10): Section 3.4 HW: Section 3.4: 2(b), 4(b), 5(a, using the divided difference approach, b).

    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.

    Week 8:

  • Monday (10/13): Sections 3.4, 3.5
    Computer project (Due: Monday Oct 27, 2014):
    Use Alg. 3.3 (Divided-Difference) to solve Exercise 10 of Section 3.4.
    Submit the code in acms40390hw/hw3 with ex3-2.cpp.
    Hint: Use the above C code of Hermite Interpolation to start.


  • Wednesday (10/15): Section 3.5 HW: Section 3.5: 2, 11, 12.
    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.

  • Friday (10/17): Section 4.1 HW: Section 4.1: 2(a), 4(a), 6(c), 8(c), 10, 13.

    Section 4.1 Numerical Differentiations



    Week 10:

  • Monday (10/27): Exam review
  • Wednesday (10/29): Exam
  • Friday (10/31): Sections 4.1, 4.3 HW: Section 4.1: 20. Section 4.3: 2(a), 4 (do 2(a)), 6(do 2(b)), 8(do 2(b)), 15, 18.

    Section 4.3 Elements of Numerical Integration

    Week 11:

  • Monday (11/03): Sections 4.3, 4.4 HW: Section 4.3: 10 (do 2(c)), 12 (do 2(c)). Section 4.4: 2(a), 4(do 2(b)), 6(do 2(c)), 9, 12, 22 (use composite Simpson's rule)

    Section 4.4 Composite Numerical Integration


  • Wednesday (11/05): Sections 4.7 HW: Section 4.7: 1(a), 2(do 1(b)), 5, 7(do n =3 case).

    Section 4.7 Gaussian Quadrature



  • Friday (11/07): Section 5.1 HW: Section 5.1: 2(a,b), 8(b).

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

  • Euler's Method Matlab code .


    Week 12:


  • Monday (11/10): Section 5.2 HW: Section 5.2: 2(a), 4 (a), 10 (b(i), c), 12.

  • Wednesday (11/12): Section 5.3 HW: Section 5.3: 2(d), 9(d(i) only need to compute w(t=1.1) in order to solve this problem).

    Section 5.3 Taylor Methods

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


  • Friday (11/14): Section 5.4 HW: Section 5.4: 2(a), 6(do 2(b)), 14(do 2(b)), 31.

  • Section 5.4 Runge-Kutta Methods

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

    Modified Euler method (C++ code)


    Week 13:

  • Monday (11/17): Section 5.6 HW: Section 5.6: 2(c) (use two-step Adams-Bashforth explicit method to compute w2, w3), 2(d) (use four-step Adams-Bashforth explicit method to compute w4, w5), 4(do 1(a), use 3-step Adams-Moulton implicit method (formula (5.38)) to compute w3, w4)
    Note: You can use the Matlab code of the 4th order RK method to compute the starting values for each of these problems respectively.

    Computer project (Due on Nov/24):
    Implement a code to solve problem 5.4.28(a) by the Runge-Kutta method of order 4. Submit the code in acms40390hw/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. The modified Euler's method C++ code is at http://www3.nd.edu/~zxu2/acms40390hw/ALG_MEuler.cpp)

  • 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/19): Section 5.10 HW: Section 5.6: 7(do 3(a), use 4th-order predictor-corrector to compute w4, w5). Section 5.10: 3 (coefficient values are in the email attachment).
    Note: You can use the Matlab code of the 4th order RK method to compute the starting values for each of these problems respectively.

    Computer project (Due on Nov/26):
    Use Alg. 5.4 (4th order predictor-corrector method) to solve Extercise 5.9. Use time step size h = 0.01. Compare with exact solution. Submit the code in acms40390/hw5 with code name ex5-2.cpp. (Hint: the base code is http://www3.nd.edu/~zxu2/acms40390hw/predict_correct.cpp)
  • Friday (11/21): Section 5.10
    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. In the difference equation, replace the approximate solution by exact values and plug these Taylor expansions into the equation. See what you have after some cancellation);
    8(instead of computing w_i to see how round-off error propagates, find the exact solution to w_i when this method is applied to solve the ODE as the example in the lecture notes does. Then tell how the error is amplified.)

  • 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 .

    Week 14:

  • Monday (11/24): Section 5.11 HW: Section 5.11: 2(a), 10, 12, 15(b).

  • 5.11 Stiff Equations

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

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


  • Stability region

    Week 15:

  • Monday (12/01): Sections 6.1 & 6.2
    HW: Section 6.2: 4(do 2(c)), 16(do 10(b)), 20 (do 10(c)), 31.


    Computer project (due on 12/08):
    Use Alg. 6.3 (Gaussian elimination with scaled partial pivoting) to solve exercise 10.d of Section 6.2 (Hint: Modify the code for Gaussian elimination with partial pivoting). Submit the code in acms40390hw/hw6 with ex6-1.cpp.

  • 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.5


  • Wednesday (12/03): Sections 7.3
    HW: Section 7.3: 2(b), 4(do 2(c)), 9(a,c).

    Computer project (due on 12/12):
    Use Alg. 7.2 (Gauss-Seidel) to solve exercise 2.d of Section 7.3 with TOL = 10^{-3} in l_2 norm. (Hint: Modify the given Gauss Seidel iteration code). Submit the code in acms40390hw/hw7 with ex7-1.cpp.

  • Section 7.3


    Jacobi iteration method C++ code .
    Gauss Seidel iteration C++ code .
  • Friday (12/05): Sections 7.4
    HW: Section 7.3: 14.
    Section 7.4: 2(b,c).

  • Section 7.4 SOR Methods & 7.5 Error Bounds


    SOR C++ code

    Week 16:

  • Monday (12/08): Sections 7.5, 8.1
    HW: Section 7.5: 4(a), 7.
    Section 8.1: 4.

  • Section 8.1 Discrete least squares approximation

    Matlab code of polynomial least squares fitting .