|
|
|
Clear Excellence |
H |
86-100 points |
|
Entirely Satisfactory |
P |
70-85 points |
|
Low Passing |
L |
50-69 points |
|
Failed |
F |
0-49 points |
The transformation to undergraduate grades is as follows:
|
|
|
B+ |
83-88 |
C+ |
65-70 |
D+ |
47-52 |
|
A |
95-100 |
B |
77-82 |
C |
59-64 |
D |
41-46 |
|
A- |
89-94 |
B- |
71-76 |
C- |
53-58 |
F |
0-40 |
Course Texts
The course will use the following textbooks:
James W. Demmel, Applied Numerical
Aslak Tveito, Ragnar Winther, Introduction to Partial Differential Equations – A Computational Approach, Springer, 1998.
David Kincaid, Ward Cheney, Numerical Analysis, Brooks-Cole, 1996.
Lecture notes on material not found in the texts will be distributed periodically through this web site.
Lecture notes
|
Lecture |
Topic |
Notes |
Lecture |
Topic |
Notes |
|
1, 1/12 |
Introduction |
21, 3/4 |
Symmetric eigenvalue problem |
||
|
2, 1/14 |
Mathematical Background |
3/7 |
Midterm review |
|
|
|
3, 1/19 |
Gauss elimination |
3/9 |
Midterm Exam |
|
|
|
4, 1/21 |
Perturbation theory |
22, 3/21 |
Two-point boundary value problems |
||
|
5, 1/24 |
Linear space geometry |
23, 3/23 |
Generating Finite Difference Formulas |
||
|
6, 1/26 |
Panel method application |
24, 3/28 |
Convergence of Finite Difference Approximation |
||
|
7, 1/28 |
Condition number estimator |
25, 3/30 |
Sturm-Liouville problem, PDE’s: First order PDE’s |
||
|
8, 1/31 |
Iterative refinement |
26, 4/1 |
Classification of PDE’s. Canonical PDE’s. |
||
|
9, 2/2 |
GEPP implementation, Strassen algorithm, BLAS, LAPACK |
27, 4/4 |
Basic iterative methods for solving linear systems |
||
|
10, 2/4 |
Cholesky algorithm |
28, 4/6 |
Eigenvalues of matrices from finite difference discretizations |
||
|
11, 2/7 |
Least Squares - Normal Eqs. |
29, 4/8 |
Homework review (see HW9 solution) |
|
|
|
12, 2/9 |
Least Squares - QR |
30, 4/11 |
Convergence of iterative methods. |
Demmel 6.5.5. |
|
|
13, 2/11 |
SVD Decomposition |
31, 4/13 |
Convergence acceleration. Chebyshev acceleration. |
Demmel 6.5.6 |
|
|
14, 2/16 |
SVD Properties, QR through Givens rotations Calling BLAS example (Linux, Intel Fortran, Intel MKL) |
32, 4/15 |
Krylov subspace methods. |
Demmel 6.6.1. |
|
|
15, 2/18 |
Householder reflections |
33, 4/18 |
Conjugate gradient method. Analysis. |
|
|
|
16, 2/21 |
Eigenvalue problem: theory |
34, 4/20 |
Fast Fourier methods & cyclic reduction. |
|
|
|
17, 2/23 |
Power method, Orthogonal iteration |
35, 4/22 |
Multigrid methods. |
|
|
|
18, 2/25 |
QR algorithm |
36, 4/25 |
Finite element methods. |
|
|
|
19, 2/28 |
|
37, 4/27 |
|
|
|
|
20, 3/2 |
|
38, 4/29 |
Final exam review |
|
Homework
Homework is given each Friday, starting Jan. 21, and is due the next Friday. Solutions are posted after the Due Date. Please download the homework from the website. The homework will involve programming, analysis of numerical experiments and theoretical work. Draft a single file with the results from your numerical experiments and your comments thereon. Include answers to any theoretical questions in this file. If you are not proficient at typesetting mathematics, you may turn in a handwritten answer to the theoretical questions included in the homework during the Friday class when the homework is due. You must include the computer code used to generate results with your homework as an e-mail attachment. Each attachment should contain the complete source code for one question, such that the attachment can be downloaded and then executed. Label the attachment as Q1.m for a Matlab code or Q1.f90 for a Fortran code for the first question, Q2 for the second question and so on. Upon completion of the homework, send it by e-mail to the course graduate grading assistant, Jie Zhou (jzhou@email.unc.edu). The grading assistant will verify correctness of codes you send in.
Notes:
1. Fortran is preferred when compiled code is called for. C, C++ is accepted. However, do not use encapsulation of basic algorithm functionality into C++ classes. I'm sure it's a very clever programming trick, but we're mostly interested in code transparency.
2. Sign the code you turn in.
|
HW # |
Due |
Problems |
Solution |
HW # |
Due |
Problems |
Solution |
|
1 |
1/28 |
7 |
3/11 |
||||
|
2 |
2/4 |
|
8 |
4/1 |
|
||
|
3 |
2/11 |
9 |
4/8 |
||||
|
4 |
2/18 |
|
10 |
4/15 |
|||
|
5 |
2/25 |
|
11 |
4/22 |
|
||
|
6 |
3/4 |
|
12 |
4/29 |
Bulletin Board
Answers to questions that arise from the homework assignments will be posted here for everyone's benefit.
|
HW # |
Question |
Answer |
Computer Code |
|
1 |
1. p-norms for HW Q4? 2. Polynomial coefficients for Q2? |
|
|
|
2 |
1. Q1. Are we allowed to use gepp.f90 that you wrote in class as a starting point for the homework assignement? 2. Q3. Following the discussion in Section 2.7.3 of Demmel, the size of width of non-zero columns in a particular row can grow to 2w+1 + w, so the result can't be stored in an n-by-(2w+1) array anymore. Does that imply that we shouldn't do *in-place* LU decomposition anymore, which is what is done for standard GEPP? 3. Q5. Is A a real square matrix? 4. Q4. Can we pick arbitary bandwidth w? 5. Q3. The use of only n-by-(2w+1) storage locations when doing pivoting might be a bit too tricky for people without some programming background, so I'll accept solutions for GE without pivoting. People who do find the correct way to implement this with partial pivoting will receive a bonus point. Hint: The problem is that row swaping seems to require some extra memory locations in the upper triangular part, but notice that some memory locations also get freed for other uses. |
1. Yes, we develop code in class with that exact purpose - to have a
starting point for assignments.
2. Try to work with GEPP as presented in
class and implemented in gepp.f90. There is a way to use just n-by-(2w+1)
storage. 3. Yes, A should be considered real & square (n by n) 4. Yes, whenever something is not explicitly specified you are at liberty to make any reasonable choice. I would interpret w=0 as unreasonable though and w=1 as marginally reasonable ;-) |
|
|
3 |
1. The homework assignment says to write Fortran 90 code for Hagar's condition number estimator. Can we do this in C++? 2. Can we use Matlab? |
1. Yes, the assignment can be carried out using C++, but there is no initial template so you'll have to start from scratch. Also, a final assignment on linear systems regarding performance issues shall be issued. In HW4 the matrix for the flow around the circle will be provided as a Fortran 90 routine. At that time, the routine will have to be translated into C/C++ or linked in. 2. No, you should use a compiled language for this homework. It will be reused next week to study performance issues. Matlab is very good at prototyping algorithms but poor in building high-performance models. One cannot claim to be a bona-fide computational scientist without a basic familiarity with at least one compiled language. |
|
|
4 |
|
|
|
|
5 |
|
|
|
|
6 |
|
|
|
|
7 |
1. Q3: Do we need to do LU factorizattion and store the multipliers? or just reduce it and forget about the multipliers? |
1. You should determine the matrices M(1),M(2),...,M(n-1) such that M(n-1)...M(1)=A. This is intended as the presentation of the mathematical algorithm. Storage questions would only arise in a computer implementation, which is not required here. |
|
|
8 |
|
|
|
|
9 |
1.
(General) i wrote down that you can set up an ivp for hyperbolic pdes but
not for parabolic or elliptic. the notes online say you can set up an ivp for hyperbolic and parabolic pdes
but not for elliptic. which is right?
in other words, can you set up an ivp for parabolic pdes? 2. (Q1)
I understand how to do this process (turning higher order PDE's into first
order systems) for equations where all the terms are the same order.
However, I am having trouble generalizing this to #1, which has a first
order and a third order term. Do you have any tips on what way I can attack
this problem compared to the other (one issue that comes up is that some
terms cannot be canceled with cross derivatives)? 3. I've read the hint about
problem #1 on the course page, and understand the introduction of the
variable v when dealing with the heat equation. What I don't understand is that
when you take p_x + Ap_t
= a for p = [u v]', p_x = [u_x v_x]' = [u_x u_xx] and p_t = [u_t v_t]' = [u_t u_xt] My question is, what do you do
with u_x and v_t = u_xt? Does the relationship between u_t
and v_x imply something about the relationship
between u_x and v_t? 4. In
the upwind algorithm I'm not sure of the significance of the choice between U^n_j-U^n_{j-1} and U^n_{j+1}-U^n_{j} based on the sign of V? since both of these are a derivative estimate just on
different sides in the grid. 5. When
actually implementing this algorithm, what would we use at the endpoints? That
is for U_j+1 if j is the largest index |
1. Parabolic PDE's are a
transition from hyperbolic to elliptic and are a bit hard to explain in
simple terms. Here's my try at it - hope it works. Hyperbolic PDE systems are
reducible to a set of separable equations of the form w_t
+ a w_x = 0 and each of these equations is
solvable along characteristics, hence we can only pose an initial value
problem. Furthermore the value of w at
any point only depends on initial value data that can be transmitted to
that point - we say that the equation has a 'finite domain of dependence' Parabolic PDE systems cannot be
reduced to diagonal form because at least two of the eigenvalues
of the system matrix correspond to colinear
eigenvectors (geometric multiplicity < algebraic multiplicity). Thus
some of the information needed for a marching procedure along the
characteristics is missing. The tricky part now is that the domain of
dependence is infinite but the influence from boundary conditions dies down
exponentially fast (this has to do with the 'fundamental solution' or Green's function for a
parabolic equation). Therefore, we can still pose an initial value problem
for parabolic equations also. Elliptic PDE systems cannot be
reduced to diagonal form because some of the eigenvalues
have non-zero imaginary parts. There is no marching procedure possible and
the domain of dependence comprises the entire definition domain and there
is no rapid decay of the influence of boundary conditions. We can only pose
boundary values problems. 2. Our
objective is to put this in a form: A q_t + B q_x = a where at least one of the matrices A or B is
invertible so we can then either obtain q_t + C q_x = c (C=inv(A)*B, c=
inv(A)*a) or q_x + D q_t
=d (D=inv(B)*A, d=inv(B)*a). Hint: When the order of derivatives in the PDE
is different, we can make an identity matrix appear as the coefficient of
the q derivative of higher order. Example: heat equation is u_t – u_xx=0. Set v=u_x. Then v_x – u_t =0. This is now a 2 by 2 system of the form q_x + A q_t=a with q=[u v]’. 3. The
two 1st order equations equivalent to the heat equation are: u_x = v v_x - u_t = 0 In
matrix form q_x + A q_t = a, with q=[u v]',
a=[v 0]' and A = [ 0 0 -1 0] We
don't need to consider the relationship v_t=u_xt, we already have a system for which we can
determine the eigenvalues. 4. The
upwind scheme is stable numerically. Math 221/2 goes into detail as to the
analysis but the qualitative reason is easy to grasp. When V>0 characteristics
propagate information about the initial condition from left to right in the
(x,t) plane as time
advances. (dU/dx)_j = (U_j
– U_{j-1}) maintains the direction of this
‘information flow’ while (dU/dx)_j=(U_{j+1}-U_j) does not. 5. This
indeed is the key question here. There are two families of characteristics,
one right-going, one left-going. The question
deliberately left the issue of what to do at the boundaries open –
you are invited to think about appropriate conditions at x=0 and x=1. One
easy approach is to impose periodic boundary conditions, i.e. U_0 = U_n, U_{n+1} = U_1, U_{-1} =
U_{n-1}. |
|
|
10 |
1. Q1:
for our matrix A as in problem 1 for instance,
should it include the separate fourth order aprox.
that we had to construct for row 1 and n-1, or are
you just looking a derivation of the eignvalues simliar to the one done in class where its more just
the interior points? |
1. Use the stencil corresponding to interior points and truncate it appropriately at the corners. |
|
|
11 |
|
|
|
|
12 |
|
|
|
General questions
Q: Questions regarding your policy on code structure.
1- We you state that classes cannot be used for "basic algorithm functionality", does this apply to something like a matrix class, which does not perform any of the algorithms we are studying but instead holds the data and performs simple (for example finding the 1-norm, or matrix addition) operations on it.
2- Does this proclude separating parts of our program into separate files for the purposes of: 1)code clarity, 2)code reusability. (Not as a separate class, just as an included file).
A: What can/cannot be included as class methods. Can be included:
· data allocation and storage
· simple functions, e.g. norm, addition
Cannot be included: the algorithm being studied in that particular HW,
e.g. QR factorization.
As to multiple files, though I understand the usefulness of properly
structuring code into files, the objective of the HW is to show that you have
understood the basics of an algorithm. I have to be able to come to a
conclusion in the few minutes I can allocate to each student. Fishing around
for an obscure definition in an .h file is not conducive to a careful
appraisal. Hence, please include everything in one C/C++ file (or better yet
in one Fortran file).
Extra Credit
Students who wish to acquire more course credit points can choose to carry out the following projects:
EC1: Present a Fortran implementation of Q4 from HW7. Test the algorithm on systems with Hilbert matrices of sizes 10,20,…,100 and compare the solutions against those from GEPP. Comment on benefits/drawbacks of this algorithm compared to GEPP. (5 points)
EC2: Present a Fortran implementation of Q6 from HW7. Test the iterative refinement technique for linear least squares on Q2 from HW4. (8 points)
EC3: Redo Q1 from HW6 for polyethylene (C2H4)m, i.e. C2H4 repeated m times. The chemical bonds to a Carbon atom are at angles corresponding to vertices of a regular tetrahedron. Present a Fortran implementation for Q2 from HW6 but using the polyethylene molecule. In particular make the implementation efficient in determining the largest eigenvalue for m+1 when that for m has been determined. (12 points).
Note: C/C++ implementations are acceptable with the course restrictions on defining C++ classes. The assignments are due on April 29.
Computing
Homework assignments will contain questions requiring numerical experimentation. The recommended programming environments are Matlab and Fortran 95, preferably on a Unix-like operating system (Linux, OS X). The UNC scientific servers (baobab, chastity) offer all the necessary tools. Students willing to manage their own systems are encouraged to install Linux on their PC-compatible machines. An installation CD-set (Suse 9.0) is available by request - send e-mail to see if it's available and then drop in during office hours to pick it up. You may keep the CD's until the next office hour session. People running Microsoft Windows might wish to install VMware workstation (an academic version is available - try Google to find a retailer) and run Linux and Windows concurrently. There are two free Fortran 95 compilers one can install under Linux:
The GNU compiler has also been packaged for OS X.
There are also various Fortran 95 implementations for Microsoft Windows.
A comprehensive list of free Fortran compilers is available.
Computing FAQ
1. Do I have to use Fortran (C)?
Yes, this course is an introduction to the world of scientific computing for which facility in the use of compiled languages is required.
Mid-term Examination
Final Examination
The Final Examination will be held on Monday, May 9 in our regular classroom, Phillips 383. The examination will concentrate on aspects presented after the Spring Break. Make-up questions for material covered in the Mid-term examination will be offered (10 points). Mid-term material can appear on final questions also (see above Part 1 synopsis)
Exam preparation 1: Questions Solution
Exam preparation 2: Questions Solution