Course 18.085: Mathematical Methods for Engineers I
(Fall 2007)

Department of Mathematics
Massachusetts Institute of Technology

[info] [homework] [notes] [class resources] [old exams] [links]

General Information

Lecturer: Gilbert Strang (gs@math.mit.edu), room 2-240

Lectures: MWF 11 room 4-370

Teaching Assistants:

There is no final exam in 18.085. The first in-class exam will be on Friday October 5.

TEXTBOOKS: The Coop has textbooks now (and probably Quantum Books too)

The associated website for this book is math.mit.edu/cse.

Here is a link to the OpenCourseWare page for this course. It has problem sets and solutions from Fall 2007.

Please look at the appendix Linear Algebra in a Nutshell.

Course outline: (ps, pdf)

Exams dates: Oct 5 Nov 2 Dec 7
All exams in 18085 will be open book and notes. Use of computers is not allowed.

Exam 1 FRIDAY October 5 solutions
Exam 2 solutions
Exam 3 solutions

Exam 2 FRIDAY November 2 in Walker Memorial (top floor) at 11
(Please sit on the west side and the back of the east side so our exam can continue if necessary at noon)

Help sessions TUESDAY 4:30 - 5:45 in 2-190 and THURSDAY 4:30 - 5:45 in 2-142

EXTRA Help sessions THURSDAY November 1 from 4:30 in 2-190 (not 2-142 or 2-242)

MATLAB cheatsheet (1-page quick-reference guide) pdf

[top]

Homework

Homework 1: due 11am WED 9/12 (pdf) (solutions)

MATLAB 1: due 11am MON 9/17 (pdf)

Homework 2: due WEDNESDAY September 19, 2007 (MATLAB #1 due Monday)

1.3 7, 8
1.4 3, 5, 7, 12
1.5 1, 4, 7, 9

Homework 3 + MATLAB: due 11am WEDNESDAY September 26

1.6 1, 3, 6, 11, 16, 19, 24, 27
1.7 12
MATLAB problem: (pdf)

MATLAB 2.2: due 11am WED Oct 10 (pdf) *** Apologies -- we just realized that the leapfrog code isn't general enough, so the homework is reduced to the first two problems only !!!

Needs these three codes also on math.mit.edu/cse

Homework 4: due WED OCT 17

2.3 7, 8, 12, 18, 22
2.4 1, 3, 8, 12, 17

MATLAB Homework 4: due FRI OCT 19 (solutions)

2.3 parts 1 and 5 of Problem 3
2.4 problems 19 and 21 (use the note on that page to construct matrices)

Three codes to help with Homework 4:

Homework 5: due WED October 24

2.7 1, 2, 3, 5
3.1 1, 5, 10, 18

Selected solutions: A B

Homework 6: before exam 2 (this homework will not be collected)

2.3 19
2.4 9, 11
3.1 11, 17
3.2 1, 6

Selected solutions: 2.4.11 and 3.1.11 and 3.1.17

Homework 7: (including MATLAB) due FRI Nov 16 *** updated due date

3.4 4, 17, 19
3.5 1, 2
3.6 5, 6 (code on cse site), 9 , 16

Partial Solutions to Problem Set 7: PDF

Solution to Problem 3.6.9:

nodevalues=[1;0;0;0;0;0;0;0];
xy=[1 0  0  0   0   0   0   0;
    1 1  0  1   0   0   0   0;
    1 0  1  0   0   1   0   0;
    1 .5 0  .25 0   0   0   0;
    1 0  .5 0   0   .25 0   0;
    1 .5 1  .25 .5  1   .25 .5;
    1 1 .5  1   .5 .25  .5  .25;
    1 1  1  1   1   1   1   1];
a=inverse(xy)*nodevalues
% a=[1,-3,-3,2,5,2,-2,-2]

% plot the function phi(x,y) in the square, with those nodevalues

[x,y] = meshgrid(0:.1:1, 0:.1:1);
z = a(1) + a(2)*x + a(3)*y + a(4)*(x.^2) + a(5)*(x.*y) + a(6)*(y.^2) +
a(7)*((x.^2).*y) + a(8)*(x.*(y.^2));

surf(x,y,z);

Thoughts and hints on Problem 3.6.6, Poisson in a triangle T:

The key steps are the lists p and t and b = boundary nodes I think you can use the squaregrid construction on page 303, m=n=5 Those triangles may go up to the right (45 degrees) and why not let T do the same (corners at 0,0 and 1,0 and 1,1). Squaregrid includes a lot of meshpoints we don't want (a whole triangle of them) but I think they can be listed carefully in b -- then the code on p.304 will remove from K to leave Kb for a triangle. (I am not expecting major programming if this fix doesn't succeed...)

