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

Department of Mathematics
Massachusetts Institute of Technology

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

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]