AMSC 460, Spring 2004: Computational Methods

News

• Final Exam: Monday, May 17, 10:30-12:30 in MTH 0101
  Bring a calculator. The exam is open book and open notes.
• Assignment 5:
  • Problem 2 and 5(b): Use fzero('G',[a,b]) with an appropriate initial interval [a,b].
  • Problem 2: Typo: It should be F(h) = /4, (not ²/4). The correct value for h should be between 0 and 1 (since h=1 means that the tank is half full).
    Write an m-file G.m so that y = G(h) computes
    y = (integral from -1 to -1 + h of ...) - /4
    using quad. Then use fzero for the function G. You must use quad and fzero in this way to get credit for this problem.
  • Problem 3(b): For the values at t=5: Print out the difference of Euler value and ode45 value, and the difference of RK2 value and ode45 value. When you use h=0.01 instead of h=0.02, how much smaller do you expect the error to be for (i) and (ii)? Compare the norms of the difference vectors.
  • Problem 4: Experiment with different values of v and plot y(t) vs. t for t in [0,40]. If v is small, y(t) converges to 0. For larger values of v, the pendulum will go over the top once before coming to rest at the bottom, and y(t) will converge to 2. Use trial and error to find the five different values for v.
  • Problem 5(a): Use y₁ := x₁, y₂ := x₁', y₃ := x₂, y₄ := x₂' .
  • Problem 5(b): How to tell ode45 to stop when a certain "event" happens: Example problem where we want to stop when y1 becomes 0:
    version for Matlab 6 , version for Matlab 5 using fev.m
    Use ode45 with odeset('Events',...) in your m-file distance.m.
    Use the m-file G.m containing
    function y = G(alpha)
    y = distance(alpha) - 0.5;
    and then use fzero for the function G.
• Euler.m: Euler method in Matlab (slightly different from version in class: first argument is now function f).
  Please try out the examples!
• How to solve ODEs in Matlab
• Matlab function Q = quad('f',a,b,tol): In the function f you must use .* ./ .^ instead of * / ^: E.g. for f(x) = 1/(1+x^2) use the following m-file f.m
  function y = f(x)
  y = 1./(1+x.^2);
• Gaussian 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 from a to b can be approximated by
  Q = (b-a)/2*(wg(1)*f(tg(1)) + ... + wg(n)*f(tg(n))
• Piecewise cubic interpolation in Matlab: Download hermite.m for Hermite interpolation, example for Hermite, complete spline, not-a-knot spline
• Faster Matlab code for membrane problem on assignment 2. This uses sparse data structures. Try it for N=40,80,160,320 (over 100,000 unknowns).
• Read Linear systems, Gaussian elimination without and with pivoting (PDF)
• If you are not familiar with Matlab:
  • Read the Matlab Primer (PDF file) and try out the commands on a computer.
  • Attend the Peer training 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.
  • A book with a gentle introduction to Matlab is "Getting Started With Matlab" by R. Pratap, ISBN: 0195150147.

Additional Course Material (Required Reading)

• Disasters caused by careless numerical computation
• 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)
• New: Linear systems, Gaussian elimination without and with pivoting (PDF)
• 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) (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)
• How to solve a linear system in Matlab:
  [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)
• Interactive demonstration of interpolation with polynomials and splines
• 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:  yi = spline(x,[yp(1);y;yp(end)],xi)
  • Not-a-knot cubic spline: yi = spline(x,y,xi)

  Example: Complete cubic spline for sin(x) using nodes 0, pi/2, pi, 3*pi/2, 2*pi:

  x = [0; pi/2; pi; 3*pi/2; 2*pi]
y = sin(x)
xi = linspace(0,2*pi,100);    % use 100 points in interval [0,2*pi] for plotting
yi = spline(x,[1;y;1],xi);    % f'(0)=1, f'(2*pi)=1
plot(x,y,'o',xi,sin(xi),'--',xi,yi)
legend('data points','exact function','complete cubic 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:

  [Q,R] = qr(A,0)
d = Q'*y
c = R\d

  Matlab shortcut (actually uses QR decomposition):

  c = A\y
  

Syllabus

Course outline handed out in the first class: Information about time & place, instructor, textbooks, syllabus, grading policy, Matlab assignments, final exam.

How to Hand in Matlab Homework

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; only print out 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 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 -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, Unix workstation on campus
• How to print text files for free