Lecture 12: Solution Methods for Frequencies and Mode Shapes

Related Resources

Study Guide (PDF)

Readings

Sections 11.1-11.6

Examples

Problems 11.1-11.27

• Download English-US transcript (PDF)

The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high quality educational resources for free. To make a donation, or view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.

PROFESSOR: Ladies and gentlemen, welcome to lecture number 12. In this lecture, I would like to discuss with you solution methods for eigenproblems. This is a very large field, and all I can do really in this lecture is to give you or discuss this you a few of the solution methods that we are using in finite element analysis. So it will be a very brief introduction, a brief survey of the methods of solution that we're using for finite element eigenvalue problems. The generalized eigenvalue problem, K phi equals lambda M phi, is a problem that we encountered in the multiple position analysis. Please recognize that I'm calling lambda here, the eigenvalue, which was our omega squared in the last lecture.

This is really the eigenproblems that we are primarily addressing ourselves to finite element analysis. The standard eigenproblems problem here, EVP stands for eigenproblem, has M equal to the identity matrix. So when we're talking about solution methods of this problem, we're really also talk about solution methods of this problem. If we have, in a [? multiple ?] position analysis, non-proportional damping, then we would have to solve this eigenproblem, which is a quadratic eigenproblem. I will not address myself to the solution of this problem, so in what follow, we are talking about a solution of this problem, and of course, the solution of this problem also, when M is equal to the identity matrix.

In particular, we are interested in the solution of large eigenvalue problems when N, is say, greater than 500. N, of course, being the order of K and M, and with the half bandwidth, I should say, being greater than 60.

The number of eigenvalues that we will be looking for, that we want to calculate, are say, 1 to 1/3 of N. 1/3 of N is already large number of eigenvalues. Generally, we are talking about 1 to and at most, 100 for large eigenproblems. When we have a small eigenvalue problem, in other words, N being small, say 50 or 100, there are a large number of different techniques that can be used to solve this eigenvalue problem, and it really does not make that much difference which technique you are using. The cost can still be quite significant, but certainly, when we're talking about very large eigenproblems, N being even larger than a 1,000 or 2,000, then it is of up most concern to select a technique that is very effective.

And I want to address myself really now towards those techniques that we feel are very effectively used in finite element analysis. In dynamic analysis, of course, with proportional damping, this was the eigenvalue problem that we wanted to solve. Now, if they are 0 frequencies present, in other words, when we have rigid body modes in the system, then we can use the following procedure to basically getting rid of these 0 frequencies. We simply add, on the left hand side, a mu M phi, and on the right hand side, a mu M phi. And taking these two together, we are getting this coefficient matrix, and a new lambda value. That lambda value now is omega squared plus mu, where omega squared is equal to lambda minus mu.

What we have been doing here, is basically, we have been applying a shift. In other words, if we were to plot the characteristic polynomial, giving us the eigenvalues, the characteristic polynomial of course, being p of lambda being determined K bar minus lambda M, where K bar is equal to K plus mu M.

If we plug this characteristic polynomial, we will see this curve here, schematically drawn, of course, where the omega squared values would start here, in other words, the lowest omega squared value is, in fact, 0. Whereas our lambda value now starts right here. This is the shift that we [? imposed. ?] So, basically, we now recognize that we have all eigenvalues being greater than 0, because this is our new starting value. All eigenvalues measured from here are greater than 0, the smallest one being mu, namely when omega squared is equal to 0. And the larger ones being here.

It is effective to perform such a shift when we have 0 eigenvalues, in other words, to operate on this eigenvalue, rather than on that one, because our techniques then are very stable and effective, and this shift here is, in general, therefore performed. So all we have to address ourselves to then from now on, is really, the solution of eigenvalues and eigenvectors, when all of the eigenvalues are greater than 0. If we start off with a problem where we have a 0 frequency, we apply the shift so that the eigenvalues now of this problem here are all greater than 0.

In buckling analysis, we have this eigenproblem. We really did not talk about-- we did not really at all about buckling analysis in the set of lectures, but you might just be interested in finding out how would we solve for these eigenvalues. Well, here is the characteristic polynomial. The smallest lambda value, of course, gives now a critical load, that is the one e would be interested in. It is effective here, to rewrite this problem in this form. All we do, is we are taking the lambda value onto the other side. In other words, we are dividing this equation here by a lambda. So that one of the lambda stands on this side, and KG phi on that side. The resulting problem is this one, where kappa is 1 over lambda and we have now here, is a K matrix.

