Programing the Finite Element Method with Matlab Jack Chessa∗ 3rd October 2002
1
Introduction
The goal of this document is to give a very brief overview and direction in the writing of finite element code using Matlab. It is assumed that the reader has a basic familiarity with the theory of the finite element method, and our attention will be mostly on the implementation. An example finite element code for analyzing static linear elastic problems written in Matlab is presented to illustrate how to program the finite element method. The example program and supporting files are available at http://www.tam.northwestern.edu/jfc795/Matlab/
1.1
Notation
For clarity we adopt the following notation in this paper; the bold italics font v denotes a vector quantity of dimension equal to the spacial dimension of the problem i.e. the displacement or velocity at a point, the bold nonitalicized font d denotes a vector or matrix which is of dimension of the number of unknowns in the discrete system i.e. a system matrix like the stiffness matrix, an uppercase subscript denotes a node number whereas a lowercase subscript in general denotes a vector component along a Cartesian unit vector. So, if d is the system vector of nodal unknowns, uI is a displacement vector of node I and uIi is the component of the displacement at node I in the i direction, or uI · ei . Often Matlab syntax will be intermixed with mathematical notation ∗
Graduate Research Assistant, Northwestern University (
[email protected])
1
which hopefully adds clarity to the explanation. The typewriter font, font, is used to indicate that Matlab syntax is being employed.
2
A Few Words on Writing Matlab Programs
The Matlab programming language is useful in illustrating how to program the finite element method due to the fact it allows one to very quickly code numerical methods and has a vast predefined mathematical library. This is also due to the fact that matrix (sparse and dense), vector and many linear algebra tools are already defined and the developer can focus entirely on the implementation of the algorithm not defining these data structures. The extensive mathematics and graphics functions further free the developer from the drudgery of developing these functions themselves or finding equivalent preexisting libraries. A simple two dimensional finite element program in Matlab need only be a few hundred lines of code whereas in Fortran or C++ one might need a few thousand. Although the Matlab programming language is very complete with respect to it’s mathematical functions there are a few finite element specific tasks that are helpful to develop as separate functions. These have been programed and are available at the previously mentioned web site. As usual there is a trade off to this ease of development. Since Matlab is an interpretive language; each line of code is interpreted by the Matlab command line interpreter and executed sequentially at run time, the run times can be much greater than that of compiled programming languages like Fortran or C++. It should be noted that the builtin Matlab functions are already compiled and are extremely efficient and should be used as much as possible. Keeping this slow down due to the interpretive nature of Matlab in mind, one programming construct that should be avoided at all costs is the for loop, especially nested for loops since these can make a Matlab programs run time orders of magnitude longer than may be needed. Often for loops can be eliminated using Matlab’s vectorized addressing. For example, the following Matlab code which sets the row and column of a matrix A to zero and puts one on the diagonal for i=1:size(A,2) A(n,i)=0; end for i=1:size(A,1) A(i,n)=0; end 2
A(n,n)=1; should never be used since the following code A(:,n)=0; A(:,n)=0; A(n,n)=0; does that same in three interpreted lines as opposed to nr +nc+1 interpreted lines, where A is a nr×nc dimensional matrix. One can easily see that this can quickly add significant overhead when dealing with large systems (as is often the case with finite element codes). Sometimes for loops are unavoidable, but it is surprising how few times this is the case. It is suggested that after developing a Matlab program, one go back and see how/if they can eliminate any of the for loops. With practice this will become second nature.
3
Sections of a Typical Finite Element Program
A typical finite element program consists of the following sections 1. Preprocessing section 2. Processing section 3. Postprocessing section In the preprocessing section the data and structures that define the particular problem statement are defined. These include the finite element discretization, material properties, solution parameters etc. . The processing section is where the finite element objects i.e. stiffness matrices, force vectors etc. are computed, boundary conditions are enforced and the system is solved. The postprocessing section is where the results from the processing section are analyzed. Here stresses may be calculated and data might be visualized. In this document we will be primarily concerned with the processing section. Many pre and postprocessing operations are already programmed in Matlab and are included in the online reference; if interested one can either look directly at the Matlab script files or type help ’function name’ at the Matlab command line to get further information on how to use these functions.
3
4
Finite Element Data Structures in Matlab
Here we discuss the data structures used in the finite element method and specifically those that are implemented in the example code. These are somewhat arbitrary in that one can imagine numerous ways to store the data for a finite element program, but we attempt to use structures that are the most flexible and conducive to Matlab. The design of these data structures may be depend on the programming language used, but usually are not significantly different than those outlined here.
4.1
Nodal Coordinate Matrix
Since we are programming the finite element method it is not unexpected that we need some way of representing the element discretization of the domain. To do so we define a set of nodes and a set of elements that connect these nodes in some way. The node coordinates are stored in the nodal coordinate matrix. This is simply a matrix of the nodal coordinates (imagine that). The dimension of this matrix is nn × sdim where nn is the number of nodes and sdim is the number of spacial dimensions of the problem. So, if we consider a nodal coordinate matrix nodes the ycoordinate of the nth node is nodes(n,2). Figure 1 shows a simple finite element discretization. For this simple mesh the nodal coordinate matrix would be as follows: 0.0 0.0 2.0 0.0 0.0 3.0 . nodes = (1) 2.0 3.0 0.0 6.0 2.0 6.0
4.2
Element Connectivity Matrix
The element definitions are stored in the element connectivity matrix. This is a matrix of node numbers where each row of the matrix contains the connectivity of an element. So if we consider the connectivity matrix elements that describes a mesh of 4node quadrilaterals the 36th element is defined by the connectivity vector elements(36,:) which for example may be [ 36 42 13 14] or that the elements connects nodes 36 → 42 → 13 → 14. So for
4
the simple mesh in Figure 1 the element connectivity matrix is 1 2 3 2 4 3 elements = 4 5 2 . 6 5 4
(2)
Note that the elements connectivities are all ordered in a counterclockwise fashion; if this is not done so some Jacobian’s will be negative and thus can cause the stiffnesses matrix to be singular (and obviously wrong!!!).
4.3
Definition of Boundaries
In the finite element method boundary conditions are used to either form force vectors (natural or Neumann boundary conditions) or to specify the value of the unknown field on a boundary (essential or Dirichlet boundary conditions). In either case a definition of the boundary is needed. The most versatile way of accomplishing this is to keep a finite element discretization of the necessary boundaries. The dimension of this mesh will be one order less that the spacial dimension of the problem (i.e. a 2D boundary mesh for a 3D problem, 1D boundary mesh for a 2D problem etc. ). Once again let’s consider the simple mesh in Figure 1. Suppose we wish to apply a boundary condition on the right edge of the mesh then the boundary mesh would be the defined by the following element connectivity matrix of 2node line elements 2 4 right Edge = . (3) 4 6 Note that the numbers in the boundary connectivity matrix refer to the same node coordinate matrix as do the numbers in the connectivity matrix of the interior elements. If we wish to apply an essential boundary conditions on this edge we need a list of the node numbers on the edge. This can be easily done in Matlab with the unique function. nodesOnBoundary = unique(rightEdge); This will set the vector nodesOnBoundary equal to [2 4 6]. If we wish to from a force vector from a natural boundary condition on this edge we simply loop over the elements and integrate the force on the edge just as we would integrate any finite element operators on the domain interior i.e. the stiffness matrix K.
5
4.4
Dof Mapping
Ultimately for all finite element programs we solve a linear algebraic system of the form Kd = f (4) for the vector d. The vector d contains the nodal unknowns for that define the finite element approximation h
u (x) =
nn X
NI (x)dI
(5)
I=1
where NI (x) are the finite element shape functions, dI are the nodal unknowns for the node I which may be scalar or vector quantities (if uh (x) is a scalar or vector) and nn is the number of nodes in the discretization. For scalar fields the location of the nodal unknowns in d is most obviously as follows dI = d(I), (6) but for vector fields the location of the nodal unknown dIi , where I refers to the node number and i refers to the component of the vector nodal unknown dI , there is some ambiguity. We need to define a mapping from the node number and vector component to the index of the nodal unknown vector d. This mapping can be written as f : {I, i} → n
(7)
where f is the mapping, I is the node number, i is the component and n is the index in d. So the location of unknown uIi in d is as follows uIi = df (I,i) .
(8)
There are two common mappings used. The first is to alternate between each spacial component in the nodal unknown vector d. With this arrangement the nodal unknown vector d is of the form u1x u1y .. . u2x u d= (9) .2y .. u nn x u nn y .. . 6
where nn is again the number of nodes in the discretization. This mapping is n = sdim(I − 1) + i. (10) With this mapping the i component of the displacement at node I is located as follows in d dIi = d( sdim*(I1) + i ). (11) The other option is to group all the like components of the nodal unknowns in a contiguous portion of d or as follows u1x u2x . . . d = unx (12) u1y u2y .. . The mapping in this case is n = (i − 1)nn + I
(13)
So for this structure the i component of the displacement at node I is located at in d dIi = d( (i1)*nn + I ) (14) For reasons that will be appreciated when we discuss the scattering of element operators into system operators we will adopt the latter dof mapping. It is important to be comfortable with these mappings since it is an operation that is performed regularly in any finite element code. Of course which ever mapping is chosen the stiffness matrix and force vectors should have the same structure.
5
Computation of Finite Element Operators
At the heart of the finite element program is the computation of finite element operators. For example in a linear static code they would be the stiffness matrix Z BT C B dΩ (15) K= Ω
7
and the external force vector f
ext
=
Z
Nt dΓ.
(16)
Γt
The global operators are evaluated by looping over the elements in the discretization, integrating the operator over the element and then to scatter the local element operator into the global operator. This procedure is written mathematically with the Assembly operator A Z K = Ae BeT C Be dΩ (17) Ωe
5.1
Quadrature
The integration of an element operator is performed with an appropriate quadrature rule which depends on the element and the function being integrated. In general a quadrature rule is as follows Z ξ=1 X f (ξ)dξ = f (ξq )Wq (18) ξ=−1
q
where f (ξ) is the function to be integrated, ξq are the quadrature points and Wq the quadrature weights. The function quadrature generates a vector of quadrature points and a vector of quadrature weights for a quadrature rule. The syntax of this function is as follows [quadWeights,quadPoints] = quadrature(integrationOrder, elementType,dimensionOfQuadrature); so an example quadrature loop to integrate the function f = x3 on a triangular element would be as follows [qPt,qWt]=quadrature(3,’TRIANGULAR’,2); for q=1:length(qWt) xi = qPt(q); % quadrature point % get the global coordinte x at the quadrature point xi % and the Jacobian at the quadrature point, jac ... f_int = f_int + x^3 * jac*qWt(q); end
8
5.2
Operator ”Scattering”
Once the element operator is computed it needs to be scattered into the global operator. An illustration of the scattering of an element force vector into a global force vector is shown in Figure 2. The scattering is dependent on the element connectivity and the dof mapping chosen. The following code performs the scatter indicated in Figure 2 elemConn = element(e,:); enn = length(elemConn); for I=1:enn; for i=1:2 Ii=nn*(i1)+sctr(I); f(Ii) = f(Ii) + f((i1)*enn+I); end end
% element connectivity % loop over element nodes % loop over spacial dimensions % dof map
but uses a nested for loop (bad bad bad). This is an even more egregious act considering the fact that it occurs within an element loop so this can really slow down the execution time of the program (by orders of magnitude in many cases). And it gets even worse when scattering a matrix operator (stiffness matrix) since we will have four nested for loops. Fortunately, Matlab allows for an easy solution; the following code performs exactly the same scattering as is done in the above code but with out any for loops, so the execution time is much improved (not to mention that it is much more concise). sctr = element(e,:); sctrVct = [ sctr sctr+nn ]; f(sctrVct) = f(sctrVct) + fe;
% element connectivity % vector scatter
To scatter an element stiffness matrix into a global stiffness matrix the following line does the trick K(sctrVct,sctrVct) = K(sctrVct,sctrVct) + ke; This terse array indexing of Matlab is a bit confusing at first but if one spends a bit of time getting used to it, it will become quite natural and useful.
5.3
Enforcement of Essential Boundary Conditions
The final issue before solving the linear algebraic system of finite element equations is the enforcement of the essential boundary conditions. Typically 9
this involves modifying the system Kd = f
(19)
so that the essential boundary condition dn = d¯n
(20)
is satisfied while retaining the original finite element equations on the unconstrained dofs. In (20) the subscript n refers to the index of the vector d not to a node number. An easy way to enforce (20) would be to modify nth row of the K matrix so that Knm = δnm
∀m ∈ {1, 2 . . . N }
(21)
where N is the dimension of K and setting fn = d¯n .
(22)
This reduces the nth equation of (19) to (20). Unfortunately, this destroys the symmetry of K which is a very important property for many efficient linear solvers. By modifying the nth column of K as follows Km,n = δnm
∀m ∈ {1, 2 . . . N }.
(23)
We can make the system symmetric. Of course this will modify every equation in (19) unless we modify the force vector f fm = Kmn d¯n .
(24)
If we write the modified k th equation in (19) Kk1 d1 + Kk2 d2 + . . . Kk(n−1) dn−1 + Kk(n+1) dn+1 + . . . + KkN dN = fk − Kkn d¯n
(25)
it can be seen that we have the same linear equations as in (19), but just with the internal force from the constrained dof. This procedure in Matlab i s as follows f = f  K(:,fixedDofs)*fixedDofValues; K(:,fixedDofs) = 0; K(fixedDofs,:) = 0; K(fixedDofs,fixedDofs) = bcwt*speye(length(fixedDofs)); f(fixedDofs) = bcwt*fixedDofValues; where fixedDofs is a vector of the indicies in d that are fixed, fixedDofValues is a vector of the values that fixedDofs are assigned to and bcwt is a weighing factor to retain the conditioning of the stiffness matrix (typically bcwt = trace(K)/N ). 10
6
Where To Go Next
Hopefully this extremely brief overview of programming simple finite element methods with Matlab has helped bridge the gap between reading the theory of the finite element method and sitting down and writing ones own finite element code. The examples in the Appendix should be looked at and run, but also I would suggest trying to write a simple 1D or 2D finite element code from scratch to really solidify the method in ones head. The examples can then be used as a reference to diminish the struggle. Good Luck!
11
A
Installation of Example Matlab Program
All the functions needed to run the example programs as well as the examples themselves can be found at http://www.tam.northwestern.edu/jfc795/Matlab/ I believe that the following files are required, but if one gets a run error about function not found chances are that I forgot to list it here but it is in one of the Matlab directories at the above web site. • MeshGenerationsquare node array.m: generates an array of nodes in 2D • MeshGenerationmake elem.m: generates elements on an array of nodes • MeshGenerationmsh2mlab.m: reads in a Gmsh file • MeshGenerationplot mesh.m: plots a finite element mesh • PostProcessingplot field.m: plots a finite element field • quadrature.m: returns various quadrature rules • lagrange basis.m: return the shape functions and gradients of the shape functions in the parent coordinate system for various elements There are many additional files that one might find useful and an interested individual can explore these on there own. These fies should be copied either the directory which contains the example script file or into a directory that is in the Matlab search path.
12
B
Example: Beam Bending Problem
The first example program solves the static bending of a linear elastic beam. The configuration of the problem is shown in Figure 3 and the program can be found at http://www.tam.northwestern.edu/jfc795/Matlab/ Examples/Static/beam.m The exact solution for this problem is as follows σ11 = −
P (L − x)y I
σ22 = 0 P 2 σ12 = (c − y 2 ) 2I Py 3 L2 − (L − x)2 + (2 + ν)(y 2 − c2 )) u1 = − 6EI Py 3 (L − x)3 − L3 − (4 + 5ν)c2 + 3L2 x + 3ν(L − x)y 2 u2 = 6EI This problem can be run with three element types; three node triangle element, a four node quadrilateral element and a nine node quadrilateral element. Also, one can choose between plane strain or plane stress assumption. % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
beam.m Solves a linear elastic 2D beam problem ( plane stress or strain ) with several element types. ^ y      > x  2c    L  with the boundary following conditions: u_x = 0 at (0,0), (0,c) and (0,c) u_y = 0 at (0,0) t_x = y t_y = P*(x^2c^2)
along the edge x=0 along the edge x=L
****************************************************************************** This file and the supporting matlab files can be found at http://www.tam.northwestern.edu/jfc795/Matlab by Jack Chessa Northwestern University
13
% % ****************************************************************************** clear colordef black state = 0; % ****************************************************************************** % *** I N P U T *** % ****************************************************************************** tic; disp(’************************************************’) disp(’*** S T A R T I N G R U N ***’) disp(’************************************************’) disp([num2str(toc),’
START’])
% MATERIAL PROPERTIES E0 = 10e7; % Young’s modulus nu0 = 0.30; % Poisson’s ratio % BEAM PROPERTIES L = 16; % length of the beam c = 2; % the distance of the outer fiber of the beam from the midline % MESH PROPERTIES elemType = ’Q9’; % % % % numy = 4; numx = 18; plotMesh = 1;
% % % %
the element type used in the FEM simulation; ’T3’ is for a three node constant strain triangular element, ’Q4’ is for a four node quadrilateral element, and ’Q9’ is for a nine node quadrilateral element. the number of elements in the xdirection (beam length) and in the ydireciton. A flag that if set to 1 plots the initial mesh (to make sure that the mesh is correct)
% TIP LOAD P = 1; % the peak magnitude of the traction at the right edge % STRESS ASSUMPTION stressState=’PLANE_STRESS’; % set to either ’PLANE_STRAIN’ or "PLANE_STRESS’ % nuff said. % ****************************************************************************** % *** P R E  P R O C E S S I N G *** % ****************************************************************************** I0=2*c^3/3; % the second polar moment of inertia of the beam crosssection. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE ELASTICITY MATRIX if ( strcmp(stressState,’PLANE_STRESS’) ) % Plane Strain case C=E0/(1nu0^2)*[ 1 nu0 0; nu0 1 0; 0 0 (1nu0)/2 ]; else % Plane Strain case C=E0/(1+nu0)/(12*nu0)*[ 1nu0 nu0 0; nu0 1nu0 0; 0 0 1/2nu0 ]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % GENERATE FINITE ELEMENT MESH %
14
% % % % % % % % % % % % % % % % % % % % %
Here we gnerate the finte element mesh (using the approriate elements). I won’t go into too much detail about how to use these functions. If one is interested one can type  help ’function name’ at the matlab comand line to find out more about it. The folowing data structures are used to describe the finite element discretization: node  is a matrix of the node coordinates, i.e. node(I,j) > x_Ij element  is a matrix of element connectivities, i.e. the connectivity of element e is given by > element(e,:) > [n1 n2 n3 ...]; To apply boundary conditions a description of the boundaries is needed. To accomplish this we use a separate finite element discretization for each boundary. For a 2D problem the boundary discretization is a set of 1D elements. rightEdge leftEdge 
a element connectivity matrix for the right edge I’ll give you three guesses
These connectivity matricies refer to the node numbers defined in the coordinate matrix node.
disp([num2str(toc),’ GENERATING MESH’]) switch elemType case ’Q4’ % here we generate the mesh of Q4 elements nnx=numx+1; nny=numy+1; node=square_node_array([0 c],[L c],[L c],[0 c],nnx,nny); inc_u=1; inc_v=nnx; node_pattern=[ 1 2 nnx+2 nnx+1 ]; element=make_elem(node_pattern,numx,numy,inc_u,inc_v); case ’Q9’ % here we generate a mehs of Q9 elements nnx=2*numx+1; nny=2*numy+1; node=square_node_array([0 c],[L c],[L c],[0 c],nnx,nny); inc_u=2; inc_v=2*nnx; node_pattern=[ 1 3 2*nnx+3 2*nnx+1 2 nnx+3 2*nnx+2 nnx+1 nnx+2 ]; element=make_elem(node_pattern,numx,numy,inc_u,inc_v); otherwise %’T3’ % and last but not least T3 elements nnx=numx+1; nny=numy+1; node=square_node_array([0 c],[L c],[L c],[0 c],nnx,nny); node_pattern1=[ 1 2 nnx+1 ]; node_pattern2=[ 2 nnx+2 nnx+1 ]; inc_u=1; inc_v=nnx; element=[make_elem(node_pattern1,numx,numy,inc_u,inc_v); make_elem(node_pattern2,numx,numy,inc_u,inc_v) ]; end % DEFINE BOUNDARIES
15
% Here we define the boundary discretizations. uln=nnx*(nny1)+1; % upper left node number urn=nnx*nny; % upper right node number lrn=nnx; % lower right node number lln=1; % lower left node number cln=nnx*(nny1)/2+1; % node number at (0,0) switch elemType case ’Q9’ rightEdge=[ lrn:2*nnx:(uln1); (lrn+2*nnx):2*nnx:urn; (lrn+nnx):2*nnx:urn ]’; leftEdge =[ uln:2*nnx:(lrn+1); (uln2*nnx):2*nnx:1; (ulnnnx):2*nnx:1 ]’; edgeElemType=’L3’; otherwise % same discretizations for Q4 and T3 meshes rightEdge=[ lrn:nnx:(uln1); (lrn+nnx):nnx:urn ]’; leftEdge =[ uln:nnx:(lrn+1); (ulnnnx):nnx:1 ]’; edgeElemType=’L2’; end % GET NODES ON DISPLACEMENT BOUNDARY % Here we get the nodes on the essential boundaries fixedNodeX=[uln lln cln]’; % a vector of the node numbers which are fixed in % the x direction fixedNodeY=[cln]’; % a vector of node numbers which are fixed in % the ydirection uFixed=zeros(size(fixedNodeX)); vFixed=zeros(size(fixedNodeY));
% a vector of the xdisplacement for the nodes % in fixedNodeX ( in this case just zeros ) % and the ydisplacements for fixedNodeY
numnode=size(node,1); % number of nodes numelem=size(element,1); % number of elements % PLOT MESH if ( plotMesh ) % if plotMesh==1 we will plot the mesh clf plot_mesh(node,element,elemType,’g.’); hold on plot_mesh(node,rightEdge,edgeElemType,’bo’); plot_mesh(node,leftEdge,edgeElemType,’bo’); plot(node(fixedNodeX,1),node(fixedNodeX,2),’r>’); plot(node(fixedNodeY,1),node(fixedNodeY,2),’r^’); axis off axis([0 L c c]) disp(’(paused)’) pause end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DEFINE SYSTEM DATA STRUCTURES % % Here we define the system data structures % U  is vector of the nodal displacements it is of length 2*numnode. The % displacements in the xdirection are in the top half of U and the % ydisplacements are in the lower half of U, for example the displacement % in the ydirection for node number I is at U(I+numnode) % f  is the nodal force vector. It’s structure is the same as U, % i.e. f(I+numnode) is the force in the y direction at node I % K  is the global stiffness matrix and is structured the same as with U and f % so that K_IiJj is at K(I+(i1)*numnode,J+(j1)*numnode) disp([num2str(toc),’ INITIALIZING DATA STRUCTURES’]) U=zeros(2*numnode,1); % nodal displacement vector
16
f=zeros(2*numnode,1); % external load vector K=sparse(2*numnode,2*numnode); % stiffness matrix % a vector of indicies that quickly address the x and y portions of the data % strtuctures so U(xs) returns U_x the nodal xdisplacements xs=1:numnode; % x portion of u and v vectors ys=(numnode+1):2*numnode; % y portion of u and v vectors % ****************************************************************************** % *** P R O C E S S I N G *** % ****************************************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE EXTERNAL FORCES % integrate the tractions on the left and right edges disp([num2str(toc),’ COMPUTING EXTERNAL LOADS’]) switch elemType % define quadrature rule case ’Q9’ [W,Q]=quadrature( 4, ’GAUSS’, 1 ); % four point quadrature otherwise [W,Q]=quadrature( 3, ’GAUSS’, 1 ); % three point quadrature end % RIGHT EDGE for e=1:size(rightEdge,1) % loop over the elements in the right edge sctr=rightEdge(e,:); sctrx=sctr; sctry=sctrx+numnode;
% scatter vector for the element % x scatter vector % y scatter vector
for q=1:size(W,1) pt=Q(q,:); wt=W(q); [N,dNdxi]=lagrange_basis(edgeElemType,pt); J0=dNdxi’*node(sctr,:); detJ0=norm(J0); yPt=N’*node(sctr,2); fyPt=P*(c^2yPt^2)/(2*I0); f(sctry)=f(sctry)+N*fyPt*detJ0*wt;
% % % % % %
quadrature loop quadrature point quadrature weight element shape functions element Jacobian determiniat of jacobian
% y coordinate at quadrature point % y traction at quadrature point % scatter force into global force vector
end % of quadrature loop end % of element loop % LEFT EDGE for e=1:size(leftEdge,1)
% loop over the elements in the left edge
sctr=rightEdge(e,:); sctrx=sctr; sctry=sctrx+numnode; for q=1:size(W,1) pt=Q(q,:); wt=W(q); [N,dNdxi]=lagrange_basis(edgeElemType,pt); J0=dNdxi’*node(sctr,:); detJ0=norm(J0); yPt=N’*node(sctr,2); fyPt=P*(c^2yPt^2)/(2*I0); fxPt=P*L*yPt/I0;
% % % % % %
quadrature loop quadrature point quadrature weight element shape functions element Jacobian determiniat of jacobian
% y traction at quadrature point % x traction at quadrature point
17
f(sctry)=f(sctry)+N*fyPt*detJ0*wt; f(sctrx)=f(sctrx)+N*fxPt*detJ0*wt; end % of quadrature loop end % of element loop % set the force at the nodes on the top and bottom edges to zero (traction free) % TOP EDGE topEdgeNodes = find(node(:,2)==c); % finds nodes on the top edge f(topEdgeNodes)=0; f(topEdgeNodes+numnode)=0; % BOTTOM EDGE bottomEdgeNodes = find(node(:,2)==c); % finds nodes on the bottom edge f(bottomEdgeNodes)=0; f(bottomEdgeNodes+numnode)=0; %%%%%%%%%%%%%%%%%%%%% COMPUTE STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp([num2str(toc),’ COMPUTING STIFFNESS MATRIX’]) switch elemType % define quadrature rule case ’Q9’ [W,Q]=quadrature( 4, ’GAUSS’, 2 ); % 4x4 Gaussian quadrature case ’Q4’ [W,Q]=quadrature( 2, ’GAUSS’, 2 ); % 2x2 Gaussian quadrature otherwise [W,Q]=quadrature( 1, ’TRIANGULAR’, 2 ); % 1 point triangural quadrature end for e=1:numelem
% start of element loop
sctr=element(e,:); % element scatter vector sctrB=[ sctr sctr+numnode ]; % vector that scatters a B matrix nn=length(sctr); for q=1:size(W,1) pt=Q(q,:); wt=W(q); [N,dNdxi]=lagrange_basis(elemType,pt); J0=node(sctr,:)’*dNdxi; invJ0=inv(J0); dNdx=dNdxi*invJ0;
% % % %
quadrature loop quadrature point quadrature weight element shape functions
% element Jacobian matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE B MATRIX % _ _ %  N_1,x N_2,x ... 0 0 ...  % B =  0 0 ... N_1,y N_2,y ...  %  N_1,y N_2,y ... N_1,x N_2,x ...  % B=zeros(3,2*nn); B(1,1:nn) = dNdx(:,1)’; B(2,nn+1:2*nn) = dNdx(:,2)’; B(3,1:nn) = dNdx(:,2)’; B(3,nn+1:2*nn) = dNdx(:,1)’; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE ELEMENT STIFFNESS AT QUADRATURE POINT K(sctrB,sctrB)=K(sctrB,sctrB)+B’*C*B*W(q)*det(J0); end
% of quadrature loop
18
end % of element loop %%%%%%%%%%%%%%%%%%% END OF STIFFNESS MATRIX COMPUTATION %%%%%%%%%%%%%%%%%%%%%% % APPLY ESSENTIAL BOUNDARY CONDITIONS disp([num2str(toc),’ APPLYING BOUNDARY CONDITIONS’]) bcwt=mean(diag(K)); % a measure of the average size of an element in K % used to keep the conditioning of the K matrix udofs=fixedNodeX; % global indecies of the fixed x displacements vdofs=fixedNodeY+numnode; % global indecies of the fixed y displacements f=fK(:,udofs)*uFixed; f=fK(:,vdofs)*vFixed; f(udofs)=uFixed; f(vdofs)=vFixed;
% modify the force vector
K(udofs,:)=0; % zero out the rows and columns of the K matrix K(vdofs,:)=0; K(:,udofs)=0; K(:,vdofs)=0; K(udofs,udofs)=bcwt*speye(length(udofs)); % put ones*bcwt on the diagonal K(vdofs,vdofs)=bcwt*speye(length(vdofs)); % SOLVE SYSTEM disp([num2str(toc),’ U=K\f;
SOLVING SYSTEM’])
%****************************************************************************** %*** P O S T  P R O C E S S I N G *** %****************************************************************************** % % Here we plot the stresses and displacements of the solution. As with the % mesh generation section we don’t go into too much detail  use help % ’function name’ to get more details. disp([num2str(toc),’
POSTPROCESSING’])
dispNorm=L/max(sqrt(U(xs).^2+U(ys).^2)); scaleFact=0.1*dispNorm; fn=1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PLOT DEFORMED DISPLACEMENT PLOT figure(fn) clf plot_field(node+scaleFact*[U(xs) U(ys)],element,elemType,U(ys)); hold on plot_mesh(node+scaleFact*[U(xs) U(ys)],element,elemType,’g.’); plot_mesh(node,element,elemType,’w’); colorbar fn=fn+1; title(’DEFORMED DISPLACEMENT IN YDIRECTION’) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE STRESS stress=zeros(numelem,size(element,2),3); switch elemType % define quadrature rule case ’Q9’ stressPoints=[1 1;1 1;1 1;1 1;0 1;1 0;0 1;1 0;0 0 ]; case ’Q4’ stressPoints=[1 1;1 1;1 1;1 1]; otherwise
19
stressPoints=[0 0;1 0;0 1]; end for e=1:numelem
% start of element loop
sctr=element(e,:); sctrB=[sctr sctr+numnode]; nn=length(sctr); for q=1:nn pt=stressPoints(q,:); [N,dNdxi]=lagrange_basis(elemType,pt); J0=node(sctr,:)’*dNdxi; invJ0=inv(J0); dNdx=dNdxi*invJ0;
% stress point % element shape functions % element Jacobian matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE B MATRIX B=zeros(3,2*nn); B(1,1:nn) = dNdx(:,1)’; B(2,nn+1:2*nn) = dNdx(:,2)’; B(3,1:nn) = dNdx(:,2)’; B(3,nn+1:2*nn) = dNdx(:,1)’; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE ELEMENT STRAIN AND STRESS AT STRESS POINT strain=B*U(sctrB); stress(e,q,:)=C*strain; end end % of element loop stressComp=1; figure(fn) clf plot_field(node+scaleFact*[U(xs) U(ys)],element,elemType,stress(:,:,stressComp)); hold on plot_mesh(node+scaleFact*[U(xs) U(ys)],element,elemType,’g.’); plot_mesh(node,element,elemType,’w’); colorbar fn=fn+1; title(’DEFORMED STRESS PLOT, BENDING COMPONENT’) %print(fn,’djpeg90’,[’beam_’,elemType,’_sigma’,num2str(stressComp),’.jpg’]) disp([num2str(toc),’ RUN FINISHED’]) % *************************************************************************** % *** E N D O F P R O G R A M *** % *************************************************************************** disp(’************************************************’) disp(’*** E N D O F R U N ***’) disp(’************************************************’)
20
C
Example: Modal Analysis of an Atomic Force Microscopy (AFM) Tip
The program presented here is found at http://www.tam.northwestern.edu/jfc795/Matlab/Examples /Static/modal afm.m In addition the mesh file afm.msh is needed. This mesh file is produced using the GPL program Gmsh which is available at http://www.geuz.org/gmsh/ This program is not needed to run this program, only the *.msh file is needed, but it is a very good program for generating finite element meshes. In this example we perform a linear modal analysis of the AFM tip shown in Figure reffig:afm. This involves computing the mass and stiffness matrix and solving the following Eigenvalue problem K − ωn2 M an = 0 (26) for the natural frequencies ωn and the corresponding mode shapes an . Here the AFM tip is modeled with eight node brick elements and we assume that the feet of the AFM tip are fixed. % modal_afm.m % % by Jack Chessa % Northwestern University % clear colordef black state = 0; %****************************************************************************** %*** I N P U T *** %****************************************************************************** tic; disp(’************************************************’) disp(’*** S T A R T I N G R U N ***’) disp(’************************************************’) disp([num2str(toc),’
START’])
% MATERIAL PROPERTIES E0 = 160; % Youngs modulus in GPa nu0 = 0.27; % Poisson ratio rho = 2.330e9; % density in 10e12 Kg/m^3 % MESH PARAMETERS quadType=’GAUSS’; quadOrder=2; % GMSH PARAMETERS fileName=’afm.msh’; domainID=50;
21
fixedID=51; topID=52; % EIGENPROBELM SOLUTION PARAMETERS numberOfModes=8; % number of modes to compute consistentMass=0; % use a consistent mass matrix fixedBC=1; % use fixed or free bcs %****************************************************************************** %*** P R E  P R O C E S S I N G *** %****************************************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % READ GMSH FILE disp([num2str(toc),’ READING GMSH FILE’]) [node,elements,elemType]=msh2mlab(fileName); [node,elements]=remove_free_nodes(node,elements); element=elements{domainID}; element=brickcheck(node,element,1); if ( fixedBC ) fixedEdge=elements{fixedID}; else fixedEdge=[]; end topSurface=elements{topID}; plot_mesh(node,element,elemType{domainID},’r’) disp([num2str(toc),’ INITIALIZING DATA STRUCTURES’]) numnode=size(node,1); % number of nodes numelem=size(element,1); % number of elements % GET NODES ON DISPLACEMENT BOUNDARY fixedNodeX=unique(fixedEdge); fixedNodeY=fixedNodeX; fixedNodeZ=fixedNodeX; uFixed=zeros(size(fixedNodeX)); vFixed=zeros(size(fixedNodeY)); wFixed=zeros(size(fixedNodeZ));
% displacement for fixed nodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE COMPLIANCE MATRIX C=zeros(6,6); C(1:3,1:3)=E0/(1+nu0)/(12*nu0)*[ 1nu0 nu0 nu0; nu0 1nu0 nu0; nu0 nu0 1nu0 ]; C(4:6,4:6)=E0/(1+nu0)*eye(3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DEFINE SYSTEM DATA STRUCTURES K=sparse(3*numnode,3*numnode); % stiffness matrix if ( consistentMass) M=sparse(3*numnode,3*numnode); % mass matrix else M=zeros(3*numnode,1); % mass vector end %****************************************************************************** %*** P R O C E S S I N G *** %******************************************************************************
22
%%%%%%%%%%%%%%%%%%%%% COMPUTE SYSTEM MATRICIES %%%%%%%%%%%%%%%%%%%%%%%%%%%% disp([num2str(toc),’ COMPUTING STIFFNESS AND MASS MATRIX’]) [W,Q]=quadrature(quadOrder,quadType,3); % define quadrature rule et=elemType{domainID}; nn=size(element,2); for e=1:numelem % start of element loop sctr=element(e,:); % element scatter vector sctrB0=[ sctr sctr+numnode sctr+2*numnode ]; % scatters a B matrix for q=1:size(W,1) pt=Q(q,:); wt=W(q); [N,dNdxi]=lagrange_basis(et,pt); J0=node(sctr,:)’*dNdxi; invJ0=inv(J0); dNdx=dNdxi*invJ0; detJ0=det(J0);
% quadrature loop % quadrature point % quadrature weight % element shape functions % element Jacobian matrix
if (detJ0 <= 0) disp([’ERROR: NEGATIVE JACOBIAN IN ELEMENT ’,num2str(e)]); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE B MATRIX B0=zeros(6,3*nn); B0(1,1:nn) = dNdx(:,1)’; B0(2,nn+1:2*nn) = dNdx(:,2)’; B0(3,2*nn+1:3*nn) = dNdx(:,3)’; B0(4,2*nn+1:3*nn) = dNdx(:,2)’; B0(4,nn+1:2*nn) = dNdx(:,3)’; B0(5,1:nn) = dNdx(:,3)’; B0(5,2*nn+1:3*nn) = dNdx(:,1)’; B0(6,nn+1:2*nn) B0(6,1:nn)
= dNdx(:,1)’; = dNdx(:,2)’;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE ELEMENT STIFFNESS AT QUADRATURE POINT K(sctrB0,sctrB0)=K(sctrB0,sctrB0)+B0’*C*B0*wt*detJ0; % COMPUTE ELEMENT MASS AT QUADRATURE POINT mQPt=N*rho*N’*wt*detJ0; if ( ~consistentMass ) mQPt=sum(mQPt)’; M(sctr) = M(sctr)+mQPt; M(sctr+numnode) = M(sctr+numnode)+mQPt; M(sctr+2*numnode) = M(sctr+2*numnode)+mQPt; else M(sctr,sctr) = M(sctr,sctr)+mQPt; M(sctr+numnode,sctr+numnode) = M(sctr+numnode,sctr+numnode)+mQPt; M(sctr+2*numnode,sctr+2*numnode) = M(sctr+2*numnode,sctr+2*numnode)+mQPt; end end % of quadrature loop end % of element loop %%%%%%%%%%%%%%%%%%% END OF SYSTEM MATRIX COMPUTATION %%%%%%%%%%%%%%%%%%%%%% % ELIMINATE FIXED DOFS FROM EIGENVALUE COMUTATION disp([num2str(toc),’ FINDING ACTIVE DOFS’])
23
activeDof=setdiff([1:numnode]’,[fixedNodeX;fixedNodeY;fixedNodeZ]); activeDof=[activeDof;activeDof+numnode;activeDof+2*numnode]; % SOLVE SYSTEM disp([num2str(toc),’ SOLVING EIGEN PROBLEM’]) if ( consistentMass ) [modeShape,freq]=eigs(K(activeDof,activeDof),M(activeDof,activeDof),... numberOfModes,0); else Minv=spdiags(1./M,0,3*numnode,3*numnode); K=Minv*K; [modeShape,freq]=eigs(K(activeDof,activeDof),numberOfModes,0); end freq=diag(freq)/(2*pi); % frequency in kHz %****************************************************************************** %*** P O S T  P R O C E S S I N G *** %****************************************************************************** disp([num2str(toc),’ POSTPROCESSING’]) disp([’THE MODE FREQUENCIES ARE:’]) for m=1:length(freq) disp([’ MODE: ’,num2str(m),’ ’,num2str(freq(m))]) % PLOT MODE SHAPE figure(m); clf; U=zeros(numnode,1); U(activeDof)=modeShape(:,m); scaleFactor=20/max(abs(U)); plot_field(node+[U(1:numnode) U(numnode+1:2*numnode) U(2*numnode+1:3*numnode)]*scaleFactor,topSurface,elemType{topID},... ones(3*numnode,1)); hold on plot_mesh(node+[U(1:numnode) U(numnode+1:2*numnode) U(2*numnode+1:3*numnode)]*scaleFactor,topSurface,elemType{topID},’k’); plot_mesh(node,topSurface,elemType{topID},’r’); title([’MODE ’,num2str(m),’, FREQUENCY = ’,num2str(freq(m)),’ [kHz]’]) view(37,36) axis off print(m, ’djpeg90’, [’afm_mode_’,num2str(m),’.jpg’]); end % ANIMATE MODE nCycles=5; % number of cycles to animate fpc=10; % frames per cycle fact=sin(linspace(0,2*pi,fpc)); m=input(’What mode would you like to animate (type 0 to exit) ’); while ( m~=0 ) U=zeros(numnode,1); U(activeDof)=modeShape(:,m); wt=20/max(abs(U)); for i=1:fpc scaleFactor=fact(i)*wt; figure(length(freq+1)); clf; plot_field(node+[U(1:numnode) U(numnode+1:2*numnode) U(2*numnode+1:3*numnode)]*scaleFactor,topSurface,elemType{topID},... ones(3*numnode,1)); hold on plot_mesh(node+[U(1:numnode) U(numnode+1:2*numnode) U(2*numnode+1:3*numnode)]*scaleFactor,topSurface,elemType{topID},’k’);
24
plot_mesh(node,topSurface,elemType{topID},’w’); hold on view(37,36) axis([70 240 30 160 10 10]) title([’MODE ’,num2str(m),’, FREQUENCY = ’,num2str(freq(m)),’ [kHz]’]) axis off film(i)=getframe; end movie(film,nCycles); m=input(’What mode would you like to animate (type 0 to exit) ’); if ( m > length(freq) ) disp([’mode must be less than ’,num2str(length(freq))]) end end disp([num2str(toc),’ RUN FINISHED’]) % *************************************************************************** % *** E N D O F P R O G R A M *** % *************************************************************************** disp(’************************************************’) disp(’*** E N D O F R U N ***’) disp(’************************************************’) % compute uexact
25
D
Common Matlab Functions
Here is a quick list of some built in Matlab functions. These discriptions are availible by using the help function in Matlab. >> help HELP topics: matlab/general matlab/ops matlab/lang matlab/elmat matlab/specmat matlab/elfun matlab/specfun matlab/matfun matlab/datafun matlab/polyfun matlab/funfun matlab/sparfun matlab/plotxy matlab/plotxyz matlab/graphics matlab/color matlab/sounds matlab/strfun matlab/iofun matlab/demos toolbox/chem toolbox/control fdident/fdident fdident/fddemos toolbox/hispec toolbox/ident toolbox/images toolbox/local toolbox/mmle3 mpc/mpccmds mpc/mpcdemos mutools/commands

General purpose commands. Operators and special characters. Language constructs and debugging. Elementary matrices and matrix manipulation. Specialized matrices. Elementary math functions. Specialized math functions. Matrix functions  numerical linear algebra. Data analysis and Fourier transform functions. Polynomial and interpolation functions. Function functions  nonlinear numerical methods. Sparse matrix functions. Two dimensional graphics. Three dimensional graphics. General purpose graphics functions. Color control and lighting model functions. Sound processing functions. Character string functions. Lowlevel file I/O functions. The MATLAB Expo and other demonstrations. Chemometrics Toolbox Control System Toolbox. Frequency Domain System Identification Toolbox Demonstrations for the FDIDENT Toolbox HiSpec Toolbox System Identification Toolbox. Image Processing Toolbox. Local function library. MMLE3 Identification Toolbox. Model Predictive Control Toolbox Model Predictive Control Toolbox MuAnalysis and Synthesis Toolbox.: Commands directory 26
mutools/subs toolbox/ncd nnet/nnet nnet/nndemos toolbox/optim toolbox/robust toolbox/signal toolbox/splines toolbox/stats toolbox/symbolic toolbox/wavbox simulink/simulink simulink/blocks simulink/simdemos toolbox/codegen
 MuAnalysis and Synthesis Toolbox  Supplement  Nonlinear Control Design Toolbox.  Neural Network Toolbox.  Neural Network Demonstrations and Applications.  Optimization Toolbox.  Robust Control Toolbox.  Signal Processing Toolbox.  Spline Toolbox.  Statistics Toolbox.  Symbolic Math Toolbox.  (No table of contents file)  SIMULINK model analysis and construction functions.  SIMULINK block library.  SIMULINK demonstrations and samples.  RealTime Workshop
For more help on directory/topic, type "help topic". >> help elmat Elementary matrices and matrix manipulation. Elementary matrices. zeros  Zeros matrix. ones  Ones matrix. eye  Identity matrix. rand  Uniformly distributed random numbers. randn  Normally distributed random numbers. linspace  Linearly spaced vector. logspace  Logarithmically spaced vector. meshgrid  X and Y arrays for 3D plots. :  Regularly spaced vector. Special variables and constants. ans  Most recent answer. eps  Floating point relative accuracy. realmax  Largest floating point number. realmin  Smallest positive floating point number. pi  3.1415926535897.... i, j  Imaginary unit. inf  Infinity. 27
NaN flops nargin nargout computer isieee isstudent why version

Time and dates. clock cputime date etime tic, toc 
NotaNumber. Count of floating point operations. Number of function input arguments. Number of function output arguments. Computer type. True for computers with IEEE arithmetic. True for the Student Edition. Succinct answer. MATLAB version number.
Wall clock. Elapsed CPU time. Calendar. Elapsed time function. Stopwatch timer functions.
Matrix manipulation. diag  Create or extract diagonals. fliplr  Flip matrix in the left/right direction. flipud  Flip matrix in the up/down direction. reshape  Change size. rot90  Rotate matrix 90 degrees. tril  Extract lower triangular part. triu  Extract upper triangular part. :  Index into matrix, rearrange matrix. >> help specmat Specialized matrices. compan gallery hadamard hankel hilb invhilb kron magic pascal rosser

Companion matrix. Several small test matrices. Hadamard matrix. Hankel matrix. Hilbert matrix. Inverse Hilbert matrix. Kronecker tensor product. Magic square. Pascal matrix. Classic symmetric eigenvalue test problem. 28
toeplitz vander wilkinson
 Toeplitz matrix.  Vandermonde matrix.  Wilkinson’s eigenvalue test matrix.
>> help elfun Elementary math functions. Trigonometric. sin sinh asin asinh cos cosh acos acosh tan tanh atan atan2 atanh sec sech asec asech csc csch acsc acsch cot coth acot acoth 
Sine. Hyperbolic sine. Inverse sine. Inverse hyperbolic sine. Cosine. Hyperbolic cosine. Inverse cosine. Inverse hyperbolic cosine. Tangent. Hyperbolic tangent. Inverse tangent. Four quadrant inverse tangent. Inverse hyperbolic tangent. Secant. Hyperbolic secant. Inverse secant. Inverse hyperbolic secant. Cosecant. Hyperbolic cosecant. Inverse cosecant. Inverse hyperbolic cosecant. Cotangent. Hyperbolic cotangent. Inverse cotangent. Inverse hyperbolic cotangent.
Exponential. exp log log10 sqrt
Exponential. Natural logarithm. Common logarithm. Square root.

29
Complex. abs angle conj imag real

Absolute value. Phase angle. Complex conjugate. Complex imaginary part. Complex real part.
Numeric. fix floor ceil round rem sign

Round towards zero. Round towards minus infinity. Round towards plus infinity. Round towards nearest integer. Remainder after division. Signum function.
>> help specfun Specialized math functions. besselj bessely besseli besselk beta betainc betaln ellipj ellipke erf erfc erfcx erfinv expint gamma gcd gammainc lcm legendre gammaln log2 pow2

Bessel function of the first kind. Bessel function of the second kind. Modified Bessel function of the first kind. Modified Bessel function of the second kind. Beta function. Incomplete beta function. Logarithm of beta function. Jacobi elliptic functions. Complete elliptic integral. Error function. Complementary error function. Scaled complementary error function. Inverse error function. Exponential integral function. Gamma function. Greatest common divisor. Incomplete gamma function. Least common multiple. Associated Legendre function. Logarithm of gamma function. Dissect floating point numbers. Scale floating point numbers. 30
rat rats cart2sph cart2pol pol2cart sph2cart

Rational approximation. Rational output. Transform from Cartesian to spherical coordinates. Transform from Cartesian to polar coordinates. Transform from polar to Cartesian coordinates. Transform from spherical to Cartesian coordinates.
>> help matfun Matrix functions  numerical linear algebra. Matrix analysis. cond  Matrix condition number. norm  Matrix or vector norm. rcond  LINPACK reciprocal condition estimator. rank  Number of linearly independent rows or columns. det  Determinant. trace  Sum of diagonal elements. null  Null space. orth  Orthogonalization. rref  Reduced row echelon form. Linear equations. \ and /  Linear equation solution; use "help slash". chol  Cholesky factorization. lu  Factors from Gaussian elimination. inv  Matrix inverse. qr  Orthogonaltriangular decomposition. qrdelete  Delete a column from the QR factorization. qrinsert  Insert a column in the QR factorization. nnls  Nonnegative leastsquares. pinv  Pseudoinverse. lscov  Least squares in the presence of known covariance. Eigenvalues and singular values. eig  Eigenvalues and eigenvectors. poly  Characteristic polynomial. polyeig  Polynomial eigenvalue problem. hess  Hessenberg form. qz  Generalized eigenvalues. rsf2csf  Real block diagonal form to complex diagonal form. 31
cdf2rdf schur balance svd

Complex diagonal form to real block diagonal form. Schur decomposition. Diagonal scaling to improve eigenvalue accuracy. Singular value decomposition.
Matrix functions. expm  Matrix exponential. expm1  Mfile implementation of expm. expm2  Matrix exponential via Taylor series. expm3  Matrix exponential via eigenvalues and eigenvectors. logm  Matrix logarithm. sqrtm  Matrix square root. funm  Evaluate general matrix function. >> help general General purpose commands. MATLAB Toolbox Version 4.2a 25Jul94 Managing commands and functions. help  Online documentation. doc  Load hypertext documentation. what  Directory listing of M, MAT and MEXfiles. type  List Mfile. lookfor  Keyword search through the HELP entries. which  Locate functions and files. demo  Run demos. path  Control MATLAB’s search path. Managing variables and the workspace. who  List current variables. whos  List current variables, long form. load  Retrieve variables from disk. save  Save workspace variables to disk. clear  Clear variables and functions from memory. pack  Consolidate workspace memory. size  Size of matrix. length  Length of vector. disp  Display matrix or text. Working with files and the operating system. 32
cd dir delete getenv ! unix diary

Change current working directory. Directory listing. Delete file. Get environment value. Execute operating system command. Execute operating system command & return result. Save text of MATLAB session.
Controlling the command window. cedit  Set command line edit/recall facility parameters. clc  Clear command window. home  Send cursor home. format  Set output format. echo  Echo commands inside script files. more  Control paged output in command window. Starting and quitting from MATLAB. quit  Terminate MATLAB. startup  Mfile executed when MATLAB is invoked. matlabrc  Master startup Mfile. General information. info  Information about MATLAB and The MathWorks, Inc. subscribe  Become subscribing user of MATLAB. hostid  MATLAB server host identification number. whatsnew  Information about new features not yet documented. ver  MATLAB, SIMULINK, and TOOLBOX version information. >> help funfun Function functions  nonlinear numerical methods. ode23 ode23p ode45 quad quad8 fmin fmins fzero fplot

Solve differential equations, low order method. Solve and plot solutions. Solve differential equations, high order method. Numerically evaluate integral, low order method. Numerically evaluate integral, high order method. Minimize function of one variable. Minimize function of several variables. Find zero of function of one variable. Plot function. 33
See also The Optimization Toolbox, which has a comprehensive set of function functions for optimizing and minimizing functions. >> help polyfun Polynomial and interpolation functions. Polynomials. roots poly polyval polyvalm residue polyfit polyder conv deconv

Find polynomial roots. Construct polynomial with specified roots. Evaluate polynomial. Evaluate polynomial with matrix argument. Partialfraction expansion (residues). Fit polynomial to data. Differentiate polynomial. Multiply polynomials. Divide polynomials.
Data interpolation. interp1  1D interpolation (1D table lookup). interp2  2D interpolation (2D table lookup). interpft  1D interpolation using FFT method. griddata  Data gridding. Spline interpolation. spline  Cubic spline data interpolation. ppval  Evaluate piecewise polynomial. >> help ops Operators and special characters. Char + * .* ^ .^
Name
HELP topic
Plus Minus Matrix multiplication Array multiplication Matrix power Array power
arith arith arith arith arith arith
34
\ / ./ kron
Backslash or left division Slash or right division Array division Kronecker tensor product
slash slash slash kron
:
Colon
colon
( ) [ ]
Parentheses Brackets
paren paren
. .. ... , ; % ! ’ =
Decimal point Parent directory Continuation Comma Semicolon Comment Exclamation point Transpose and quote Assignment
punct punct punct punct punct punct punct punct punct
== Equality <,> Relational operators & Logical AND  Logical OR ~ Logical NOT xor Logical EXCLUSIVE OR
relop relop relop relop relop xor
Logical characteristics. exist  Check if variables or functions are defined. any  True if any element of vector is true. all  True if all elements of vector are true. find  Find indices of nonzero elements. isnan  True for NotANumber. isinf  True for infinite elements. finite  True for finite elements. isempty  True for empty matrix. isreal  True for real matrix. issparse  True for sparse matrix. isstr  True for text string. isglobal  True for global variables. 35
>> help lang Language constructs and debugging. MATLAB as a programming language. script  About MATLAB scripts and Mfiles. function  Add new function. eval  Execute string with MATLAB expression. feval  Execute function specified by string. global  Define global variable. nargchk  Validate number of input arguments. lasterr  Last error message. Control flow. if else elseif end for while break return error

Conditionally execute statements. Used with IF. Used with IF. Terminate the scope of FOR, WHILE and IF statements. Repeat statements a specific number of times. Repeat statements an indefinite number of times. Terminate execution of loop. Return to invoking function. Display message and abort function.
Interactive input. input  Prompt for user input. keyboard  Invoke keyboard as if it were a Scriptfile. menu  Generate menu of choices for user input. pause  Wait for user response. uimenu  Create user interface menu. uicontrol  Create user interface control. Debugging commands. dbstop  Set breakpoint. dbclear  Remove breakpoint. dbcont  Resume execution. dbdown  Change local workspace context. dbstack  List who called whom. dbstatus  List all breakpoints. dbstep  Execute one or more lines. 36
dbtype dbup dbquit mexdebug

List Mfile with line numbers. Change local workspace context. Quit debug mode. Debug MEXfiles.
>> help plotxy Two dimensional graphics. Elementary XY graphs. plot  Linear plot. loglog  Loglog scale plot. semilogx  Semilog scale plot. semilogy  Semilog scale plot. fill  Draw filled 2D polygons. Specialized polar bar stem stairs errorbar hist rose compass feather fplot comet
XY graphs.  Polar coordinate plot.  Bar graph.  Discrete sequence or "stem" plot.  Stairstep plot.  Error bar plot.  Histogram plot.  Angle histogram plot.  Compass plot.  Feather plot.  Plot function.  Cometlike trajectory.
Graph annotation. title  Graph title. xlabel  Xaxis label. ylabel  Yaxis label. text  Text annotation. gtext  Mouse placement of text. grid  Grid lines. See also PLOTXYZ, GRAPHICS. >> help plotxyz Three dimensional graphics. 37
Line and area fill commands. plot3  Plot lines and points in 3D space. fill3  Draw filled 3D polygons in 3D space. comet3  3D cometlike trajectories. Contour and other 2D plots of 3D data. contour  Contour plot. contour3  3D contour plot. clabel  Contour plot elevation labels. contourc  Contour plot computation (used by contour). pcolor  Pseudocolor (checkerboard) plot. quiver  Quiver plot. Surface and mesh plots. mesh  3D mesh surface. meshc  Combination mesh/contour plot. meshz  3D Mesh with zero plane. surf  3D shaded surface. surfc  Combination surf/contour plot. surfl  3D shaded surface with lighting. waterfall  Waterfall plot. Volume visualization. slice  Volumetric visualization plots. Graph appearance. view  3D graph viewpoint specification. viewmtx  View transformation matrices. hidden  Mesh hidden line removal mode. shading  Color shading mode. axis  Axis scaling and appearance. caxis  Pseudocolor axis scaling. colormap  Color lookup table. Graph annotation. title  Graph title. xlabel  Xaxis label. ylabel  Yaxis label. zlabel  Zaxis label for 3D plots. text  Text annotation. 38
gtext grid
 Mouse placement of text.  Grid lines.
3D objects. cylinder  Generate cylinder. sphere  Generate sphere. See also COLOR, PLOTXY, GRAPHICS. >> help strfun Character string functions. General. strings abs setstr isstr blanks deblank str2mat eval

About character strings in MATLAB. Convert string to numeric values. Convert numeric values to string. True for string. String of blanks. Remove trailing blanks. Form text matrix from individual strings. Execute string with MATLAB expression.
String comparison. strcmp  Compare strings. findstr  Find one string within another. upper  Convert string to uppercase. lower  Convert string to lowercase. isletter  True for letters of the alphabet. isspace  True for white space characters. strrep  Replace a string with another. strtok  Find a token in a string. String to number conversion. num2str  Convert number to string. int2str  Convert integer to string. str2num  Convert string to number. mat2str  Convert matrix to string. sprintf  Convert number to string under format control. sscanf  Convert string to number under format control.
39
Hexadecimal to number conversion. hex2num  Convert hex string to IEEE floating point number. hex2dec  Convert hex string to decimal integer. dec2hex  Convert decimal integer to hex string.
Also the MathWorks web site has a lot of good tutorials, examples and reference documentation. http://www.mathworks.com A good tutorial is at http://www.mathworks.com/access/helpdesk/help/techdoc/ learn_matlab/learn_matlab.shtml
40
List of Figures 1 2 3
4
A simple finite element mesh of triangular elements . . . . . An example of a element force vector f e scattered into a global force vector f . . . . . . . . . . . . . . . . . . . . . . . . . . Diagram of beam used in beam bending example. The following displacement boundary conditions are applied: ux = 0 at the points (0, ±c) and (0,0), uy = 0 at (0, 0). The following traction boundary conditions are used tx = y on x = 0 and ty = P (x2 − c2 ) on x = L. . . . . . . . . . . . . . . . . . . . AFM tip modeled in modal analysis example . . . . . . . . .
41
. 42 . 43
. 44 . 45
(2,6)
6
5
4
3 3
4 2
1
2
(0,0) 1
Figure 1: A simple finite element mesh of triangular elements
42
1x
4
2x e
3x 4x
2
5x 1x 6x 2x 1y 1y 2y 2y 3y
f
4y
e
5y 6y
f Figure 2: An example of a element force vector f e scattered into a global force vector f
43
y
x
L
2c
Figure 3: Diagram of beam used in beam bending example. The following displacement boundary conditions are applied: ux = 0 at the points (0, ±c) and (0,0), uy = 0 at (0, 0). The following traction boundary conditions are used tx = y on x = 0 and ty = P (x2 − c2 ) on x = L.
44
Figure 4: AFM tip modeled in modal analysis example
45