Determining the Characteristic Polynomial of the Companion Matrix by Use of a Large Matrix

Most proofs of the characteristic polynomial of the companion matrix–an important specific case–proceed by induction, and start with a 2\times2 matrix.  It strikes me that an inductive proof has more force (or at least makes more sense) if a larger matrix is used.  In this case we will use a “large” (numerical analysts will laugh at this characterisation) 10\times10 matrix.

Let us begin by making a notation change. Consider the general polynomial


For this to be monic (one of the requirements for the polynomial in question) we should divide by the last coefficient, thus


Our object is thus to prove that this (or a variation of this, as we will see) is the characteristic polynomial of


The characteristic polynomial of this is the determinant of the following:


(For another application of the characteristic polynomial and the companion matrix, click here.)

To find the determinant, we expand along the first row. But then we discover that only two minors that matter: the one in the upper left corner and the one in the upper right. Breaking this up into minors and cofactors yields the following:


The second matrix, however, is an upper triangular matrix with ones for all of its diagonal entries. Its determinant, therefore, is unity. Also rewriting the coefficient of the second term, we have




Repeating this process for the next set of minors and cofactors yields


Note carefully the inclusion of -\lambda in the second term. We can also write this as


Repeating this process until the end, it is easy to see that


or more generally


where n is the degree of the polynomial (and the size of the companion matrix.) If we drop the terms we used to make the polynomial monic, we have at last



“I Need to Find the Way”

When numerical methods are combined with fluid mechanics, you ultimately get Computational Fluid Dynamics (CFD.)  CFD has been a test of numerical methods for a long time, especially with the solution of the Navier-Stokes equations.  The non-linear nature of these equations has made mating them with CFD a challenging problem.

When you’re first introduced to numerical methods you’re told that there are ways one can predict the stability and consistency of these methods with a given set of physics.  These methods are helpful but things go wrong.  On top of that coding these solutions into a language like C++ or Fortran for parallel processing (or any other computer application) is a challenging business; the endless parade of bugs we see exploited in internet security can be also seen in modelling physics such as fluid flow, just different bugs.  Probably the most important lesson one learns in mastering such an art is that, just because it comes out of a computer, it can’t always be trusted.

That situation confronted a Turkish friend of mine who was working on a particularly sticky CFD coding problem.  One frustrating day he exclaimed, “I need to find the way,” i.e., out of his CFD coding problem.

Finding the way has been a key issue in the Middle East for a long time.  The last night before he was crucified, Jesus Christ had this exchange with his disciples:

“We do not know where you are going, Master,” said Thomas; “so how can we know the way?”

Jesus answered: “I am the Way, and the Truth, and the Life; no one ever comes to the Father except through me. (John 14:5-6 TCNT)

This, then, is a key difference between Jesus Christ and the others.  They knew (or they thought they knew) the way, or they have an idea of what the way is.  But Jesus Christ is the way, something he showed when he rose from the dead shortly thereafter.

My Turkish friend was doubly blessed; not only did he find the solution to his CFD problem, but he also found the way in Jesus Christ, and has been baptised.  You can see a video of another Turk who found the way here.

And for the rest of us…

For more information click here.

Mohr’s Circle Analysis Using Linear Algebra and Numerical Methods


Mohr’s Circle-or more generally the stress equilibrium in solids-is a well known method to analyse the stress state of a two- or three-dimensional solid. Most engineers are exposed to its derivation using ordinary algebra, especially as it relates to the determination of principal stresses and invariants for comparison with failure criteria. In this presentation, an approach using linear algebra will be set forth, which is interesting not only from the standpoint of the stress state but from linear algebra considerations, as the problem makes an excellent illustration of linear algebra concepts from a real geometric system. A numerical implementation of the solution using Danilevsky’s Method will be detailed.

1. Introduction

The analysis of the stress state of a solid at a given point is a basic part of mechanics of materials.Although such analysis is generally associated with the theory of elasticity, in fact it is based on static equilibrium, and is also valid for the plastic region as well. In this way it is used in non-linear finite element analysis, among other applications. The usual objective of such an analysis is to determine the principal stresses at a point, which in turn can be compared to failure criteria to determine local deformation.In this approach the governing equations will be cast in a linear algebra form and the problem solved in this way, as opposed to other types of solutions given in various textbooks. Doing it in this way can have three results:

  1. It allows the abstract concepts of linear algebra to be illustrated well in a physical problem.
  2. It allows the physics of the determination of principal stresses to be seen in a different way with the mathematics.
  3. It opens up the problem to numerical solution, as opposed to the complicated closed-form solutions usually encountered, of the invariants, principal stresses or direction cosines.

2. Two-Dimensional Analysis

2.1. Eigenvalues, Eigenvectors and Principal Stresses

The simplest way to illustrate this is to use two-dimensional analysis. Even with this simplest case, the algebra can become very difficult very quickly, and the concepts themselves obscured. Consider first the stress state shown in Figure 1, with the notation which will be used in the rest of the article.


Figure 1: Stress State Coordinates (modified from Verruijt and van Bars [8])

The theory behind this figure is described in many mechanics of materials textbooks; for this case the presentation of Jaeger and Cook [4] was used. The element is in static equilibrium along both axes. In order for the summation of moments to be zero,

\tau_{xy}=\tau_{yx} (1)

The angles are done a little differently than usual; this is to allow an easier transition when three-dimensional analysis is considered. The direction cosines based on these angles are as follows:

l=cos\alpha\ (2)

m=cos\beta\ (3)

Now consider the components p_{x},\,p_{y} of the stress vector p on their respective axes. Putting these into matrix form, they are computed as follows:

\left[\begin{array}{cc} \sigma_{{x}} & \tau_{{\it xy}}\\ \noalign{\medskip}\tau_{{\it xy}} & \sigma_{{y}} \end{array}\right]\left[\begin{array}{c} l\\ \noalign{\medskip}m \end{array}\right]=\left[\begin{array}{c} p_{{x}}\\ \noalign{\medskip}p_{{y}} \end{array}\right] (4)

which is a classic Ax = b type of problem. At this point there are some things about the matrix in Equation 4 that need to be noted (DeRusso et al. [2]):

  1. It is square.
  2. It is real and symmetric. Because of these two properties:
    1. The eigenvalues (and thus the principal stresses, as will be shown) are real. Since for two-dimensional space the equations for the principal stresses are quadratic, this is not a “given” from pure algebra alone.
    2. The eigenvectors form an orthogonal set, which is important in the diagonalization process.
  3. The sum of the diagonal entries of the matrix, referred to as the trace, is equal to the sum of the values of the eigenvalues. As will be seen, this means that, as we rotate the coordinate axes, the trace remains invariant.

At this point there are two related questions that need to be asked. The first is whether static equilibrium will hold if the coordinate axes are rotated. The obvious answer is “yes,” otherwise there would be no real static equilibrium. The stresses will obviously change in the process of rotation. These values can be found using a rotation matrix and multiplying the original matrix as follows (Strang [7]):

\left[\begin{array}{cc} cos\alpha & -sin\alpha\\ sin\alpha & cos\alpha \end{array}\right]\left[\begin{array}{cc} \sigma_{{x}} & \tau_{{\it xy}}\\ \noalign{\medskip}\tau_{{\it xy}} & \sigma_{{y}} \end{array}\right]=\left[\begin{array}{cc} \sigma'_{{x}} & \tau'_{{\it xy}}\\ \noalign{\medskip}\tau'_{{\it xy}} & \sigma'_{{y}} \end{array}\right] (5)