K is positive definite, and now it is more effective to operate on this equation. How do we obtain now the largest kappa? Because the largest kappa is the one that we are interested in, which corresponds to the smallest lambda, the smallest buckling load. Well, we would now perform a shift to the right, subtracting mu K from this side, from this side, and from that side, to obtain this eigenvalue problem. That means we are really shifting from here up to there, and now we are interested in solving kappa N by operating on this eigenproblems problem. It is this eigenvalue problem that we can operate on directly, using the techniques that I would be discussing later, to determine search method and the subspace iteration method.

But before we get into those techniques, I would like to discuss with you some preliminary considerations that are important when we look at an eigenvalue solution. The important first point that I want to make is that there is a traditional approach for the solution of this eigenvalue problem here. Let's look at this eigenvalue problem. Where of course, when we're talking about a buckling analysis, we would have different matrices here, but basically, we still solve this eigenvalue problem.

The traditional approach lies in factorizing the M matrix into Cholesky factors. We talked about the Cholesky factor when we talked about solution of equations. So here it comes up again. We would factorize the M matrix into Cholesky factors, as shown here. We define now a new eigenvector, phi curl being equal as shown here, L curl transpose times 5. We substitute from here into there, and we directly obtain this eigenvalue problem here. Now we have the standard eigenproblem.

Now once we have a standard eigenproblem, there are a large number of techniques that can be used. For example, the QR method, the standard Jacobi method, et cetera. So this is, what we might call. a traditional approach to solve the generalized eigenproblem, this is the one here that we are really interested in. But it involves this transformation, and this transformation here can involve a very significant cost, because we have to find these Cholesky factors and then, of course, we have to calculate this K curl. Another way of performing a transformation of the generalized item problem to a standard form would be to solve the eigenproblem of M first. W here and W transpose [? store ?] the eigenvectors of the M matrix, and D squared is a diagonal matrix of the eigenvalues of M. And then having factored M in this form, an equation similar to this one here can be written. Where now, the transformation, of course, involves the eigenvectors, W.

Well, this approach is sometimes still quite effective, but in general, it cannot be recommended because there's too much cost involved in the factorizing of the M matrix, and then particular, in the transformation here, and the solution of this problem and the, of course, to get from this problem, the phi vector, we have to go back to this equation. A direct solution of the generalized eigenvalue problem that we want to solve is really more effective. And now, we are ordering the lambdas as shown here. Remember, our smallest lambda now is greater than 0, the eigonvectors are stored as shown here, and for most of a position analysis, we are really interested in solving only for a certain number of eigen pairs, i equals 1 to p, as discussed in the last lecture, or i equals r to s. If we have a vibration excitation problem, we really only want to solve for specific frequencies and corresponding mode shapes in a certain interval, as shown here.

The solution procedures on this generalized eigenvalue problem can be categorized into the following forms. The first solution procedure, or a whole class of solution procedures, that I should really say, which we might call the vector iteration techniques. The idea is the following-- if we want to solve this problem, then why not solve it in the following way? Start off with a vector on this side, that we assume, calculate a new vector on this side, then take that new vector, plug it back on the right hand side, and cycle through this way.

This is indeed inverse iteration. What do we do here? We are putting a vector on the right hand side, x1, when K is equal to 1, that vector we have to assume. Then we are multiplying out the right hand side. Now we have a set of simultaneous equations that we would solve just in the same way as we solve the equations in static analysis. We can think of this as a load vector. We are calculating this displacement vector, and then having got this displacement vector, and not taking care of this lambda, or having, basically, taken this lambda out of this equation. Which means that this vector will grow in lengths, or collapse in lengths. Therefore, want to scale this vector down to a good length-- a length that does not hit overflow, so to say, or underflow in the machine. And that scaling is achieved as shown here.

So, the final result then, this is just a scalar here, as you can see, is a vector here, which we can put on to the right hand side and keep on iterating around. In fact, we find as that iteration procedure converges to the lowest eigenvector, provided the starting vector is not orthogonal, is not M orthogonal to phi 1. This is one first iteration technique. Maybe I should put a side note down. A side note that many of you might have in your mind. Many of you might have the following question in your mind-- why he iterate, why does he not solve directly for the eigenvalues? Just like we saw a set of simultaneous equations. Well, there's one very important point, and this is the point-- if we are solving an eigenvalue problem, and shown here, then really that means the following-- we are saying that determinant, K minus lambda M shall be equal to 0.

Why? Because if we write this equation in the following form-- K minus lambda M times phi equals 0, then we have here a set of simultaneous equations with a 0 vector on the right hand side. And there is an important theorem in linear algebra that says, only a solution to this set of equations when this coefficient matrix is singular. Which means that the determinant of this coefficient matrix must be 0. So, what do we want? We really want to solve p of lambda determinant, K minus lambda M. We want to solve for the 0's of this polynomial. That's really what we want.