Graphical notes on the pagoda function

Homework 8: for WEDNESDAY NOV 29 Updated

4.1 3, 10, 13, 18
4.2 1, 2, 20
4.3 2, 6, 8, 10, 15

Solutions

Homework 9 (the last!) on convolution and Fourier integrals can be turned in with the exam on FRI Dec 7 (or it can be turned in MON Dec 10).

4.4 2, 8, 9 and also problem 4.6.1
4.5 1, 2, 10, 12, 22

The exam on FRI DEC 7 in Walker will cover homeworks 7-8-9

Gradient-Divergence-Laplace (powers of x+iy, 5-point difference matrix,...)
Fourier series, N-point DFT, convolution, Fourier integral (NOT 4.2)
Review in class WED DEC 5 (bring questions) and help sessions Tues/Thursday

[top]

Important Notes

Here is a finite element code with many comments -- I thought you might like to see how it works and try some examples by varying boundary conds. Homeworks to learn from, not to turn in! How does error drop from n to 2n? Moving on to Fourier after discussion Monday of finite elements.

This FEM code solves Poisson with f=1 on the square with standard mesh. You will be able to run the code as is/ boundary conditions easy to change HOMEWORK: USE SEVERAL m=n to find u at the center node (U=0 on boundary). HOMEWORK 2: Change boundary conditions (nodes in b) to make U=0 on *3* sides Solve again for U at the center // NOTICE how boundary conditions come last

2D Finite Element Code

