Most proofs of the characteristic polynomial of the companion matrix–an important specific case–proceed by induction, and start with a 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) 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

or

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

Note carefully the inclusion of 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 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

Abstract

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.

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,

(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:

(2)

(3)

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

(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]):

(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:^

(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 :

(7)

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

(8)

If we define

(9)

and

(10)

Equation 8 can be rewritten as

(11)

The quantities 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

(12)

and the second

(13)

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

  (14)

and the second is

  (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

(16)

and

(17)

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

(18)

where

(19)

This is referred to as a similarity or collinearity transformation. Similar matrices are matrices with the same eigenvalues. The matrix , although simpler in form than , 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 , 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

(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,

(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.

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   .By direct substitution (and carrying the results to precision unjustified in geotechnical engineering) we obtain the following:

• (Equations 12 and 20.)
• (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.
• (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.
• (Equations 2, 3 and 6.)
• In both cases the trace of A and \Lambda is the same, namely (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,

(22)

where 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:

(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

(24)

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

(25)

where

(26)

(27)

(28)

It is easy to show the following:

• Equation 26 reduces to Equation 9 if .
• Equation 27 reduces to Equation 10 if additionally .
• For all three conditions, Equation 25 reduces to 11.

From the previous derivation, we can thus expect

(29)
(30)
(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:

(32)

(33)

(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:

(35)

(36)

(37)

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

(38)

and

(39)

The characteristic equation of the matrix in Equation 39 is

(40)
where

(41)

(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 ( ) 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

(43)

(44)

(45)

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

  (46)

(47)

(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 format, but it can be expanded to any size square matrix. The idea is, through a series of similarity transformations, to transform a matrix (such as is shown in Equation 23) to the Frobenius normal form, which is

(49)

From here we take the determinant of , thus

(50)

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

(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 by carrying the third row into the second, dividing all the elements by the element in the last row by the element , then subtracting the second column, multiplied by 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 by the matrix

(52)

The result will not be similar to , but if we premultiply that result by ,

(53)

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

(54)

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

(55)

(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

(57)

From this Equation 49 is easily formed.
The coefficients and 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

(58)

Premultiplying it by (see Equation 50) yields

(59)

For y to be non-trivial,

(60)

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

(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

(62)

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

(63)

The second transformation, following Equation 57,

(64)

From this the characteristic equation is, from Equation 51,

(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

(66)

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

(67)

(68)

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

(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.

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
sigmax=120.
sigmay=55.
sigmaz=-85.
tauxy=-55.
tauyz=33.
tauxz=-75.
! 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)
CLOSE (7)
STOP
END

! Subroutine danil -- Danilevsky's Method specifically written
! for 3 x 3 matrices and stress-state solution
SUBROUTINE danil (sigmax, sigmay, sigmaz, tauxy, tauyz, tauxz,
&ityp)
REAL p(3),x(3,3),sigma(3),i1,lambda
! Convert Normal Stresses to Deviatoric Stresses if ityp = 1
IF (ityp.eq.1) THEN
i1=0.3333333333*(sigmax+sigmay+sigmaz)
sigmax=sigmax-i1
sigmay=sigmay-i1
sigmaz=sigmaz-i1
END IF
! Compute Invariants
IF (ityp.eq.0) p(1)=sigmax+sigmay+sigmaz
IF (ityp.eq.1) p(1)=0.0
p(2)=-sigmax*sigmay-sigmax*sigmaz+tauxy**2+tauxz**2-sigmay*sigmaz+
&tauyz**2
p(3)=sigmax*sigmay*sigmaz-sigmax*tauyz**2+2*tauxy*tauxz*tauyz-
&tauxy**2*sigmaz-tauxz**2*sigmay
! 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
fun=lambda**2*p(1)-lambda**3+p(2)*lambda+p(3)
der=2*lambda*p(1)-3*lambda**2+p(2)
lambda=lambda-fun/der
delta=abs(fun/lambda)
IF (delta.lt.1.0e-5) GO TO 60
50 CONTINUE
60 sigma(j)=lambda
! Determine Eigenvectors
DO 80 i=1,3,1
lambda=sigma(i)
x(1,i)=-1/(-tauxz*sigmax*tauyz+tauxy*tauxz**2-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)
x(2,i)=-tauxz/tauyz*(-1/(-tauxz*sigmax*tauyz+tauxy*tauxz**2-
& 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
x(3,i)=1.0
! Normalise Eigenvectors
xnorm=1.0/sqrt(x(1,i)**2+x(2,i)**2+x(3,i)**2)
DO 70 j=1,3,1
70 x(j,i)=x(j,i)*xnorm
80 CONTINUE
! 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)
RETURN
END

References

[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.

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 4 matrix.

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

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

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 be normalized, or

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

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

Let us make the last “” or

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

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 vector, invert the former using the Gauss-Jordan technique and then multiply this by the vector before performing the Newton iteration.

Now let us consider the case of the tridiagonal matrix

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.

There are two different quantities tracked here:

• The residual, i.e., the Euclidean norm of the entire 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 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 () 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 1010.

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.

FORTRAN Code

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 1010 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),
 &xt(nplus1),resid(mcycle),eigerr(mcycle)
 call tstamp
 c Define original matrix a
 do 1 i=1,nn,1
 do 1 j=1,nn,1
 a(i,j)=0.0
 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
 vstart=float(mm)
 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
 eigold=x(nplus1)
 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
 fj(nplus1,nplus1)=0.0
 goto 6
 endif
 if(i.eq.j) then
 fj(i,i)=a(i,i)-x(nplus1)
 goto 6
 endif
 if(j.eq.nplus1) then
 fj(i,nplus1)=-x(i)
 goto 6
 endif
 if(i.eq.nplus1) then
 fj(nplus1,j)=2.0*x(j)
 goto 6
 endif
 fj(i,j)=a(i,j)
 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
 if(i-nplus1)21,22,22
 21 f(i)=-x(nplus1)*x(i)
 do 23 j=1,nn,1
 f(i)=f(i)+a(i,j)*x(j)
 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
 x(k)=x(k)-xt(k)
 7 write(2,*)x(k)
 c Compute norm of residual and error
 eigerr(m)=abs(x(nplus1)-eigold)
 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'
 write(2,101)mm,mm
 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
 close(unit=2)
 stop
 end

My first semester at Texas A&M University, I was required to take Analytical Geometry. My teacher was a former seminary student (come to think of it, so was my Calculus teacher!) Generally it wasn’t a difficult course but it had its moments.

One of them came with a particularly difficult problem we had for homework. I really didn’t know how to solve it, so I bluffed my way through it the best I could. At the top of the page I placed the following:

Evidently he was quite impressed with this show of Greek, so, in handing the papers back, he wrote my heading out on the board, pointing out that it appeared on my paper, and informing the class that he had in fact shown no mercy. One of the students, obviously unaware that any Aggie would know Greek, asked, “How did he manage to write that?”

“It was said by a very famous man,” the teacher replied.

That “very famous man,” of course, is Jesus Christ, and the passage in English is “Blessed are the merciful, for they shall obtain mercy.” (Matthew 5:7.) It is one of the Beatitudes, which open His Sermon on the Mount. The theme of being merciful and of forgiveness runs through the Gospel; in the next chapter, after the Lord’s Prayer He says again, “For and if ye shall forgive other men their trespasses, your father in heaven shall also forgive you. But and ye will not forgive men their trespasses, no more shall, your father forgive your trespasses.” (Matthew 6:14-15, Tyndale) Forgiveness is not just a nice thing to do; for the Christian, forgiveness is mandatory for eternal life.

And forgiveness is in short supply these days.  Today we live in a bitter, divided society with a record incarceration rate and creeping euthanasia. Schools call the police for acts that, a generation or two ago, would occasion a call to the parents. Makes one think of Thomas Hobbes’ characterisation of life as “brutish and short.”

God’s standard for forgiveness–from Him and from us–has not changed. If we do not forgive, we are tormented in this life and the life to come. There are certainly earthly consequences for the things that people do, but these should not be confused with our response of forgiveness. And what others do should never obscure the need for us to seek forgiveness of our own sins from God. Our “zero-tolerance” society teaches us that no one is free from mistakes. But our God sent His Son to eradicate those mistakes and make us a way to eternal life.

For more information click here.

The object of this is to solve the differential equation

for the following boundary conditions and parameters:

• 
• 
• 
• 

Conventional wisdom would indicate that, because of the high order of the derivatives, this problem cannot be solved using a scalar implementation of simple shooting. However, this is not the case.

To see why this is so, let us begin by implementing the following notation:

If we substitute these into the differential equation and solve for the third derivative, we have

Now we can construct a series of first-order differential equations for our Runge-Kutta integration scheme as follows:

Now let us consider things at our first boundary point, namely

Making the appropriate substitutions yields

We thus see that, if we use as our independent variable for shooting purposes, we also assume a value of as well. Put another way, since we are given the dependent variable and its first derivative of the ODE at the first boundary, if we guess the second derivative the third derivative automatically follows in spite of the non-linearity of the problem. So we can use a scalar shooting scheme for this problem as well. The dependent variable for root-finding purposes is

which goes to zero as the far boundary condition is met.

To implement this we used a same simple shooting routine with a regula falsi root-finding technique. One thing that was varied was the number of integration steps in the interval; we wanted to see how this affected the convergence.

The results of this are summarized below.

The number of shooting iterations is fairly constant with the variation in integration steps; however, the final value of seems to be refining itself until around 100 integration steps, at which point the accumulation of numerical error begins to affect the precision of the result.

A plot of the results is shown below.

FORTRAN 77 code is below.  Some of the “includes”, subroutines and functions are as follows:

• matsize.for is a brief snippet of code which includes a parameter statement which sizes the variable size arrays.
• tstamp calls an OPEN WATCOM routine that returns the date and time and places it in the output.
• scrsho is a “line plotter” routine, very old school.

c Solution of Third-Order Ordinary Differential Equation
 c Using Simple Shooting Method and Regula Falsi root-finding technique
 c Solution based on Carnahan, Luther and Wilkes (1969)
 include 'matsize.for'
 parameter(neqs=3,nx=10)
 character*40 vinput(nmax)
 character*16 voutpt(nmax)
 c Plot title
 character*60 titlep
 dimension y(neqs),dy(neqs),xyplot(0:nx,2)
 c Define third derivative function
 y4(y1,y2,y3,xmu)=xmu*y3-xmu*y3*y1**2-2.*xmu*y1*y2**2-y2
 c Define second boundary condition function
 bfin(y1,y1fin)=y1-y1fin
 c Enter Boundary Conditions
 data y1strt/0.0/,y2strt/0.5/,y1fin/1.0/xstrt/0.0/,xend/2.0/
 c Enter initial upper and lower bounds for y(3) (Second Derivative)
 data y3left/2.0/,y3rite/0.0/
 c Enter number of shots
 data n/50/
 c Enter friction coefficient
 data xmu/0.5/
 c Enter nomenclature describing initial data as data statements
 data vinput(1)/'Left Initial Bracket'/
 data vinput(2)/'Right Initial Bracket'/
 data vinput(4)/'Delta x increment'/
 data vinput(5)/'First Boundary Condition'/
 data vinput(6)/'Second Boundary Condition'/
 data vinput(7)/'mu'/
 data vinput(8)/'Number of Shots'/
 c Enter nomenclature describing output for individual case
 data vinput(11)/'Shot Number'/
 data vinput(12)/'Value of Guess for Second Derivative'/
 data voutpt(1)/'x'/
 data voutpt(2)/'y(x)'/
 data voutpt(3)/'y`(x)'/
 data voutpt(4)/'y``(x)'/
 data voutpt(5)/'y```(x)'/
 data titlep/'Plot of Final Boundary Condition vs. x'/
 c Description of variables below the comment:
 write(*,*)'Math 5610 Spring 2012 Final Exam'
 write(*,*)'Simple Shooting Method for Third-Order ODE'
 call tstamp
 c Compute value of integration step size for dx
 dx = (xend-xstrt)/float(nx)
 write(*,*)
 c Interval of uncertainty (left)
 1 write(*,*)'Initial Parameters for Problem:'
 write(*,200)vinput(1),y3left
 c Interval of uncertainty (right)
 write(*,200)vinput(2),y3rite
 write(*,200)vinput(4),dx
 write(*,200)vinput(5),xstrt
 write(*,200)vinput(6),xend
 write(*,200)vinput(7),xmu
 write(*,200)vinput(8),n
 do 21 iter=1,n,1
 c ..... Set and print initial conditions
 c ..... Since regula falsi isn't "self starting" it is necessary
 c to compute the left and right values for y(3) (regula
 c falsi independent variable) and compute the corresponding
 c result for bfnact (regula falsi dependent variable)
 c
 c This is done in the first two iterations
 c Iteration 1: Right Value
 c Iteration 2: Left Value
 c Susequent iterations adjust boundaries according to the
 c method
 x = xstrt
 y(1)=y1strt
 y(2)=y2strt
 if(iter.eq.1)then
 y(3) = y3rite
 write(*,*)'Iteration for Right Estimate'
 elseif(iter.eq.2)then
 y(3) = y3left
 write(*,*)'Iteration for Left Estimate'
 else
 y(3) = (y3left*bfinr-y3rite*bfinl)/(bfinr-bfinl)
 endif
 bfnold =bfnact
 y3zero = y(3)
 c Set print index to zero; print index changed at the first, middle
 c and last iteration and printing/plotting takes place
 write(*,*)
 write(*,210)iter
 write(*,200)vinput(1),y3left
 write(*,200)vinput(12),y3zero
 write(*,200)vinput(2),y3rite
 write(*,*)
 write(*,212)(voutpt(np),np=1,5,1)
 write(*,*)
 bfnact=bfin(y(1),y1fin)
 write(*,202)x,y(1),y(2),y(3),bfnact
 c ..... Set Selected variable for plot routine
 xyplot(0,1)=x
 xyplot(0,2)=bfnact
 do 17 icycle=1,nx,1
 c ..... Call Runge-Kutta Subroutine .....
 8 if(irunge(neqs,y,dy,x,dx).ne.1)goto 10
 dy(1)=y(2)
 dy(2)=y(3)
 dy(3)=y4(y(1),y(2),y(3),xmu)
 goto 8
 c ..... Print Solutions, Plot bfnact vs. x values .....
 10 bfnact=bfin(y(1),y1fin)
 write(*,202)x,y(1),y(2),y(3),bfnact
 c ..... Set Selected variable for plot routine
 nplot = icycle
 xyplot(icycle,1)=x
 xyplot(icycle,2)=bfnact
 if(icycle.eq.nx)call scrsho(xyplot,nx,titlep)
 17 continue
 c ..... Ending routine if difference between current and past
 c values of the boundary condition function are small enough
 if(iter.eq.1)then
 bfinr=bfnact
 elseif(iter.eq.2)then
 bfinl=bfnact
 else
 err = abs((bfnold-bfnact)/xmu)
 if(err.lt.1.0e-06)goto 23
 if(bfnact*bfinl.gt.0)then
 y3left = y3zero
 bfinl = bfnact
 else
 y3rite = y3zero
 bfinr = bfnact
 endif
 endif
 21 continue
 23 continue
 200 format(A40,2h ,g15.7)
 202 format(1x,f7.4,2f16.7,2f16.8)
 210 format(35hResults for Regula Falsi Iteration ,i2)
 212 format(1x,5a16)
 stop
 end

function irunge(n,y,f,x,h)
 c Fourth-order Runge-Kutta integration routine
 c Taken from Carnahan, Luther and Wilkes (1969)
 c
 c The function irunge employs the fourth-order Runge-Kutta method
 c with Kutta's coefficients to integrate a system of n simultaneous
 c first-order ordinary differential equations f(j)=dy(j)/dx,
 c (j=1,2,...,n), across one step of length h in the independent
 c variable x, subject to initial conditions y(j), (j=1,2,...,n).
 c Each f(j), the derivative of y(j), must be computed four times
 c per integration step by the calling program. The function must
 c be called five times per step (pass(1)...pass(5)) so that the
 c independent variable value (x) and the solution values (y(1)...y(n))
 c can be updated using the Runge-Kutta algorithm. M is the pass
 c counter. Irunge returns as its value 1 to signal that all derivatives
 c (the f(j)) be evaluated or 0 to signal that the integration process
 c for the current step is finished. Savey(j) is used to save the
 c initial value of y(j) and phi(j) is the increment function for the
 c j(th) equation.
 c
 c The size of the savey and phi arrays depends upon the value of the parameter
 c which is included.
 include 'matsize.for'
 dimension phi(nmax),savey(nmax),y(n),f(n)
 data m/0/
 m=m+1
 goto(1,2,3,4,5), m
 c ..... Pass 1 .....
 1 irunge=1
 return
 c ..... Pass 2 .....
 2 do 22 j=1,n,1
 savey(j)=y(j)
 phi(j)=f(j)
 22 y(j)=savey(j)+h*f(j)/2.0
 x=x+h/2.0
 irunge=1
 return
 c ..... Pass 3 .....
 3 do 33 j=1,n,1
 phi(j)=phi(j)+2.0*f(j)
 33 y(j) = savey(j)+0.5*h*f(j)
 irunge=1
 return
 c ..... Pass 4 .....
 4 do 44 j=1,n,1
 phi(j)=phi(j)+2.0*f(j)
 44 y(j)=savey(j)+h*f(j)
 x=x+0.5*h
 irunge=1
 return
 c ..... Pass 5 .....
 5 do 55 j=1,n,1
 55 y(j)=savey(j)+(phi(j)+f(j))*h/6.0
 m=0
 irunge=0
 return
 end