Now, if we want to solve for the 0's of this polynomial, then we also have to understand that when k is of order 4, there will be four 0's. The polynomial itself will be of order 4. If K is of order, N, then there will be N 0's. So the polynomial is order N. However, if we want to find the solution of a polynomial, or order N, then we very quickly find that, in fact, we can only look up the solution to such polynomials if they are of low enough order, lower than say, 4.

Above for larger order polynomials, there's no way you can look up the solution. You have to iterate to obtain the solution of the 0's of this polynomial. However, since this problem here, is completely equivalent to this problem. To solve this problem, you have to iterate. And this is the important conclusion. It is of up most importance that you understand that to solve an eigenvalue problem of order larger than say, 4, you will have to iterate. Because you're solving for the 0's of a polynomial of order larger than 4. Only for very small order polynomials, we can look up the solutions in handbooks.

So having grasped the fact that we have to iterate to solve for the 0's of this polynomial, to solve for the eigenvalues of this problem, we are stuck with iteration, and since we're stuck with iteration, we now want to come up with techniques to iterate. And the first technique to iterate is this vector iteration technique.

And the basic idea, once more, is if you want to solve this equation, well, then why not iterate in the following way? Put some starting value in here, this number is just a scalar, let's neglect that effect for the moment. And having put a starting value in here, calculate this one and keep on cycling around this way. And that's what we call inverse iteration. Inverse iteration here, a vector inverse iteration, because we have to solve a set of equations here, with K as the coefficient matrix.

There are other iteration techniques. There's also forward iteration. And basically, what we're saying there is-- well if we are starting off here and going around in this circle, why not start off here and go around in that circle, the red line. And that, in fact, gives us forward iteration.

One can prove that that forward iteration converges to the largest eigenvalue. We are usually interested in the smallest frequencies, the smallest eigenvalues, and therefore, we are using inverse iteration quite commonly.

Then there is Rayleigh Quotient iteration. The Rayleigh Quotient iteration is again, an inverse iteration, however, with an acceleration scheme in it, a very effective acceleration scheme. This Rayleigh Quotient iteration can be employed to calculate one eigenvalue and vector, and then, of course, we have to deflate afterwards the matrix, the coefficient matrix for the eigenvector that we have calculated already, and we then could go on towards the calculation of additional eigenvalues and vectors. The same deflation also would be used in inverse iteration alone, and in forward iteration alone. In a Rayleigh Quotient iteration, we converge to an eigen pair, but which one is not guaranteed.

So this is a class of techniques, iteration methods that we are using.

Another class of techniques, again iteration methods, we are stuck with iteration, as I mentioned earlier, is the polynomial iteration methods. Here, we are saying that this problem is equivalent to that, we want to solve for these 0's. Therefore, for the 0's of this polynomial, sketched out here, and one way is to use Newton iteration. Here we have p, here we have p prime. Well, of course, if we have p here, we have to also recognize we should have the coefficients here. Well, these coefficients are given here, A0 to a n. And we could rewrite this polynomial in this form, once we know the zeroes. Now, if we don't know the zeroes, of course, we would be stuck on iterating on this matrix, where we would have to calculate these coefficients.

However, there are now really, two different techniques. The first technique is explicit polynomial iteration. In the explicit polynomial iteration, we expand the polynomial and iterate for the 0's, by that I mean, we actually expand this polynomial. However, this technique is not suitable for larger problems because first of all, there's much work to obtain the ai's, and then it's an unstable process. Unstable because small errors in these coefficients here, give large errors in the zeroes.

So this is not a suitable technique, and we are better off using implicit polynomial iteration. In the implicit polynomial iteration, we evaluate the determinant all this matrix here. In others, p at mu i, via the L D L transpose factorization, which we discussed in the Gauss elimination solution, when we talked about Gauss elimination solution. To obtain the determinant of this matrix here then, we need simply to multiply all the diagonal elements. This technique is accurate, provided we do not encounter large multipliers L matrix. And if we do, we have to shift away, we have to change our mu i. Also we directly solve for lambda 1 and so on. And for the more we can use effectively a secant iteration, instead of calculating, in other words, the p prime at mu. This is here approximately equal to p prime, at mu i.

Instead of calculating this derivative, we are putting this value here, which gives us an approximation to that derivative, this is what we call the secant iteration. Once we have calculated one eigenvalue, of course, we will have to deflate to calculate the next eigenvalue, and this is the process that he would pursue. Pictorially, we have here, our p lambda, this is the polynomial. Let me show it to you here.

