2
>> sigma=[40 20 -18; 20 28 12; -18 12 14];
>> e = MaxEig(Sigma)
e =
56.4898
1
5.11 Write a user-defined MATLAB function that determines the smallest eigenvalue of an matrix
by using the inverse power method. For the function name and argument use e = MinEig(A), where A is
the matrix and e is the value of the smallest eigenvalue. Inside MinEig MATLAB’s built-in function inv
for calculating the inverse of the matrix A. Use the function MinEig for calculating the smallest eigen-
value of the matrix of Problem 5.8. Check the answer by using MATLAB’s built-in function for finding the
eigenvalues of a matrix.
Solution
The listing of the user-defined function MinEig is:
function e = MinEig(A)
Ainv = Inverse(A);
one_over_e = MaxEig(Ainv);
e = 1/one_over_e;
nn×()
1
5.12 Write a user-defined MATLAB function that determines all the eigenvalues of an matrix by
using the QR factorization and iteration method. For the function name and argument use
e=AllEig(A), where A is the matrix and e is a vector whose elements are the eigenvalues. Use the func-
tion AllEig for calculating the eigenvalues of the matrix of Problem 5.8. Check the answer by using
MATLAB’s built-in function for finding the eigenvalues of a matrix.
Solution
The listing of the user-defined function AllEig is (the user-defined function QRFactorization is a
subfunction inside AllEig):
nn×()
2
Q = I;
for j = 1:n-1
c = R(:,j);
c(1:j-1) = 0;
e(1:n,1)=0;
The user defined function AllEig is used in the Command Window for calculating the smallest eigen-
value of the matrix of Problem 5.8:
σij
40 20 18–
20 28 12
18–1214
=
1
5.13 Write a user-defined MATLAB function that determines the principal stresses and the directions of
the principal stresses for a given three-dimensional state of stress. For the function name and arguments,
use [Ps Pd] = PrinplStre(S), where S is a matrix with the values of the stress tensor, Ps is
a column vector with the values of the principal stresses, and Pd is a matrix in which each row lists
a unit vector in a principal direction. Use MATLAB built-in functions.
Use the function for determining the principal stresses and principal directions for the state of stress
given in Problem 5.8: ksi.
Solution
The principal stresses are the eigenvalues of the matrix and the directions of the principal stresses are the
eigenvalues. The user-defined function PrinplStre first calculates the eigenvalues and then the eigen-
vectors. The eigenvalues are calculated with a user-defined function named AllEig (see listing below)
which uses the QR factorization and iteration method. The eigenvectors are calculated by using the follow-
ing procedure:
33×()
33×()
σij
40 20 18–
20 28 12
18–1214
=
2
The third eigenvalue and eigenvector are: (known) and (unknown). They are
related by:
By setting , the following two equations can be written:
σ11 λ2
–σ12 σ13
σ21 σ22 λ2
–σ23
σ31 σ32 σ33 λ2
–
u12()
u22()
u32()
0
0
0
=
λ3
u3() u13()u23()u33()
,,[]=
σ11 λ3
–σ12 σ13
σ21 σ22 λ3
–σ23
σ31 σ32 σ33 λ3
–
u13()
u23()
u33()
0
0
0
=
σ11 λ3
–()u13() σ12u23()
+σ13
–=
σ11 λ2
+σ12
–=
–=
3
% direction.
I=eye(3);
Ps = AllEig(S);
% Once the eigenvalues are found, the eigenvectors are u(:,1) corresponding
% to the first eigenvalue, u(:,2) corresponding to the second, and u(:,2)
% corresponding to the second.
u([1:2],3)=Third\rhs;
% Assign the eigenvectors to Pd.
Pd=[u(:,1) u(:,2) u(:,3)];
% Normalize the eigenvectors so they are unit vectors:
Pd(:,1)=Pd(:,1)./sqrt(sum(Pd(:,1).^2));
Pd(:,2)=Pd(:,2)./sqrt(sum(Pd(:,2).^2));
Pd(:,3)=Pd(:,3)./sqrt(sum(Pd(:,3).^2));
4
% Input argument:
% A A matrix
% Output argument:
% e A column vector with the values of the eigenvalues.
end
if con == 0
break
end
end
The listing of the user-defined function QRFactorization that is used inside AllEig (was written in
Example 5-4 (Program 5-1)) is:
function [Q R] = QRFactorization(R)
% The function factors a matrix [A] into an orthogonal matrix [Q]
5
nmatrix = size(R);
n = nmatrix(1);
I = eye(n);
Q = I;
for j = 1:n-1
c = R(:,j);
c(1:j-1) = 0;
e(1:n,1)=0;
if c(j) > 0
e(j) = 1;
else
e(j) = -1;
end
6
Pd =
0.8396 -0.2406 0.4870
0.4989 0.6961 -0.5163
-0.2148 0.6765 0.7045
>> % Checking the answer with MATLAB’s built-in function eig:
>> [v d]=eig(sigma)
v =
1
5.14 The structure of the C2H2 (acetylene) molecule may be ide-
alized as four masses connected by two springs (see discussion in
Problem 5.6). By applying the equation of motion, the following
system of equations can be written for the amplitudes of vibration
of each atom:
where is the frequency, kg/s2 and kg/s2 are the restoring force spring
constants representing the C–H and C–C bonds, respectively, and and are the
masses of the atoms ( kg).
(a) Determine the eigenvalues (frequencies) and the corresponding wavelengths (where
m/s is the speed of light).
(b) Determine the eigenvectors corresponding to the eigenvalues found in part (a). From the eigenvectors,
deduce the relative motion of the atoms (i.e., are they moving toward or away from each other?).
Solution
The following script file solves this problem:
kCH
mH
——––ω2
–kCH
mH
——––
–00
kCH
mC
——––
–kCH kCC
+()
mC
—————-————–ω2
–kCC
mC
——-–
–0
0kCC
mC
——-–
–kCH kCC
+()
mC
—————–———––ω2
–kCH
mC
——––
–
00 kCH
mH
——––
–kCH
mH
——––ω2
–
A1
A2
A3
A4
0
0
0
0
=
ω
kCH 5.92 102
×=
kCC 15.8 102
×=
mH1amu=
mC12amu=
1amu 1.6605 10 27–
×=
ω
λ2πc
ω
——––
=
c310
8
×=
2
temp=sqrt(eigenvalues); icount=1;
% Looking for the real eigenvalues and corresponding eigenvectors.
end
disp(‘The real eigenvalues are:’)
omega
disp(‘The corresponding wavelengths (except for the zero eigenvalue) in meters
are:’)
for i=2:length(omega)
lambda=2*pi*c./omega;
end
lambda
disp(‘The eigenvectors corresponding to the eigenvalues are:’)
eigenvectors
When the script file is executed, the following output is displayed in the Command Window:
3
6.0400e-001 -7.0057e-001 7.0466e-001
(b) The columns of the eigenvector matrix are the individual eigenvectors themselves. They represent
amplitudes of the motion or vibration of each atom in the molecule. The first column represents the
eigenvector for the eigenvalue associated with the symmetric stretch of the molecule where the hydrogen
C2H2