Taking the Last Voyage with Newton and Pascal

He’s not widely known outside of the fields he specialised in, but Adhémar Jean Claude Barré de Saint-Venant (1797-1886, usually known in the Anglophone world as simply Saint-Venant) was one of the premier scientists, engineers and mathematicians of the nineteenth century.  His accomplishments were many and include the following:

  • saint-venant-torsionSuccessful derivation of the Navier-Stokes Equations for a viscous flow before Stokes; these equations are the basis of computational fluid dynamics and the analysis of things that fly.
  • Systematisation and development of methods in the theory of elasticity of solids, including his semi-inverse methods for torsion, important in things such as automobile crank shafts.
  • Methods for the analysis of wave mechanics in bars, which we see in many places, from musical instruments to driven foundation piles.

Saint-Venant was born into a royalist, aristocratic, traditionally Roman Catholic family at a time when it was not safe to be any of these: the French Revolution, at that point stumbling from the Reign of Terror to control of France–and most of Europe–by Napoleon Bonaparte.  It was about the latter where Saint-Venant made a statement about himself that got him into trouble with the “new” Europe.  As described in S. Timoshenko’s History of Strength of Materials:

The political events of 1814 had a great effect on Saint-Venant’s career.  In March of that year, the armies of the allies were approaching Paris and the students of the École Polytechnique were mobilized.  On March 30, 1814, they were moving their guns to the Paris fortification when Saint-Venant, who was the first sargeant of the detachment, stepped out from the ranks with the exclamation: “My conscience forbids me to fight for an usurper…” His schoolmates resented that action very much and Saint-Venant was proclaimed a deserter and never allowed to resume his study at the École Polytechnique.

Saint-Venant’s statement of conscience was at once a political and religious statement, and “progressives” of his day didn’t miss either.  The French, then and now innocent of anti-discrimination legislation or sentiments, made his life miserable. The École Polytechnique was and is France’s premier technical institute of higher learning; getting kicked out of it was the equivalent of, say, being expelled from Princeton or MIT.  He worked in the powder industry for nine years, then was admitted to the École des Ponts et Chausées, where his fellow students shunned him.  He graduated first in his class anyway and began his illustrious career in technical things both theoretical and practical.

In spite of his difficulties within France, his reputation outside of her was another matter.  When François Napoleon Moigno wrote his book on statics, he discovered the following:

He (Moigno) wanted the portion on the statics of elastic bodies to be written by an expert in the theory of elasticity, but every time he asked for the collaboration of an English or a German scientist, he was given the same answer: “You have there, close to you, the authority par excellence, M. de Saint-Venant, consult him, listen to him, follow him.” One of them, M. Ettingshausen, added: “Your Academy of Sciences makes a mistake, a great mistake when it does not open its doors to a mathematician who is so highly placed in the opinion of the most competent judges.” In conclusion Moigno observes: “Fatally belittled in France of which he is the purest mathematical glory, M. de Saint-Venant enjoys a reputation in foreign countries which we dare to call grandiose.”

The French finally broke down and admitted Saint-Venant into the Academy of Sciences in 1868.  He continued his work, much of it from his home, up until the time of his death.  When the President of the Academy announced that passing, he made the following statement:

Old age was kind to our great colleague.  He died, advanced in years, without infirmities, occupied up to the last hour with problems which were dear to him and supported in the great passage by the hopes which had supported Pascal and Newton.

Europeans of the time would not have missed the import of the last statement: Pascal and Newton were Christians, and Saint-Venant was being identified with them as one also.  It was also a statement that Saint-Venant, for all of his achievements and interests which have enriched the world, also had an eternal goal as well.

There’s no evidence that Saint-Venant was ostentatious in his faith walk; descriptions of his life show the contrary.  And–shock to today’s atheist–there’s no evidence that it ever impeded the progress of his research or his thought.  As the statistician and eugenicist Karl Pearson, no friend of Christianity, noted:

The more I studied Saint-Venant’s work, the more new directions it seemed to me to open up for original investigation of the most valuable kind. It suggested innumerable unsolved problems in atomic physics, in impact, in plasticity and in a variety of other branches of elasticity, which do not seem beyond solution, and the solution of which if obtained would be of extreme importance. I felt convinced that a study of Saint-Venant’s researches would be a most valuable directive to the several young scientists, whose recent memoirs shew their interest in elasticity as well as their mathematical capacity. Many of the problems raised by Saint-Venant’s suggestive memoirs were quite beyond my powers of analysis, and I recognised that the most useful task I could undertake, was by a careful account of the memoirs themselves to lead the more competent on to their solution.


The biggest impediment he had to face was the blowback from his stand at the École Polytechnique, and that came from his secularist colleagues.  But, when the end came, all of his colleagues knew where he stood, in this life and the next one.

I spend a lot of time on this site and others talking about sea (and sometimes air) voyages.  And I’ve spent most of my career (and all the academic part of same) in the applied sciences.  But when I take my last voyage into eternity, I want to do it in the same hope of Newton and Pascal–and Saint-Venant and Euler for that matter–namely that which comes from following Jesus Christ out of the grave and into eternal life.

And you should too.

Note: I have a more complete account of Saint Venant’s story with a list of references here.


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 et.al. (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 (delta.lt.1.0e-5) 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.