Jacob Fish and Ted Belytschko
Preface This book is intended as a textbook for junior and senior level undergraduate students in science and engineering. The text material evolved from over 50 years of combined teaching experience of the authors in undergraduate finite element courses. The book deals with a formulation and application of the finite element method and it differs from the related finite element textbooks in the following three respects: 1. It is introductory introductory and self-contained self-contained . Only a modest background covered by most engineering and science curricula in the first two years is required. The relevant topics in mathematics, such as matrix algebra, differential and integral equations, and in physics, such as conservation laws and constitutive equations, are reviewed as needed. reinforced. 2. It is generic generic. While most introductory finite element textbooks are app lication specific, such as being based on energy principles for linear elasticity, the finite element method in this book is formulated as a general purpose numerical procedure for solving partial differential equations. Consequently, students from various engineering and science disciplines will benefit eq ually from the exposition of the subject. 3. It is a hands-on experience . The book integrates finite element theory, finite element code development and application of commercial software packages. Finite element code development is introduced through MATLAB exercises and a MATLAB program, whereas ABAQUS is chosen for demonstrating the use of the finite element method in the commercial arena. The material in the book can be covered in a single semester. The course material is organized in three chronological units one month month each: (1) the finite element formulation formulation for one-dimensional problems, (2) the finite element formulation for scalar field problems in two dimensions, (3) finite element programming and a pplication to scalar field problems; and finite element formulation for vector field problems in two dimensions and beams. The following course outline is recommended: Part 1: Finite element formulation for trusses and one-dimensional problems Chapter 1: Introduction Chapter 2: Direct Approach Chapter 3: Strong and weak forms for one-dimensional problems Chapter 4: Trial and weight functions for one-dimensional problems
2
Chapter 5: Finite element formulation for one-dimensional problems ******* Quiz 1********** Part 2: Finite element formulation for for scalar field problems in two dimensions Chapter 6: Strong and weak forms for scalar field problems in two dimensions Chapter 7: Trial and Weight Functions in two dimensions Chapter 8: Finite element formulation for two-dimensional heat flow *******Quiz 2********** Part 3: Finite element program programming programming and application to scalar field problems problems Chapter 11: Finite element programming with MATLAB. *********Project assignment***** Chapter 12: Application of general purpose finite element software – ABAQUS tutorial for linear heat conduction and elasticity problems Chapter 9: Finite element formulation for vector field problems in twodimensions Chapter 10: Finite element formulation for beams *******Project Due**********
3
Chapter
1
Introduction
1.1 Background Many physical phenomena in engineering and science can be described in terms of partial differential equations. In general, solving these these equations by classical analytical methods for arbitrary shapes is is almost impossible. impossible. The finite element method method is a numerical approach by which these partial differential equations can be solved approximately. From an engineering standpoint, the finite element method is a method for solving engineering problems such as stress analysis, heat transfer, fluid flow and electromagnetics by computer simulation.
Figure 1.1: Geometry, loads and finite element meshes
Millions of engineers and scientists worldwide are using the Finite Element Method (FEM) to predict the behavior of structural, mechanical, thermal, electrical and chemical systems in both design and performance analysis. Its popularity can be gleaned by the fact that over $1 billion is spent annually in the United States on FEM software and computer time. More than 250 books and over 20,000 articles have been published on this subject. To explain the basic approach of the finite element method, consider a plate with a hole as shown shown in Figure 1.1. Suppose that our objective objective is to find the temperature distribution in the the plate. In principle, it is possible to write a heat balance equation at each point in the plate. However, the solution of the partial differential equations that governs energy balance for a complicated geometry such as an engine block becomes almost impossible by classical methods like like separation of variables. Numerical methods such as finite difference methods also are quite awkward for arbitrary shapes, software developers have not marketed finite difference programs that can deal with the arbitrary geometries that are commonplace in engineering. Similarly, in stress analysis problems, equilibrium is governed by complex partial differential equations that are very difficult to solve except for very simple geometries, like rectangles, and engineering problems seldom have such simple shapes. The basic idea behind the FE method is to divide the body or region into a finite number of finite elements, connected by nodes, and obtain an approximate solution in terms of the temperatures at these nodes. These elements may be one, two or threedimensional. A popular two-dimensional element is the triangular element. When a two dimensional region is divided into non-overlapping triangles, we can see that essentially any planar geometry can be easily represented by a union (in the language of mathematical sets) of these elements. Instead of determining the temperature at every point in the plate, let us consider determining the temperature at only a finite number of points, the vertices (or nodes) of the triangles. The finite element method method provides a systematic systematic methodology whereby the temperatures at the the nodes (nodal temperatures) can be determined. Interpolation can then be used to obtain the temperature throughout the problem domain. For linear problems, we will need to solve a system of linear equations with the number of unknowns (nodal (nodal temperatures) equal to the number of nodes. To obtain a reasonably accurate solution, thousands of nodes are usually needed, so computers are essential in solving these equations. Generally, the accuracy of the solution improves as the number of elements (nodes) increases, but the computer time, and hence the cost, also increases. The computer program determines the temperature at each node and the heat flow through each element. The results are usually presented as computer visualizations, such as contour plots, although selected results are often printed. This information is then used in the engineering design process. Although we have just presented the concept in the context of heat transfer, the same basic concept is applicable in other areas. In stress analysis, the field variables are the displacements, in chemical systems the variables are material concentrations, and in electromagnetics, the potential field. This course is limited limited to linear finite element analysis. The preponderance of finite element analyses in engineering design is today still linear FEM. In heat conduction, linearity requires requires that the conductance be independent independent of temperature. In
5
stress analysis, linear FEM is only applicable if the material behavior is linear elastic and the displacements are small. Additional discussions of the assumptions underlying linear analysis are given later. The finite element method was developed in the 1950’s in the aerospace industry. The major players were Boeing and Bell Aerospace (long vanished) in the U.S. and Rolls Royce in the U.K. M. J. Turner, R. W. Clough, H. C. Martin, Martin, and L. J. Topp published one of first papers that laid out the major ideas in 1956. It established the procedures of element matrix assembly and element formulations that you will learn in this course, but did not use the term “finite elements.” The second author in this paper, Ray Clough, was a professor at Berkeley, who was at Boeing for a summer summer job. Subsequently he wrote a paper that first used the name finite elements, and he is given much credit as one of the founders of the method. He only worked in finite elements for for a few more more years, and then turned to experimental methods, but his work instigated a tremendous effort at Berkeley, led by the younger professors, primarily E. Wilson and R.L. Taylor and graduate students such such as T.J.R. Hughes, C. Felippa, and K.J. Bathe. These established Berkeley as the center of linear finite element research research for many years. years. Their work coincided with the rapid growth of computer power, and the method quickly became widely used in the nuclear power, defense, automotive and aeronautics industries. Much of the academic community first viewed finite element methods very skeptically, and some journals refused to publish papers on FEM-the typical resistance of mankind (and particularly academic communities) to the new. Nevertheless, several capable researchers early recognized its potential, most notably O.C. Zienkiewicz (then at Northwestern) and R.H. Gallagher (at Cornell). O.C. Zienkiwiecz built a renowned group at Swansea in Wales that included B. Iron, R. Owen and many others who pioneered concepts like the isoparametric isoparametric element and nonlinear analysis methds. Other important early contributors were J.H. Argyris and J.T. Oden. Subsequently mathematicians discovered a 1943 paper by R. Courant, in which he used triangular elements with variational principles to solve vibration problems. Consequently, many mathematicians have claimed that this was the original discovery of the method (though it is somewhat reminiscent of the claim that the Vikings discovered America instead of Columbus). It is interesting that for many years the finite element method lacked a theoretical basis, i.e. any logical underpinning that finite element solutions gave the right answer, and this original work of Courant’s was duplicated by the engineering community to provide a mathematical justification for for the method. In the late 1960’s, the field field became of great interest interest to many mathematicians, who showed that that for linear problems, such as we will deal with in this course, finite element solutions converge to the correct solution of the partial differential equation (provided that certain aspects of the problem are sufficiently smooth). smooth). In other words, it has been shown that as the number of elements is increased the solutions will improve and tend in the limit to the exact solution of the partial differential equations. The improvements in personal computers since the 1980’s have had a major impact on the finite finite element applications, and hence engineering. Until about the early 1990’s, personal computers were too feeble to handle the many equations needed for practical finite element analysis. However, around 1990, personal computers crossed that threshold, and PCs rapidly became capable of treating finite element analyses of linear problems required in engineering practice. practice. The software that is most widely used by
6
engineers in stress analysis, such as NASTRAN, ANSYS, ABAQUS, and LSDYNA, were adapted to PCs. Today a large part of finite element analysis in industry is done by PCs, identical to yours.
1.2 Matrix Algebra with MATLAB MATLAB is an interactive software system for numerical computations and graphics. As the name suggests, MATLAB is especially designed for matrix computations. In addition, it has a variety of graphical and visualization capabilities, and can be extended through programs written in its own programming language. Here, we introduce only some basic ideas so that you can perform essential matrix operations and basic programming needed for understanding and development of the finite element program.
1.2.1 Definition of matrices A matrix is an mxn array of numbers or variables arranged in m rows and n columns; such a matrix is said to have dimension mxn as shown below ⎡ a11 ⎢ ⎢ a 21 a = ⎢⎢ ⎢ ⎢a ⎢⎣ m 1
a12 a 22
a 1n ⎤
⎥ ⎥ ⎥ ⎥ ⎥ a mn ⎥⎥ ⎦
Bold letters will denote matrices or vectors. The elements of a matrix matrix a are denoted by a ij , where i is the row number and j is the column number. Note that in both describing the dimension of the matrix and in the subscripts identifying the row and column number, the row is always placed first. An example of a 3x3 matrix is: ⎡1 ⎢ ⎢ a = ⎢4 ⎢ ⎢7 ⎢⎣
2 3⎤⎥ ⎥
5 6⎥ ⎥
8 0⎥⎥ ⎦
The above matrix a is is an example of a square matrix since the number of rows and columns are equal. The following commands show how to enter matrices in MATLAB (>> is the MATLAB prompt; it may be different with different computers or different versions of MATLAB.)
7
>> a = [1 2 3; 4 5 6; 7 8 0] a = 1
2
3
4
5
6
7
8
0
Notice that rows of a matrix are separated by semicolons, while the entries on a row are separated by spaces (or commas). The order of o f matrix a can be determined from
size(a ) The transpose of any matrix matrix is obtained by interchanging rows and columns. So for example, the transpose of a is: ⎡1 ⎢ ⎢ T a = ⎢2 ⎢ ⎢3 ⎢⎣
4 7⎤⎥ ⎥ 5 8⎥ ⎥
6 0⎥⎥ ⎦
In MATLAB the transpose of a matrix is denoted by an apostrophe (‘). Square matrices may be symmetric or diagonal. If aT = a we say a is symmetric . A matrix is called a column matrix or a vector if n=1, e.g. ⎡b ⎤ ⎢ 1⎥ ⎢ ⎥ b = ⎢b2 ⎥ ⎢ ⎥ ⎢b ⎥ ⎢⎣ 3 ⎥⎦
A column vector in MATLAB is defined as
>> b = [1 2 3]' b = 1 2 3 The components of the vector b are b1, b2 ,b3 . The transpose of b is given by a row vector
8
bT = ⎡⎢b1 b2 ⎣
b3 ⎤⎥ ⎦
or in MATLAB
>> b = [1 2 3] b = 1
2
3
A matrix is called a diagonal matrix if only the diagonal components are nonzero, i.e., aij = 0, i ≠ j . For example, the matrix below is a diagonal matrix:
⎡1 0 0 ⎤ ⎢ ⎥ a= 0 5 0 ⎢ ⎥ ⎢⎣0 0 6 ⎥⎦ Diagonal matrix in MATLAB is constructed by first defining a row vector b = [1 5 6], then
>> b = [1 5 6]; >> a = diag (b) a = 1
0
0
0
5
0
0
0
6
A diagonal matrix where all diagonal components are equal to one is called an identity or unit matrix and is denoted by I. For example, 2 × 2 identity matrix is given by ⎡1 I = ⎢⎢ ⎢⎣0
0⎤ ⎥ 1⎥⎥ ⎦
The MATLAB expression for an order n unit matrix is
eye(n ) (2 ) gives the matrix above. Thus, the MATLAB expression I = eye(2) A matrix in which all components are zero is called a zero matrix and is denoted by 0. In MATLAB, B = zeros (m, n) creates m×n matrix B of zeros. A random m×n matrix can be created by rand (m,n).
9
In finite element method matrices are often sparse, i.e., they contain many zero terms. MATLAB has the ability to store and manipulate sparse matrices, which greatly increases its usefulness for realistic problems. The command sparse (m, n) creates an m × n zero matrix in sparse format. The nonzero entries can then be added one-by-one or in a loop.
>> a = sparse (3,2) a = All zero sparse: 3-by-2 >> a(1,2)=1; >> a(3,1)=4; >> a(3,2)=-1; >> a a = (3,1)
4
(1,2)
1
(3,2)
-1
Notice that the display can be suppressed by ending the line with a semicolon. The inverse of a square matrix is defined such that a −1a = aa −1 = I
provided that matrix a is not singular. The MATLAB expression is inv(a ) . The inverse can also be obtained using backslash operator as shown in Section 1.2.3. Matrix a is nonsingular if the determinant, denoted by det (a) , is not equal to zero. A determinant of order 2 is defined as ⎡a11 a=⎢ ⎢a22 ⎣
a 12 ⎤ ⎥ a 22 ⎥ ⎦
det (a) = a11a22 − a12a21
The MATLAB expression for determinant is
det(a ) For example,
10
>> a = [1 3; 4 2]; >> det (a) ans = -10
1.2.2 Operation with matrices 1.2.2.1 Addition and Subtraction ⎡ a 11 ⎢ ⎢a ⎢ 21 c=a±b= ⎢ ⎢ ⎢ ⎢a ⎣⎢ m 1
± b11
a12 ± b12
± b21
a22 ± b22
a1n ± b1n ⎤
± bm1
a mn
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ± bmn ⎥⎥ ⎦
In example of matrix addition in MATLAB is given below:
>> a = [1 2 3;4 5 6;7 8 9]; >> a = [1 1 1;2 2 2;3 3 3]; >> c = [1 2;3 4;5 6]; >> a+b ans = 2
3
4
6
7
8
10
11
12
>> a+c ??? Error using ==> + Matrix dimensions must agree 1.2.2.2 Multiplication Multiplication 1. Scalar-matrix (vector) multiplication ⎡ a11 ⎢ ⎢ a 21 ca = c ⎢⎢ ⎢ ⎢a ⎢⎣ m 1
a12 12 a22 22
a1n ⎤
⎡ ca11 ⎥ ⎢ ⎥ ⎢ ca21 ⎥=⎢ ⎥ ⎢ ⎥ ⎢ amn ⎥⎥ ⎢⎢cam1 ⎦ ⎣
11
ca c a12 ca22
ca1n ⎤
⎥ ⎥ ⎥ ⎥ ⎥ camn ⎥⎥ ⎦
2. Scalar product
aT b = ⎡⎢a1 ⎣
a2
⎡b1 ⎤ ⎢ ⎥ ⎢b ⎥ ⎢ 2⎥ … an ⎤⎥ ⎢ ⎥ = ⎦⎢⎥ ⎢ ⎥ ⎢b ⎥ ⎣⎢ n ⎦⎥
n
∑ aibi i =1
In MATLAB the scalar product as defined above is either a * b ′ or dot(a, b) . The length of a vector a is denoted by |a| and is given by 2 a = a12 + a22 + + a n
3. Product of two matrices
The product of two matrices a (m × k ) and b (k × n ) is defined as ⎡ k ⎢ ∑ a1 jbj 1 ⎢ j =1 ⎢ ⎢ k ⎢ a b ⎢ ∑ 2 j j 1 c = ab = ⎢ j =1 ⎢ ⎢ ⎢ k ⎢ ⎢ ∑ amjbj 1 ⎢ j =1 ⎣
k
k
∑ a1 jb j 2
j =1
⎤
∑ a1 jbjn ⎥⎥ j =1
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ k ⎥ a b ∑ mj jn ⎥⎥ j =1 ⎦
k
∑ a2 jbj 2 j =1
The components of matrix c are given as n
cij = ∑ aikbkj
k =1
The product of two matrices a and b can be taken only if the number of columns in a equals the number of rows in a. In other words, if a is an (m × k ) matrix, then b must be an (k × n ) matrix, where k is arbitrary. The product c will then have the same number of rows as a and the same number of columns as b, i.e. it will be an m × n matrix. An important fact to remember is that matrix multiplication is not commutative, i.e. in general ab ≠ ba . The MATLAB expression for matrix multiplication is
12
c = a *b Consider the same matrices a and c as before. An example of matrix multiplication with MATLAB is:
>> a*c ans = 22
28
49
64
76 100 >> c*c ??? Error using ==> * Inner matrix dimensions must agree.
4. Other matrix operations T
a) Transpose of product: (ab) = bT aT b) Product with identity matrix: aI = a c) Product with zero matrix: a0 = 0
1.2.3 Solution of system of linear equations Consider the following system of n equations with n unknowns, d k , k = 1, 2, , n : ⎧K11d1 + K 2d2 + + K1ndn = f1 ⎪ ⎪ ⎪ ⎪ ⎪K 21d1 + K 22d2 + + K 2ndn = f2 ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ K d + K n 2d2 + + Knndn = fn ⎪ ⎪ n1 1 ⎩
We can rewrite this system of equations in matrix notation as follows: Kd = f
where
13
⎡ K 11 ⎢ ⎢K ⎢ 21 K=⎢ ⎢ ⎢ ⎢K ⎢⎣ n 1
K 12 K 22
K1n ⎤
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ K nn ⎥⎥ ⎦
⎡ f1 ⎤ ⎢ ⎥ ⎢f ⎥ ⎢ 2⎥ f =⎢ ⎥ ⎢⎥ ⎢ ⎥ ⎢f ⎥ ⎢⎣ n ⎥⎦
⎡d1 ⎤ ⎢ ⎥ ⎢d ⎥ ⎢ 2⎥ d=⎢ ⎥ ⎢⎥ ⎢ ⎥ ⎢d ⎥ ⎢⎣ n ⎥⎦
The symbolic solution of the above system of equation can be found by multiplying both sides with inverse of K, which yields − d = K 1f
MATLAB expression for solving the system of equations is
d = K \ f or
d = inv (K ) * f An example of solution of system of equations equ ations with MATLAB is given below:
>> A = rand (3,3) A = 0.2190
0.6793
0.5194
0.0470
0.9347
0.8310
0.6789
0.3835
0.0346
>> b = rand (3,1) b = 0.0535 0.5297 0.6711 >> x = A\b x = -159.3380 314.8625 -344.5078
14
Remarks a) In practical FEM analysis, the number of equations n can be very large. PCs can today solve thousands of equations in a matter of minutes if they are sparse (as they are in FEM analysis-you will learn about this later) but sometimes millions of equations are needed, as for an aircraft carrier or a full model of an aircraft; parallel computers are then ne eded. b) Efficient solution techniques that take advantage of the sparsity and other advantageous properties of FEM equations are essential for treating even moderately large systems. The issue of how to efficiently solve large systems will not be considered in this course. c) In this course, we will see that • The system matrix for a finite element model K is non-singular (often called regular), i.e., K−1 exists if the correct boundary conditions are prescribed and the elements elements are properly formulated. Furthermore, for good models it is usually well-conditioned, which means it is not very sensitive to roundoff errors. • K is symmetric, i.e. KT = K . • K is positive definite, i.e., xT Kx > 0 ∀x (meaning for any value of x) Alternatively, K is said to be positive definite if all the eigenvalues are strictly positive. The eigenvalue problem consists of finding nonzero eigenvectors y and the corresponding eigenvalues λ satisfying Ky = λy
The MATLAB expression for the eigenvalues problem is:
>> K=[2 -2;-2 4]; >> [y, lamda]=eig(K) y = 0.8507 -0.5257 -0.5257
0.8507
lamda = 0.7639 0
15
0 5.2361
1.3 Mathematical Preliminaries 1.3.1 Scalar Peoduct Theorem Consider scalar product of two vetors w and p equal to zero for any vector w denoted as ∀w wT p = 0 ∀w or w1p1 + w2 p2 + … + wn pn = 0 ∀w . Since the vector component wi are arbitrary it is convenient to make the following choices:
w1 = 1 w2 = w3 = … = wn = 0 ⇒ p1 = 0 w2 = 1 w1 = w 3 = … = wn = 0 ⇒ p2 = 0
wn = 1 w1 = w 2 = … = wn −1 = 0 ⇒ pn = 0 which yields p = 0 .
1.3.2 Taylor’s Taylor’s formula with remainder remainder and the mean mean value theorem theorem Taylor’s formula with remainder plays a key role in understanding the behavior of the finite element method discussed in Chapter 5.7. A function f (x ) defined on interval 0 ≤ x ≤ l can be expanded around a point 0 ≤ x 0 ≤ l as follows:
1 2 f (x ) = f (x 0 ) + (x − x 0 ) f,x (x 0 ) + (x − x 0 ) f,xx (x 0 ) 2 1 + (x − x 0 )k f,x ...x (c) k ! k times times
(1.1)
where c is some point in the interval 0 ≤ c ≤ x 0 . The mean value theorem is obtained as a special case for k=1
f (x ) = f (x 0 ) + (x − x 0 ) f,x (c )
(1.2)
Consider a special case of x 0 = 0 and x = l then
f,x (c ) =
f (l ) − f (0) l
16
(1.3)
which means that for any differentiable function f (x ) there is a point in the interval 0 ≤ c ≤ l which has a slope defined by Eq. (1.3). See problem 1.3 for a specific choice of function f (x ) .
1.4 Programming with MATLAB MATLAB is very convenient to write simple finite element programs written in its own programming language. It provides the standard constructs, such as loops and conditionals; these constructs can be used interactively to reduce the tedium of repetitive tasks, or collected in programs stored in ''m-files'' (nothing more than a text file with extension ``.m'').
1.4.1 Conditional and Loops MATLAB has a standard if-elseif-else conditional. The general form An example >> t = 0.76; if expression1 >> if t > 0.75 statements1 s = 0; elseif expression2 elseif t < 0.25 statements2 s = 1; … else … s = 1-2*(t-0.25); … end else >> s statements s= end 0 MATLAB provides two types of loops, a for-loop (comparable to a Fortran do-loop or a C for-loop) and a while-loop. A for-loop repeats the statements in the loop as the loop index takes on the values in a given row vector; the while-loop repeats as long as the given expression expression is true (nonzero): The general form for index = start:increment:end
statements end
while expression
statements end
Examples >> for i=1:1:3 disp(i^2) end 1 4 9 >> x=1; >> while 1+x > 1 x = x/2; end
17
>> x x= 1.1102e-16
1.4.2 Functions Functions allow the user to create new MATLAB commands. A function is defined in an m-file that begins with a line of the following form: function [output1,ou [output1,output2,...] tput2,...] = cmd_name(inp cmd_name(input1,input2 ut1,input2,...) ,...)
The rest of the m-file consists of ordinary MATLAB commands computing the values of the outputs and performing other desired actions. Here is a simple example of a function; func tion; 2 it computes the quadratic function f (x ) = x − 3x − 1 . The following commands should be stored in the file fcn.m (the name of the function within MATLAB is the name of the m-file, without the extension):
function y = fcn( x ) y=x^2-3*x-1; Then type command: >> fcn(0.1) ans = -1.2900
Problems Chapter 1 Problem 1.1: Write a MATLAB program program to generate generate a set of linear algebraic equations Ax=b where A is an n × n matrix given by
⎧-1 if i = j − 1 or i = j + 1 ⎪ Aij = ⎨ 2 if i = j ⎪ 0 otherwise ⎩ −1
−1
and compute A . Then check how closely A A corresponds to I . Do this for n = 5,10,1000 and a larger value you choose. The accuracy of the computed product can be compared to the correct results by computing a norm of the error given by
1 err = 2 n
n
n
∑∑(B i =1 j =1
Repeat the above with the matrix A defined by
Aij =
18
1 i+ j
2
ij − I ij )
Repeat the above with n = 3,4,5 error is greater than one, since ,4,5,,6,7... . Stop when the error the solution is then meaningless. The top equation is of the form we will see in finite element equations. They can be accurately solved even for very large system for they are not susceptible to roundoff error. They are known as well-conditioned. The second matrix is called a Hilbert matrix. It is an example of an estremely ill-conditioned matrix. Problem 1.2: Consider a system of linear equations
⎧ ⎪ 8x 1 + x 2 + 6x 3 = 7.5 ⎪ ⎪ ⎪ ⎪ 3x + 5x + 7x = 4 ⎨ 1 2 3 ⎪ ⎪ ⎪ 4x + 9x 2 + 2x 3 = 12 ⎪ ⎪ 1 ⎩ a) Write the system system of equations in matrix notation Ax = b and solve for the unknown x using MATLAB b) Suppose we impose and additional constraint on the solution new g (x) = x 1 + x 2 + x − 1 = 0 . Using MATLAB find a new vector x so that it will satisfy exactly the constraint equation
g (xnew )
and will minimize the error
T
err (xnew ) = (Axnew − b) (Axnew − b) x 2 (3l − 2x ) Problem 1.3: Consider a cubic function f (x ) = defined on the interval l 3 0 ≤ x ≤ l and its linaer approximation f(x ) defined such that f(0) = f (0) and f(l ) = f (l ) . Verify Eq. (1.3) and find a point c, where f,x (c) = f ,x f(x)
1
f(x ) x l Figure 1.2: Cubic function f(x) and its linear approximation g(x) Problem 1.4: Consider the following symmetrix matrices K :
⎡k1 + k2 ⎢ ⎢ −k 2 ⎢⎣
where
k, k1, k 2
−k2 ⎤ ⎥; k2 ⎥⎥ ⎦
⎡k 0⎤ ⎢ ⎥ ⎢ 0 k ⎥ ; ⎢⎣ ⎥⎦
⎡k ⎢ 1 ⎢ ⎢−k1 ⎢ ⎢ 0 ⎢⎣
−k 1
0 ⎤⎥ ⎥
k1 + k2 −k2 ⎥ ; ⎥
−k2
k 2 ⎥⎥ ⎦
are positive constants.
a) Check if the above three matrices matrices are positive definite. Recall that if for any vector
x ≠ 0 we have xT Kx > 0 then matrix K is called Symmetric Positive Definite
19
T
(SPD). If, on the other hand, x Kx ≥ 0 for any vector x ≠ 0 then the matrix K is symmetric semi-positive definite. Choose one of the semi-positive definite matrices shown above and show that for any right hand side vector, f , the system of equations Kd = f has no unique solution. b) Verify your results by computing the eigenvalues for the above three three matrices
20
