18.085 SUMMER 2012



You FINISHED 18.085
!!! CONGRATULATIONS !!!

!! Good luck with your future endeavours !!

HmWk4 MATLAB

FRIDAY is LAST lecture. Review of topics we covered in course.
A chance for you to ask questions about anything in the whole course, or homeworks ...
and give feedback (e.g. how can we make this course more useful to grad students in engineering?) about what we could do better next time.

COMMENT ON HmWk7, Q2: you have to change a boundary condition. This can be done by changing one line (for you to find and change). HINT: the one line is already in the code, and you simply have to uncomment it; i will ask in class how people are going with this on WEDNESDAY. Depending on how everyone is going, there may need to be more help with this?
COMMENT ON HmWk7, Q5: You may have to convert x,y coordinates to polar coordinates.

Homeworks 4 and 5 due on Monday? Is that okay?
Remember that more homeworks are coming, and all grades have to be in to central admin by a fixed date in the third week of August, and i need to grade them first!


Will reorganize this page soon, so that we can actually find things!

Homework 8 (LAST homework):
Question 1. VERY approximately, what two frequencies occur at maximum power in the whale recording?
Question 2. What is the first note the Cello plays?
Look up a table of frequencies on the web, maybe try
http://www.phy.mtu.edu/~suits/notefreqs.html.

FFTDemo2.m
whale1.wav
Cello.wav


Homework 7:

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. You may have to convert x,y coordinates to polar coordinates. 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.
---- End Homework 7 ----

Quiz2Q4.m

7 votes for Fourier and 3 votes for Trusess. Fourier is the winner.
Now CONFIRMED: QUIZ 2 in lecture time Friday Aug 10. Same rules as Quiz 1.

