ACMS 40390: Fall 2015 - Numerical Analysis

Course Information

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

Office phone number: (574)631-3423

Class time: M, W, F, 11:30am - 12:20pm.

Classroom: 117 Haggar Hall

Teaching assistant: Michael Machen, mmachen@nd.edu, 215 Hayes-Healy

Office hours:

Xu: T 3:00pm - 5:00pm or by appointment and drop-in (when available), 242 Hayes-Healy
Michael Machen: T 3:30pm - 4:30pm, 229 Hayes-Healy

Textbook: Richard L. Burden and J. Douglas Faires, Numerical Analysis, 10th 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/21 in class
First in class test Wednesday, 09/23 in class
Second test review Monday, 10/26 in class
Second in class test Wednesday, 10/28 in class
Final test reviews Saturday (12/12), Sunday (12/13): both at 3:30pm - 5:00pm 127 Hayes-Healy
Final exam Monday, 12/14. 4:15pm - 6:15pm 117 Haggar Hall

Tutorial Session
Computer Lab 1 Time: TBA Place: TBA

Week 1:

  • Wednesday (08/26): Sections 1.1 HW:

    Section 1.1 Review of Calculus

    Computer Lab (Bring your own laptop): When: TBA. Where: TBA.
    Increase AFS space quota: email oithelp@nd.edu to ask for 200 MB AFS disk space if you are NOT engineering students.
  • Friday (08/28): Section 1.1 HW: Section 1.1: 2(c), 4(a), 6(a), 12
    Matlab code for finding extreme values of functions .

    Week 2:

  • Monday (08/31): Sections 1.1, 1.2 HW: Section 1.1: 19, 25(b). Section 1.2: 19(b).

    Section 1.2 Arithmetic


    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)


  • Wednesday (09/02): Section 1.2 HW: Section 1.2: 3(b), 4(a), 6(c), 14(a, b, d(only find relative error for value obtained in (b))).
  • Friday (09/04): Section 1.2 HW: Section 1.2: 16(a), 21, 26, 29.

    Week 3:

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

  • Monday (09/07): Sections 1.2, 1.3
    HW: Section 1.2: 25; Section 1.3: 6(a,c), 8, discussion question 6(a,b)
  • Wednesday (09/09): Section 2.1
    HW: Section 1.3: 8; Section 2.1: 1,3,11.

    Section 2.1 The Bisection Method

    Matlab code for bisection method .
    The sample code for the bisection method is ALG021_bisection.cpp .
  • Friday (09/11): Section 2.2
    HW: Section 2.2: 1(a,b),
    4(a,b) (Only derive the fixed point method, i.e., find necessary algebraic manipulations for transforming f(x)=0 into x = g(x)),
    10(Use Thm 2.3 to show that g(x) has a unique fixed-point on [1/3, 1]) .

    Computer project (due on 09/18):
    Use Alg. 2.1 (the bisection method) to solve Exercise 2.1.16.
    Submit the code under folder acms40390hw/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/14): Sections 2.2, 2.3
    HW: Section 2.2:
    2(a) (Do functions g1 and g2 in exercise 1),
    10(Use Corollary 2.5 to estimate the number of iterations required to achieve 10^-4 accuracy.),
    20, 23(a).

    Computer project (due on 09/25):
    Use Alg. 2.2 (the fixed-point iteration method) to solve Exercise 2.2.8.
    Submit the code under folder acms40390hw/hw2 with file name 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


  • Wednesday (09/16): Sections 2.3, 2.4
    HW:
    Section 2.3: 2, 4(a), 19 (Only need to derive the equation of iteration of Newton's method).
    Section 2.4: 6, 8(a).

    Computer project (due on 09/30):
    Use Alg. 2.3 (Newton's method) and Alg. 2.4 (Secant method) to solve Exercise 2.3.22, respectively.
    (Hint: Modify the demo C++ codes for these two methods. The demo codes are available on the class web site. )
    Use |pn -pn-1|/|pn|<10-4 as the criterion to stop iteration.
    Submit the code under folder acms40390hw/hw3 with file name ex3-1.cpp (for Newton's method) and ex3-2.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


  • Friday (09/18): Section 2.4
    HW:
    Section 2.4: 4(Do problem 2(a), use p0=0.5 and compute till p4), 10.

    Week 5:

  • Monday (09/21): Review
  • Wednesday (09/23): Exam
  • Friday (09/25): Discussion
    HW: Section :

    Week 6:

  • Monday (09/28): Sections 2.5, 3.1.
    HW:
    Section 2.5: 4, 6.

  • Section 2.5 Accelerating Convergence

    ALG026_ALG026_Steffensen.cpp
  • Wednesday (09/30): Section 3.1.
    HW:
    Section 3.1: 2(c)(Construct interpolation polynomial of degree at most 2 to approximate f(1.4) and find the absolute error), 9.

  • Section 3.1 Interpolation and the Lagrange Polynomial


    Using Matlab to plot solution of exercise 2(a)
    syms x
    xi=[1.0, 1.25, 1.6]
    fi = sin(3.14159*xi)
    xx=1.0:0.02:1.6
    plot(xx,sin(3.14159*xx), '--g', xi, fi, 'or')
    hold on
    ezplot('sin(3.14159*1.25)*(x-1)*(x-1.6)*80/(-7)+sin(1.6*3.14159)*100/21*(x-1)*(x-1.25)',[1,1.6])


    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/02): Section 3.1.
    HW:
    Section 3.1: 5(a), 14(b), 17.

    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.

    Week 7:

  • Monday (10/05): Section 3.3.
    HW:
    Section 3.3: 2(a), 14.
  • Wednesday (10/07): Sections 3.3, 3.4.
    HW:
    Section 3.3: 11, 16, 20.

  • Friday (10/09): Section 3.4.
    HW:
    Section 3.4: 2(a Use both Theorem 3.9 and divided difference method to construct Hermite polynomial, respectively), 4(a), 7(Find error bound for H3(x) approximating f(x)).

    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/12): Section 3.5.
    HW:
    Section 3.5: 11, 13.

    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.

  • Wednesday (10/14): Section 3.5.
    HW:
    Section 3.5: 4(c), 6(c), 21(c).

    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

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

    Section 4.1 Numerical Differentiations



    Week 10:


  • Monday (10/26): Exam Review

  • Wednesday (10/28): Exam
  • Friday (10/30): Sections 4.1, 4.2
    HW:
    Section 4.1: 20, 26.
    Section 4.3: 2(b), 4(do exercise 2(b)).

    Section 4.3 Elements of Numerical Integration


    Week 11:


  • Monday (11/2): Section 4.3
    HW:
    Section 4.3: 6 (do exercise 2(b)), 8(do exercise 2(b)), 13, 21.
  • Wednesday (11/4): Sections 4.3, 4.4
    HW:
    Section 4.3: 16(c, d. Only use formula 4.32 to do numerical integration). Section 4.4: 2(a), 4(do exercise 2(b)), 6(do exercise 2(a)), 9, 11(b).

    Section 4.4 Composite Numerical Integration


  • Friday (11/6): Sections 4.4, 4.7
    HW:
    Section 4.7: 2(b), 4(do exercise 2b), 13( do n = 2 case).

    Section 4.7 Gaussian Quadrature


    Week 12:

  • Monday (11/9): Sections 5.1, 5.2
    HW:
    Section 5.1: 4, 5(a. Only show the given equation implicitly defines a solution).

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

  • Euler's Method Matlab code .
  • Euler's Method C++ code .

  • Wednesday (11/11): Sections 5.2, 5.3
    HW: Section 5.2: 2(d), 4(do exercise 2d), 9(a, b(ii)).

  • Friday (11/13): Sections 5.3, 5.4
    HW: Section 5.3: 2(a), 4(do exericse 2a), 9(b, d) (do (ii) for both. Use the fact that w5 = 3.910985, w6 = 5.643081).

    Section 5.3 Taylor Methods

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


    Week 13:

  • Monday (11/16): Section 5.4
    HW: Section 5.4: 2(a), 6(do ex 2(b)), 14(do ex 2(b)), 32.

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



  • Wednesday (11/18): Section 5.6
    HW:
    Section 5.6: 4(a. Use AB 2-step method to compute w2 and w3; Use AB 4-step method to compute w4 and w5, respectively. Use the exact solution to get necessary starting values for calculations. Compare results to exact solution, respectively.).
    5(do exercise 1(a). Use AM 2-step method to compute w2 and w3; Use AM 4-step method to compute w4 and w5, respectively. Use the exact solution to get necessary starting values for calculations. Compare results to exact solution, respectively.)


    Computer project (Due on Nov/30):
    Implement a C++ code to solve exercise 5.4.27 by the Runge-Kutta method of order 4. use time step h = 0.001. 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 .

  • Friday (11/20): Sections 5.6, 5.10
    HW:
    Section 5.6: 10 (Do exercise 4(d). Use Adams 4th order predictor-corrector method to compute solutions w4, w5. Use the exact solution to get necessary starting values for calculations. Compare results to exact solution, respectively.).
    Section 5.10: 3, 4(a,b. Only compute w2, w3).

    Week 14:


  • Monday (11/23): Sections 5.10
    HW:
    Section 5.10: 5(a) (Hint: for analyzing consisteny by local truncation error, do 3rd order Taylor expansion for y_{i-2}, y_{i-1} 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).

  • 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 15:


  • Monday (11/30): Sections 5.10, 5.11
    HW:
    Section 5.10: 8 (Instead of computing wi exactly for i = 2,...,6, compute the zeros of the associated characteristic polynomial. Then solve for the general solution to wi, and determine the stability of the method).
    Section 5.11: 10.

  • 5.11 Stiff Equations

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

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


  • Stability region
  • Wednesday (12/02): Sections 5.11, 6.1
    HW:
    Section 5.11: Discussion Question 1 (Hint: Use theorem 5.24).
    Section 6.1: 6(b).

  • Section 6.1 & Section 6.2


  • Friday (12/04): Sections 6.2
    HW:
    Section 6.2: 16 (do exercise 10.b), 20(do exercise 10.c).

  • Section 7.3 Jacobi and Gauss-Siedel methods


    Jacobi iteration method C++ code .
    Gauss Seidel iteration C++ code .
    Computer project test: 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 .