Welcome!

18.085 Fall 2011

Announcements

• Exam 3 Solutions: Exam 3 Solutions Fall 2011 [pdf]
  
  
• Exam 3 will be returned in the LAST CLASS MONDAY. The average was 90 with sigma = 7


• Homework 11 Solutions: Solutions [pdf]
  
  
• Homework 11 will NOT BE COLLECTED / it is practice for the exam We plan to have Monday help at 6:00 in 2-190 if we can Then review in class WED, exam on FRI, no Matlab sessions
  
  The exam in Walker will have 4 questions on Chapter 4
  1. Fourier series
  2. Discrete Fourier Transform
  3. Convolution
  4. Fourier integrals
  
  Suggested practice problems are 4.5 1, 2, 8, 11, 12, 22, 23, 24, 25
  
  
• Homework marks
  
  18085_hw_1_7_anonymous.xls
  
  
• Homework 10 due WEDNESDAY NOV 30
  
  3.5: 2, 6
  4.3: 8
  4.4: 1, 6, 7, 8, 9
  
  
• MATLAB Help session will be 5 p.m. Friday in the same room as lectures, as usual


exam 2 solutions Fall 2011

sols.pdf





• Homework for MONDAY NOV 21
  4.1: 1, 3, 6, 10, 13, 18
  4.2: 1
  4.3: 2, 6
  
  
  MATLAB Homework
  The homework is about the important "Gibbs phenomenon". Create a figure like Fig. 4.2 showing the partial sums of the Fourier series for the SQUARE WAVE in Exercise 4.1.2. It jumps from -1 to +1 at x = 0.
  Can you determine the amount of the overshoot? In what way does this same overshoot number appear in Fig. 4.2 for the delta function? What is the size of the overshoot if the step fcn jumps from 0 to A ?
  Find a relation between the WIDTH of the overshoot (the first lobe above +1) and the NUMBER of terms in the partial sum for the square wave.
  New question (to me): Beyond the overshoot in Fig. 4.3 it looks to me that there is an UNDERSHOOT below +1. Is this true in your graph of partial sums? What depth of undershoot? I hope you can blow up the figure near the jump.
  Extra: Also graph and blow up the Fourier series for the ramp function f(x) = max(0,x) near the corner (not a jump!) at x = 0.
  



• Questions about the grading of the quiz - 1st problem to Annie / 2nd problem to Sasha / 3rd problem to Vinoth / 4th problem to Inna



• Four questions in Wednesday's exam (11:00 in Walker) will be on
  
  1. Trusses
  2. Finite elements and weak form for 1-D equations
  3. div grad curl u(x,y) s(x,y) Green's formula Divergence thm
  4. Solving Laplace's eqn by f(z) z=x+iy Summary page 275
  (not 'mapping' in 3.4) (not formulas for the Fourier coeffs a b)
  