QUIZ 2 covers four sections of textbook:
section 3.1 (finite elements in 1D with hat functions), examples from class
section 3.3 and 3.4 (gradient, divergence, streamlines, equipotentials, Laplace's equation and complex variables)
and section 4.1 Fourier series, good becuase of orthogonal basis functions, Gibbs phenomena, and the matlab demo in class of the power at different frequencies (but only a very simple question)
Some past exams and solutions that we looked at (also available on the 18.085 site but here too in case it helps)
Note we only looked at certain questions, e.g. Finite Elements in 1D, Laplace's equation, Gradient Divergence, vector field tests.
quiz2f07.pdf
quiz2f07solutions.pdf
quiz2_18085_f06.pdf
q2_sol_18085_f06.pdf
quiz2_18085f02.pdf
f02q2sol.pdf

2 MATLAB homeworks coming on: 1. Finite Elements in 2D, Laplace and Poisson equations, in circle and maybe L-shape domain.
2. Fourier Transform and FFT of signals.


Help session Friday after class in SAME room as lecture. Will do one example of a truss and write the associated matrix explicitly. Also, discussion of MATLAB or homeworks, or other questions people might have.



Small experiments with lecture format: our lecture time is very long. So planning to break it up a bit. At least one 5 min break to stretch, get a drink etc. And maybe the last part of the lecture can be a change of pace.


! Quiz 1 was harder than i intended. Apologies! Grading will take this into account. Overall, it was a very good performance on quiz 1. Congratulations!


First quiz coming soon, perhaps near July 20. Exact date still to be decided.
NOW CONFIRMED: QUIZ on Friday July 20. Be there.
OPEN BOOK. Bring your textbook, notes, etc.
Calculator okay.
But no smart electronic gadgets.
Quiz is in what would normally be a lecture, so you have the normal time alloted for lecture (approx 1.5 hours) for your quiz. And probably more time if you need it, so long as the room continues to be free.



! Help session Thursday 12.30 in SAME room as lecture (assuming we don't get booted out) !



Tentative proposal: second (and last) quiz be in the second last week of lectures, which is the week of August 10
Now CONFIRMED: QUIZ 2 in lecture time Friday Aug 10. Same rules as Quiz 1.

Very okay scores on HomeWork2: class average was 88%
Lowest homework mark is dropped for final grades, so if you did not do so well this time, there are still more chances.



!!! NEW DUE DATE for HomeWork 2: Friday July 6 !!!



VERY OPTIONAL review of material so far,
This TUESDAY (July 3) 3-5 in the SAME room as lectures.
Extended office hours / Help session. This is a chance to ask questions about MATLAB
I will also plan some MATLAB examples, and some more examples from class
This is TENTATIVE (still awaiting confirmation that we can have the room)... NOW CONFIRMED.

M W F
9.30 - 11
2-131
Dr Shev Macnamara
Room 2-371, office hours: Tuesday 3-5


Lecture 1: section 1.1 of CSE textbook.
Special matrices K,T,B,C
symmetric, tridiagonal, invertible or singular

First homework: problem set 1.1, q26, remember your NAME!!
DUE MONDAY June 18


Lecture 2: section 1.2 of CSE textbook
Second differences from 1,-2,1
-u''=f(x) becomes Ku=f
boundary conditions, first and last rows, K, fixed-fixed, T, free-fixed



HomeWork2 Due Monday July 2 is an OLD date
!!! NOW EXTENDED !!!
!!! NEW DUE DATE: Friday July 6 !!!
I had hoped to go over the homework in class but ran out of time. Keeping the momentum from the last good lecture was the top priority. So Monday's lecture i will go over the homework briefly in class.

Here is a short MATLAB code that may help to get started on question 6:
HmWk2Q6.m

Here is MATLAB code for solution to question 6:
Here is a MATLAB code solution to question 7:
Here is a MATLAB code solution to question 8:
ComplexExpPlot.pdf
Use the MATLAB code DeltaExample.m (below) to make matrices like K and T in the homework.

Here is another short MATLAB code, that may help to get started on question 9:
HmWk2Q9.m
Here is a MATLAB code solution to question 9:

Remember to check your answers in MATLAB!
For example, try:
ezplot('real(exp((-1 + 20i)*x))',[0 2])

Here is another short MATLAB code - an example of a LogLog plot:
LogLogPlot.m
Another short MATLAB code - an example of a "for loop":
LoopExample.m



Lecture 3: section 1.3 of CSE textbook
Elimination, LU factorization, K=LU, K=LDL^T and symmetry


Lecture 4: Reviewed material from first three lectures
Examples in MATLAB
Lorenz equations in MATLAB


Lecture 5: section 1.4 of CSE textbook
Inverses, Delta functions, Green functions
Continuous and discrete case
Here is an example MATLAB script from class to solve T U = F for a Delta function on the right hand side

DeltaExample.m


Lecture 6: section 1.5 of CSE textbook
Eigenvalues and eigenvectors
Eigenvalue decomposition of square matrix [V D] = eig(A)
Symmetric matrices are nice: A = V D V^T
Examples in MATLAB
eigshow, eigenwalker!


Lecture 7: More on eigenvectors and eigenvalues
Given eigenvalue decomposition A = V D V^-1, we found
exp(A) = V exp(D) V^-1, and solved linear differential equations of time.
Started positive definite matrices: x^T K x > 0
and section 1.6 of the CSE textbook




Lecture 8: !!! SPECIAL GUEST LECTURE !!!
section 2.4 of CSE textbook
Framework for applied mathematics: A^T C A .
! MAIN THEME OF COURSE ! IMPORTANT !
Examples (!IMPORTANT EXAMPLES!):
*Masses and springs
*Networks/graphs, voltages and currents
*Least squares: NORMAL equations A^T A u = A^T b


Lecture 9: Continued the A^T C A framework from the previous special guest lecture
Started masses and springs examples of the A^T C A framework, in section 2.1 of CSE textbook.
displacements u; elongations e=Au; internal forces in springs w=Ce, balance with external forces f = A^T w.
More examples of positive definite matrices.


POSITIVE DEFINITE DEMO
A MATLAB code to draw the energy function P(u) = u^T A u

PositiveDefiniteDemo.m

Lecture 10:Finished the masses and springs example and of the A^T C A framework, (section 2.1 of CSE textbook).
Started SINGULAR VALUE DECOMPOSITION (SVD).
A = U S V^T
Given an SVD, we can:
1. read off the rank of (A), and
2. write the PSUEDOINVERSE A+ = V S+ U^T.
3. Least squares problem (A x = b, overdetermined by tall and thin A) has a neat formula: x = A+ b.

PRACTICE QUIZ
PracticeQuiz.pdf
PRACTICE QUIZ SOLUTIONS
PracticeQuizSols.pdf
MATLAB code to experiment with the matrix in the practice quiz
PracticeQuizQuestionExample.m

No lecture July 4

HomeWork3 Due (new date) MONDAY July 16, and HELP SESSION Friday after class, so long as no one kicks us out of the room



Lecture 11: (Section 2.2 of CSE textbook) Masses and springs again, but now masses moving with TIME, so displacement u = u(t). Still same K=A^T C A framework, and now two new equations:
M u'' + K u = 0, and
K x = omega^2 M x
Second is an eigenvector equation, solve in MATLAB with [V, D] = eig(K,M). M is the mass matrix, and omega is the natural frequancy.

*Also started Least Squares* (Section 2.3 of CSE textbook).
Want to solve A x = b, when A is tall and thin.
Usually, there is no exact solution. Too many equations, not enough unkowns. Instead, we find the approximation that minimizes the error, e = Ax-b, in a least squares sense: find x to minimize ||e||^2 = ||Ax-b||^2.
How do we know which x to choose?
Answer: solve NORMAL equations A^T A u = A^T b
We stated this important result without justification, and did not complete the example. So that is coming next lecture.


Lecture 12: (Section 2.3 of CSE textbook) Least Squares
We saw NORMAL equations A^T A u = A^T b as example of projection, and saw a picture of the column space of A, and the error vector was orthogonal to this space


Here is a MATLAB code for the Least Squares Example in the lecture
It includes an example of how to solve this with the SVD (Very good to know.)
LeastSquaresExample.m LeastSquaresExample2.m

Here is a MATLAB code to learn about the Householder method of computing the QR factorization
QRHouseholder.m

Lecture 13: (Section 2.3 of CSE textbook) Finished Least Squares
Another brief review of section 2.4, A^TCA, Kirchhoff's law, incidence matrix A, degree matrix D, adjacency matrix W, graph laplacian A^T A = D-W.

Lecture 14: (Section 2.3) Looked at factorizing A = QR.
Mentioned Gram-schmidt very briefly. Talked about Householder, saw original paper, and have our own Householder MATLAB code to play with.
Practice quiz on least squares.
Started trusses, A^T C A, is it stable or are there collapse mechanisms? Is there a nontrivial solution to Au = 0? to be continued...


Hoping these homeworks are fast, and maybe even fun?! Good to look at HomeWork 4 NOW.
HomeWork4.pdf

HomeWork5.pdf

Lecture 15: (Section 2.7) Continued with trusses.


Lecture 16: (Section 3.1) Started FINITE ELEMENTS.
Talked about WEAK FORM ! IMPORTANT !
Agreed we could get from the weak form to KU = F.
But did not succeed in a small example with hat functions. So that is coming next lecture.


HomeWork6.pdf

Lecture 17 : (Section 3.1) Continued Finite Elements.
Got most of the way through a small 1D example with hat functions. Will finish that example next lecture, and spend more time on Finite Elements. Also saw some neat applications of A^T C A, and of least squares, turning out to be useful in imaging for example.

Lecture 18 : (Section 3.1) Continued Finite Elements.
Small 1D example with hat functions done ad nauseam. Two ways to make the same stiffness matrix K: (i) use the entry-by-entry formula for Kij, or (ii) assemble K from "element matrices": K = K1+ K2 + K3.

Lecture 19 : (Section 3.3) Started Gradient and Divergence.
When is a vector field the gradient of a potential? Did an example of the zero vorticity test. Continuous analogues of Kirchhoff's Laws.

Lecture 20 : (Section 3.4) Laplace's equation: u_xx + u_yy = 0
Solutions come in pairs of a potential function u(x,y) and a stream function s(x,y). They are the real and imaginary parts of a complex function. Easy to see the pattern in polar coordiantes (r, theta). Euler's formula is useful: exp(i theta) = cos theta + i sin theta. Equipotentials perpendicular to streamlines, via Cauchy-Rieman equations. Found the general solution u(r, theta) to Laplace's equation inside the unit circle, given a boundary condition on the unit circle. We still have not done Fourier series, but we had to use a tiny bit of Fourier for the last example.

Lecture 21 : (Very Optional: Fun with Moebius transofrmations and conformal maps for Laplace's equation.)
(Section 3.5) Started finite differences for Poisson's equation in 2D.
Talked about how to make matrix K2D. Still to think about how to solve K2D U = F.

Lecture 22 : (Section 3.6) Started finite elements for Poisson's equation in 2D. Will write the weak form in 2D on the board again, Gauss-Green formula is integration by parts in 2D
(Finished Section 3.5 MATLAB demo in class of the matrix K2D. Solve K2D U = F with the eigenvector decomposition, and the FFT.)

Lecture 23 : (Section 4.1) Started Fourier series. Sine series (odd functions). Cosine series (even functions). Complex fourier series.
Gibbs phenomena near jumps, and looked at MATLAB demo of partial sums trying to converge to square wave, delta function and ramp
Differentiation is like mulitplication by ik in Fourier space: functions get smoother when you integrate and rougher when you take derivatives. Rate of decay of coefficients is linked to smoothness e.g. smooth functions: coefficients decay fast! e.g. 1/k for step function and still to say more.

Lecture 24 : (Section 4.3) Started Fast Fourier Transform. Saw an example in MATLAB.Still to say more.
Continued section 4.1: Important fact: Energy in function = energy in physical space = energy in frequency space. Will continue both of these in Monday's lecture.
Lecture 25 : Review of weak form. Boundary conditions in weak form.
MATLAB examples of Finite Elements in 2D for Laplace equation in circle geometry, and Poisson equation on L shaped geometry.
Did an example of rolling two dice: the sum of two random variables, has a probability distribution that is the CONVOLUTION
section 4.4 is convolution. I don't know how far we will get with this important topic, plus we still have to tidy up a few Fourier FFT topics.
Lecture 26 : Finished convolution.