Here, we have the polynomial, here we have lambda 1, here we have lambda 2. In the secant iteration, we start off with the value of mu i minus 1 here, mu i, these two values are given. Usually, in a practical analysis, mu 0 would be this value here, and mu 1 would be a lower bound on lambda 1. And plugging these two values in, we get mu i plus 1, and you can see that this process would then converge to this eigenvalue. Once we have lambda 1, we can repeat the iteration. However, working now on p lambda divided by lambda minus lambda 1. And that means we are operating on this polynomial here, and polynomial then gives us the iteration on that polynomial then gives us the second eigenvalue. So we are deflating the actual p lambda by the lambda 1 value in order to obtain a convergence to the second value.

Well, this is the one class of methods. Another class of methods, the third one, are the Sturm sequence methods. And here, the basic idea is the following-- if we look at this eigenproblem, and I use a 4x4 matrix here, then if we take a shift on this eigenproblem, mu S, and this is the shift here. Then if we factorize K minus mu S M, equal to L D L transpose, again, we're using our Gauss elimination procedures basically. Then the number of negative elements in D is equal to the number eigenvalues smaller than mu S, a very important point.

Let us look at this example here of the 4x4 matrix. Here we have plotted characteristic polynomials. The characteristic polynomial of this 4x4 matrix, with that M matrix, of course, is drawn along here, four eigenvalues. And I have put mu S here as the red line. What this theorem says is if I factorize K minus mu S, M, this being the K matrix, the total K matrix, and the 4x4 M matrix goes in here. Then there would be one negative element in D. One only because mu S is large than lambda 1. If mu S were on this side here, there would be two negative elements in D. If mu S were on this side, there would be three negative elements in D-- very important point because qw can use this procedure to directly calculate eigenvalues the way I will be discussing it just now.

However, let us look a little bit more at why this hole is here. Well, it derives from the associated constraint problems with the original eigenproblem. The original eigenproblem corresponds to this beam problem, which we considered earlier already in the discussion of Gauss elimination. If I constrain this degree of freedom here, I get a new polynomial corresponding to this problem, and the important point is that the lowest root, the lowest value of this problem here lies in between these two-- lambda 1 and lambda 2.

The second one here, lies between lambda 2 and lambda 3, and the third one here lies between lambda 3 and lambda 4. The second associated constraint problem has two eigenvalues, and this one lies between these two, this one lies between these two. And the final eigenvalue problem here, the third associate constraint problem obtained by knocking out these degrees of freedom here, has one eigenvalue, and that one lies in between these two values here. It's this important fact that is used to derive this final result. Now, we have not time to go through the actual details, but this is the important fact here. In other words, that the [? i's ?] eigenvalue of this associated constraint problem, for example, the second eigenvalue of this associated constraint problem here, bisects or lies in between, not bisect-- lies in between lambda 2 and lambda 3 of this eigenproblem here, and so on.

That result is used to derive this final result, which we use effectively then in the eigen solution. And the class of techniques then that we are referring to using this fact are the Sturm sequence methods. So basically, what we are saying is the following-- take this coefficient matrix, now I make the mu x available. We factorize that coefficient matrix into L D L transposed, and the number of negative elements in D is equal to the number of eigenvalues smaller than mu S, I. How do we use this fact?

Well, if this is here, the polynomial, we could start off by putting mu S 1 into this point here, and we would find out the number of eigenvalues smaller than mu S 1. Then we are taking, mu S 2, and we find out the number of eigenvalues smaller than mu S 2. Of course, we would find that they are one more eigenvalue smaller than mu S 2. Then mu S 1, namely this one. So we know there's one eigenvalue in between mu S 1 and mu S 2. Now we could bisect this interval, and we would find that, in fact, there is still one more eigenvalue smaller than mu S 3, then there was smaller than mu S 1. Because this eigenvalue still lies on the left side of mu S 3.

So we bisect again, we get mu S 4. Notice this distance here, is equal to that distance here. And now, we would find that the eigenvalue still lies in between these two lines, between mu S 3 and mu S 4. And like that. we could continue bisecting the interval until we finally find a very small interval in which this eigenvalue must lie.

We, of course, have to take care in the L D L transpose factorization, because we are not working on a positive definite matrix. The multipliers must be small, the L, I, J elements must remain small in the factorization. But provided that is the case, the round of errors are small, and the procedure is very stable.