• MATLAB Homework 8
  Due Monday Nov 14
  No book homework on exam week
  Finite Element Method in 2D
  Two MATLAB codes solve two equations:
  *Poisson's equation on a L-shape geometry,
  *Laplace's equation on a circle geometry.
  These codes need distmesh, which is a great matlab code for meshing. The MATLAB codes are here:
  
  PoissonLFEM.m
  LaplaceCircleFEM.m
  
  fd_L.m
  refine.m
  distmesh2d.m
  drectangle.m
  dunion.m
  huniform.m
  
  The two codes for you to run are PoissonLFEM.m and LaplaceCircleFEM.m
  Two sample outputs are:
  
  SampleOutputFEMPoissonL.pdf
  SampleOutputFEMLaplaceCircle.pdf
  
  
  *The Poisson equation on the L-shape domain*
  
  The MATLAB code is PoissonLFEM.m
  1. The FEM code solves Poisson's equation,
  
        -u_xx - u_yy = 1
  
  on the L-shape domain, with essential boundary conditions U=0 on the boundary, and plots the output.
  You could experiment by making the parameter h smaller or by refining the mesh in other simple ways. With h=0.1, please print and hand in one plot of the output.
  
  2. Change the boundary conditions to be U=0 everywhere on the boundary, except for the y-axis, where the condition is now the natural boundary condition:
  
  du/dn=0 on the y-axis.
  
  There are two points on the boundary that are on both the y-axis and another boundary, they are the corners (0,0),and (0,1). The boundary condition on these two points is still zero. This can be done by changing one line (for you to find and change) of the MATLAB code.
  Please print and hand in:
  *The one new line of code
  *A plot of the new solution (so the y axis can be seen)
  
  
  *The Laplace equation in a circular domain*
  
  3. Give the weak form of Laplace's equation (in 2D) with essential boundary conditions U=Ub on the boundary. (The boundary condition on the test functions is still v=0.) This can be written by hand.
  
  4. The LaplaceCircleFEM.m code solves Laplace's equation,
  
        -u_xx - u_yy = 0
  
  in a circle, with prescribed boundary conditions U=U(b) on the boundary. The boundary conditions are U=+1 on the top half of the circle and U=-1 on the bottom half of the circle.
  
  Please print and hand in one plot of the solution.
  
  A comment about boundary conditions: This code imposes boundary conditions in a different way to the previous code. Suppose we have formed the free stiffness matrix K for Laplace's equation, as in the MATLAB code provided, and now we want to impose essential boundary conditions so that U=U(b) on the boundary. This time, U(b) may not be zero. Let b be a vector such that p(b,:) is the list of boundary nodes of the mesh, and U(b) is the prescribed essential boundary conditions at these nodes. In 1 line of MATLAB code we can move the known values to the right side of KU=F (with minus sign) and compute the load vector F as: F=-K(:,b)*u(b).
  
  5. The exact analytic solution is given in the textbook (p273, section 3.4). Use the formula to evaluate the exact solution Uexact = U(0.5,0.5) at the point x=0.5 and y=0.5. Please evaluate to 5 decimal places and hand this number in.
  
  We can use this to see how accurate our FEM code is. Use the code to compute the solution, with each of the following values of h:
  
  h=[0.5 0.25 0.1 0.05].
  
  For each value of h, record the FEM solution at the same point at which we evaluated the exact solution, i.e. UFEM=U(0.5,0.5). Estimate the error at this one point point as
  
  error = |Uexact - UFEM|.
  
  Plot the error (vertical axis) as a function of h (horizonatal axis), as four dots on a loglog plot. Estimate the slope of the line, to the nearest integer. The power of h in the error is approximately the order of accuracy (it might be VERY approximate in this example and that is okay).
  
  6. Please hand in the loglog plot, together with your (integer) estimate of the slope of the line.



• Homework 7 for FRIDAY November 4
  3.3: 6, 8, 10, 11
  3.4: 2, 3, 4, 6, 8
  3.5: 1, 3
  Planning to return these on Monday the 7th before the quiz on Wednesday 9th



• Monday October 31st - Help session will be 5-6 p.m. (only on this date!)


• Homework 6 for Friday OCT 28
  3.1: 5, 10, 12, 13, 14, 18
  3.2: 1, 6, 7, 17-18 (those are Matlab)
  



• Quiz 1 average 83 and standard deviation 11.4
  Almost everybody lost 9 or more points on 1d (which was harder than I intended, apologies !)
  For questions on 1, 2, 3 Inna can help / 4 goes to Vinosh
  Inna will do lead the help session at 6:00 in 2-190 MONDAY
  Quiz solutions will be posted soon
  Very OK results on Quiz 1 *****
  
  Solutions
  
  q1.pdf
  q2.pdf
  q3.pdf
  q4.pdf
  
  
  
  


