NARRATOR: The following content is
provided by MIT OpenCourseWare under a Creative Commons
license.
Additional information about our license
and MIT OpenCourseWare in general
is available at ocw.mit.edu.
PROFESSOR: So I'm thinking of large sparse matrices.
So the point about A -- the kind of matrices that we've been
writing down, maybe two-dimensional,
maybe three-dimensional difference matrices.
So it might be five or seven non-zeros on a typical row
and that means that you don't want to invert such a matrix.
These are big matrices.
The order could easily be a million.
What you can do fast is multiply A by a vector,
because multiplying a sparse matrix by a vector,
you only have maybe five non-zeros times the n --
the vector of length n -- so say 5n operations, so that's fast.
I'm thinking that we are at a size where elimination,
even with a good ordering, is too expensive in time
or in storage.
So I'm going from direct elimination methods
to iterative methods.
Iterative means that I never actually get
to the perfect answer x, but I get close and what I want to do
is get close quickly.
So a lot of thought has gone into creating
iterative methods.
I'll write down one key word here.
Part of the decision is to choose a good preconditioner.
I'll call that matrix P. So it's a matrix that's supposed to be
-- that we have some freedom to choose and it very much speeds
up -- so the decisions are, what preconditioner to choose?
A preconditioner that's somehow like the matrix A,
but hopefully a lot simpler.
There are a bunch of iterative methods
and that's what today's lecture is about.
Multigrid is an idea that just keeps developing.
It really is producing answers to very large problems very
fast, so that will be the following lectures.
And then come these methods -- you may not be familiar with
that word Krylov.
Let me mention that in this family, the best known method
it is called conjugate gradients.
The conjugate gradient method.
So I'll just write that down, conjugate gradient.
So that will come next week.
That's a method that's a giant success
for symmetric, positive definite problem.
So I don't assume in general that A is
symmetric or positive definite.
Often it is if it comes from a finite difference
or finite element approximation.
If it is, then conjugate gradients
is highly recommended.
I'll start with the general picture.
Let me put the general picture of what
an iteration looks like.
So an iteration computes the new x.
Let me call it x-new or x_(k+1) from the known x.
So we start with x_0, any x_0, In principle.
It's not too important for x_0 to be close.
In linear methods, you could start even at 0.
In nonlinear methods, Newton methods,
you do want to be close, and a lot of attention
goes into a good start.
Here the start isn't too important;
it's the steps that you taking.
Let me -- so I guess I'm going to --
I'm trying to solve A*x equal b.
I'm going to split that equation into --
I'm just going to write the same equation this way.
I could write it as x equals I minus A x plus b.
That would be the same equation, right?
Because x and x cancel.
A*x comes over here and I have A*x equal b.
I've split up the equation and now actually, let me bring --
so that's without preconditioner.
You don't see a matrix P. That's just an ordinary iteration,
no preconditioner.
The preconditioner will come -- I'll have P and it'll go here
too.
Then, of course, I have to change that to a P.
So that's still the same equation, right?
P*x is on both sides, cancels -- and I have A*x equal b.
So that's the same equation, but it suggests the iteration
that I want to analyze.
On the right-hand side is the known approximation.
On the left side is the next approximation, x_(k+1).
I multiply by P so I -- the idea is, P should be close to A.
Of course, if it was exactly A -- if P equaled A,
then this wouldn't be here, and I would just be solving,
in one step, A*x equal b.
So P equal A is one extreme.
I mean, it's totally, perfectly preconditioned,
but the point is that if P equals A,
we don't want to solve that system,
we're trying to escape from solving that system.
But we're willing to solve a related system that is simpler.
Why simpler?
Here are some possible preconditioners.
P equals the diagonal.
Just take the diagonal of A and that choice
is named after Jacobi.
That's Jacobi's method.
Actually, it's already a good idea,
because the effect of that is somehow to properly scale
the equations.
When I had the identity here, I had no idea
whether the identity and A are in the right scaling.
A may be divided by delta x squared or something.
It could be totally different units, but P --
if I take the diagonal of A, at least they're comparable.
Then you see what -- we'll study Jacobi's method.
So that's a very simple choice.
Takes no more work than the identity, really.
Second choice would be -- so that would be P for Jacobi.
A second choice, that you would think of pretty fast,
is the diagonal and the lower triangular part.
So it's the -- I'll say the triangular --
I'm going to say lower triangular,
triangular part of A, including the diagonal.
This is associated with Gauss's great name and Seidel.
So that's the Gauss-Seidel choice.
So why triangular?
We'll analyze the effect of this choice.
If P is triangular, then we can solve --
I'm never going to write out an inverse matrix.
That's understood, but I can compute the new x
from the old x just by back substitution.
I find the first component of x_(k+1), then the second,
then the third.
It's just because the first equation -- if P is triangular,
the first equation will only have one term,
one new term on the left side and everything else will be
over there.
The second equation will have a diagonal and something
that uses the number I just found for the first component
and onward.
So in other words, triangular matrices are good.
Then people thought about combinations of these.
A big lot of activity went in something called
overrelaxation, where you did a --
you had a weighting of these and it turned out to be if you
adjusted the weight, you could speed up the method,
sometimes called SOR-- successive overrelaxation.
So when I was in graduate school, that was a very --
that was sort of the beginning of progress beyond Jacobi
and Gauss-Seidel.
Then after that came a very natural --
ILU is the next one that I want -- the I stands for incomplete,
incomplete LU.
So what are L and U?
That's the letters I always use as a signal for the lower
triangular and upper triangular factors of A, the ones
that exact elimination would produce.
But we don't want to do exact elimination.
Why?
Because non-zeros fill in.
Non-zeros appear in those spaces where A itself is 0.
So the factors are much -- have many more non-zeros than
the original A. They're not as sparse.
So the idea of incomplete LU is -- and I'll come back to it --
the idea of incomplete LU, you might say,
should have occurred to people earlier --
only keep the non-zeros that are non-zeros in the original A.
Don't allow fill-in or allow a certain amount of fill-in.
You could have a tolerance and you would include a non-zero
if it was big enough, but the many, many small non-zeros
that just fill in and make the elimination long,
you throw them out.
So P is -- in this case, P is L --
is an approximate L times an approximate U.
So it's close to A, but because we've refused some
of the fill-in, it's not exactly A.
And you have a tolerance there where if you set the tolerance
high, then these are exact, but you've got all that fill-in;
if you set the tolerance to 0, you've got the other extreme.
So that would give you an idea of preconditioners,
of reasonable choices.
Now I'm going to think about -- first, how do --
how to decide whether the x's approach the correct answer?
I need an equation for the error and because my equations are
linear, I'm just going to subtract --
I'll subtract the upper equation from the lower one
and I'll call the -- so the error e_k will be x minus
the true answer.
This is x exact.
Minus x, my k-th approximation.
So when I subtract that from that,
I have e_(k+1) and here I have P minus A --
when I subtract that from that, I have e_k,
and when I subtract that from that, 0.
So very simple error equation.
So this is an equation, this would be the error equation.
You see the b has disappeared, because I'm only
looking at differences now.
Let me write that over on this board,
even slightly differently.
Now I will multiply through by P inverse.
That really shows it just so clearly.
So the new error is, if I multiply by P inverse,
P inverse times P is I, and I have P inverse A, e_k.
So what does that say?
That says that at every iteration,
I multiply the error by that matrix.
Of course, if P equals A, if I pick
the perfect preconditioner, then P inverse A is the identity,
this is the zero matrix, and the error in the first step is 0.
I'm solving exactly.
I don't plan to do that.
I plan to have a P that's close to A, so in some way,
this should be close to I. This whole matrix should
be close to 0, close in some sense.
Let me give a name for that matrix.
Let's call that matrix M. So that's
the iteration matrix, which we're never going to display.
I mean, it would be madness to display.
We don't want to compute P inverse explicitly or display
it, but to understand it -- to understand convergence or not
convergence, it's M that controls everything.
So this is what I would call pure iteration or maybe
another word that's used instead of pure
is stationary iteration.
What does stationary mean?
It means that you do the same thing all the time.
At every step, you just use the same equation, where --
these methods will be smarter.
They'll adapt.
They'll use more information.
They'll converge faster, but this is the simple one --
and this one plays a part in those.
So what's the story on convergence now?
What about that matrix M?
If I take a starting vector, e_0 --
my initial error could be anything.
I multiply by M to get e_1.
Then I multiply by M again to get e_2.
Every step, I multiply by M. So after k steps,
I've multiplied k times by M. I want to know --
so the key question is, does that go to 0?
So the question, does it go to 0, and if it does, how fast?
The two parts, does it go to 0 and how fast,
are pretty much controlled by the same property of M. So
what's the answer?
When does -- if I take -- I don't know anything about e_0,
so what property of M is going to decide whether its powers
get small, go to 0, so that the error goes to 0?
That's convergence.
Or, maybe blow up, or maybe stay away from 0.
What controls it?
One word is the answer.
The eigenvalues.
It's the eigenvalues of M, and in particular, it's
the biggest one.
Because if I have an eigenvalue of M, as far as I know,
e_naught could be the eigenvector,
and then every step would multiply
by that eigenvalue lambda.
So the condition is that all eigenvalues of M
have to be smaller than 1, in magnitude of course,
because they could be negative, they could be complex.
So that's the total condition, that's
exactly the requirement on M. How do
we say how fast they go to 0?
How fast -- if all the eigenvalues are below 1,
then the slowest convergence is going to come
for the eigenvalue that's the closest to 1.
So really it will be -- and this is called the spectral radius
and it's usually denoted by a Greek rho of M,
which is the maximum of these guys.
It's the max of the eigenvalues.
And that has to be below 1.
So it's this number, that number,
because it's that number which really says what I'm
multiplying by at every step.
Because most likely, e_naught, my initial unknown error,
has some component in the direction of this biggest
eigenvalue, in the direction of its eigenvector.
Then that component will multiply at every step
by lambda, but that one will be the biggest guy, rho.
So it's the maximum eigenvalue.
That's really what it comes to.
What is the maximum eigenvalue?
Let me take Jacobi.
Can we figure out the maximum eigenvalue for Jacobi?
Of course, if I don't know anything about the matrix,
I certainly don't want to compute eigenvalues.
I could run Jacobi's method.
So again, Jacobi's method, I'm taking
P to be the diagonal part.
Let me make that explicit and let me take the matrix
that we know the best.
I'll call it K. So k -- or A, this is the A --
is my friend with 2's on the diagonal, minus 1's above.
I apologize -- I don't really apologize,
but I'll pretend to apologize -- for bringing this matrix
in so often.
Do we remember its eigenvalues?
We computed them last semester, but that's another world.
I'll say what they are.
The eigenvalues of A -- this isn't M now --
the eigenvalues of this matrix A -- I see 2's on the diagonal.
So that's a 2.
That accounts for the diagonal part, 2 times the identity.
I have a minus and then the part that comes from these 1's,
which are forward and back.
What would von Neumann say?
What would von Neumann do with that matrix?
Of course, he would test it on exponentials and he would get 2
minus e to the i*k minus e to the minus i*k.
He would put the two exponentials together
into cosines and that's the answer.
It's 2 minus 2 cosine of whatever angle theta.
Those are the eigenvalues.
The j-th eigenvalue is -- the cosine is at some angle
theta_j.
It's worth since -- maybe just to take again a minute on this
matrix.
So the eigenvalues are between 0 and 4.
The eigenvalues never go negative, right?
Because the cosine can't be bigger than 1.
Actually, for this matrix the cosine doesn't reach 1.
So the smallest eigenvalue is 2 minus 2 times
that cosine close to 1.
It's too close for comfort, but it is below 1,
so the eigenvalues are positive.
And they're between 0 and 4.
At the other extreme, theta is near pi,
the cosine is near minus 1, and we have something near 4.
So the eigenvalues of that -- if you look at that matrix,
you should see.
Special matrix eigenvalues in this interval, 0 to 4.
The key point is, how close to 0?
Now what about M -- what about P first?
What's P?
For Jacobi, so P for Jacobi is the diagonal part,
which is exactly 2I.
That's a very simple diagonal part, very simple P.
It produces the matrix M. You remember,
M is the identity matrix minus P inverse A. That's beautiful.
So P is just 2I.
So I'm just dividing A by 2, which puts a 1 there,
and then subtracting from the identity.
So it has 0's on the diagonal.
That's what we expect.
It will be 0's on the diagonal and off the diagonal, what
do I have?
Minus -- with that minus, P is 2, so it's 1/2.
It's 1/2 off the diagonal and otherwise, all 0's.
That's the nice matrix that we get for M in one dimension.
I have to say, of course, as you know,
we wouldn't be doing any of this for this one-dimensional matrix
A. It's tridiagonal.
We can't ask for more than that.
Direct elimination would be top speed,
but the 2D case is what we understand that it has a big
bandwidth because it's two-dimensional,
three-dimensional even bigger -- and the eigenvalues will come
directly from the eigenvalues of this 1D case.
So if we understand the 1D case, we'll
knock out the two- and three-dimensional easily.
So I'll stay with this one-dimensional case,
where the matrix M -- and I'll even,
maybe should write down exactly what the iteration does --
but if I'm looking now at convergence,
do I have or don't I have convergence, and how fast?
What are the eigenvalues and what's the biggest one?
The eigenvalues of A are this and P is just --
P inverse is just 1/2 times the identity.
So what are the eigenvalues?
Again, I mean, I know the eigenvalues of A,
so I divide by 2.
Let me write down what I'm going to get now
for the eigenvalues of M, remembering this connection.
So the eigenvalues of M are 1 minus 1/2 times the eigenvalues
of A. Better write that down.
1 minus 1/2 times the eigenvalues of A.
Of course, it's very special that we
have this terrific example matrix where
we know the eigenvalues and really can see what's happening
and we know the matrix.
So what happens?
I take half of that -- that makes that a 1.
When I subtract it from that 1, it's gone.
1/2 of that, that 2 is cancelled -- I just get cos theta.
Theta sub -- whatever that angle was.
It's a cosine and it's less than 1.
So convergence is going to happen and convergence is going
to be slow or fast according as this is very near 1 or not --
and unfortunately, it is pretty near 1, because the first,
the largest one -- I could tell you what the angles are.
That's -- to complete this information,
I should have put in what -- and von Neumann got them right.
That's j -- is there maybe a pi over N plus 1?
I think so.
Those are the angles that come in,
just the equally spaced angles in the finite case.
So we have N eigenvalues.
And the biggest one is going to be when j is 1.
So the biggest one, rho, the maximum of these guys,
of these lambdas, will be the cosine of, when j is 1,
pi over N plus 1.
That's the one that's slowing us down.
Notice the eigenvector that's slowing us down,
because it's the eigenvector that
goes with that eigenvalue that's the biggest eigenvalue that's
going to hang around the longest, decay the slowest.
That eigenvector is very smooth.
It's the low frequency.
It's frequency 1.
So this is j equal to 1.
This is the low frequency, j equal 1, low frequency.
The eigenvector is just samples of the sine,
of a discrete sine.
Lowest frequency that the grid can sustain.
Why do I emphasize which eigenvector it is?
Because when you know which eigenvector is slowing you down
-- and here it's the eigenvector you would think would be
the simplest one to handle -- this is the easiest
eigenvector, higher-frequency eigenvectors have eigenvalues
that start dropping.
Cosine of 2*pi, cosine of 3*pi over N plus 1 is going to be
further away from 1.
Now I have to admit that the very highest frequency,
when j is capital N -- so when j is capital N I have to admit it
happens to come back with a minus sign,
but the magnitude is still just as bad.
The cosine of N*pi over N plus 1 is just the negative of that
one.
So actually, we have two eigenvalues
that are both hitting the spectral radius
and both of them are the worst, the slowest converging
eigenvectors.
This eigenvector looks very smooth.
The eigenvector for that high-frequency one
is the fastest oscillation.
If I graph the eigenvectors, eigenvalues, I will.
I'll graph the eigenvalues.
You'll see that it's the first and the last that
are nearest in magnitude to 1.
I was just going to say, can I anticipate
multigrid for a minute?
So multigrid does -- is going to try to get these guys to go
faster.
Now this one can be handled, you'll see,
by just a simple adjustment of Jacobi.
I just weight Jacobi a little bit.
Instead of 2I, I take 3I or 4I.
3I would be very good for the preconditioner.
That will have the effect on this high frequency one
that it'll be well down now, below 1.
This will still be the slowpoke in converging.
That's where multigrid says, that's low frequency.
That's too low.
That's not converging quickly.
Take a different grid on which that same thing
looks higher frequency.
We'll see that tomorrow.
So if I stay with Jacobi, what have I learned?
I've learned that this is rho.
For this matrix, the spectral radius
is that number and let's just see how big it is.
How big is -- that's a cosine of a small number.
So the cosine of a small number, theta near 0,
the cosine of theta is 1 minus -- so it's below 1,
of course -- 1/2 theta squared.
Theta is pi over N plus 1 squared.
All these simple estimates really
tell you what happens when you do Jacobi iteration.
You'll see it on the output from -- I mean, you can --
I hope will -- code up Jacobi's method,
taking P to be a multiple of the identity, so very fast,
and you'll see the convergence rate, you can graph.
Here's an interesting point and I'll make it again.
What should you graph in seeing --
to display convergence and rate of convergence?
Let me draw a possible graph here.
Let me draw two possible graphs.
I could graph -- this is K. As I take more steps,
I could graph the error, well, the size of the error.
It's a vector, so I take its length.
I don't know where I start -- at e_naught.
This is k equals 0.
This is the initial guess.
If I use -- I'm going to draw a graph,
which you'll draw much better by actually getting MATLAB to do
it.
Before I draw it, let me mention the other graph we could draw.
This tells us the error in -- this is the error
in the solution.
If I think that's the natural thing to draw,
it is, but we could also draw the error
in the equation, the residual error.
So I'll just put up here what I mean by that.
If I put in -- so A*x equals b exactly.
A*x_k is close.
The residual is -- let me call it r, as everybody does --
is the difference A*x minus A*x_k It's how close we are
to b, so it's b minus A*x_k.
The true A*x gets b exactly right.
The wrong A*x gets b nearly right.
So you see, because we have this linear problem,
it's A times the error. x minus x_k is e_k.
So this is the residual and I could graph that.
That's how close my answer comes to getting b right,
how close the equation is.
Where this is how close the x_k is.
So that would be the other thing I could graph.
r_k, the size of the residual, which is A*e_k.
So now I have to remember what I think these graphs look like,
but of course, it's a little absurd for me to draw them when
-- I think what we see is that the residual drops pretty fast.
The residual drops pretty fast and keeps going,
for, let's say, Jacobi's method, or, in a minute or maybe next
time, another iterative method --
but the error in the solution starts down fast
and what's happening here is the high-frequency components are
getting killed, but those low-frequency components are
not disappearing, because we figured out that they give
the spectral radius, the largest eigenvalue.
So somehow it slopes off and it's --
the price of such a simple method is that you don't get
great convergence.
Of course, if I put this on a log-log plot,
the slope would be exactly connected
to the spectral radius.
That's the rate of decay.
It's the -- e_k is approximately --
the size of e_k is approximately that worst eigenvalue
to the k-th power.
Do you see how this can happen?
You see, this can be quite small,
the residual quite small.
I think it's very important to see the two
different quantities.
The residual can be quite small, but then, from the residual,
if I wanted to get back to the error,
I'd have to multiply by A inverse
and A inverse is not small.
Our matrix A is pretty -- has eigenvalues close to 0.
So A inverse is pretty big.
A inverse times A*e_k -- this is small but A inverse picks up
those smallest eigenvalues of A --
picks up those low frequencies and it's bigger,
slower to go to 0 as we see in this picture.
So it's this picture that is our problem.
That shows where Jacobi is not good enough.
So what more to say about Jacobi?
For this matrix K, we know its eigenvalues.
The eigenvalues of M we know exactly.
The matrix M we know exactly.
Do you see in that matrix -- now if I had said,
look at that matrix and tell me something about its
eigenvalues, what would you say?
You would say, it's symmetric, so the eigenvalues are real
and you would say that every eigenvalue is smaller than 1.
Why?
Because every eigenvalue is -- if you like to remember this
business of silly circles, every eigenvalue is in a circle
centered at this number -- in other words, centered at 0 --
and its radius is the sum of those, which is 1.
So we know that every eigenvalue is in the unit circle
and actually, it comes inside the unit circle,
but not very much -- only by 1 over N squared.
What do I learn from this 1 over N squared?
How many iterations do I have to take
to reduce the error by, let's say, 10 or e or whatever?
If the eigenvalues are 1 minus a constant over N squared,
I'll have to take the N squared power.
So if I have a matrix of order 100,
the number of iterations I'm going to need
is 100 squared, 10,000 steps.
Every step is fast, especially with the Jacobi choice.
Here's my iteration -- and P is just a multiple of the identity
here.
So every step is certainly fast.
Every step involves a matrix multiplication
and that's why I began the lecture by saying,
matrix multiplications are fast.
A times x_k -- that's the work and it's not much.
It's maybe five non-zeros in A -- so 5N.
That's fast, but if I have to do it 10,000 times
just to reduce the error by some constant factor,
that's not good.
So that's why people think, okay, Jacobi gives us the idea.
It does knock off the high frequencies quite quickly,
but it leaves those low frequencies.
Well, I was going to mention weighed Jacobi.
Weighted Jacobi is to take, instead of p --
put in a weight.
So weighted Jacobi -- let me put this here and we'll come back
to it.
So weighted Jacobi simply takes P
to be the diagonal over a factor omega.
For example, omega equals 2/3 is a good choice.
Let me draw the -- you can imagine that if I --
this is 2 times the identity, so it's just 2 over omega I --
I'm just choosing a different constant.
I'm just choosing a different constant in P
and it has a simple affect on M. I'll draw the picture.
The eigenvalues of M are these cosines.
So I'm going to draw those cosines, the eigenvalues of M.
So they start -- I'll just draw the picture of cos theta if I
can.
What does a graph of cos theta look like?
From 0 to pi?
So I'm going to see -- this is theta equals 0.
This is theta equal pi.
If I just draw cos theta -- please tell me if this
isn't it.
Is it something like that?
Now where are these theta_1 -- the bad one -- theta_2,
theta_3, theta_4, theta_5, equally bad, down here.
So there's cos -- there's the biggest one.
That's rho and you see it's pretty near 1.
Here I'm going to hit -- at this frequency,
I'm going to hit minus rho.
It'll be, again, near 1.
So the spectral radius is that, that's rho.
The cosine of that smallest angle, theta over N plus 1.
What happens when I bring in these weights?
It turns out that -- so, now the eigenvalues of M,
with the weights, will be 1 -- it turns out --
1 minus omega plus omega cos theta.
You'll see it in the notes, no problem.
So it's simply influenced by omega.
If omega is 2/3 -- if omega is 2/3,
then I have 1 minus omega is 1/3 plus 2/3 cos theta.
Instead of graphing cos theta, which I've done --
let me draw -- let me complete this graph for cos theta.
The eigenvalues unweighted, with weight omega equal 1.
Now it's this thing I draw.
So what happens?
What's going on here?
For a small theta, cosine is near 1.
I'm again near 1.
In fact, I'm probably even little worse, probably up here.
I think it goes down like that.
So it starts at near 1, but it ends when theta is pi.
This is 1/3 minus 2/3.
It ends at minus 1/3.
You see, that's a lot smarter.
This one ended at minus 1, that one ended at near minus 1,
but this one will end with omega -- this is the omega equal 2/3,
the weighted one.
All you've done is fix the high frequency.
When I say, oh, that was a good thing to do.
The high frequencies are no longer a problem.
The low frequencies still are.
So the only way we'll get out of that problem
is moving to multigrid.
Can I, in the remaining minute, report on Gauss-Seidel
and then I'll come back to it?
But quick report on Gauss-Seidel.
So what's the difference in Gauss-Seidel?
What does Gauss-Seidel do?
It uses the latest information.
As soon as you compute an x, you use it.
Jacobi computed all -- had to keep all these x's until it
found all the new x's, but the Gauss-Seidel idea is,
as soon as you have a better x, put it in the system.
So the storage is cut in half, because you're not
saving a big old vector while you compute the new one.
I'll write down Gauss-Seidel again.
So it's more up to date and the effect
is that the Jacobi eigenvalues get squared.
So you're squaring numbers that are less than 1.
So Gauss-Seidel, this gives the Jacobi eigenvalues squared.
The result is that if I draw the same curve for Gauss-Seidel,
I'll be below, right?
Essentially, a Gauss-Seidel step is worth two Jacobi steps,
because eigenvalues get squared as they
would do in two Jacobi steps.
If I do Jacobi twice, I square its eigenvalues.
So it seems like absolutely a smart move.
It's not brilliant, because this becomes cosine squared.
I lose the 1/2, but I'm still very, very close to 1.
It still takes order N squared iterations to get anywhere.
So Jacobi is, you could say, replaced
by Gauss-Seidel as being one step better,
but it's not good enough.
And that's why the subject kept going.
This gives rho as 1 minus order of 1
over N, first power, where these were 1 minus 1 over N squared.
So this is definitely further from 1
and then this can be way further from 1.
So you can see some of the numerical experiments
and analysis that's coming.
I'll see you Wednesday to finish this and move into multigrid,
which uses these guys.