However, convergence can be very slow, can be extremely slow because you have to cut down the intervals further and further until you finally get a very small interval, in which you now know the eigenvalue does lie. Of course, at some point, you might want to switch to another iteration scheme. Once you know there is an eigenvalue in between say mu S 3 and mu S 4, you might want to switch to a iteration scheme that converges to this particular value much faster.

Finally, a whole class of methods that I would like to briefly mention to you are the methods of transformations, where we are operating on this eigenvalue problem, recognizing that phi transpose K phi is equal to lambda. And phi transpose M phi is equal to I. Where phi is defined as shown here, lambda is defined as shown here. This fact, of course, we use already earlier in the [? multiple ?] position analysis. Now how can we get phi and lambda? That is our objective. Well, if we know that this matrix diagonalizes the K matrix, and it also diagonalizes the M matrix, then why not try to construct it by iteration? And this is the basic idea in the transformation methods.

What we do is the following-- we take the original K matrix, we are operating with P1 transpose K P1. On K, and on M, such as to make this resultant matrix here closer to a diagonal form. Ultimately, we want to get a diagonal matrix. Well, one step in the iteration is to make this P1 transpose K P1 matrix, and P1 transpose M P1 matrix, both of these matrices a little closer to diagonal form. Then K and M were by themselves.

So this is the first step. Now, we have here, let's call this a K2 matrix, and this here an M2 matrix. Now, we are operating with P2 transpose on K2, and also there's a P2 here. So let me take a different color to show that now we are operating with this part here, and that part there. And this, of course, we are doing in such a way as to have that P2 transposed, K2, P2 is a little bit closer to diagonal form than was K2. I hope you can see these lines-- let me put them in once more. What I'm saying is we have K2 after having applied P1 transpose on P1, on K, and similar for M.

And now we want to make K2 still a little bit closer to diagonal form. Well, we are now applying P2 transposed on K2 and the P2 on K2, this one here has been replaced by K2, this one here has been replaced by M2. And what we're saying is that this final result, in this final result should be closer to diagonal form. In fact, this should be closer to I matrix, and this should be closer to the lambda matrix, ideally.

Well, this is the basic process that is used in the generalized Jacobi method. We can show that if we choose specific matrices, P I matrices. That zero always one [? of ?] diagonal element. In other words, the P I matrix, zeroes a particular off diagonal element in this current coefficient matrix, similarly for the M2, for the current M coefficient matrix. Then we know in the generalized Jacobi method that the process finally converges. It converges and we get a diagonal matrix here.

With proper scaling, it will be the identity matrix. And here, we are getting also a diagonal matrix, and if this is an identity matrix, this will actually be the matrix storing the eigenvalues. Well, this is the generalized Jacobi method here. They are other techniques that can also be used, of course, this is a very common approach taken when we have M being equal to the identity matrix.

However, the disadvantage of this approach is that we are calculating all eigen pairs simultaneously. In finite element analysis, we are only interested in specific eigenvalues, so lowest P, or a certain number of eigenvalues in an interval. Here, we have to calculate all of the eigenvalues stored in the lambda matrix, and all of the eigenvectors. And the eigenvectors here, would be these vectors here. These are the the matrix phi, which is P1 times P2, up to P K, after K steps of iteration. This gives us the matrix phi, storing all of the eigenvectors.

So the disadvantage is that we are calculating all of that eigenvectors, which in finite element analysis is hardly ever required. All we want is the total number of eigenvalues and eigenvectors. Remember, also if K is of [? order ?] 1,000, this would be a tremendous amount of work involved and, in fact, it would be beyond the current state of computing capabilities to use a generalized Jacobi method on a K matrix 1,000x1,000 and an M matrix 1,000x1,000. You will certainly stretch the current computing capabilities tremendously by using the generalized Jacobi method on large eigenproblems.

For large eigenproblems, it is really best to combine the basic approaches that I just mentioned. Let's recall the basic approaches. There were iteration methods, vector iteration methods, there were Sturm sequence methods, there was the transformation method. And then, of course, there were polynomial iteration methods. So really, four classes of methods that we are talking about.

Each of them, we find, have some advantages. Well, as engineer, surely we want to now combine the advantages of each of these methods to come up with optimum techniques for our specific problems. Well, there are a number of thoughts that we would go through. First of all, a determinant search, a polynomial iteration method would be effective to get near a root. Vector iteration then could be used to calculate an eigenvalue and an eigenvector. The transformation method would be good to orthogonalize the current iteration vectors if we iterate with more than just one vector. And the Sturm sequence method is extremely powerful to assure that we have actually calculated the required eigenvalues and eigenvectors.

