AMSC 460, Spring 2006: Computational Methods

News

Syllabus

includes information about textbook, office hours, grading policy

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 (note that you have to put the files where Matlab can find them)
• Solving linear systems using Gaussian elimination (PDF)
• Errors for linear systems (PDF)
• How to solve a linear system in Matlab: x = A\b or
  [L1,U] = lu(A); % compute factorization A=L₁ U using Gaussian elimination with pivoting, choosing pivot candidates with largest absolute value
  y = L1\b;   % "forward substitution"
  x = U\y;    % back susbstitution
  condition = condestdec(L1,U) % estimate condition number of A (for 1-norm)
• vector and matrix norms in Matlab: norm(x,1) , norm(x,2) , norm(x,inf)
  condition numbers in Matlab: cond(A,1), cond(A,2), cond(A,3) (in all cases inv(A) is computed)
  estimate for cond₁(A): condest(A) (actually computes LU decomposition of A)
  If decomposition A=L₁ U is already computed, use condestdec(L1,U)
  (download condestdec.m , for Matlab 5 use condestdec5.m instead)
• Interactive demonstration of interpolation with polynomials and splines
• Polynomial interpolation:
  The 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.

  Then

  |(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
• How to solve linear least squares problem ||Ac-y|| = min in Matlab:
  using normal equations:

  B = A'*A; z = A'*y; c = B\z

  using the QR decomposition ("economy version"):

  [Q0,R0] = qr(A,0)
b0 = Q0'*y
c = R0\b0

  Matlab shortcut (actually uses QR decomposition):

  c = A\y

• Solving a nonlinear equation with Matlab: xs = fzero('f',[a,b]) for function in m-file f.m
  Example with inline function: f=inline('sin(x)','x'); xs=fzero(f,[3,4])
• Gauss-Legendre integration: Download findgauss.m.
  [xg,wg] = findgauss(n) gives n Gauss nodes xg(1),...,xg(n) and weights wg(1),...,wg(n) . With tg = (a+b)/2 + xg*(b-a)/2 the integral of f(x) from a to b can be approximated by
  Q = (b-a)/2*(wg(1)*f(tg(1)) + ... + wg(n)*f(tg(n))
  The routine findgauss.m implements the method for finding the Gauss weights and nodes as explained in class. However, a more efficient and numerically stable way to compute this is given by gauleg.m.
• Gauss-Jacobi integration: Download gaussjac.m.
  [xg,wg] = gaussjac(n,alpha,beta) gives n Gauss nodes xg(1),...,xg(n) and weights wg(1),...,wg(n) . With tg = (a+b)/2 + xg*(b-a)/2 the integral of f(x)*(x-a)^alpha*(b-x)^beta from a to b can be approximated by
  Q = ((b-a)/2)^(1+alpha+beta)*(wg(1)*f(tg(1)) + ... + wg(n)*f(tg(n)).
  Example: gaussj_ex.m
• Numerical integration with Matlab: Q = quad(f,a,b,tol) and Q = quadl(f,a,b,tol) (high order method). [Q,N] = quadl(f,a,b,tol) also returns the number of function evaluations.
  f may be string f='x.^2', inline function f=inline('x.^2','x'), or anonymous function f=@(x) x.^2 . If the function is in an m-file f.m use quad('f',...) or quad(@f,...). Note: You MUST use .* ./ .^ instead of * / ^ in the definition of f.
• Euler method in Matlab: Euler.m
• Solving ODEs with ode45
• Specifying the tolerance for ode45: Default for RelTol is 1e-3, default for AbsTol is 1e-6

  opts = odeset('RelTol',1e-4,'AbsTol',1e-7);
[t,y] = ode45('f',[t0,T],y0,opts);

• Matlab documentation for ODE solvers ode*, error tolerance, event location
• Matlab ODE solvers:
  ode45: try this first (uses methods of order 4,5)
  ode23: sometimes better for lower accuracy (uses methods of order 2,3)
  ode113: sometimes better for higher accuracy (uses variable order)
  ode15s: for "stiff" problems, high accuracy (uses variable order)
  ode23s: for "stiff" problems, low accuracy (uses order 2)

How to Hand in Matlab Homework

You will only receive credit for your homework if you follow the following rules:

• For each problem, all the commands which produce the numerical results and graphics have to be put in an m-file. Printouts of interactive input at the Matlab prompt will not be accepted!
• All output and graphics must be clearly labeled (see below for instructions)
• Print out only those values which were asked for in the problem!

For each problem you must hand in

1. a printout of the m-file (and all m-files called from it) containing the Matlab commands.
   • Include comments to explain what you are doing.
   • Print out only values asked for in the problem. Use ``;'' at the end of lines to suppress the output of intermediate results (especially if these are large arrays).
   • Clearly label all output and graphics. Use the disp or fprintf commands to for text output, see examples. Use the title, xlabel, ylabel, legend commands to label your graphs, see example.
2. a printout of the exact output of your m-file (where all results are properly labeled).
3. a printout of all graphics (properly labeled) (The command print -dps figXY.ps in Matlab produces a postscript file figXY.ps of the current figure which you can then send to the printer.)
4. separate text pages on which all the questions asked in the problem are answered (can be handwritten).

If you use Matlab 7 you can also use the publish command to generate an html file and print this from your web browser. This printout replaces 1.-3. from above.

Matlab Information

• Matlab is available in many computer labs on campus:
  • On Windows machines select Matlab from the Start menu to start the integrated enviroment (with editor and debugger).
  • On WAM or Glue Unix machines type tap matlab and then matlab & . This starts the integrated enviroment (with editor and debugger).
    In order to start Matlab in a terminal window type matlab -nodesktop. In this way you can also use Matlab remotely (but you won't see any graphs unless you have an X-server running on your local machine). You can use your favorite editor (nedit, pico, vi, emacs) for writing m-files.
• For an introduction to Matlab read the Matlab Primer (PDF file).
  Another Introduction to Matlab.
• Common Problems with Matlab
• Matlab documentation
  Matlab command:  

Using Computers on Campus

• Computer labs where Matlab is available
• Currently available WAM/Glue Unix machines
• How to get a WAM/Glue account, Printing account
• Text editors available on WAM, Glue (for writing m-files): nedit, pico, vi, emacs
• How to use a PC, Macintosh computer, Unix workstation in the campus labs
• How to print text files for free