% [p,t,b] = squaregrid(m,n) % create grid of N=mn nodes to be listed in p
% generate mesh of T=2(m-1)(n-1) right triangles in unit square
m=11; n=11; % includes boundary nodes, mesh spacing 1/(m-1) and 1/(n-1)
[x,y]=ndgrid((0:m-1)/(m-1),(0:n-1)/(n-1)); % matlab forms x and y lists
p=[x(:),y(:)]; % N by 2 matrix listing x,y coordinates of all N=mn nodes
t=[1,2,m+2; 1,m+2,m+1]; % 3 node numbers for two triangles in first square
t=kron(t,ones(m-1,1))+kron(ones(size(t)),(0:m-2)');
% now t lists 3 node numbers of 2(m-1) triangles in the first mesh row
t=kron(t,ones(n-1,1))+kron(ones(size(t)),(0:n-2)'*m);
% final t lists 3 node numbers of all triangles in T by 3 matrix 
b=[1:m,m+1:m:m*n,2*m:m:m*n,m*n-m+2:m*n-1]; % bottom, left, right, top 
% b = numbers of all 2m+2n **boundary nodes** preparing for U(b)=0

% [K,F] = assemble(p,t) % K and F for any mesh of triangles: linear phi's
N=size(p,1);T=size(t,1); % number of nodes, number of triangles
% p lists x,y coordinates of N nodes, t lists triangles by 3 node numbers
K=sparse(N,N); % zero matrix in sparse format: zeros(N) would be "dense"
F=zeros(N,1); % load vector F to hold integrals of phi's times load f(x,y)

for e=1:T  % integration over one triangular element at a time
  nodes=t(e,:); % row of t = node numbers of the 3 corners of triangle e
  Pe=[ones(3,1),p(nodes,:)]; % 3 by 3 matrix with rows=[1 xcorner ycorner] 
  Area=abs(det(Pe))/2; % area of triangle e = half of parallelogram area
  C=inv(Pe); % columns of C are coeffs in a+bx+cy to give phi=1,0,0 at nodes
  % now compute 3 by 3 Ke and 3 by 1 Fe for element e
  grad=C(2:3,:);Ke=Area*grad'*grad; % element matrix from slopes b,c in grad
  Fe=Area/3; % integral of phi over triangle is volume of pyramid: f(x,y)=1
  % multiply Fe by f at centroid for load f(x,y): one-point quadrature!
  % centroid would be mean(p(nodes,:)) = average of 3 node coordinates
  K(nodes,nodes)=K(nodes,nodes)+Ke; % add Ke to 9 entries of global K
  F(nodes)=F(nodes)+Fe; % add Fe to 3 components of load vector F
end   % all T element matrices and vectors now assembled into K and F

% [Kb,Fb] = dirichlet(K,F,b) % assembled K was singular! K*ones(N,1)=0
% Implement Dirichlet boundary conditions U(b)=0 at nodes in list b
K(b,:)=0; K(:,b)=0; F(b)=0; % put zeros in boundary rows/columns of K and F 
K(b,b)=speye(length(b),length(b)); % put I into boundary submatrix of K
Kb=K; Fb=F; % Stiffness matrix Kb (sparse format) and load vector Fb

% Solving for the vector U will produce U(b)=0 at boundary nodes
U=Kb\Fb;  % The FEM approximation is U_1 phi_1 + ... + U_N phi_N

% Plot the FEM approximation U(x,y) with values U_1 to U_N at the nodes 
trisurf(t,p(:,1),p(:,2),0*p(:,1),U,'edgecolor','k','facecolor','interp');
view(2),axis equal,colorbar

Gradient and Divergence / Parallel Table ( ps | pdf )

normalmodescode for Mu'' + Ku = 0

     % inputs M, K, uzero, vzero, t
     [vectors,values] = eig(K,M); eigen = diag(values); % solve Kx = (lambda)Mx
     A = vectors\uzero; B = (vectors*sqrt(values))\vzero;
     coeffs = A.*cos(t*sqrt(eigen)) + B.*sin(t*sqrt(eigen));
     u = vectors*coeffs;  % solution at time t to Mu'' + Ku = 0

improved GRID2Dcode from Professor Strang

     N=3; % N*N nodes and 2N*N - 2N edges in a square grid
     col=[-1; zeros(N*N - 1,1)]; % -1 on diagonal of AH
     row=[-1 1 zeros(1,N*N -2)]; % +1 next to the diagonal
     AH=toeplitz(col,row); % Incidence matrix/Horizontal edges
     AH(N:N:N^2,:)=[]; % Remove N nonexistent edges/end nodes
     ROW=[-1 zeros(1,N-1) 1 zeros(1,N*N-N-1)]; % off-diagonal 1
     COL=[-1; zeros(N*N-N-1,1)]; % -1 on diagonal of AV
     AV=toeplitz(COL,ROW); % Incidence matrix/Vertical edges
     A=[AH;AV]; % Combine horizontal and vertical edges into A
     ATA=A'*A; % Conductance matrix (singular) of order N*N
     norm(ATA*ones(N*N,1)) % Check that ATA(1;...;1)=0

     B=toeplitz([2 -1 zeros(1,N-2)]); B(1,1)=1; B(N,N)=1;
     fastATA=kron(B,eye(N)) + kron(eye(N),B); % 2D from 1D
     norm(ATA - fastATA) % Check that both ATA's are correct
     % Voltages 0 and 1 at nodes k and j/can remove columns j,k from A
     % Easier way ! Create a current source f between nodes j and k
     % Ground a node (which can be k) and find u=ATA\f and u(j)
     % This is the voltage needed at j for unit current from j to k
     ATA(:,k)=[]; ATA(k,:)=[]; % Ground node k to make ATA invertible
     f=zeros(N*N - 1); f(j)=1; u=ATA\f; u(j) % Expect 1/2 for neighbors

Short paper on infinite networks

Another code:

     A=zeros(2*N*(N-1),N*N);
     for i=1:N-1, % first row horizontal edges
         A(i,i)=-1;
         A(i,i+1)=1;
     end
     for i=1:N, % ith vertical edges
         for j=1:N-1,
             A(i+j*(2*N-1)-N,i+(j-1)*N)=-1;
             A(i+j*(2*N-1)-N,i+j*N)=1;
         end
     end
     for i=1:N-1, % horizontal edges from second row to last row
         for j=1:N-1,
             A(i+j*(2*N-1),i+j*N)=-1;
             A(i+j*(2*N-1),i+1+j*N)=1;
         end
     end

[top]