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 9.
In this lecture, I would like to discuss with you the
solution of equilibrium equations in static analysis.
In the earlier lectures, we derived equations KU equals R.
And we now have to solve these equations for the
displacements stored in the vector U.
This is an important phase of the finite element analysis,
because much of the computational effort goes into
the solution of these equations.
In the earliest finite element analyses, iterative methods
were used, among them the Gauss-Seidel method.
Now, the basic disadvantage of using an iterative method is
that we do not know how many iterations have to be
performed for solution, and thus it is very difficult to
estimate the total solution effort to solve the equations
KU equals R. Now we use almost exclusively direct methods, or
computer programs that are in industry, in abundance use,
are basically using direct methods.
And these direct methods are basically variations of Gauss
elimination, a method that Gauss some 100 years ago
proposed for the solution of equations.
Now we are talking about static condensation,
substructuring, frontal solution, LDL transpose
factorization, Cholesky decomposition.
We are also talking about the crude method, the column
reduction skyline solver.
But all of these methods are indeed variations of the basic
Gauss elimination procedure.
These are direct methods because the number of
operations can be counted for the solution of the equations.
In fact, the number of operations for a system with n
equations, a half bandwidth m, is 1/2 nm squared, where we of
course use a half bandwidth m that is
constant for the system.
If the bandwidth varied, some adjustments have to be made to
that variable m.
A mean value has to be used.
However, the important point is that we actually can count
the number of operations to obtain the solution, which was
not the case in the iterative method, because we did not
know how many iterations have to be performed.
What I would like to do now in this lecture is to very
briefly review with you how the basic Gauss elimination is
performed, and then I would like to show to you how this
method is indeed the basis of static condensation
substructuring, et cetera.
We will finally also consider the column reduction skyline
solver, which is very effectively used, and that is
in fact the equation solver that is used in Sab and the
ADINA program.
Well, I have prepared, as before, some view graphs here.
And the first view graph here shows a system of equations.
This is the K matrix here.
This is a displacement vector.
This a load vector.
And I will show you later on a physical system which in fact
gives this stiffness matrix K. But let's look now at the
system simply as an example for solving equations.
In the basic Gauss elimination procedure,
we proceed as follows.
We subtract a multiple of this equation here from the second
equation and the third equation in order to produce
here a 0, and to produce here a 0.
Well, to obtain a 0 here, we have to take 1/5 of this
equation, 1/5 of that equation, and subtract it from
this equation.
1/5 because then we obtain a 1 here, and that 1 knocks that 1
out to obtain a 0.
We have to subtract minus 4/5 of this equation from the
second equation, and we obtain a 0 here, too.
While this process is shown on the next view graph, the
result of that process is this set of equations.
We now have produced our 0 here, our 0 here.
Of course we had already a 0 here.
We will see later on that the matrix in here, the 3-by-3
matrix, bounded by the dash lines on the side, in fact
represents a stiffness matrix.
That's important, because a stiffness matrix has always
positive elements on the diagonal.
Otherwise it is not a stiffness matrix.
It does not govern a stable structure.
See, here we have diagonal elements that are all
positive, and if we can show that in the Gauss elimination
procedure, we are still always producing new stiffness
matrices, here a 3-by-3 matrix, then we directly can
conclude that the diagonal elements must remain positive.
Well, in the Gauss elimination procedure, the same process
that we use to obtain from the 4-by-4 or 3-by-3 system is now
carried on to obtain, from the 3-by-3 system,
a new 2-by-2 system.
This means that we have to subtract a multiple of the
third equation from the fourth equation and the fifth
equation in order to obtain a 0 here and a 0 here.
For example, 5/14 of the second equation is subtracted
from the fourth equation to obtain a 0 here.
The result of that operation is this matrix here.
Now we have a 2-by-2 matrix here.
Again, positive elements on the diagonal.
In fact, this is a 2-by-2 stiffness matrix, as I will be
showing to you just now.
Finally we proceed again, once more, with the same step, to
produce a 0 here, which means that we have to subtract from
this fourth equation a minus 20/7 divided by
15/7 of this equation.
And the result, then, is this set of equations here.
Notice that we now have a single element here, and in
fact, this is a stiffness matrix, in this case, a 1-by-1
stiffness matrix, just one element, that governs the
behavior of a single spring.
The solution off the equations is now completed by solving
for u4 from this relation down here.
This is a 1 degree of freedom system.
Knowing u4, we can go back to this equation here, substitute
u4 in, and the only unknown variable is u3.
We solve, therefore, u3 from this equation.
Now we know u4 and u3.
We now go to this equation.
Knowing u3 u4, the only unknown in this
equation here is now u2.
We can calculate u2, and repeating the process, we can
calculate u1.
A very simple, straightforward procedure,
but extremely powerful.
The result of that solution is shown here.
Notice u1 is simply obtained by this operation here.
Of course, here we're having u2 and u3 in that we
calculated already.
u2 is calculated here. u4 is calculated here.
It might be easier to actually look at it the way we
discussed it just now.
We call that the back substitution, because we are
coming from below upwards, and in the back substitution, we
have u4 first.
We then calculate u3.
u4 is known now, remember.
We can substitute directly u4 into here, and we get u3.
Knowing now u3 and u4, we can get u2, as shown here.
And finally, we get u2.
This is the basic Gauss elimination procedure, an
extremely powerful process.
And I would like to show you know how this procedure
relates to the static condensation and relates to
the other techniques, and how we then actually use this
basic Gauss elimination procedure in an actual finite
element program.
Well, static condensation is the first topic that I'd like
to talk about.
The basic process is here that we are partitioning the total
matrix into the following parts.
A Kaa, Kac, Kca, and Kcc part.
The cc part goes with Uc, the displacements Uc, which we
want to condense statically out.
Notice I'm putting the degrees of freedom that I want to
condense out at the bottom of the U vector.
Well, by static condensation we mean that we express Uc in
terms of Ua, which this equation here is obtained from
the bottom equation here.
We then substitute Uc in the top equation here, and the
result is this equation here.
Notice that the results can be remembered quite easily
because we have here an aa, Kaa, then we have an a here, a
c, a cc, ca.
Notice that one might say, just for remembering it, that
this c knocks out that c, this c knocks out that c, and what
we are left with are aa's.
And that means an aa here, an aa here.
So the dimensions are compatible of the matrices,
and the result is, I call it here, K bar aa.
Of course the right-hand side has also changed, because the
components must have an effect on the equation that give Ua.
Well, let's look at this process, called the static
condensation process, for our simple example.
And in this example now, I deliberately
put the Kcc up here.
In other words, the Uc degree of freedom is now the U1.
This is all our Uc vector only containing now one component.
The static condensation was originally proposed by
considering or statically condensing out of the last
equation, and that is how it's usually used also in finite
element analysis.
As an example, if we have incompatible modes, we might
condense out the incompatible degrees of freedom, just as
shown here.
Or if we have internal degrees of freedom in finite elements
that only that only couple with the bound the degrees of
freedom of a finite element, we could condense them out
this way, and then we would have this as the last degree
of freedom.
This would be Uc degrees of freedom, and we would condense
them out, just as shown here.
However, in order to show you how the static condensation
relates to Gauss elimination, I want to look at the example
that we looked at earlier already, and this the
partitioning I have used.
U1 now being the UC, U2 U4 being the
Ua degrees of freedom.
So I want to now statically condense out Uc.
The result of the K bar aa is shown here.
Notice we are taking this part, this part here, this
part here goes right there, this vector here is this
vector here, 1/5 is the inverse of Kcc, and this
vector here is that vector there.
Well, performing these operations, the
result is this one here.
That is the result.
And if we go back to our earlier Gauss elimination
procedure, we in fact notice that this matrix, K bar aa, is
nothing else than the matrix that we derived from the
4-by-4 system when we eliminated the non-zero
elements here.
In other words, we subtracted in the Gauss elimination from
the second equation and the third equation a multiple of
the first equation, we got these zeros here.
The result was this 3-by-3 matrix.
And that 3-by-3 matrix is exactly this matrix.
So what have we been doing then in the Gauss elimination?
We have really statically condensed out the first degree
of freedom, and we have obtained a stiffness matrix, K
bar aa, in that Gauss elimination procedure.
We can now also look at the physics that we have been
going through, and this is, in my opinion, a very
important part also.
We start off with this matrix here.
Notice the following00 this now is the stiffness matrix of
this beam for the U1 to the U4 degrees of freedom.
This stiffness matrix actually has been derived by finite
differences.
It is a nice example to show demonstratively what happens
in the Gauss elimination.
Of course, we could've also used a finite element system,
but I chose here a finite difference matrix, which is
still a stiffness matrix, and which I can use very nicely to
show to you the basic procedures of the Gauss
elimination.
So this stiffness matrix governed this system here.
Notice the following.
The first column, this column here, represents the forces.
When I displace this node, if you think of this as a note,
by unity, and keep all the other degrees
of freedom to 0--
in other words, if I put this deflected
shape onto the beam--
0 displacements right here, but a unit displacement here.
That is a unit displacement.
Now notice, these are the forces that have to act at
these degrees of freedom in order to give
this deflected shape.
The force here, if I put them here in
blue, is 5, as an example.
The force here is minus 4.
The force here is 1, and the force here turns out to be 0.
That element, or that element, of course, we
have a symmetric matrix.
And of course, in the same way, we can interpret other
elements in this matrix.
1 column represents always forces that have to act onto
the system in order to obtain unit displacements at 1 degree
of freedom.
If you look at the I's column, then it's the
I's degree of freedom.
And all the other degrees of freedom are set to 0.
Well, if we perform our Gauss elimination or static
condensation-- we have shown that it is the same--
on this system, in the first step, we
obtain a 3-by-3 matrix.
Physically, what we have obtained is the stiffness
matrix of this system now.
And the following important point is now to be noted.
If I look at this column just in the same way as I looked at
the first column of the 4-by-4 system, then this column
represents now the forces that are required at these degrees
of freedom to impose a unit displacement here at this
degree of freedom, this now being unity, with 0
displacements here.
And the force required here is 14/5.
That is the force that I have to apply here to get a unit
displacement here.
At the same time, I have to apply here a force of minus
16/5, and I have to apply here a force of 1.
So I'm talking again about forces acting onto the system
in this matrix here.
And these forces are stored in the columns or the rows off
the matrix, and in the I's column or row, I have the
forces required to impose a unit displacement at the I's
degree of freedom with all the degrees of freedom being 0.
And the important point is that having performed now the
Gauss elimination, here I do not have
anymore a degree of freedom.
The system is free to go down here.
And that's what we mean by static condensation, really.
We have taken a degree of freedom out of the system.
Proceeding in the same way, we obtain this 2 by 2 system.
Now we only have 2 degrees of freedom.
I have condensed out this degree of freedom and that
degree of freedom, which were there.
The beam is now free to move down here, so the deflected
shape corresponding to U3, this one here, would be like
that, with 0 displacement here, and unit
displacement there.
The forces required at these degrees of freedom are stored
in this column here.
Now, we immediately can observe that since we are
talking about a physical system into which we have to
put energy to deform it, we can conclude that these
elements here must be positive elements.
In other words, the diagonal elements in the stiffness
matrices that we derive must be positive.
They must be positive because if I push down the U3 degree
of freedom to obtain a unit displacement, surely the force
that I have to apply here must be positive.
In other words, it must be a force into the direction of
the degree of freedom U3.
And that means this must be a positive element.
So if, in a Gauss elimination procedure, we obtain, all of a
sudden, a negative diagonal element, this one or that one,
as an example, then we would know that the structure we are
analyzing is unstable.
In the final step of the Gauss elimination, we obtain a
single degree of freedom system.
And that single degree of freedom system
behaves like that.
Notice no more degree of freedoms here.
And the equation--
I have to apologize that this matrix has been written a
little bit down.
It should be a little bit higher.
They should, of course, be aligned these three.
But basically, we have 5/6 times U4 equals 7/6.
And this here is really nothing else than the
stiffness of a spring.
A spring K having us a number 5/6.
In other words, a spring stiffness of 5/6 here.
This is U4, the U4 degree of freedom.
And the force applied to that spring is 7/6.
We have a single degree of freedom system.
But notice that this single degree of freedom system, this
stiffness now, K equals 5/6, contains really all the
physical properties of that beam.
The effect of the other degrees of freedom, U1, U2,
U3, have been carried over into this spring stiffness.
Similarly, 7/6 is the force which is a result off the
other forces that were applied to the
other degrees of freedom.
And now from the single degree of freedom system, of course,
we can solve for U4, and then we go back to the 2 degree of
freedom system to solve for U3, and so on.
By that back substitution process, we obtain thus the
displacement that the beam is actually undergoing.
So the important point is that when we do perform Gauss
elimination, we really operate always on physical systems.
We operate always from one stiffness matrix on to another
stiffness matrix, and the next stiffness matrix is basically
obtained by static condensation.
We are statically condensing out 1 degree of
freedom after the next.
Well, if that is very well understood, and I urge you to
read up in the book on this procedure, then I think you
can very well appreciate how the substructuring process,
frontal solution method, directly relates, again, to
Gauss elimination.
And here I summarize once the basic procedure that we use in
substructuring analysis.
We really go through static condensation on the internal
degrees of freedom of a substructure.
Here I'm showing a substructure, which is an
assemblage, in this case, of the plane stress 4-node
elements, just for illustration.
Of course, there could be beam elements internally here or
externally here.
There could be shell elements.
Any kind of elements.
But the important point is that we have an assemblage of
elements, which by itself, of course, gives
us a structure already.
We use static condensation on the internal degrees of
freedom of the substructure.
Well, the result, then, is a new stiffness matrix of the
substructure involving boundary
degrees of freedom only.
Let us look, then, at this example.
Here we have 2 degrees of freedom at each node.
At the boundary nodes, which are shown solid, here, and the
hollow points here, the circles here are
the internal nodes.
2 degrees of freedom at each.
The result is that we have a 50-by-50 stiffness matrix,
because we have 1, 2, 3, 4, 5 nodes along this length, and 5
layers this way.
So 25 nodes times 2 degrees of freedom per node gives us a
50-by-50 stiffness matrix.
Now what we do in the substructuring analysis, we
statically condense out the internal degrees of freedom,
just in the way I have shown how we condense out in that
beam element, in that beam analysis, U1 to U3, the same
way we are condensing out the internal degrees of freedom.
Notice that we do not need to perform Gauss elimination from
U1 to U3 or U4 onwards.
We could also first deal with U2, then with U4, then with
U1, and then solve for U3.
So we have really a lot of flexibility in the Gauss
elimination procedure.
The important point is that now we're using it to condense
these internal degrees of freedom out.
The result, then, is a stiffness matrix that only
involves the boundary degrees of freedom.
Now in this case, we have 1, 2, 3, 4, 5 nodes here.
Another 5 is 10.
3 here and 3 here means 16, times 2 means it's going to be
is 32 by 32 stiffness matrix.
Notice that the internal degrees of freedom have been
statically condensed out, which means that if I impose,
for example, a unit displacement here, keeping all
the other displacements 0 at the other nodes, and I keep
this one also 0, then in this case, I would have internal
displacements at these nodes.
Whereas in this case, if I put a unit displacement only
corresponding to this degree of freedom and keep all other
displacements equal to 0, I would have 0 displacements at
these nodes, at these points.
Here these nodes are free to go wherever they want to go in
the finite element solution.
And if I put a unit displacement on here, they
will be taking on a particular value.
These nodal displacements will take on a particular value.
What we really have derived, and this is important, is a
new stiffness matrix.
A new, we might call it, element stiffness matrix.
A superelement, or substructure stiffness matrix.
These are all equivalent words that, however, only involve
boundary degrees of freedom.
This is now a new finite element that can be used just
as any ordinary finite element in the assemblage process to
obtain the stiffness matrix of the complete structure.
Let me show to you a simple example.
Here is a very simple structure, a truss element
that has only 3 degrees of freedom.
U1, U2, U3.
The stiffness matrix of that truss element is
shown here as 3-by-3.
Matrix here, we see U1, U2, U3, and the forces
corresponding to these degrees of freedom, of course, are
listed here.
The objective in this example is to statically condense out
this degree of freedom, and thus obtain a new element
stiffness matrix that only involve the boundary degrees
of freedom.
Well, the first step that I've pursued here was to rearrange
the equation so as to have U1 and U3 as the top equations,
and U2 as the bottom equation.
So we want to do static condensation now from the
back, or from the bottom, the way it's usually done.
And it's usually done that way because it is numerically
effective to proceed that way.
The result, then, is shown here.
We write down just to identify, to
use the various parts.
This matrix here is this 2-by-2 matrix.
This minus 20, minus 28 is that one here.
Here we get the Kcc inverse, which is this part here, and
then of course here we have this vector, shown here.
If we perform these multiplications here, we
obtain this stiffness matrix here.
A 1 minus 1 minus 11 with a constant in front.
Of course, the load vector also has changed.
There is a carry over from the 2 degree of freedom in to the
1 and 3 degrees of freedom.
The U2 degree of freedom then can be recovered via this
equation here, going back to the original equation.
The important point is that we have now obtained the
stiffness matrix of the truss element, this stiffness matrix
of the trust element, involving U1 and U3 only.
U2 is not anymore present.
It has been statically condensed out.
The U2 will adjust itself to a particular value, and that
value, of course, is given right here.
But the stiffness matrix that we have now here only
corresponds to U1 and U3.
Well, that process can then be repeated, and we can go
through a process that we call multilevel substructuring.
In other words, we don't just go once through a
substructuring process, but we go repeatedly through
substructuring processes.
And basically, what we are doing is we are solving the
system of equations KU equals R of the complete system in a
very effective way.
Let me show you here an example.
Here we have a simple truss assemblage.
We have 1, 2, 3, 4 elements.
The important point is here that the
elements are all similar.
The area is varying 2A1, 4A1, 8A1, 16A1.
And what I want to do is, I want to deal now, via
substructuring with the equations to obtain a very
effective solution.
The first step is to look at this element, which involves
U1 and U3 and U2.
And here you see, once more, a sketch of that element.
That is the element here.
Now I condense out U2, and that's what we just did in the
previous example.
The only degrees of freedom that we are left with,
then, are U1 U3.
So now we have an element stiffness matrix, this element
stiffness matrix, but in terms of U1 and U3 only.
The next step, then, is to say, let
us repeat this process.
Let us now say, since we have this element stiffness matrix
here in terms of U1 and U3 only, and since this adamant
is similar to this element here, there is only a scalar
involved because the area is varying as shown, from A1 to
2A1 to 4A1 and so on, we can use this stiffness matrix
here, this element stiffness matrix, which, however, is
already a substructure the way I've defined it, to assemble
the stiffness matrix corresponding to
both of these elements.
And that process is schematically shown here.
This stiffness matrix now will involve U1, U3, and U5.
Notice U2 is gone, and the important point U4 is gone
also, because I've dealt with U2, I have statically
condensed out U2, and by that, of course, I also have
statically condensed out U4 simultaneously, because I'm
using this element stiffness matrix here again.
So U4 is also gone already.
And now all I need to do is condense out U3, and I obtain
a second level substructure, which I can use again now to
assemble this element here with that element.
And that is shown on the next view graph here.
I repeat that process using this element here now, which
is this part here, involving only U1 and U5, and this
element here, which is similar to that element.
Again, an area variation there, of course.
And I obtain an element stiffness matrix, or a
structural stiffness matrix, whichever way you want to look
at it, that involves U1, U5, and U9 only.
I can now statically condense out U5.
The only unknowns I'm left with are U1 and U9, and of
course then I can simply solve for U1 and U9.
Of course, we would have to have also one of the
displacement being prescribed to
actually perform a solution.
Otherwise we have a rigid body mode in the system.
Assuming that, say, U1 is 0, we now would impose U1 being
0, solve for U9, go back, solve for U5, and proceed
backwards to solve for all the displacements.
The important point is that via this substructuring
process, a multilevel substructuring process, we
have saved operations.
We have saved operations because in this particular
case, you see, I have condensed out these degrees of
freedom already, and I use this fact in this assemblage
process and solution in such a way as that I do not need to
condense out these degrees of freedom, that I do not need to
deal with these degrees of freedom again, because I've
dealt with them in effect already by condensing out
these degrees of freedom here.
There are also other advantages using
substructuring, and one important advantage that I
should mention is that in the input definition phase, we
only need to define the elements, nodal points, and so
on of a typical substructure, and we can use that definition
in a computer program, then, again and again, to assemble
the complete structure.
In other words, if I have defined, typically for a
high-rise building, one typical floor, and I have
assembled the floor stiffness matrix with columns, beams,
and so on, that pertain to that floor, if I then
statically condense out the degrees of freedom that we
want to condense out for that part of the substructure, I
can use that floor substructure, so to say, again
and again, to assemble the complete building structure.
And this means I only have to define really the particular
substructure input, and I deal, then, with that
substructure numerically, as I have
shown here in the example.
The next important procedure that I like to discuss very
briefly with you is the frontal solution.
And the frontal solution really, again, only represents
a variation of the basic Gauss elimination procedure.
It represents, really, a static condensation procedure.
And that static condensation, however, is performed in a
particular way.
And this is really where the advantages of the frontal
solution can lie.
It's performed via looking at the elements.
Let me describe to you the basic process.
Here we have, as an example, a two-dimensional plane stress
element mesh.
We have 4-node elements, as shown here, element 1, element
2, element 3, element 4, q, q plus 1, and so on.
In other words, there's a whole lay of elements this
way, and that way, and so on.
This mesh could represent, for example, a cantilever plate
that acts in plane stress.
In the frontal solution, the solution of the KU equals R
equation is performed by elements.
So we're putting in element 1, which couples into node 1,
node 2, node m, and node n plus 1.
Since the degrees of freedom at node 1 only obtain
stiffness from element 1, and no other element, the
assemblage of this element into the global structure
stiffness matrix gives the only stiffness
contribution at node 1.
Therefore, we have the final equations at node 1, and we
can deal with them already.
And that's how the frontal solution proceeds now.
It statically condenses out these degrees of freedom.
Then we move on to the next node of element
1, this node here.
However, this node gets stiffness from element
1 and element 2.
So in order to obtain the final equations here, we have
to add the element 2 stiffness matrix to the element 1
stiffness matrix.
And that means in this way, then, of course, we obtain the
final equation, corresponding to node 2.
We now can statically condense out the degrees of
freedom at node 2.
We then move on to the next node, involving element 1 and
element 2, which is node 3, say.
And now we have to, however, assemble the stiffness matrix
of elements 3, also into the equations that we have
assembled already.
Having done so, we can now use static condensation to
condense out these degrees of freedom.
We proceed, therefore, via elements, and we are talking
about a wave front--
which is this one for node 1, this is the wave front for
node 2, this is the wave front--
in other words, what we are saying is, the wave front
really embodies or contains all the elements that we need
to have assembled before we can condense
out degrees of freedom.
To condense out the degrees of freedom at node 2, we have to
have a wave front such as shown here.
We have to have assembled element 1 and 2, because
element 1 and element 2 give stiffness to node 2.
The important point is that we are proceeding via static
condensation from node 1 to node 2 and so on, and this
means really that we are performing
Gauss elimination again.
The process is summarized here once more.
The frontal solution consists of successive static
condensation of nodal degrees for freedom.
The solution is performed in the order
of the element numbering.
That's important.
Same number of operations are performed in the frontal
solution as in the skyline solution, if the element
numbering in the wave front solution corresponds to the
nodal point numbering in the skyline solution.
So if the element numbering in the frontal solution
corresponds to the nodal point numbering in the skyline
solution, we perform in the skyline solution, which I have
not really talked about yet, and that's what I want to get
on now, the same number of operations as
in the frontal solution.
By corresponding, I mean in this particular case here, we
would condense out in the frontal solution the degrees
of freedom at node 1, then 2, then 3, then 4, and so on.
And of course, in the skyline solution, or in the Gauss
elimination procedure, we also proceed via node 1, 2,
3, 4, and so on.
So in this particular example, the frontal solution and the
skyline solution would give the same number of operations.
Well, let me go on to the skyline solution, then.
And the basis of that solution is LDL transpose
factorization.
Basically again, Gauss elimination.
The basic step here is the following.
We are having our stiffness matrix, K here, and now I'm
going back to that example which I discussed earlier, the
beam example.
We're talking about the 4-by-4 stiffness matrix again.
We're using this stiffness matrix, and we premultiply it
by an L1 inverse matrix, this matrix here.
In fact we never perform formally this
premultiplication via matrix, because what this does is, it
subtracts a multiple of this first row from the rows below
it to update a 0 here and a 0 there.
That's all this premultiplication does.
And the result is this matrix here, which I've shown to you
earlier already.
The next step is then to operate on this
matrix in the same way.
However, before doing so, let me show you here a view graph
in which I have written out specifically the L1 inverse.
And the important point, L1.
Notice that L1, the inverse of this matrix, is obtained by
trivial operation.
Namely, simply changing the plus here to a minus sign, and
the minus sign to a plus here.
Changing the signs of the off-diagonal elements gives us
directly the inverse of that matrix.
Well, having obtained, therefore, via this operation,
or having carried out via this operation, I should say, the
first step of the basic Gauss elimination procedure, namely
subtracting the top row from the rows below it to obtain
zeroes in the first column below the diagonal element, we
proceed in the same way.
and via this procedure here--
now, of course, we have an L2 inverse, and so on-- we obtain
what we call an upper triangular matrix.
S, an upper triangular matrix.
Zeroes all below the diagonal, and only non-zero elements
right here.
Of course, if the matrix has a bandwidth originally, that
bandwidth would be preserved, and we would have no non-zero
elements outside that band.
In other words, all elements outside the band would, in
fact, be also zero elements.
This upper triangular matrix, then, obtained--
this is the S matrix, here--
will be used in the back substitution.
This is here the element that gives us the single degree of
freedom stiffness, which I talked about earlier, from
which we can solve Un.
Here I'm talking about a general case with
n degrees of freedom.
However, looking at this equation, let us first
identify that we can put all of these inverse matrices onto
the right-hand side, and we obtain this equation here.
Now notice that this L1 is obtained by trivial operation
from L1 inverse.
The same holds for L2 and so on.
So we never really had to invert a matrix to
obtain this part here.
Furthermore, this L1 was an identity matrix with only
non-zero, off-diagonal elements in the first column.
This one is a matrix which looks like the identity
matrix, but has also non-zero, off-diagonal elements in the
second column.
And so on.
And the product of these matrices here, then, is very
simply obtained by writing the 1's into the diagonal
positions, and the off-diagonal elements here are
those of L1 in the first column, those of L2 in the
second column, and so on.
So even this multiplication is obtained in a very simple way.
We only need to write the first column of L1 into here,
the second column of L2 into here, the third
column of L3 into here.
And the result, all of that, I call the L matrix.
Now notice that once again, we do not perform any matrix
multiplication.
We have n minus 1 matrices here, but we do not perform a
matrix multiplication, because all we do is, we take the
first column of L1 and put it into the first column of L.
The second column of L2 is the second column in L. The n
minus second column in L is nothing else than the n minus
second column in L n minus 2.
And so on.
The result, therefore, is that we have written k equal to Lf,
L being equal to this product here.
Now because K is symmetric, we can, in fact, write this F
matrix here as D times L transposed, where d is a
diagonal matrix, dii being equal to sii.
So now we have factorized the K matrix into LDL transposed.
The d elements, or the dii elements, must all be positive
if our structure is stable, as I pointed out earlier.
The Cholesky factorization, just by remark, uses this
form, K being L curl times L curl transpose, where we
should note that L curl is really nothing else than L
times d to the 1/2, where d to the 1/2 is obtained by using
the D matrix and taking square roots of
the diagonal elements.
In other words, the ith elemental of d to the 1/2 is
nothing else than square root out of dii.
Well, having performed the LDL transpose factorization, we
now are ready to solve the equations, and they are solved
in the following steps.
We have LV equals R, where V is defined as shown here.
Notice that if I were to put this part into there, I really
obtain nothing else than LDL transpose times U equals R,
where this is, of course, nothing else than K times U
equals R, if I were to substitute
from here into there.
However, we separate these two steps out, and we call this
the forward reduction of R, which is really obtained by
subtracting specific elements of R from other elements.
We never really perform a matrix product here.
We write it that way, but really, it means nothing else
than taking, say, the ith element off R and subtracting
that element, or multiple of that element, from the
elements below it.
That is a typical operation performed when
we form the V vector.
Having obtained V, we go back to this equation here, and we
perform this step.
Notice that this step is nothing else than the division
off the V elements by the D elements.
In other words, here we take Vi and divide it by dii.
The final step, then, is the back substitution, which gives
us U. The back substitution is performed
as I discussed earlier.
Already we're going from UN, from the last degree of
freedom, upwards to calculate the values of
the degrees of freedom.
Let us now look at the actual implementation.
The implementation has to be performed very
effectively, of course.
All I have been giving to you now so far were the
mathematical equations that we are basically operating on.
The LDL transpose factorization is, of course, a
way of writing a set of steps of solution that we are
performing.
However, we want to perform them as effectively as
possible, and the scheme that we're using is related,
really, to the crude scheme--
we call it the column reduction scheme--
that we are going through for the solution of the equation.
And the procedure that we're using is as follows.
Here we have the stiffness matrix.
Notice that we only store the upper part.
Since we have a asymmetric matrix, we can directly obtain
the lower part elements from the upper part elements.
We do not need to store them in the computer program.
We also perform the Gauss elimination in a very, you
might say, strange way, when you look at it first.
And the way we're performing it is the following.
That we are reducing this column first.
In other words, we are going column-wise.
We are reducing this column first.
The result is that one.
In the next step, then, we are going through the reduction of
the next column.
See here, we are now reducing this column, and the
result is that one.
The dashed line always is on the right-hand side of that
part of the matrix, which has already been reduced.
And now we are looking here at the third column, and this is
the reduction off the third column The final step is then
that we are reducing the fourth column,
and this is the result.
Now notice that having performed these operations
column-wise--
and I have not given you the details.
We do not have time in this lecture
to look at the details.
The details are in the book, if you want to read them up.
All I'd like to convey to you now are the procedures that
we're using to show you the effectiveness.
Notice that the final result here in this operation is the
D matrix on the diagonal and the L transpose in the
off-diagonal elements.
And we have done it column-wise.
Well, what does this mean, then, in practice?
In practice it means the following.
Notice that I have now a system which has different
column heights for different equations.
The column heights are being defined by the first non-zero
element in a column.
So this is the column height for this last column here,
because this is the first non-zero
element in the column.
The zero elements below this element here, in general,
become non-zero elements.
In general, not always.
But because they in general become non-zero elements, we
have defined our column height this way.
And what we are now doing in the actual solution is that we
are reducing, column by column, just the way I have
shown it, to obtain finally the elements of the D matrix
on the diagonal, and the elements of the L transpose
matrix in the off-diagonal elements here.
And we can do so--
this is important--
by simply replacing elements.
If you go back to an earlier step that we looked at, what
we have done here is looked at this matrix here, and we, in
the reduction process, simply replaced this element by 1/5,
replaced this element by the minus 4 element, by minus 8/7,
over and so on.
So we don't need additional storage.
We can simply replace elements.
In the same way, we proceed now here, as shown, on the
next view graph.
Here we have filled in the zero element, and the
reduction has been performed column after column.
The effectiveness of this whole procedure is really seen
on this view graph here, when we talk about large systems.
I want to denote here this as a large system, because the
high speed core storage is not sufficient to store the
whole matrix in.
And so what we have been doing is, we are storing
the matrix in blocks.
This is here the first block, shown by the dashed line.
This is the second block.
This is the third block, and that's the fourth block.
What I'm saying is that we have divided the total
stiffness matrix into blocks, and in fact, we divided the
blocks in such a way as to be able to take always at least
two blocks into core.
In other words, in this solution process, we have to
be able to take this block into core,
that block into core.
Two blocks.
We also would have to be able to take this block into core,
that block into core, and so on.
Of course, the way this is actually automatized in a
computer program, is that we say we have so
much storage available.
Say that we have 40,000 storage locations available.
Then for the solution of the equations, or for the storing
of the stiffness matrix, then 20,000 locations
would be used per block.
And what the solution routine simply does-- it counts these
elements and puts this dashed line in such a way that 20,000
elements are in here, 20,000 elements are in here, 20,000
elements are in here, and so on.
The solution procedure, or the solution of the equation, is
now performed in the following way,
using the block structure.
We proceed with the reduction off this column, of that
column, of that column and that column.
Now, I should point out the following important point.
When we do we do reduce this column, this column coupled
into the first equation.
So we have to have these elements in core.
Of course, since this is 1 block, and since we can store
that block in core anyway, there is no problem.
We can reduce this column directly.
These elements that we already have obtained
are indeed in core.
Now let's go on to the second block.
The second block coupled into the first block, because these
elements here couple into these equations.
So in order to reduce the second block column by column,
we have to have the first block in core.
And that's what we do.
We take this block into core and reduce the second block.
Having finished the second block reduction, we go on to
the fourth block.
And the fourth block now is reduced by first taking the
first block into core, because the first block also couples
into the third block.
Then taking the second block into core together with this
one, and reducing this one, and then of course, finally,
we reduce this one.
So the general procedure, therefore, is that we take a
particular block into core that we want to reduce, we
have reduced already all the previous blocks, we have
storage available also for one of the previous blocks, and we
sequentially take the previous blocks into core to take into
account all the coupling that these previous blocks have
onto the current block that we want to reduce.
In this example, we're taking this block into core.
We then take first also this block into core.
That means we have two blocks in core.
We take the effect of this block into
account onto that block.
Then we throw this book out.
We take this block into core.
We take the effect of this block into account onto this
block, and finally, we reduce the block by itself.
In the same way, we would proceed with this block, just
as another example.
We take this block into a core.
This one couples into all previous blocks, so what we
have to do is take this into core, and first take this into
core, take the effect of this coupling into account onto
this block, throw this out, take this one into core, put
this block coupling onto that block, throw this out, take
this one into core, and finally take the effect of
this block coupling onto that one into account, and finally
then, of course, reduce this block by itself.
This, then, completes the complete
solution of a large system.
I only was able to give you the basic principles, the
basic ideas.
Some of the details are described in the textbook, and
you can read up there further on it.
Thank you very much for your attention.
This was all what I wanted to say in this lecture.