The rotation matrix is normally associated with Wallace Givens, who taught at the University of Tennessee at Knoxville. The primed values represent the stresses in the rotated coordinate system. We can rewrite the rotation matrix as follows, to correspond with the notation given above:^

\left[\begin{array}{cc} l & -m\\ m & l \end{array}\right]\left[\begin{array}{cc} \sigma_{{x}} & \tau_{{\it xy}}\\ \noalign{\medskip}\tau_{{\it xy}} & \sigma_{{y}} \end{array}\right]=\left[\begin{array}{cc} \sigma'_{{x}} & \tau'_{{\it xy}}\\ \noalign{\medskip}\tau'_{{\it xy}} & \sigma'_{{y}} \end{array}\right]\ (6)

2.1 Eigenvalues, Eigenvectors and Principal Stresses

The second question is this: is there an angle (or set of direction cosines) where the shear stresses would go away, leaving only normal stresses? Because of the properties of the matrix, the answer to this is also “yes,” and involves the process of diagonalizing the matrix. The diagonalized matrix (the matrix with only non-zero values along the diagonal) can be found if the eigenvalues of the matrix can be found, i.e., if the following equation can be solved for \lambda :

\left[\begin{array}{cc} \sigma_{{x}}-\lambda & \tau_{{\it xy}}\\ \noalign{\medskip}\tau_{{\it xy}} & \sigma_{{y}}-\lambda \end{array}\right]=0\ (7)

To accomplish this, we take the determinant of the left hand side of Equation 7, namely

{\lambda}^{2}-\left(\sigma_{{x}}+\sigma_{{y}}\right)\lambda-\left({\tau_{{\it xy}}}^{2}-\sigma_{x}\sigma_{y}\right)=0 (8)

If we define

J_{1}=\sigma_{{x}}+\sigma_{{y}} (9)


J{}_{2}={\tau_{{\it xy}}}^{2}-\sigma_{x}\sigma_{y} (10)

Equation 8 can be rewritten as

\lambda^{2}-J_{1}\lambda-J_{2}=0 (11)

The quantities J_{1},\,J_{2} are referred to as the invariants; they do not change as the axes are rotated. The first invariant is the trace of the system, which was predicted to be invariant earlier. These are very important,

especially in finite element analysis, where failure criteria are frequently computed relative to the invariants and not the principal stresses in any combination. This is discussed at length in Owen and Hinton [6].

This is the characteristic polynomial of the matrix of Equation 7. Most people generally don’t associate matrices and polynomials, but every matrix has an associated characteristic polynomial, and conversely a polynomial can have one or more corresponding matrices.The solution of Equation 8 produces the two eigenvalues of the matrix. The first eigenvalue of the matrix is

\lambda_{1}=1/2\,\sigma_{{y}}+1/2\,\sigma_{{x}}+1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}  (12)

and the second

\lambda_{1}=1/2\,\sigma_{{y}}+1/2\,\sigma_{{x}}-1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}} (13)

We can also determine the eigenvectors from this. Without going into the process of determining these,the first eigenvector is

\bar{x}_{1}=\left[\begin{array}{c} {\frac{-1/2\,\sigma_{{y}}+1/2\,\sigma_{{x}}+1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}}{\tau_{{\it xy}}}}\\ 1 \end{array}\right]\   (14)

and the second is

\bar{x}_{2}=\left[\begin{array}{c} {\frac{-1/2\,\sigma_{{y}}+1/2\,\sigma_{{x}}-1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}}{\tau_{{\it xy}}}}\\ 1 \end{array}\right]   (15)

One thing we can do with these eigenvectors is to normalize them, i.e., have it so that the sum of the squares of the entries is unity. Doing this yields