Well, the first method then that I would like to discuss with you is the determinant search method. In this method, we have combined some of these techniques that I pointed out to you earlier in an optimum way. Here, we have the polynomial P lambda, we want to, of course, solve for the roots of that polynomial-- that is our objective. Here is lambda 1, there is lambda 2. I also show a point, A, because I will refer to it just now.

The basic scheme in the determinant search method is to look for the zeroes of the determinant, of the determinant K minus mu I M. Of course, these zeroes give us R equal to the eigenvalues. That's what we're interested in calculating. Well, we are calculating the determinant by the factorization of this matrix into L D L transpose, and we notice that this determinant is simply the product of the [? D I I's, ?] as I mentioned earlier.

Now, having a certain starting value, mu 0, say, and mu 1. We can directly use the secant iteration, that I mentioned to you earlier to get close to an eigenvalue. Now, I also slipped in here an acceleration factor. This acceleration factor would be equal to 1 in the normal secant iteration that I showed to you earlier already.

However, since we want to use this polynomial iteration technique only to get close to an eigenvalue, and we might actually obtain with eta equal to 1, sometimes slow convergence. What we do is we simply set it larger to accelerate the process. And this is what we're doing here. If we do jump beyond an eigenvalue because we have accelerated the secant iteration. In other words, if we do a jump across here, then since we are factorizing the matrix here, we can also look at the number of negative elements in the D matrix. These would tell us that we, in fact, have jumped across an unknown eigenvalue.

So, in this particular case, of course, there would be one negative element in the D matrix. And we would know now we have jumped across. The jumping across an unknown eigenvalue is not possible when eta is equal to 1. However, it is possible when eta is greater than 1. And however, that jumping does not hurt us. In fact, all we will find is that the Sturm sequence count will tell us that we have jumped across an eigenvalue, and then, of course, we have to go to another strategy of solving for the eigenvalue more accurately.

But the important point is that all we want to do via this procedure is to get into the vicinity of the eigenvalue. We do not want to calculate it accurately, all we want to do is get into the vicinity of an unknown eigenvalue. First, of course, lambda 1, then lambda 2, and so on.

There's also a concern, of course, that we do not want to jump beyond A, because if we, in our process of looking for lambda 1, we are jumping too close to lambda 2, then via the other solution techniques that we are using later on then to calculate the eigenvalue more accurately, we would actually converge to lambda 2. And, of course, that we don't want, we want to calculate lambda 1.

So, the process is then to start off with eta equal to 1, iterate, see how we are progressing towards an unknown eigenvalue. If that progressing is very slow, we are choosing eta larger. There are specific strategies that we're using to choose the right amount all eta. And until we find that we are close to an eigenvalue, either from the left or we have jumped across it already, but you should never jumped to far. And we are making sure that we never jump beyond that point A. and having jumped, we find that the Sturm sequence count tells us about it.

And now, we are going on to another strategy, and that strategy is the vector iteration, vector inverse iteration method that I discussed with you earlier. What we are doing now, is we have a new coefficient matrix, that's the one. We are sitting here, say, and now I have assumed that we have, in fact, jumped. It's not necessary. We might have also approached from this side and come very close to lambda J. And now, we are using vector iteration to get this eigenvalue accurately.

This coefficient here, or this value here corresponds to that length. And if you are adding this value here to this shift, we are getting an approximation to lambda J. And that, of course, is our objective. Once we have lambda J, we have calculated this eigenvalue-- the eigenvector comes out, in x k plus 1 as k increases, so we have calculated this eigen pair, and now, of course, we would use deflation method, just the way I showed it to your earlier, for the polynomial iteration method. To deflate the polynomial of this value, implicitly-- implicitly, that's important-- and we would continue our iteration process towards that eigenvalue next.

In this process, we also want to make sure that the iteration vector is deflated of the previously eigenvectors, for example, Gram-Schmidt orthogonalization. If the convergence is slow in this particular iteration here, we're using also the really Rayleigh-quotient iteration to speed up iteration itself. The advantage of the method-- well, it calculates only the eigen pairs that are actually required., there's no prior transformation of the eigenprobelm necessary. The disadvantage of the method is that we have to perform many triangular factorizations. Therefore, the method is really only effective for small banded systems.

If we want to calculate eigenvalues of large banded systems, we really need an algorithm with less factorizations, and more vector iterations. And that is the subspace iteration method. In fact, we have now extended the subspace iteration method to be also most effective for small banded systems. And I will just refer to it very briefly later.

The basic steps of the subspace iteration method are the following. Here, we have the first equation used. We have on the right hand side, the current iteration vectors. Now, notice that in this particular case, we're iterating with q vectors, where q ip larger than p simultaneously. This gives us here, q right hand side vector. We are solving for q vectors here. This is an inverse iteration step on q vector simultaneously. Then we are calculating the projections of these k and m matrices on these vectors.

Notice that this pat here, of course, is equal to this right hand side. So we would not calculate this here. We simply use the right hand side, which we have evaluated already and plug it right in there. But having calculated these two matrices, which are now of order q by q, we are solving eigenvalue problem all of these q by q matrices. And the solution of that eigenvalue problem is shown here. Notice this matrix here is a q by q matrix. We're using q vectors. This is a q by q matrix storing the eigenvectors of this eigen problem here, of the eigenproblem corresponding to these matrices.

This here is the diagonal matrix listing the eigenvalues of the eigenproblem corresponding to these matrices. Having calculated these matrices here, and of course, that one too, we are performing this operation to get a new set of iteration vectors, which then are plugged back into here. Notice that the solution of this eigenproblem, q by q, might be a 50x50 matrix here, 100x100 matrix, maybe a 200x200, but that already would be a large number of vectors if we're using-- a very large number of vectors. This solution of this eigenproblem can be effectively achieved via the generalized Jacobi method, which I referred to briefly earlier.

Of course here, we are using our solution algorithm for static analysis. There's L D L transpose factorization involved here. And here, this is a simple matrix multiplication that has to be carried out. This part here, is the eigenproblem solution, and there's another vector multiplication carried out here. Now, under specific conditions, this eigensolution algorithm converges and the vectors stored in x k plus 1, they are q vectors now here, with appropriate ordering, converge to the phi matrix. This is now an n by q matrix, storing the lowest eigenvectorr, or I should say, the q eigenvectors corresponding to the lowest eigenvalues. And here, we have it written out, and of course, this is a diagonal matrix. The condition that we have to satisfy here is that the starting subspace spanned by these vectors here must not be orthogonal to the subspace spanned by the p eigenvectors that we are interested in.

Of course, here, I'm talking about q eigenvectors. I really am only interested in p eigenvectors, in the lowest p eigenvectors. So in general, what happens is that I will only be interested in these eigenvectors here, and the lowest eigenvalues. One comes out here, I'm not too much concerned about. In fact, we are carrying these along only to accelerate the iteration scheme. And I will show you just know what the convergence rate is, and why we are carrying them along.

So since we only interested in the lowest p eigenvalues and corresponding eigenvectors, we want to get convergence towards those, and what that means is that our starting vectors in here must not be m orthogonal to the eigenvectors that we want to calculate, phi 1 to phi p. If that condition is satisfied, one can show, theoretically, and of course, we see it in practice also, that this iteration scheme will always converge to the lowest eigenvectors and corresponding eigenvalues.

One important point here, once we have calculated, say, a certain number of iteration vectors, and we have tested for convergence, we are saying that convergence is reached, say, when the current iterates on the eigenvalues don't change very much anymore. This tol might be 10 to the minus 6, as an example. Then once this condition is satisfied, we stop the iteration, and we really want to make then sure that we really converge to the eigenvalues of interest. And how is that accomplished? Well, if these are the p eigenvalues, we are calculating error bounds on these eigenvalues, indicated here by the blue lines here. And then just to the right of the error bound. Or we can really use that error bound, we apply a Sturm sequence shift, mu S. What this then says that k minus mu S factorize into L D L transpose. And we must now have in the D matrix, p negative elements. And only p, not more, not less, just p negative elements must be occurring in the D matrix.

This check is absolutely sure. If this check a satisfied, we can be absolutely sure that we have calculated the eigenvalues of interest, the lowest p eigenvalues, in this particular case. So I repeat, after we have satisfied this convergence rate, or this conversions limit, I should say-- typically 10 to the minus 6, or 10 to the minus 8, depending on what accuracy you want. We calculating eigenvalue bounds, as shown here. We are putting a mu S Sturm sequence shift just to the right of the last eigenvalue bound, and if this is the p's eigenvalue, we must have here, exactly p negative elements indeed.

The convergence rate of the iteration, we find is given by lambda i divided by lambda q plus 1. And the convergence rate to the eigenvalues of interest are just this value squared. Notice that we have here, q plus 1. So if we are interested in the p's eigenvalue,== the black pen doesn't work too well, let me take the blue one here-- if you're interested in the p's eigenvalue as the largest eigenvalue, notice that the convergence rate is that lambda p plus-- lambda p divided by lambda q plus 1. That is the worst convergence rate that we are reaching because the other eigenvalues below lambda p are smaller, giving us better convergence rate. So this is the worst that we can reach, if we are only interested in p eigenvalues and corresponding eigenvectors.

Notice however, that lambda q plus 1 must be greater than lambda p. If lambda q plus 1 is equal to lambda p, we would never converge. So this is the reason why we are choosing q to be larger than p, substantially large than p. For example, we might choose q to be typically, as an example, q being the maximum of 2p and p plus 8. This is a formula that we have used with quite a lot of success. In other words, when p is equal to 4, we would use-- no, sorry, I should say the minimum here. the In other words, when p is equal to 4, we would have 8 here. When p is equal to 8, it's the same. Because these two values give us the same. And when p is equal to 100, then we just have 8 more. Q is equal to 108. This 8 here is a safeguard, so that this cannot be equal to 1 in practice. It would only be equal to 1 if we would have, basically. 9 time the same root, and of course, that is a case which we would hardly encounter in practice.

Now, regarding the choice of the starting values. For the vectors, there are two choices. x1 is simply a vector with 1 only. xj would be a vector with zeroes everywhere, but with a 1 somewhere. And that one is in the j's entry-- sorry, the k's entry. The k's entry for ek, ek has in the case, entry of 1. And otherwise, everywhere zeroes. This would be the second to the q minus first vector. And the last vector, we take a random vector. That is one choice of starting vectors. The Lanczos method can also be used to generate starting vectors very effectively.

Once we have selected the starting vector, we go through the iteration schemes the way I have been displaying it. And there, it is important, once again, to perform after the convergence- the Sturm system sequence check, and also we might want to calculate error bounds on the eigenvalues and eigenvectors. These error bounds would be calculated by simply saying, well, if we want to satisfy-- or if we want to solve for lambdas and phis that satisfies this equation, then why not see how close we will satisfy. Make that a minus and put here an error, say, R, on the right hand side.

Substituting then the last values for lambda and phi that we just kind calculated, the last iterative values, we get this top line here. We're taking a norm the way I've been talking about earlier. And we are also taking a norm at the bottom to get, basically, a relative value, top to bottom. And then relative value gives us epsilon i. So epsilon i should be of order, 10 to the minus 3, if we want to have three-digit accuracy in the eigenvalues and eigenvectors. In fact, the eigenvalues will generally be more accurate than the eigenvectors. That is, of course, a very low accuracy. Ideally, we have much more accuracy after the iteration schemes. We would have to tighten the tolerance, the tol value that I referred to earlier sufficiently to make epsilon i small enough.

Now, I like to refer, very briefly also, to an accelerated subspace iteration method that we have been developing during the last years, which is a very substantial extension of this subspace iteration method. The idea here is, basically, that we wanted to combine the sum of the effective techniques that we are using in the polynomial iteration method, with the standard subspace iteration method that I have discussed with you here. Basically, in the accelerated subspace iteration method, we are not keeping a constant coefficient matrix, as we have done here. The k matrix never changed. Here, we are shifting the coefficient matrix, we are shifting through the spectrum. Much in the same way as we are pursuing it in the determined search method. And that is number one.

And number two-- we are also not iterating with vectors that are always larger than the number of eigenvectors that we want to calculated. In fact, if we want to calculate, say, 100 eigenvalues, we might just use four iteration vectors at a time, and shift through the spectrum, always with four iteration vectors, capturing the eigenvalues values and eigenvectors that we are really interested in, in bundles. So this is a way how we have very substantially accelerated the standard subspace iteration method that I have been talking to you about just now.

If you are interested in seeing some solution times, please we look at this reference. Ladies and gentlemen, this then brings me to the conclusion of this lecture. And in fact, it also brings us to the conclusion of the whole set of lectures that I have been presenting to you. These set up lectures have been brief and very compact. They've been a survey introduction to the finite element method-- a very brief and compact introduction and survey. Brief because there are only 12 lectures. Compact because I wanted to present to you as much as possible in these 12 lectures. And there is, of course, a lot of information that we could have talked about.

However, I hope that the information that I had given in the lectures has been of interest to you, and it will also be valuable to you in your work and general decision making process of whether you want to use a finite element method, how you want to use it, and so on. I like to take the opportunity at this point to thank the Center of Advanced Engineering Studies for the terrific effort they've made in preparing the tapes. I'm very happy that they have been cooperative with me in the way they did, in the preparation of the tapes. I also would like to thank you as listeners for your patience of listening to the tapes and, once again, I hope you found the tapes interesting and have obtained some value information on them. Best of luck in your further studies of the finite element method.