AMSC/CMSC 466, Fall 2012: Introduction to Numerical Analysis

News

• Final Exam: Friday, Dec. 14, 8:00-10:00am in MTH0409
  You can bring a cheat sheet (letter size, use front & back, must be handwritten)
• Notes about Numerical Integration
  For the final exam you should know the boxed formulas.
  Practice Problems for Numerical Integration
  Solutions
• Practice problems for exam 2
  Solutions
• Exam 2 will be on Thursday, Dec. 6
• Assignment 4, due date extended to Wednesday, Dec. 5, 3:00pm.
  Solution of Assignment 4
  Problem 2(b): How to plot contours where two functions are zero:
  Example: To plot the contours where x1 + x2 -x1*x2 - 2 = 0 and x1 - x2 + x1^3 + 3 = 0 in the rectangle [-4,4]×[-5,5] use
  

    [X1,X2] = meshgrid(linspace(-4,4,100),linspace(-5,5,100));
    contour(X1,X2, X1 + X2 - X1.*X2 - 2,[0 0],'b-');  hold on 
    contour(X1,X2, X1 - X2 + X1.^3 + 3,[0 0],'r:'); hold off 
  

  Note that you must use .* ./ .^ instead of * / ^.
  Problem 3: Read the notes about the fixed point iteration below, please download the latest version.
  Problem 3(c): Use the a-priori estimate (equation (2) in the notes).
• Assignment 3 was handed out on Nov.1, due date extended to Nov. 13.
  Solution of Assignment 3
  Problem 2(b): download membr.m. You have to modify a_diag, a_left, a_right, a_down, a_up.
  Measuring computation time in Matlab: Use tic; (your Matlab code goes here) time=toc;
  Problem 2(c): download spmembr.m. You have to modify a_diag, a_left, a_right, a_down, a_up.
  Problem 3: Note that the error bounds obtained from the error formula may be much larger than the actual error.
  Problem 4(a): See the divided difference algorithm below. The m-file evnewt.m should also work if xt is a vector: Initialize v as a vector: v = d(1)*ones(size(xt)); and use ".*" instead of "*".
  Problem 4(b): For (C) use hermite.m, see the material about piecewise cubic interpolation below and the example code.
• Practice problems for Exam 1
  solutions
• Exam 1 was on Tuesday, Oct. 23. It will cover error propagation, roundoff error and linear systems.
• Assignment 2, due date Tuesday, Oct. 9
  Solution
• Assignment 1, due date extended to Thursday, September 13
  Solution
  Typo: In 4(b) it should say: "Check in Matlab that fl(2+yhat)>2" (instead of 4).
• Instructions for Assignments
  • Handing in homework
    When you slide your homework under my office door you MUST write the current time and date on the homework. Otherwise it will count as one additional day late. Writing an incorrect time/date will be considered academic dishonesty.
  • Matlab
    All computations have to be done in Matlab. Matlab homework will not be accepted unless it is in the following format:
    • You must write the commands for each problem into an m-file (e.g. problem2a.m). Just entering a sequence of commands at the Matlab prompt (>>...) and printing this out is not allowed.
    • For each problem you hand in (1) all m-files, (2) the output from running the m-files, (3)  all graphs produced by the m-file, (4) additional text (can be handwritten) which answers all questions asked in the problem.
      If you want you can use "publish" in Matlab for printing (1)-(3).
    • All output (printed numbers and graphs) must be clearly labeled (example for numbers, example for graphs).
    • Only print out values explicitely asked for in the problem! Make sure that you have ";" at the end of each line to suppress printing of intermediate output.
• If you have not used Matlab before or would like to review the basics:
  Read the Matlab Primer (PDF file) and (important!) try out the commands on a computer.

Additional Course Material (Required Reading)

• Disasters caused by careless numerical computation
• Introduction: Errors (PDF)
• Error propagation, forward error analysis, numerical stability (PDF)
• How to approximate a real number by a machine number.
• IEEE 754 machine numbers and machine arithmetic
• IEEE 754 double precision machine numbers in Matlab
  Download the files num2bin.m and bin2num.m (put them in the same directory/folder as your other m-files.)
• Solving Linear Systems
• Errors for Linear Systems
• How to solve a linear system in Matlab:
  • [L,U,p] = lu(A,'vector'); % Gaussian elimination choosing the pivot candidate with largest absolute value
    y = L\b(p)                % forward substitution with permuted vector b
    x = U\y                   % back substitution
  • or use
    x = A\b
    which does the above steps automatically (but without the option to save the result of the LU-decomposition).
• Estimating the condition number if L,U are known: Download invnormest.m.
  For the 1-norm condition number use

  c = norm(A,1)*invnormest(L,U)

  For the infinity-norm condition number use

  c = norm(A,Inf)*invnormest(U',L')

  This works also for a sparse matrix A (see right below).
  The basic algorithm for the norm estimation is invnormest_simple.m