\bar{x}_{1}=\left[\begin{array}{c} {\frac{-1/2\,\sigma_{{y}}{\tau_{{\it xy}}}^{-1}+1/2\,\sigma_{{x}}{\tau_{{\it xy}}}^{-1}+1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}{\tau_{{\it xy}}}^{-1}}{\sqrt{1/2\,{\frac{{\sigma_{{y}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-{\frac{\sigma_{{x}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}-1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{{\sigma_{{x}}}^{2}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{x}}}{{\tau_{{\it xy}}}^{2}}}+2}}}\\ {\frac{1}{\sqrt{1/2\,{\frac{{\sigma_{{y}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-{\frac{\sigma_{{x}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}-1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{{\sigma_{{x}}}^{2}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{x}}}{{\tau_{{\it xy}}}^{2}}}+2}}} \end{array}\right] (16)


\bar{x}_{2}=\left[\begin{array}{c} {\frac{-1/2\,\sigma_{{y}}{\tau_{{\it xy}}}^{-1}+1/2\,\sigma_{{x}}{\tau_{{\it xy}}}^{-1}-1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}{\tau_{{\it xy}}}^{-1}}{\sqrt{1/2\,{\frac{{\sigma_{{y}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-{\frac{\sigma_{{x}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{{\sigma_{{x}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{x}}}{{\tau_{{\it xy}}}^{2}}}+2}}}\\ {\frac{1}{\sqrt{1/2\,{\frac{{\sigma_{{y}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-{\frac{\sigma_{{x}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{{\sigma_{{x}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{x}}}{{\tau_{{\it xy}}}^{2}}}+2}}} \end{array}\right] (17)

The eigenvectors, complicated as the algebra is, are useful in that we can diagonalize the matrix as follows:

S^{-1}AS=\Lambda (18)


S=\left[\begin{array}{cc} \bar{x}_{1_{1}} & \bar{x}_{2_{1}}\\ \bar{x}_{1_{2}} & \bar{x}_{2_{2}} \end{array}\right] (19)

This is referred to as a similarity or collinearity transformation. Similar matrices are matrices with the same eigenvalues. The matrix \Lambda , although simpler in form than A , has the same eigenvalues as the “original”matrix.

Either the original or the normalized forms of the eigenvectors can be used; the result will be the same as the scalar multiples will cancel in the matrix inversion. The normalization was done to illustrate the
relationship between the eigenvectors, the diagonalization matrix, and the Givens rotation matrix, since for the last sin^{2}\alpha+cos^{2}\alpha=1, an automatically normalized relationship. It is thus possible to extract the angle of the principal stresses from the eigenvectors.

Inverting the result of Equation 19 and multiplying through Equation 18 yields at last

\Lambda=\left[\begin{array}{cc} \lambda_{1} & 0\\ 0 & \lambda_{2} \end{array}\right]=\left[\begin{array}{cc} \sigma_{1} & 0\\ 0 & \sigma_{3} \end{array}\right] (20)

The two eigenvalues from Equations 12 and 13 are the principal stresses at the stress point in question.(The use of different subscripts for the eigenvalues and principal stresses comes from too many years dealing with Mohr-Coulomb failure theory.) One practical result of Equations 9 and 20 and the underlying theory is that, for any coordinate orientation,

\sigma_{x}+\sigma_{y}=\sigma_{1}+\sigma_{3} (21)
This is a very handy check when working problems such as this, especially since the algebra is so involved.

2.2. Two-Dimensional Example

To see all of this “in action” consider the soil stress states shown in Figure 2.


Figure 2: Stress State Example (from Navy [5])

Consider the stress state “A.” The stresses are expressed as a ratio of a pressure P at the surface; furthermore, following geotechnical engineering convention, compressive stresses are positive and the ordinate is the “z” axis.. For this study the convention of Figure 1 will be adopted and thus \sigma_{x}=-0.77,\,\sigma_{y}=-0.98,\,\tau_{xy}=0.05   .By direct substitution (and carrying the results to precision unjustified in geotechnical engineering) we obtain the following:

  • \lambda_{1}=\sigma_{1}=-.7587029665P (Equations 12 and 20.)
  • \lambda_{2}=\sigma_{3}=-.9912970335P (Equations 12 and 20.) Note that the principal stresses are reversed; this is because the sign convention is reversed, and thus what was formerly the smaller of the stresses is now the larger. Also, the axis of the first principal stress changes because the first principal stress itself has changed.
  • S=\left[\begin{array}{cc} .9754128670 & -.220385433\\ .2203854365 & .9754128502 \end{array}\right] (Equation 19.) Note that the normalized values of the eigenvectors (Equations 16 and 17) are used. Also note the similarity between this matrix and the Givens rotation matrix of Equation 5.  The diagonalization process is simply a rotation process, as the physics of the problem suggest. (The diagonal terms should be equal to each other and the absolute values ofthe off-diagonal terms should be also; they are not because of numerical errors in Maple where they were computed.
  • \alpha=12.73167230^{\circ};\,\beta=77.26832747^{\circ} (Equations 2, 3 and 6.)
  • In both cases the trace of A and \Lambda is the same, namely 1.75P (Equation 21.)

The results are thus the same. The precision is obviously greater than the “slide rule era” Navy [5], but the results illustrate that numerical errors can creep in, even with digital computation.

3. Three-Dimensional Analysis

3.1. Regular Principal Stresses

It is easy to see that, although the concepts are relatively simple, the algebra is involved.It is for this reason that the two-dimensional state was illustrated. With three dimensional stresses, the principles are the same, but the algebra is even more difficult.To begin, one more direction cosine must be defined,

n = cos\gamma (22)

where \gamma is of course the angle of the direction of the rotated coordinate system relative to the original one. The rotated system does not have to be the principal axis system, although that one is of most interest.
For three dimensions, Equation 4 should be written as follows:

\left[\begin{array}{ccc} \sigma_{x} & \tau_{xy} & \tau_{xz}\\ \tau_{xy} & \sigma_{y} & \tau_{yz}\\ \tau_{xz} & \tau_{yz} & \sigma_{z} \end{array}\right]\left[\begin{array}{c} l\\ m\\ n \end{array}\right]=\left[\begin{array}{c} p_{x}\\ p_{y}\\ p_{z} \end{array}\right] (23)

All of the moment summations of the shear stresses-which result in the symmetry of the matrix-have been included. The eigenvalues of the matrix are thus the solution of

\left[\begin{array}{ccc} \sigma_{x}-\lambda & \tau_{xy} & \tau_{xz}\\ \tau_{xy} & \sigma_{y}-\lambda & \tau_{yz}\\ \tau_{xz} & \tau_{yz} & \sigma_{z}-\lambda \end{array}\right]=0 (24)

The characteristic equation of this matrix-and thus the equation that solves for the eigenvalues-is

\lambda^{3}-J_{1}\lambda^{2}-J_{2}\lambda-J_{3}=0 (25)


J_1=\sigma_{x}+\sigma_{y}+\sigma_{z} (26)

J_2=\tau_{xy}^{2}+\tau_{xz}^{2}+\tau_{yz}^{2}-\left(\sigma_{x}\sigma_{y}+\sigma_{x}\sigma_{z}+\sigma_{y}\sigma_{z}\right) (27)

J_3=\sigma_{x}\sigma_{y}\sigma_{z}+2\tau_{xy}\tau_{xz}\tau_{yz}-\left(\sigma_{x}\tau_{yz}^{2}+\sigma_{y}\tau_{xz}^{2}+\sigma_{z}\tau_{xy}^{2}\right) (28)

It is easy to show the following:

  • Equation 26 reduces to Equation 9 if \sigma_{z}=0 .
  • Equation 27 reduces to Equation 10 if additionally \tau_{xz}=\tau_{yz}=0 .
  • For all three conditions, Equation 25 reduces to 11.

From the previous derivation, we can thus expect

 \lambda_{1}=\sigma_{1} (29)
 \lambda_{2}=\sigma_{2} (30)
\lambda_{3}=\sigma_{3} (31)
Solving a cubic equation such as this can be accomplished using the Cardan Formula. Generally, it is possible that we will end up with complex roots, but because of the properties of the matrix we are promised real roots. The eigenvalues are as follows:

\lambda_{1}  =  1/6\,\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}\nonumber \\ -6\,{\frac{1/3\,J_{{2}}-1/9\,{J_{{1}}}^{2}}{\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}}}-1/3\,J_{{1}} (32)

\lambda_{2}  =  -1/12\,\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}\nonumber \\ +3\,{\frac{1/3\,J_{{2}}-1/9\,{J_{{1}}}^{2}}{\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}}}\nonumber \\ -1/3\,J_{{1}}+1/12\,\sqrt{-3}\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}\nonumber \\ +2\sqrt{-3}\,{\frac{1/3\,J_{{2}}-1/9\,{J_{{1}}}^{2}}{\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}}} (33)

\lambda_{3} =  -1/12\,\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}\nonumber \\ +3\,{\frac{1/3\,J_{{2}}-1/9\,{J_{{1}}}^{2}}{\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}}}\nonumber \\ -1/3\,J_{{1}}-1/12\,\sqrt{-3}\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}\nonumber \\ -2\sqrt{-3}\,{\frac{1/3\,J_{{2}}-1/9\,{J_{{1}}}^{2}}{\sqrt[3]{36\,J_{{2}}J_{{1}}-108\,J_{{3}}-8\,{J_{{1}}}^{3}+12\,\sqrt{12\,{J_{{2}}}^{3}-3\,{J_{{2}}}^{2}{J_{{1}}}^{2}-54\,J_{{2}}J_{{1}}J_{{3}}+81\,{J_{{3}}}^{2}+12\,J_{{3}}{J_{{1}}}^{3}}}}} (34)

3.2. Deviatoric Principal Stresses

The algebra of Equations 32, 33 and 34 is rather involved. Is there a way to simplify it? The answer is”yes” if we use a common reduction of the Cardan Formula. Consider the following change of variables:

 \sigma_{x}'=\sigma_{x}-\frac{J_{1}}{3} (35)

 \sigma_{y}'=\sigma_{y}-\frac{J{}_{1}}{3} (36)

\sigma_{z}'=\sigma_{z}-\frac{J{}_{1}}{3}  (37)

Since the change is the same in all directions, we can write

\left[\begin{array}{ccc} \sigma'_{x} & \tau_{xy} & \tau_{xz}\\ \tau_{xy} & \sigma'_{y} & \tau_{yz}\\ \tau_{xz} & \tau_{yz} & \sigma'_{z} \end{array}\right]\left[\begin{array}{c} l\\ m\\ n \end{array}\right]=\left[\begin{array}{c} p'_{x}\\ p'_{y}\\ p'_{z} \end{array}\right] (38)


\left[\begin{array}{ccc} \sigma'_{x}-\lambda' & \tau_{xy} & \tau_{xz}\\ \tau_{xy} & \sigma'_{y}-\lambda' & \tau_{yz}\\ \tau_{xz} & \tau_{yz} & \sigma'_{z}-\lambda' \end{array}\right]=0 (39)

The characteristic equation of the matrix in Equation 39 is

\lambda'^{3}-J'_{2}\lambda'-J'_{3}=0 (40)

J'_{2}=\tau_{xy}^{2}+\tau_{xz}^{2}+\tau_{yz}^{2}-\left(\sigma'_{x}\sigma'_{y}+\sigma'_{x}\sigma'_{z}+\sigma'_{y}\sigma'_{z}\right) (41)

J'_{3}=\sigma'_{x}\sigma'_{y}\sigma'_{z}+2\tau_{xy}\tau_{xz}\tau_{yz}-\left(\sigma'_{x}\tau_{yz}^{2}+\sigma'_{y}\tau_{xz}^{2}+\sigma'_{z}\tau_{xy}^{2}\right) (42)

Although our motivation for this was to simplify the characteristic equation, the deviatoric stress infact has physical significance. As Jaeger and Cook [4] point out, “…essentially (\frac{I_{1}}{3} ) determines uniform compression or dilation, while the stress deviation determines distortion. Since many criteria of failure are concerned primarily with distortion, and since they must be invariant with respect to rotation of axes, it appears that the invariants of stress deviation will be involved.”

The one thing to be careful about is that the eigenvalues from Equation 24 are not the same as those from Equation 39. The latter are in fact “deviatoric eigenvalues” and are equivalent to deviatoric principal stresses, thus

\lambda'_1 = \sigma'_1 (43)

\lambda'_2 = \sigma'_2 (44)

\lambda'_3 = \sigma'_3 (45)

These can be converted back to full stresses using Equations 35, 36 and 37. The eigenvalues for the deviatoric stresses are

\lambda'_{1}=1/6\,\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}+2\,{\frac{J'_{{2}}}{\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}}}   (46)

\lambda'_{2} = -1/12\,\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}-{\frac{J'_{{2}}}{\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}}}\nonumber \\ +1/2\,\sqrt{-1}\sqrt{3}\left(1/6\,\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}-2\,{\frac{J'_{{2}}}{\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}}}\right) (47)

\lambda'_{3} = -1/12\,\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}-{\frac{J'_{{2}}}{\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}}}\nonumber \\ -1/2\,\sqrt{-1}\sqrt{3}\left(1/6\,\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}-2\,{\frac{J'_{{2}}}{\sqrt[3]{108\,J'_{{3}}+12\,\sqrt{-12\,{J'_{{2}}}^{3}+81\,{J'_{{3}}}^{2}}}}}\right) (48)

Solving Equations 24 and 39 will yield eigenvalues and principal stresses of one kind or another. The eigenvectors and direction cosines can be determined using the same methods used for the two-dimensional problem. It should be readily evident, however, that the explicit solution of these eivenvectors and diagonalizing rotational matrices is very involved, although using deviatoric stresses simplifies the algebra considerably. If the material is isotropic, and its properties are the same in all directions, than the direction of the principal stresses can frequently be neglected. For anisotropic materials, this is not the case.

4. Finding the Eigenvalues and Eigenvectors Using the Method of Danilevsky

For the simple 3 * 3 matrix, computing the determinant, and from it the the invariants, is not a major task. With larger matrices, it is far more difficult to find the determinant, let alone the characteristic polynomial.There have been many solutions to this problem over the years. For larger matrices the most common are Householder reflections and/or Givens rotations. We discussed the latter earlier. By doing multiple reflections and/or rotations, the matrix can be diagonalized and the eigenvalues “magically appear.” The whole problem of solving the characteristic equation is avoided, although the iterative cost in getting to the result can be considerable.

The method we’ll explore here is that of Danilevsky, first published in 1937, as presented in Faddeeva [3].

4.1. Theory

The math will be presented in 3\times3 format, but it can be expanded to any size square matrix. The idea is, through a series of similarity transformations, to transform a 3\times3 matrix (such as is shown in Equation 23) to the Frobenius normal form, which is

D=\left[\begin{array}{ccc} p_{{1}} & p_{{2}} & p_{{3}}\\ \noalign{\medskip}1 & 0 & 0\\ \noalign{\medskip}0 & 1 & 0 \end{array}\right] (49)

From here we take the determinant of D-\lambda I , thus

D\left(\lambda\right)=\left|\begin{array}{ccc} p_{1}-\lambda & p_{2} & p_{3}\\ 1 & -\lambda & 0\\ 0 & 1 & -\lambda \end{array}\right| (50)

The determinant-and thus the characteristic polynomial-is easy to compute, and setting it equal to zero,

{\lambda}^{3}-{\lambda}^{2}p_{{1}}-p_{{2}}\lambda-p_{{3}}=0 (51)

This is obviously the same as Equation 25, but now we can read the invariants straight out of the similar matrix.

So how do we get to Equation 49? The simple answer is through a series of row reductions. We start for the matrix A by carrying the third row into the second, dividing all the elements by the element in the last row by the element a_{3,2} , then subtracting the second column, multiplied by a_{3,1},\,a_{3,2},\,a_{3,3} respectively from all the rest of the columns. Row and column manipulation is common with computational matrix routines, but a more straightforward way is to postmultiply A by the matrix

M=\left[\begin{array}{ccc} 1 & 0 & 0\\ \noalign{\medskip}-{\frac{a_{{3,1}}}{a_{{3,2}}}} & {a_{{3,2}}}^{-1} & -{\frac{a_{{3,3}}}{a_{{3,2}}}}\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right] (52)

The result will not be similar to A , but if we premultiply that result by M^{-1} ,

M^{-1}=\left[\begin{array}{ccc} 1 & 0 & 0\\ \noalign{\medskip}a_{{3,1}} & a_{{3,2}} & a_{{3,3}}\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right] (53)

that result will be by Equation 19, although M does not diagonalize the result the way S did. Performing this series of matrix multplications yields

B=M^{-1}AM=\left[\begin{array}{ccc} a_{{1,1}}-{\frac{a_{{1,2}}a_{{3,1}}}{a_{{3,2}}}} & {\frac{a_{{1,2}}}{a_{{3,2}}}} & -{\frac{a_{{1,2}}a_{{3,3}}}{a_{{3,2}}}}+a_{{1,3}}\\ \noalign{\medskip}a_{{3,1}}\left(a_{{1,1}}-{\frac{a_{{1,2}}a_{{3,1}}}{a_{{3,2}}}}\right)+a_{{3,2}}\left(a_{{2,1}}-{\frac{a_{{2,2}}a_{{3,1}}}{a_{{3,2}}}}\right) & {\frac{a_{{1,2}}a_{{3,1}}}{a_{{3,2}}}}+a_{{2,2}}+a_{{3,3}} & a_{{3,1}}\left(-{\frac{a_{{1,2}}a_{{3,3}}}{a_{{3,2}}}}+a_{{1,3}}\right)+a_{{3,2}}\left(-{\frac{a_{{2,2}}a_{{3,3}}}{a_{{3,2}}}}+a_{{2,3}}\right)\\ \noalign{\medskip}0 & 1 & 0 \end{array}\right] (54)

We continue this process by moving up a row, thus the new similarity transformation matrices are

N  =  \left[\begin{array}{ccc} {b_{{2,1}}}^{-1} & -{\frac{b_{{2,2}}}{b_{{2,1}}}} & -{\frac{b_{{2,3}}}{b_{{2,1}}}}\\ \noalign{\medskip}0 & 1 & 0\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right] (55)

N^{-1}  =  \left[\begin{array}{ccc} b_{{2,1}} & b_{{2,2}} & b_{{2,3}}\\ \noalign{\medskip}0 & 1 & 0\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right] (56)

Although we formed the M and N matrices first, from a conceptual and computational standpoint itis easier to form the “inverse” matrices first and then invert them. One drawback to Danilevsky’s method is that there are many opportunities for division by zero, which need to be taken into consideration when employing the method.

In any case the next step yields

C=N^{-1}BN=\left[\begin{array}{ccc} b_{{1,1}}+b_{{2,2}} & b_{{2,1}}\left(-{\frac{b_{{1,1}}b_{{2,2}}}{b_{{2,1}}}}+b_{{1,2}}\right)+b_{{2,3}} & b_{{2,1}}\left(-{\frac{b_{{1,1}}b_{{2,3}}}{b_{{2,1}}}}+b_{{1,3}}\right)\\ \noalign{\medskip}1 & 0 & 0\\ \noalign{\medskip}0 & 1 & 0 \end{array}\right] (57)

From this Equation 49 is easily formed.
The coefficients p_1,\,p_2 and p_3 are simply the invariants, which have already been established earlier for both standard and deviatoric normal stresses. From that standpoint the method seems to be overkill.

The advantage comes when we compute the eigenvectors/direction cosines. Consider the vector

y=\left[\begin{array}{c} y_{1}\\ y_{2}\\ y_{3} \end{array}\right] (58)

Premultiplying it by D-\lambda I (see Equation 50) yields

\left[\begin{array}{ccc} p_{1}-\lambda & p_{2} & p_{3}\\ 1 & -\lambda & 0\\ 0 & 1 & -\lambda \end{array}\right]\left[\begin{array}{c} y_{1}\\ y_{2}\\ y_{3} \end{array}\right]=\left[\begin{array}{c} \left(p_{{1}}-\lambda\right)y_{{1,1}}+p_{{2}}y_{{2,1}}+p_{{3}}y_{{3,1}}\\ \noalign{\medskip}y_{{1,1}}-\lambda\,y_{{2,1}}\\ \noalign{\medskip}y_{{2,1}}-\lambda\,y_{{3,1}} \end{array}\right]=0 (59)

For y to be non-trivial,

y=\left[\begin{array}{c} \lambda^{2}\\ \lambda\\ 1 \end{array}\right] (60)

It can be shown (Faddeeva [3]) that the eigenvectors x_{n} can be found as follows (successive rotations:)

\bar{x}_{n}=MNy=\left[\begin{array}{c} {\frac{{\lambda_{n}}^{2}}{b_{{2,1}}}}-{\frac{b_{{2,2}}\lambda_{n}}{b_{{2,1}}}}-{\frac{b_{{2,3}}}{b_{{2,1}}}}\\ \noalign{\medskip}-a_{{3,1}}\left({\frac{{\lambda_{n}}^{2}}{b_{{2,1}}}}-{\frac{b_{{2,2}}\lambda_{n}}{b_{{2,1}}}}-{\frac{b_{{2,3}}}{b_{{2,1}}}}\right){a_{{3,2}}}^{-1}+{\frac{\lambda_{n}}{a_{{3,2}}}}-{\frac{a_{{3,3}}}{a_{{3,2}}}}\\ \noalign{\medskip}1 \end{array}\right] (61)

These can, as was shown earlier, be normalized to direction cosines and can form the diagonalizing matrices.

4.2. Three-Dimensional Example

An example of this is drawn from Boresi et al. [1]. It involves a stress state, which referring to Equation 23 can be written (all stresses in kPa) as

A=\left[\begin{array}{ccc} 120 & -55 & -75\\ \noalign{\medskip}-55 & 55 & 33\\ \noalign{\medskip}-75 & 33 & -85 \end{array}\right] (62)

For the first transformation, we substitute these values into Equation 54 and thus

B=\left[\begin{array}{ccc} 1 & 0 & 0\\ \noalign{\medskip}-75 & 33 & -85\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right]\left[\begin{array}{ccc} 120 & -55 & -75\\ \noalign{\medskip}-55 & 55 & 33\\ \noalign{\medskip}-75 & 33 & -85 \end{array}\right]\left[\begin{array}{ccc} 1 & 0 & 0\\ \noalign{\medskip}{\frac{25}{11}} & 1/33 & {\frac{85}{33}}\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right]=\left[\begin{array}{ccc} -5 & -5/3 & -{\frac{650}{3}}\\ \noalign{\medskip}2685 & 95 & 22014\\ \noalign{\medskip}0 & 1 & 0  \end{array}\right] (63)

The second transformation, following Equation 57,

C=\left[\begin{array}{ccc} 2685 & 95 & 22014\\ \noalign{\medskip}0 & 1 & 0\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right]\left[\begin{array}{ccc} -5 & -5/3 & -{\frac{650}{3}}\\ \noalign{\medskip}2685 & 95 & 22014\\ \noalign{\medskip}0 & 1 & 0 \end{array}\right]\left[\begin{array}{ccc} {\frac{1}{2685}} & -{\frac{19}{537}} & -{\frac{7338}{895}}\\ \noalign{\medskip}0 & 1 & 0\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right]=\left[\begin{array}{ccc} 90 & 18014 & -471680\\ \noalign{\medskip}1 & 0 & 0\\ \noalign{\medskip}0 & 1 & 0 \end{array}\right] (64)

From this the characteristic equation is, from Equation 51,

{\lambda}^{3}-90\,{\lambda}^{2}-18014\,\lambda+471680=0 (65)

which yields eigenvalues/principal stresses of 176.8, 24.06 and -110.86 kPa.

So how does this look with deviatoric stresses? Applying Equations 26, 35, 36 and 37 yields

A'=\left[\begin{array}{ccc} 90 & -55 & -75\\ \noalign{\medskip}-55 & 25 & 33\\ \noalign{\medskip}-75 & 33 & -115 \end{array}\right] (66)

It’s worth noting that, for the deviatoric stress matrix, the trace is always zero. The two similarity transformations are thus:

B'  =  \left[\begin{array}{ccc} 1 & 0 & 0\\ \noalign{\medskip}-75 & 33 & -115\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right]\left[\begin{array}{ccc} 90 & -55 & -75\\ \noalign{\medskip}-55 & 25 & 33\\ \noalign{\medskip}-75 & 33 & -115 \end{array}\right]\left[\begin{array}{ccc} 1 & 0 & 0\\ \noalign{\medskip}{\frac{25}{11}} & 1/33 & {\frac{115}{33}}\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right]=\left[\begin{array}{ccc} -35 & -5/3 & -{\frac{800}{3}}\\ \noalign{\medskip}2685 & 35 & 23964\\ \noalign{\medskip}0 & 1 & 0 \end{array}\right] (67)

C' =  \left[\begin{array}{ccc} 2685 & 35 & 23964\\ \noalign{\medskip}0 & 1 & 0\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right]\left[\begin{array}{ccc} -35 & -5/3 & -{\frac{800}{3}}\\ \noalign{\medskip}2685 & 35 & 23964\\ \noalign{\medskip}0 & 1 & 0 \end{array}\right]\left[\begin{array}{ccc} {\frac{1}{2685}} & -{\frac{7}{537}} & -{\frac{7988}{895}}\\ \noalign{\medskip}0 & 1 & 0\\ \noalign{\medskip}0 & 0 & 1 \end{array}\right]=\left[\begin{array}{ccc} 0 & 20714 & 122740\\ \noalign{\medskip}1 & 0 & 0\\ \noalign{\medskip}0 & 1 & 0 \end{array}\right] (68)

and the characteristic polynomial (missing the squared term, as expected)

\lambda^3-20714\,\lambda-122740=0 (69)

which yields eigenvalues/principal deviatoric stresses of 146.8, -5.94 and -140.86 kPa.

4.3. Numerical Implementation

Even with a “simple” problem such as this, the algebra of the solution can become very involved. Moreover, as is frequently the case with numerical linear algebra, there are two “loose ends” that need to be tied up: the division by zero problem in Danilevsky’s Method, and solving for the eigenvalues/principal stresses. A numerical implementation, shown below, will address both issues, if not completely.

The FORTRAN 77 code is in the appendix. It is capable of solving the problem both for standard and deviatoric stresses. The example problem given above is used. The steps used are as follows:

  1. If deviatoric stresses are specified, the normal stresses are accordingly converted as shown earlier.
  2. The invariants are computed. For a more general routine of Danilevsky’s Method, the matrices associated with the similarity transformations-and their inverses-would have to be computed and the matrix multiplications done. In this case the results of these multiplications-the invariants themselves-are “hard coded” into the routine, which saves a great deal of matrix multiplication. It also eliminated problems with division by zero up to this point.
  3. The eigenvalues/principal stresses are computed using Newton’s Method. Obviously they could be computed analytically as shown earlier. For this type of problem, this is possible; for characteristic polynomials beyond degree four, another type of solution is necessary. The major problem with New-ton’s Method is picking initial values so as to yield distinct eigenvalues. The nature of the problem suggests that the initial normal, diagonal stresses can be used as starting values. For the sample problem this worked out very well; it was not tested on a variety of stress states.
  4. Taking the three resulting eigenvalues, the eigenvectors were computed using Equation 61. Again the matrix multiplications were hard-coded into the routine. The problem of division by zero of certain values of shear stresses was not addressed. In some cases Equation 61 will result in division by zero.The eigenvectors were then normalized, which yielded three sets of direction cosines.

The results for the regular stresses are shown in Figure 3 and those for deviatoric stresses are shown in Figure 4.


Figure 3: Regular Normal Stress Results for Test Case


Figure 4: Deviatoric Normal Stress Results for Test Case

The principal stress results agree with the “hand” calculations earlier. The direction cosines-and thus the diagonalizing matrices-are the same for both cases.

5. Conclusion

We have explored the analysis of Mohr’s Circle for stress states in both two and three dimensions.The equations can be derived from strictly linear algebra considerations, the principal stresses being the eigenvalues and the direction cosines being the normalized eigenvectors. It was also shown that numericalmethods can be employed to analyze these results and the invariants, in this case using the computationally efficient Danilevsky’s Method.

Note: since we first published this on a companion site two years ago, we have made several improvements to the methods described here.  These can be found in this post.

Appendix: Source Code for Numerical Implementation

! Routine to Determine Stress State from Danilevsky's Method
! Taken from Faddeeva (1959) with Example Case from Boresi (1993)
! Define Variables for Test Case
CHARACTER*12 filnam
! Set Variable for Eigenvalue Solution
! ityp = 0 Regular
! ityp = 1 Deviatoric
WRITE (6,*) 'Standard (0) or Deviatoric (1) Stress:'
READ (5,*) ityp
IF (ityp.eq.0) filnam='danil0.txt'
IF (ityp.eq.1) filnam='danil1.txt'
! Open Output File
OPEN (7,file=filnam)
CALL danil (sigmax, sigmay, sigmaz, tauxy, tauyz, tauxz, ityp)

! Subroutine danil -- Danilevsky's Method specifically written
! for 3 x 3 matrices and stress-state solution
SUBROUTINE danil (sigmax, sigmay, sigmaz, tauxy, tauyz, tauxz,
REAL p(3),x(3,3),sigma(3),i1,lambda
! Convert Normal Stresses to Deviatoric Stresses if ityp = 1
IF (ityp.eq.1) THEN
! Compute Invariants
IF (ityp.eq.0) p(1)=sigmax+sigmay+sigmaz
IF (ityp.eq.1) p(1)=0.0
! Use Newton's Method to Determine Eigenvalues
DO 60 j=1,3,1
GO TO (10,20,30),j
10 lambda=sigmax
GO TO 40
20 lambda=sigmay
GO TO 40
30 lambda=sigmaz
40 DO 50 i=1,10,1
IF ( GO TO 60
60 sigma(j)=lambda
! Determine Eigenvectors
DO 80 i=1,3,1
& tauyz*sigmay*tauxz)*tauyz*lambda**2+(tauxy*tauxz+sigmay*tauyz+
& sigmaz*tauyz)/(-tauxz*sigmax*tauyz+tauxy*tauxz**2-tauxy*tauyz**
& 2+tauyz*sigmay*tauxz)*lambda-(tauxz*tauxy*sigmaz-tauyz*tauxz**
& 2+tauyz*sigmay*sigmaz-tauyz**3)/(-tauxz*sigmax*tauyz+tauxy*
& tauxz**2-tauxy*tauyz**2+tauyz*sigmay*tauxz)
& tauxy*tauyz**2+tauyz*sigmay*tauxz)*tauyz*lambda**2+(tauxy*
& tauxz+sigmay*tauyz+sigmaz*tauyz)/(-tauxz*sigmax*tauyz+tauxy*
& tauxz**2-tauxy*tauyz**2+tauyz*sigmay*tauxz)*lambda-(tauxz*
& tauxy*sigmaz-tauyz*tauxz**2+tauyz*sigmay*sigmaz-tauyz**3)/(-
& tauxz*sigmax*tauyz+tauxy*tauxz**2-tauxy*tauyz**2+tauyz*sigmay*
& tauxz))+1/tauyz*lambda-sigmaz/tauyz
! Normalise Eigenvectors
DO 70 j=1,3,1
70 x(j,i)=x(j,i)*xnorm
! Write Results
WRITE (7,*) 'Stress State Using Danilevskys Method:'
IF (ityp.eq.0) WRITE (7,*) 'Regular Normal Stresses'
IF (ityp.eq.1) WRITE (7,*) 'Deviatoric Normal Stresses'
WRITE (7,*) 'Original Stress Matrix/Tensor:'
WRITE (7,100) sigmax,tauxy,tauxz
WRITE (7,100) tauxy,sigmay,tauyz
WRITE (7,100) tauxz,tauyz,sigmaz
WRITE (7,*) 'Invariants:'
WRITE (7,100) (p(j),j=1,3,1)
WRITE (7,*) 'Principal Stresses:'
WRITE (7,100) (sigma(j),j=1,3,1)
WRITE (7,*) 'Direction Cosines:'
DO 90 j=1,3,1
90 WRITE (7,110) (x(j,i),i=1,3,1)
100 FORMAT (3f15.2)
110 FORMAT (3f15.3)


[1] Boresi, A. P., Schmidt, R. J., Sidebottom, O. M., 1993. Advanced Mechanics of Materials, 5th Edition.John Wiley & Sons, Inc.

[2] DeRusso, P., Roy, R., Close, C., 1965. State Variables for Engineers. John Wiley & Sons, Inc., New York,NY.

[3] Faddeeva, V., 1959. Computational Methods of Linear Algebra. Dover Publications, Inc.

[4] Jaeger, J., Cook, N., 1979. Fundamentals of Rock Mechanics, 3rd Edition. Chapman and Hall.

[5] Navy, U., 1986. Soil mechanics. Tech. Rep. DM 7.01.

[6] Owen, D., Hinton, E., 1980. Finite Elements in Plasticity: Theory and Practice. Pineridge Press, Swansea,Wales.

[7] Strang, G., 1993. Introduction to Linear Algebra. Wellesley-Cambridge Press.

[8] Verruijt, A., van Bars, S., 2007. Soil Mechanics. VSSD, The Netherlands.

When the Pathfinder Gets Lost

Most people are familiar with the America’s Cup, the sailboat race which has attracted prominent competitors such as Sir Thomas Lipton and Ted Turner.  What most people don’t know is that there is a fresh water counterpart to that race, the Canada’s Cup, sailed on the Great Lakes.  The genesis of that race, how the Americans lost it the first time and gained it back the second, is one of the most interesting in yachting history.

In 1896 the Lincoln Park Yacht Club in Chicago challenged the Royal Canadian Yacht Club in Toronto to a cross-Great Lakes race.  The first race was held that year in Toledo, OH.  The Canadians were victorious, thus they named the race the Canada’s Cup.


The Pathfinder at its birthplace, the Racine Boat Mfg. Co.

In the same year, another member of the Lincoln Park club, Chicago rubber magnate Fred W. Morgan, commissioned the building of his new yacht, the Pathfinder.  Built in Racine, WI, its maiden voyage to Chicago was “the event of the season,” and people gathered on the docks to watch it glide into the harbour.

The Pathfinder was 140′ long (a sizeable yacht then or now) and resembled the U.S. Navy’s battleships of the day.   It had telephones and electric lighting in an era when the vast majority of American homes had neither.


The Pathfinder, underway in the Great Lakes.

Two years later the Chicago Yacht Club, with Fred Morgan as its Commodore and the Pathfinder the club’s flagship, issued a challenge to the RCYC for the Canada’s Cup.  The following year the Pathfinder steamed from Chicago to Toronto for the rematch.

The first day of racing was 21 August 1899.  The object was to start near the RCYC, round two buoys, and return to the club.  The Canadians’ yacht, the Beaver, had an accident right at the start and was out of the race.   The American yacht, the Genessee, had a chance for a default for the first race (it was two out of three to win.)

But such was not to be.  The Pathfinder was acting as a kind of “pace yacht” for the race, but in the haze on Lake Ontario itself got lost and missed the first buoy.  Behind it was not only the Genessee but steam yachts Siren and Canada!  The entire entourage mistook another buoy for the first official one, where the judges’ ship waited in vain.

Because of this the Americans missed the chance for a forfeit.  We eventually (and when I say “we” I mean not only the Americans but the Chicago Yacht Club, where my great-grandfather was Commodore the following year) won the race, but only by 31 seconds on a five hour course.

In those days it was easy for everyone to follow the largest, most magnificent yacht on the Great Lakes.  But they still got lost.  Unfortunately things haven’t changed as much as we would like to think.

Our country and our world is directed by those who are supposedly the most educated and enlightened among us.  Yet we still experienced the crash of 2008 and its aftermath.  Many of us were participants in that by virtue of taking out loans on our houses and just about everything else we owned (and didn’t in some cases) because were were informed by “knowledgeable” people that we could afford it.  That’s little comfort now that those homes and possessions are going away in foreclosure and repossession, especially since many of us are without income.

We have also been directed for the last forty years by a “knowledge class” that assured us that our “old ways” were hopelessly passé and that we would be happier in our new family structures.  In a country with high divorce and incarceration rates and single-parent poverty, their self-confident assurances aren’t much comfort either.

In the end the only one we can trust for the truth is he that is the truth, Jesus Christ.   Through him all things were created, so he is most knowledgeable as to what we need.  And, of course, he is the real “pathfinder” from the ultimate challenge of life, death itself, as we walked out of the tomb after those who didn’t care for his challenge of their authority had him executed.

It’s time for all of us who have followed what looked to be the biggest and most magnificent thing “on the water” to turn to he who actually walked on it.

For more information click here.

Use of Netwon’s Method to Determine Matrix Eigenvalues

The problem here is to develop a routine that will determine one or more eigenvalues of a matrix using Newton’s method and considering the eigenvalue problem to be that of a nonlinear equation solution problem.

The simplest way to illustrate the problem and its solution is to use a 4 \times4 matrix.

Let us consider a square matrix A. The definition of an eigenvalue requires that

Ax=\lambda x

For our test case writing out the solutions for the left and right hand sides yields

\left[\begin{array}{c} a_{{1,1}}x_{{1}}+a_{{1,2}}x_{{2}}+a_{{1,3}}x_{{3}}+a_{{1,4}}x_{{4}}\\ a_{{2,1}}x_{{1}}+a_{{2,2}}x_{{2}}+a_{{2,3}}x_{{3}}+a_{{2,4}}x_{{4}}\\ a_{{3,1}}x_{{1}}+a_{{3,2}}x_{{2}}+a_{{3,3}}x_{{3}}+a_{{3,4}}x_{{4}}\\ a_{{4,1}}x_{{1}}+a_{{4,2}}x_{{2}}+a_{{4,3}}x_{{3}}+a_{{4,4}}x_{{4}} \end{array}\right]=\left[\begin{array}{c} \lambda\, x_{{1}}\\ \lambda\, x_{{2}}\\ \lambda\, x_{{3}}\\ \lambda\, x_{{4}} \end{array}\right]

The problem here is that, with the intrusion of the eigenvalue, we
have one more unknown than we have equations for. We can remedy this
by insisting that the values of x be normalized, or


Let us now construct a vector such that we have a closed system, all on the left hand side:

\left[\begin{array}{c} a_{{1,1}}x_{{1}}+a_{{1,2}}x_{{2}}+a_{{1,3}}x_{{3}}+a_{{1,4}}x_{{4}}-\lambda\, x_{{1}}\\ a_{{2,1}}x_{{1}}+a_{{2,2}}x_{{2}}+a_{{2,3}}x_{{3}}+a_{{2,4}}x_{{4}}-\lambda\, x_{{2}}\\ a_{{3,1}}x_{{1}}+a_{{3,2}}x_{{2}}+a_{{3,3}}x_{{3}}+a_{{3,4}}x_{{4}}-\lambda\, x_{{3}}\\ a_{{4,1}}x_{{1}}+a_{{4,2}}x_{{2}}+a_{{4,3}}x_{{3}}+a_{{4,4}}x_{{4}}-\lambda\, x_{{4}}\\ {x_{{1}}}^{2}+{x_{{2}}}^{2}+{x_{{3}}}^{2}+{x_{{4}}}^{2}-1 \end{array}\right]=0

Since the object of Newton’s method is to arrive at this, we need an intermediate vector F, or

F=\left[\begin{array}{c} a_{{1,1}}x_{{1}}+a_{{1,2}}x_{{2}}+a_{{1,3}}x_{{3}}+a_{{1,4}}x_{{4}}-\lambda\, x_{{1}}\\ a_{{2,1}}x_{{1}}+a_{{2,2}}x_{{2}}+a_{{2,3}}x_{{3}}+a_{{2,4}}x_{{4}}-\lambda\, x_{{2}}\\ a_{{3,1}}x_{{1}}+a_{{3,2}}x_{{2}}+a_{{3,3}}x_{{3}}+a_{{3,4}}x_{{4}}-\lambda\, x_{{3}}\\ a_{{4,1}}x_{{1}}+a_{{4,2}}x_{{2}}+a_{{4,3}}x_{{3}}+a_{{4,4}}x_{{4}}-\lambda\, x_{{4}}\\ {x_{{1}}}^{2}+{x_{{2}}}^{2}+{x_{{3}}}^{2}+{x_{{4}}}^{2}-1 \end{array}\right]

Let us make \lambda the last “x” or


The Jacobian of this vector with respect to the expanded x vector is

J\left(F,x\right)=\left[\begin{array}{ccccc} a_{{1,1}}-\lambda & a_{{1,2}} & a_{{1,3}} & a_{{1,4}} & -x_{{1}}\\ \noalign{\medskip}a_{{2,1}} & a_{{2,2}}-\lambda & a_{{2,3}} & a_{{2,4}} & -x_{{2}}\\ \noalign{\medskip}a_{{3,1}} & a_{{3,2}} & a_{{3,3}}-\lambda & a_{{3,4}} & -x_{{3}}\\ \noalign{\medskip}a_{{4,1}} & a_{{4,2}} & a_{{4,3}} & a_{{4,4}}-\lambda & -x_{{4}}\\ \noalign{\medskip}2\, x_{{1}} & 2\, x_{{2}} & 2\, x_{{3}} & 2\, x_{{4}} & 0 \end{array}\right]

Now we can solve the Newton’s method equation iteratively:


It is certainly possible to invert the Jacobian symbolically, but the algebra becomes very complicated, even for a relatively small matrix such as this one. A more sensible solution is to obtain numerical values for both Jacobian and F vector, invert the former using the Gauss-Jordan technique and then multiply this by the F vector before performing the Newton iteration.

Now let us consider the case of the tridiagonal matrix

A=\left[\begin{array}{cccc} 2 & -1 & 0 & 0\\ \noalign{\medskip}-1 & 2 & -1 & 0\\ \noalign{\medskip}0 & -1 & 2 & -1\\ \noalign{\medskip}0 & 0 & -1 & 2 \end{array}\right]

The eigenvalues of this matrix are


or numerically 3.618033989, 1.381966011, 2.618033989, and .381966011. The FORTRAN 77 code used for this is given at the end of the problem. Although the matrix is “hard coded” into the program, it is only necessary to change one parameter to change the matrix size.


Summary Results of Newton’s Method, 4×4 Matrix

There are two different quantities tracked here:

  • The residual, i.e., the Euclidean norm of the entire x vector.
  • The error, i.e., the change in the eigenvalue from one step to the next.

Both these are tracked for three different starting places. For convenience,
all of the x vector entries were initialized at the same value.

Some comments on the solution are as follows:

  • As mentioned earlier, there were four eigenvalues to be determined. The method only returns one eigenvalue for each starting point, and that eigenvalue depends upon the starting point. The eigenvalues determined were 0.381966 (Start = 1 and Start = 2) and 2.6180339 (Start = 3.) Like the power method, this solution determines both an eigenvalue and eigenvector. Unlike the power method, it does not necessarily return the largest and/or smallest eigenvalue or even one easily predictable by the choice of starting value. With Newton’s Method, this is unsurprising; with multiple roots of an equation, which root is found very much depends upon the starting point. However, in the case where the number of roots is (in simple terms) the basic size of the matrix, this correspondence can become very complicated, and in some cases it is possible that certain eigenvalues will be missed altogether.
  • As a corollary to the previous point, with some starting values the convergence almost comes apart, especially where Start = 3. With larger values than this, the program generates “nan” values and terminates abnormally. This result is a part of the method; poor initial estimates can led successive iterations “off the track” and thus fail to converge. Thus it is necessary to know the range of eigenvalues before one actually determines them, which to a great extent defeats the purpose of the method.
  • The method is better at finding eigenvalues than finding eigenvectors. The residuals of the vector (which include the eigenvector plus the eigenvalue) converge much more slowly than the error. Those runs that were not limited by the termination criterion of the eigenvalue (|\Delta\lambda|<1.0\times10^{-6}) showed a very slow convergence continue with the eigenvectors. For the eigenvectors, the convergence rate is reasonable.

The general impression of this method, therefore, is that under proper circumstances it is capable of producing eigenvalues and (to a lesser extent) eigenvectors, but that it is necessary to a) have some idea of the values of the eigenvectors going in and b) be prepared for the method to collapse with the wrong initial guesses. One possible application (given the convergence limitations discussed above) is to use it in conjunction with, say, the Householder-Givens Method to determine the eigenvectors, which do not come out as a result of this method. The known eigenvalues give excellent starting values.

These results are confirmed when we expand the matrix to 10\times10.


Summary Results of Newton’s Method, 10×10 Matrix

We see that all of the errors and residuals go up and down from the initial guesses until convergence was achieved. The eigenvalues determined were 0.6902785 (Start = 1,) 1.7153704 (Start = 2,) and 1.7153703 (Start = 3.) Thus as before with three initial starting values only two eigenvalues were determined. The convergence was more lengthy than with the smaller matrix but not by much. This means that the method may be, to some extent, insensitive to matrix size, which would make it useful with larger matrices.


The code calls several “standard” routines and includes whose significance is as follows:

  • matsize.for is an include which includes a parameter statement for the actual matrix size. For the 10\times10 matrix, it read as follows:
    • parameter (nn=10, nplus1=nn+1,mcycle=100)
  • xmtinv is a subroutine to invert a matrix using Gauss-Jordan elimination.
  • matvec performs matrix-vector multiplication (in that order).
  • veclen computes the Euclidean norm of a vector.
 c Determination of Eigenvalue(s) of Tridiagonal Matrix using
 c Newton's Method
 include 'matsize.for'
 dimension a(nn,nn),x(nplus1),f(nplus1),fj(nplus1,nplus1),
 call tstamp
 c Define original matrix a
 do 1 i=1,nn,1
 do 1 j=1,nn,1
 if (i.eq.j) a(i,j) = 2.0
 if (abs(i-j).eq.1) a(i,j) = -1.0
 1 continue
 open (unit=2,file='m5610p2d.csv')
 write(2,*)'Original Matrix'
 do 2 i=1,nn,1
 2 write(2,3)(a(i,j),j=1,nn,1)
 3 format(100(g10.3,1h,))
 c Make initial guesses of eigenvalues and values of x
 do 100 mm=1,3,1
 write(2,*)'Eigenvalue Results for vstart =',vstart
 do 4 i=1,nplus1,1
 4 x(i)=vstart
 c Cycle Newton's method
 do 5 m=1,mcycle,1
 c Compute Jacobian for current step
 do 6 i=1,nplus1,1
 do 6 j=1,nplus1,1
 if(i.eq.nplus1.and.j.eq.nplus1) then
 goto 6
 if(i.eq.j) then
 goto 6
 if(j.eq.nplus1) then
 goto 6
 if(i.eq.nplus1) then
 goto 6
 6 continue
 c Output Jacobian
 write(2,*)'Jacobian Matrix'
 do 10 i=1,nplus1,1
 10 write(2,3)(fj(i,j),j=1,nplus1,1)
 c Compute Newton "F" vector
 do 20 i=1,nplus1,1
 21 f(i)=-x(nplus1)*x(i)
 do 23 j=1,nn,1
 23 continue
 goto 20
 22 f(i)=-1.0
 do 24 j=1,nn,1
 24 f(i)=f(i)+x(i)**2
 20 continue
 c Invert Jacobian
 call xmtinv(fj,nplus1,nplus1)
 c Postmultiply Jacobian by f vector
 call matvec(fj,f,xt,nplus1,nplus1)
 c Update Value of x vector and output the result
 write(2,*)'Result Vector for m = ',m
 do 7 k=1,nplus1,1
 7 write(2,*)x(k)
 c Compute norm of residual and error
 call veclen(xt,resid(m),nplus1)
 if(resid(m).lt.1.0e-07.or.eigerr(m).lt.1.0e-07)goto 30
 5 continue
 c Output residual plot
 30 write(2,*)'Residuals and Errors'
 101 format(27hIteration,Residual (Start =,i2,
 &16h),Error (Start = ,i2,1h))
 do 31 i=1,m-1,1
 31 write(2,32)i,resid(i),eigerr(i)
 32 format(i3,2(1h,,g10.3))
 100 continue