AMSC 460, Fall 2004: Computational Methods

News

Final Exam is on Saturday, Dec. 18, 8:00 in MTH 0405
- You can bring a calculator and a handwritten 8.5 by 11 sheet with formulas (front & back)
- This covers the whole semester, including numerical integration. Some topics which will not be on the final exam: splines, Gaussian quadrature. For ODEs you should know Euler, Improved Euler, backward Euler, but not RK4.
- You have to know how to use the Matlab commands "\" (for linear systems and least squares), lu, fzero, ode45, quad
- Old exam for practicing, solutions using Matlab
Practice problems for exam 2, Solutions for 1-3(c), Solution 3(d)
Practice problems for exam 1, Solutions
Projects:
- TA Stacey Nicholls, son@math.umd.edu
  office hours: Wed Fri 11-12 in CSS 4100 or by appointment (contact her by e-mail)
- Projects are due Fri Dec. 10 10:00pm.
  You should hand in m-files, output, graphs on paper, together with your observations for each part of the problem.
  You should also e-mail me m-files (as attachments) problem_a.m, problem_b.m etc. for each part of the problem, so that running e.g. problem_a generates all results, graphs, animations for part (a), properly labeled, without any user input.
- Projects:
  1. Sound waves in ocean
  2. Predator-prey models
    Investigating critical points of an autonomous system
  3. Spacecraft in earth-moon system
  4. Double Pendulum
    Derivation, Matlab function dpanim.m for animating results
  5. Cannon ball
    Correction: The ODE should be
    x₁'' = -v·(x₁' - w)
    x₂'' = -1 - v·x₂'
    where v = [(x₁' - w)² + (x₂')²]^1/2.
    Use w=-0.1 in (d) (not -0.3).
  6. Roller coaster
    coastertrack.m
  7. Putting a golf ball
    Hint for (a)
- General Matlab instructions:
  - Plotting curve using plot3 and comet3, plotting graph of function f(x,y) using surf
    First download nice3d.m
  - Evaluating splines and their derivatives with ppval, ppval1, ppval2
  - Solving one nonlinear equation with fzero
  - Minimizing a function of one variable with fminbnd
  - Solving ODEs with ode45
  - Adjusting the accuracy of Matlab ODE solvers
  - Stopping ode45 at an "Event"
  - Plotting the vector field for autonomous system of two first order ODEs
  - You can pass parameters from main program to function f using global: pendulum example
  - Animation using comet and comet3: You need output with uniformly spaced t-values by using e.g. [t,y]=ode45('f',0:.01:10,y0): pendulum example
This web page can also be accessed under
www.wam.umd.edu/~petersd/460
Use this if there is a problem with the normal address www.glue.umd.edu/~tvp/460
TEXTBOOK: "Numerical Computing with MATLAB" by Cleve Moler,
available
- free download as PDF files
- for $42.50 available at Campus Bookstore, Book Exchange soon
If you are not familiar with Matlab:
- Read the Matlab Primer (PDF file) and try out the commands on a computer.
- A book with a gentle introduction to Matlab is "Getting Started With Matlab" by R. Pratap, ISBN: 0195150147.

Syllabus

office hours of instructor and grader, 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)
Linear systems: Gaussian elimination without and with pivoting (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
Using ppval, ppval1, ppval2 for evaluating splines: First download ppval1.m, ppval2.m :
Example: Interpolating the sine function using a not-a-knot spline
x = linspace(0,2*pi,5); y = sin(x); % use 5 points pp = spline(x,y) % find coefficients pp of not-a-knot spline p(x) xe = linspace(0,2*pi,9); % evaluate at 9 points xe ye = ppval(pp,xe) % values of spline p(x) ye1 = ppval1(pp,xe) % values of p'(x) ye2 = ppval2(pp,xe) % values of p''(x)
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: 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])
Minimizing a function of one variable: xs = fminbnd('f',a,b) for function in m-file f.m
Example with inline function: f=inline('sin(x)','x'); xs=fminbnd(f,1,7)
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)
Stopping ode45 at an "Event": Example problem where we want to stop when y1 becomes 0:
version for Matlab 6 (note: you MUST use ode45(@f,...) even if f is in a separate file f.m), version for Matlab 5 using fev.m

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

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 uninteresting 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.
a printout of the exact output of your m-file (where all results are properly labeled).
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.)
separate text pages on which all the questions asked in the problem are answered (can be handwritten).

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 Macintosh OS X computers: Click on the face icon on the bottom left to open the Finder. Select "Applications" on the left in the Finder, then double-click on "Matlab" on the right in the Finder.
- On WAM or Glue Unix machines typetap matlaband then matlab & . This starts the integrated enviroment (with editor and debugger).
  In order to start Matlab in a terminal window typematlab -nojvm. 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).
Common Problems with Matlab
Matlab documentation
Matlab command:
Matlab Peer Training: If you are new to Matlab you might consider attending the class Introduction to Matlab (see schedule). This class costs $10. It is recommended to register early (WAM account required). Take a look at the class notes.

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