• Homework 5 due Friday 21 OCT
  
  There is some book homework and some MATLAB homework, both due Friday.
  
  Book homework
  2.6: 2, 4, 6
  2.7: 1, 2, 3, 5, 7
  3.1: 8, 9
  
  
  MATLAB homework
  Guidelines for what to hand in for MATLAB homework: There are five questions (q1-q5). They require about one line each to answer. q3: You should answer yes or no, plus give an example.
  
  Answers:
  q1: The boundary conditions are u(0)=1 and u(1)=0.
  q2: For large t.
  q3: Yes.
  Example: The sample output shows Newton's method for two different initial conditions, one takes many iterations (green dashed line), whereas the other takes only a few iterations (red line).
  q4: Differentiate the reaction term: -A*exp(-E./(u-p)) + A*(1-u)*E.*exp(-E./(u-p))./(u-p).^2
  q5: It solves the steady state equation: 0 = v u_x + D u_xx + A (1-u) exp(-E/(u-p))
  
  
  This week we look at finite differences for advection-diffusion-reaction problems.
  
  u_t = v u_x + D u_xx + A (1-u) exp(-E/(u-p))
  
  The reactions are nonlinear. A, E and p are real parameters characterizing the reactions.
  
  Four MATLAB codes are needed:
  
  ADRFiniteDiffs.m
  adr_dt.m
  FPIteration.m
  NewtonIteration.m
  
  The main code to run is ADRFiniteDiffs.m
  
  Here, just as in the previous finite difference homework, you can:
  *Try varying the parameters for: diffusion (D) and advection (v)
  *Try increasing the number of meshpoints N (thus reducing the spatial discretization size h).
  
  See the following three sample outputs.
  For A=0:
  OutPutAzero.pdf
  For A>0:
  OutPutApositive1.pdf
  OutPutApositive2.pdf
  
  
  Try the following two exercises.
  
  
  Part (i) Set A=0 in the MATLAB code ADRFiniteDiffs.m
  
  This turns “off” the reactions and you are left with the LINEAR advection-diffusion problem, that is almost the same as the last homework exercise, except for two things:
  
  1. Now the boundary conditions are different. q1: Can you see what the boundary conditions are this time?
  
  2. Now the problem is time-dependent. We reach steady state when u_t =0, i.e. when
  
  0 = v u_x + D u_xx
  
  The code solves for the steady state in two different ways. One method is to use the linear solver, via “backslash”, that you already tried in a previous homework exercise. A different method is to use an ODE-solver that is part of MATLABs built-in functions, and integrates from an initial condition, u(0), to a solution, u(t), at a later time t, for a large value of t. You can try experimenting with small values of t and large values of t, and see how the solutions from the two methods compare when you plot both in the same figure. Try t=0.01 and t=10, for example. q2: Is the solution from the ODE-solver a better approximation to the solution from the linear solve (via backslash) for small t or for large t?
  
  
  Part (ii): Now try A>0.
  
  This keeps the reactions in the equation and they are NONLINEAR. At equilibrium we have
  
  0 = v u_x + D u_xx + A (1-u) exp(-E/(u-p))
  
  So we can not simply use backslash, as we did when A=0. Try solving the problem for equilibrium, in THREE different ways:
  
  1. MATLAB's ODE solver, which starts with an initial condition and integrates in time until some large value of t.
  
  2. Newton's method.
  Try at least two different choices for u(0). q3: Does it matter where you start with Newton's method?
  
  For Newton's method the Jacobian is important. In this example the Jacobian is J=D+A+R, where D is for diffusion, A is for advection and R is for the reactions. In this example, R is diagonal. q4: Can you give the jth term on the diagonal of R?
  
  3. Fixed point iteration.
  Let
  
  Q ( u ) = v u_x + D u_xx
  R ( u ) = A (1-u) exp(-E/(u-p)).
  
  Suppose we start with an initial guess u_1 and then for n>=2 we compute the next term u_{n+1} from the current term, u_n, by
  
  Set b = R( u_n ). Solve Q u_{n+1}= - b
  
  Suppose we are lucky, and we happen to reach a fixed point u* such that u*=-Q^{-1}R(u*). q5: What equation does u* solve?
  
  Notice that: Q(u) is linear so it is easy to compute Qu, or to solve Qu=b for u, (via MATLAB's “backslash”, say). R(u) is nonlinear so we can easily compute R(u) explicitly but it is not so easy to solve R(u)=b for u. This is why the fixed point iteration method above first computes R, and then solves with Q (and not the other way around).
  
  
  



• EXAM next Wednesday 11-12:30 in Walker: 4 questions
  1. Finite difference matrices
  2. Least squares and projection of b
  3. Line of springs and oscillation
  4. Small network
  
  


• Tuesday Sept. 27th - Inna's office hours will be 2-3pm.


• HELP SESSIONS WILL BE MONDAY at 6:00 IN 2-190
  Inna will hold the first help session on Monday Sept 19 at 6 in 2-190
  Ready for questions about Chapter 1 and positive definite matrices
  Probably we will not do Sections 1.7-1.8 of the book in full detail (although the SVD = Singular Value Decomposition is neat)



• Homework 4 for FRIDAY OCTOBER 7
  2.3: 7, 8, 10, 12, 18, 22
  2.4: 1, 3, 8, 17
  
  There will be a REVIEW IN CLASS October 7 of all we have done.
  


• Homework 3 is due FRIDAY Sept 30
  1.7: 12
  2.1: 3, 6, 7, 8
  2.2: 6, 7, 8
  
  MATLAB Question : This uses the codes just added to the website for advection-diffusion -Du_xx + vu_x = 1 in the earlier hwk problem 1.2.19. Those codes use finite differences (u_x is centered or upwind) and compare with the true u(x).
  
  Part 1: Derive the true solution given in the code. A particular solution is u=x/v. Add the general solution to -D u_xx + v u_x = 0 with constants A and B, then find A and B.
  Part 2: Decide the accuracy (both centered and upwind) using the new code Ch1Q19FiniteDiffsConvergence.m that is also on the website The max error decreases with what power of h as h goes to 0 ?
  Part 3: Continue the experiment for very small h. We think the error may start to increase and we don't know why. Why ?
  

Notes on Homework 2
- [L,U] = lu(A) includes any row exchanges in L (so it may not be triangular) [L,U,P] = lu(A) will separate out that permutation matrix
- The product of pivots in U is MINUS the determinant of A after an odd number of row exchanges
- derivative of delta function : Integrate by parts to see that it picks out -g'(0)
- to avoid row exchanges for positive definite matrices K, use chol(K) where chol = Cholesky



• Homework 2 is due FRIDAY Sept 23 (no class Wednesday)
  Your class number will be on Homework 1 returning Monday
  Can you put that number on all homeworks in large clear letters?
  1.3 11, 12, 13
  1.4 1, 4, 15
  1.5 1
  1.6 3, 9, 14, 15
  
  
• Homework 1 is due Wednesday September 14.
  Please staple and write clearly / could you PRINT your name.
  Discussing problems with other students is allowed / this is NOT a test.
  Grading will be fairly quick / important to do homeworks!
  MATLAB results can always be printed out / not retyped.
  Problems will come partly from the CSE textbook and sometimes new.

• Section 1.1 of the book is on math.mit.edu/cse (along with many codes)
  This first time I will scan in the problems from 1.2: Problems from Ch. 1.2
  
  1.1: 2, 5, 9, 14, 26
  1.2: 11, 19

	Ch1Q19FiniteDiffs.m
	output1.pdf
	output2.pdf
	output3.pdf
	Ch1Q19FiniteDiffsConvergence.m
	This is a code for Problem 1.2.19: Finite differences for the linear advection-diffusion equation - D * u_xx + v * u_x = 1 in Homework 1 [1.2.19] You could test this code with different parameters D, v, h as suggested below. The code solves and then plots the solutions. It compares the true analytic solution with the solution from finite differences, using (1) centered differences and (2) upwind differences. See sample outputs. It is easy to * see the matrices K (diffusion) and Del0 or Del_+ (centered/forward) for small N * see the boundary conditions u(0) = 0 and u(1) = 0 * see whether centered or upwind is more accurate * make h smaller and see the approximations converge (increase the number of mesh points, N, to make h smaller) What is the order of convergence (try the second code for this)? * try different values of the parameters D, v (i) v=D=1 (ii) v << D (and see the numerical approximation is good, even with large h) (iii) v >> D (and see a boundary layer because diffusion is so weak)

Class Stuff

• Time/Location: MWF 11-12, Room 4-370
• Lecturer: Gilbert Strang
• Office: 2-240
• e-mail: gs@math.mit.edu
• TAs
  • Inna Entova-Aizenbud
    e-mail: innaento@math.mit.edu
    Office Hours: Tuesday, 13:00-14:00 at 2-251, or by appointment
  • Vinoth Nandakumar
    e-mail: vinothn@math.mit.edu
    Office Hours: Wednesday, 3:30-4:30 at 2-090, or by appointment
• Graders
  • Annie Huang
    e-mail: a_huang@mit.edu
  • Sasha Tsymbaliuk
    email: sasha_ts@mit.edu

Course Topics

• Applied Linear Algebra
• Applied Differential Equations
• Fourier Methods
• Algorithms
• Course outline 2008

Additional Information

• Quiz dates 11-12:30 in Walker: October 12, November 9, December 9
• Goals for the Course: See applications of calculus, ODE, linear algebra, and discrete methods without going into too much proof.
• Textbook: Computational Science and Engineering (Wellesley-Cambridge, 2007).
• Grades: Homework 40%, 3 quizzes 60%, no final exam. Please email Prof. Strang about conflicts with quiz dates.
• Homework: Due Wednesdays. Please write neatly, staple, and include your class number (to be assigned after Homework 1).
• Use of MATLAB for tedious calculations is encouraged, however you need to know how to do the basic algorithms taught in the course by hand (at least for small matrices) for the quizzes.

Class Resources

Movie of elimination   moe.m   (also need realmmd.m )

Code to create K,T,B,C as sparse matrices

MATLAB's backslash command to solve Ax = b   (ps, pdf)

Getting started with Matlab: http://ocw.mit.edu/OcwWeb/Mathematics/18-06Spring-2005/RelatedResources/