• How to solve sparse linear systems in Matlab:
  A sparse matrix is a matrix where most elements are zero. In this case it is much more efficient to use the special sparse data structures in Matlab. All operations like *, \, lu have special efficient algorithms for matrices in sparse format. In particular the backslash command x=A\b uses clever row and column permutations to minimize "fill-in".
  • initializing a sparse matrix: A sparse matrix can be described by giving a list of the nonzero elements, as vectors iv (contains i-values), jv (contains j-values), av (contains matrix entries). Then a sparse matrix is initialized by

    A = sparse(iv,jv,av);

    Example: The full matrix given given by A = [0 7 0; 0 0 8; 9 0 0] can be defined as a sparse matrix by

    A = sparse([1 2 3],[2 3 1],[7 8 9])

  • Solving a linear system with a sparse matrix: use row and column permutations to minimize "fill-in":
    • [L,U,p,q] = lu(A,'vector'); % Gaussian elimination with row and column permutations so that L*U = A(p,q)
      y = L\b(p)                  % forward substitution with permuted vector b
      x(q,1) = U\y                % back substitution, yielding permuted vector x
    • or use
      x = A\b
      which does the above steps automatically (but without the option to save the result of the LU-decomposition).
• linear systems with symmetric positive definite matrix:
  • C = chol(A);
    y = C'\b; x = C\y;
  • or use
    x = A\b
    which does the above steps automatically (but without the option to save the result of the LU-decomposition).
• Interactive demonstration of interpolation with polynomials and splines
• Divided difference algorithm: Given are the distinct nodes x₁,...,x_n and the function values y₁,...,y_n.

  divdiff computes the divided differences d₁=f[x₁,...,x_n], d₂=f[x₂,...,x_n],...,d_n=f[x_n]
  Algorithm d=divdiff(x,y)
  d := y
  For k=1 to n-1
     For i=1 to n-k
        d_i := (d_i+1-d_i)/(x_i+k-x_i)

  evnewt evaluates the interpolating polynomial at the point t
  Algorithm v=evnewt(d,x,t)
  v := d₁
  For i=2 to n
     v := v·(t-x_i) + d_i

• Polynomial interpolation:
  Equidistant nodes on the interval [a,b] are given by

  h = (b-a)/(n-1);    x_j = a + (j-1)h for j = 1,...,n

  In Matlab you can use x = linspace(a,b,n) to define a vector of equidistant nodes.
  The node polynomial has the bound

  |(x-x₁)...(x-x_n)| ≤ hⁿ(n-1)!/4 for x in [a,b]

  Chebyshev nodes on the interval [a,b] are given by

  x_j = (a+b)/2 + cos(π(j-1/2)/n) (b-a)/2 for j = 1,...,n.

  The node polynomial has the bound

  |(x-x₁)...(x-x_n)| ≤ 2 ((b-a)/4)ⁿ for x in [a,b]

• Piecewise cubic interpolation in Matlab: The column vector x contains the x-coordinates of the nodes, the column vector y contains the function values at the nodes, the vector yp contains the derivatives at the nodes. The vector xi contains the points where we want to evaluate the interpolating function, the vector yi contains the corresponding values of the interpolating function.
  • Cubic Hermite interpolation: Download hermite.m . Then use yi = hermite(x,y,yp,xi)
  • Complete cubic spline: points are given by column vectors x, y. Slopes at left and right endpoint are sl and sr.
    yi = spline(x,[sl;y;sr],xi)
  • Not-a-knot cubic spline: yi = spline(x,y,xi)
  Example: example for Hermite, complete spline, not-a-knot spline
  In all three cases these we can first compute the spline as pp (a special structure for representing piecewise polynomials), and later evaluate the spline: E.g., for the not-a-knot spline:

  pp = spline(x,y);
  yi = ppval(pp,xi);

• Fixed Point Iteration and Contraction Mapping Theorem (latest version: Dec. 3, 6pm)
• Solving a nonlinear equation with Matlab:
  We want to find a zero of the function f in the interval [a,b] where f(a)f(b)<0:

  xs = fzero(@f,[a,b])

  where the function f is given by the m-file f.m.
  Example with @-function: f=@(x)sin(x); xs=fzero(f,[3,4])
• Newton method for systems:
  given: function y=f(x) (where x,y ∈ Rⁿ)
              initial guess x⁽⁰⁾
  first find Jacobian f′(x) (n×n matrix of 1st order partial derivatives)
  for i = 0,1,2,3,...
  1. evaluate b := f(x⁽ⁱ⁾)
  2. evaluate A := f′(x⁽ⁱ⁾)
  3. solve linear system A d = -b
  4. x⁽ⁱ⁺¹⁾ := x⁽ⁱ⁾ + d

  Matlab code for example from class

• Solving a nonlinear system in Matlab using fsolve:
  • Without using the Jacobian:
    Write an m-file f.m containing the function f:

    function y = f(x)
    y = ...

    Then use x = fsolve(@f,x0) where x0 is the initial guess.
  • With using the Jacobian:
    Write an m-file f.m containing the function and the Jacobian:

    function [y,A] = f(x)
    y = ...
    A = ...

    Then use

    opt = optimset('Jacobian','on');
    x = fsolve(@f,x0,opt);

    where x0 is the initial guess.
  • You can specify additional options using optimset, e.g. the tolerance for the size of the function value and the tolerance for the change in x:

    opt = optimset('TolFun',1e-5,'TolX',1e-5);
    x = fsolve(@f,x0,opt);

    There are many other options, type doc fsolve for more information.

Matlab Information

• How to hand in Matlab homework
• Matlab is available in many computer labs on campus
• For an introduction to Matlab read the Matlab Primer (PDF file).
  Another Introduction to Matlab: Part I, Part II
• Common Problems with Matlab
• Matlab documentation
  Matlab command: