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 to view additional materials
from hundreds of MIT courses, visit MIT OpenCourseWare
at ocw.mit.edu.
PROFESSOR: So today we're going to talk
about partial differential equations.
Let me just give a quick review where we are.
So we're going to have a PDE lecture today.
Then Professor Swan is going to do a review on Friday.
There's no lecture on Monday, but you have the quiz, too,
Monday night.
Wednesday, next week is Veterans Day so it's a holiday.
There's no classes at MIT.
So the next time I'll be lecturing to you
will be a week from Friday.
The homework following the quiz will be posted
and so you guys can get started on that.
That homework involves a COMSOL problem
which, for those of you who had trouble with COMSOL,
might take some extra time, since it's
the first time you're doing it.
So I would suggest you at least give it a try early.
I'm sure the TAs would be very happy to help you
to force COMSOL to cooperate.
Bend COMSOL to your will.
All right, so we're going to talk about partial differential
equations.
And, as it sounds, what's different
about a partial differential equation is,
it has partial derivatives in it.
And so we have a lot of equations like this.
So for example, I work a lot on this equation.
It's a poor choice, isn't it?
OK, so this is the conservation equation for a chemical species
in a reacting flow.
So the change in the mass fraction of the species z,
at a point, xyz at time, t.
That differential is given by a diffusion term, a convection
term, and a reaction term.
So this is a pretty famous equation.
This is part of the Navier-Stokes equations
for the whole reactive flow system, which
I hope you guys have seen already,
and for sure you will see, if you haven't seen them yet.
And there's a lot of other partial differential equations.
I work a lot, also, on the Schrodinger equation.
That's another partial differential equation.
And what's special about the partial differential
equations is that, in this case, this partial derivative
is respect to time, holding all the spatial coordinates fixed.
And these partial derivatives are respect
to space, holding time fixed.
A very important thing about partial derivatives,
which you probably have encountered in 1040 already,
is, it really depends what you hold fixed.
And you get different results if you
hold different things fixed.
Now the implication, when you write an equation like this,
is that whatever partials you see--
When I wrote this, if I have this
and it has a term that's like d squared,
dx squared, is one of the terms in this thing.
This says t.
This says x.
The convention is, oh boy, I better hold x fixed.
And here must be, I must hold t fixed.
Because there's another one in the same equation.
Now, you can change coordinates however you want.
Just really watch out that you carefully
follow the rules about what you hold
fixed, when you change things.
Otherwise you can cause all kinds of craziness.
Now you guys probably took multivariable calculus
at some point in your distant past.
And they told you very carefully what the rules were.
So follow them.
And I suppose, when you're doing thermo right now?
Have you started to do all this partial derivative stuff?
Yeah, so you've seen how completely confusing
it can be with negative signs showing up all over the place.
And other things like this.
I don't know if you've encountered this yet?
At least, I thought it was confusing when I took it.
So just be aware.
It's the same thing.
Now, often you do want to change the equations.
So you can write down an equation like this.
But, for example, if you have a special symmetry
of the problem, like it's cylindrical, for example,
then you might want to change the cylindrical coordinates.
And then you have to just be careful
that you write that correct Laplacian and cylindrical
coordinates.
Make sure it doesn't mess up anything about what
you're holding fixed.
Sometimes you might want to be crafty and use
some weirdo coordinates.
For example, if you wanted to solve the Schrodinger equation
for h2 plus, it turns out there's
a special coordinate system called
the elliptic coordinate system.
And in that coordinate system, the partial differential
equation is separable.
And so, if you do it in that coordinate system,
you can actually solve analytical solution,
which is pretty unusual for the Schroedinger equation.
But it's a pretty goofball coordinate system
and you really have to watch out, when you convert,
to make sure you do it just right.
But it's just the rules, chain rule and the rules
about how to keep track of what's being held constant.
And when you change, will tell them, how,
being held constant, how to change it.
Now in principle, this problem is
like no different than the ODE BBP
problems we were just doing.
You just have a few more terms.
And so, you can use finite element.
You can use finite differences.
You can use basis set methods.
And it's really no different.
So all those methods basically look the same.
You just get equations that have more unknowns in them.
Because, for example, if you're using finite element methods
or basis set methods with a local basis set,
you need to put points both in the x direction and the y
direction, and in this case, the z direction,
and then maybe also in the time direction.
So now you have a lot of mesh points
in all different directions.
And so you get a lot of mesh points.
And so fundamentally, there's no problem,
but in practice, there's a big problem
if the number of unknowns gets to be so large.
That's the main problem.
So think about it.
You have some kind of equation like this.
It's in three spatial dimensions and in time.
How many points do you think you need
to discretize your solution in the spatial coordinate?
There might be somewhere, might be 100, I don't know,
points might be enough.
And so now I have x and y and z.
So I'll need, well, just in COMSOL, you saw this on Monday.
Suppose I just have two dimensions.
Suppose I just have x and z.
And I start putting points in.
Well, I want to have some points along this way.
And I want to have some points along this way.
And so, I actually have an unknown here and an unknown
here, and one there and one there,
and one there and one there.
There's quite a few of them.
I think of how many points I have.
Each of those is a point where I have a little basis function.
And each one of those has a magnitude of how much weight
I have on that basis function, things
they called the d's before in the basis functions.
And I have a lot of them.
Now, it's going to be, if I have 100 points this way
and 100 points this way, I have 10,000 of these little points.
That means I have 10,000 unknowns I have to solve for.
10,000 unknowns are a lot harder to solve for than 100 unknowns.
and if you saw, on some of the problems we did with the ODE
BBP problems, even with 100 unknowns,
we were having a lot of trouble, sometimes,
to get the real solution.
All right?
So you can just see how it just gets to be, like,
a big numerical problem.
So that's two dimensions.
The third dimension, you'll get another 100 times
as many points.
You might have a million unknowns.
And once you get up to a million unknowns,
we're starting to talk serious math, now.
OK?
Even for your good computers that you have.
In addition, it's not just you have the number of mesh points,
but it's the number of mesh points
times the number of variables at each mesh point.
So in the Navier-Stokes equations, what variables
do you have?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Louder.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Velocity and pressure, OK.
So you have pressure.
We have the different components of velocity.
What else we got?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Sorry, what?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yeah, though, yeah.
Density or temperature.
Yeah, temperature maybe, so.
So maybe temperature would be what you might use.
And then you'd have a state function maybe with these guys.
What else would you have?
So then you'd have all these species, right,?
The mesh fractions are all the species.
So you have, like, y, species one.
Mesh fraction, species two.
However many you got.
So that's a pretty big set of variables
and you need to know all these numbers at each point.
So this might be--
how many have we've got?
One, two, three, four, five, plus however many species
we have?
So this is n species plus 5 times the mesh.
And I told you the mesh was what, a million?
So maybe I have 10 million unknowns.
So the issue, to a large extent, has to do with just,
can I solve it?
Can I solve a system of 10 million unknowns?
10 million equations, 10 million unknowns.
All right, so that's going to be one big issue.
Now the other issues, we still have the same problems
we had in ordinary differential equations,
that we have to worry about numerical stability.
We have to worry about actually, is the whole problem
even stable to begin with?
Physically?
Because if it's not stable physically, then probably
we're not going to get a stable solution.
So you have all those kinds of problems
that we had before in the ODE case,
and now it's just sort of amplified by the fact,
now we have boundary conditions in more than one dimension.
So we have multiple chances to screw it up.
And we're trying to figure out how to set the boundary
conditions correctly.
And we have to worry about stability in all
the different directions.
So those are the different kinds of things
we have to worry about.
But it's all the same problems you had before, just amplified.
Now, because of this problem of the very large number
of unknowns, currently, if you have
a problem that has two dimensions,
you can kind of very routinely solve it.
And things like COMSOL on your laptop,
you can usually do those.
When you get up to three dimensions,
you might be able to solve it or might not be able to solve it.
And it could be, like, a really big challenge sometimes
to solve it.
And then in reality, we have problems like the Navier-Stokes
is four dimensions, right?
There's time, as well.
So then that gets to be really hard.
And then the Schrodinger equation
has three times the number of electrons dimensions in it.
So that gets impossibly hard.
And so, specialized methods are developed
to handle those systems with a lot of dimensionality.
So, like, for the Schrodinger equation,
people almost always do it with basis set methods.
And there's been a huge effort, over decades,
to figure out really clever basis sets that
are really close to the real solutions,
so that you don't need to use too many basis functions to get
the good solution.
So that's the way to keep the number of variables down.
In the Navier-Stokes world, I'll show you
what people do to try to deal with that.
There's many other problems.
I don't know all the problems in the world.
But in general, when you have three or more dimensions
in your PDE system, you're looking up special tricks.
You're looking up, how do people in this field
deal with this particular PDE?
Special ways that are good for that kind of PDE.
All right, so I said, one problem
was that you could have, that the PDE system, itself,
can be intrinsically unstable.
So one example I know like that is, detonation problems.
So if I have a system that has some explosive in it,
or a mixture of hydrogen and oxygen even,
and I have it in a tube or something,
if I change my initial conditions,
I can go from the case where it's stable, that it just
sits there, a hydrogen and oxygen.
Or it could turn into a regular burning.
Or if I change the conditions a little bit differently,
it could change into a detonation, where it actually
sends a shockwave down the tube faster
than the speed of sound, totally different situation.
So those kind of systems can be really
sensitive to the physical initial conditions.
A very small spark can make a really big difference
in a system like that.
And also, in your numerical solution,
if the thing is physically sensitive,
it'll typically be sensitive also to numerical noise.
If you get numerical noise at some point, that
might be enough to kick it over from, say, not
igniting at all to having a detonation,
which is a completely different solution.
But another case, like in biology,
if you're trying to study the growth of a cell culture
or the growth of a tumor, as you know,
if somebody has a tumor and one of the cancer cells mutates
and does something different, it can
have a really different outcome of what
happens to the patient, what happens to the tumor.
The whole morphology can change completely.
Everything about it can change.
So that's the kind of system, also, that
can be very sensitive to small details.
A small change in one cell, for example,
could be a really big difference.
So this is not just in combustion problems
you have these kinds of instabilities.
All kinds of problems have this kind of thing.
You can have just regular chemical reactors,
like in multiple steady states, little tiny fluctuations
can make it jump from one steady state to another one.
That's another very famous kind of problem.
There's another class of PDEs where
it's stable in some sense, but it causes a problem.
We call those hyperbolic.
And those are like wave equations.
The solutions are propagating waves.
So the equations of acoustics, the equations
of electromagnetism, are wave equations
where, basically, a wave propagates more or less
without any damping.
So what happens, then, is you have some numerical noise that
will make a little wavelet that will propagate, more or less,
without any damping, too.
And it'll keep bouncing back and forth inside your domain
numerically.
And if you keep introducing numerical noise,
eventually the whole thing is just full of noise.
And it doesn't have much relation
to the real physical solution.
So hyperbolic systems of differential equations
need special solution methods.
That in the numerical method, it's
carefully trying to damp out that kind of propagating wave
that's coming from the noise.
And they set them up a special way.
I'm not going to talk about how you solve them in this class,
but there's a special, whole group of solution methods
that people use who are solving acoustics equations
and solving electromagnetic equations.
If you get into solving shockwave equations--
I actually work for Professor [? Besant. ?]
He has, like, electrochemical shocks in some of his systems.
You have the same kind of thing.
You have to use special mathematics to correctly handle
that, special numerical tricks.
But we're not going to get into that in this class.
If you get into those things, you
should take the PDE class for hyperbolic equations
and they'll tell you just all about the solvers
for those kinds of things and special tricks.
And there's journal papers all about it.
Here, what we're going to focus primarily
on two kinds of problems: elliptic problems
and parabolic problems.
And those problems are ones where
there's some dissipation that's very significant.
And so, in an elliptic problem, you
introduce some numerical noise and as you solve it,
the numerical noise kind of goes away.
Sort of like intrinsically stable all throughout.
Now the real definition of these things,
parabolic, hyperbolic, elliptic, has
to do with flow of information.
So in a hyperbolic system, information--
there's regions of the domain that you can make a change here
and it doesn't make any effect on some domain part over there.
And so that causes-- because that
is the real situation, that also propagates
into the numerical solution method.
So for example, if I'm modeling a shockwave that's
moving faster than the speed of sound,
I can introduce any pressure fluctuation
I want behind the shockwave, and it
has no effect on the shockwave.
Because the pressure waves from that pressure fluctuation
are traveling at the speed of sound.
They won't catch up to the shock.
So they will have no effect whatsoever.
So if I was judging what was happening
by looking at what's happening in the shock,
I can introduce any kind of random noise I want back over,
downstream of the shock, I guess.
Yeah?
And it won't make any difference at all.
And so, that's part of the reason why I
need a special kind of solver.
If I have elliptic equations, every point
affects every other point.
A famous case for that is, like, this steady state heat
equation.
You have some temperature source here, a heat source here,
and you have some coolness source over here.
And there's some heat flowing from the hot to the cold.
And the heat flow is kind of slow enough
that everything sort of affects everything else, since you
actually have, the heat is trying to flow by,
sort of, every path.
And there's some insulation that's resisting its flow.
And it has been in the connectivity
of flow a certain way.
And those kinds of equations are called elliptic.
And the methods that we use for the relaxation mesh
that we showed are really good for those kinds of methods,
for those kind of problems.
And then another kind of problem that we get into a lot
is like this one.
This is called parabolic.
And typically, in the cases we see that are parabolic,
we have time as part of the problem.
Always what we think should happen
is that what happens in the future
depends on what happened in the past, but not vice versa.
So you would expect that when you're
computing what's happening at the next time point,
it's going to depend on what happened at earlier time
points.
But if you try to do it the other way around,
it would be kind of weird, right?
You wouldn't expect that what in the future
really affected stuff in the past.
So it naturally leads us to sort of pose it as an ODE IVP,
or IVP kind of problem.
And we know initially, some time zero something,
and we're trying to calculate what
happens at some later time.
And so those kinds of problems we call parabolic.
And the methods that we're going to talk about mostly
are ones for solving parabolic and elliptic problems.
And you can read, in the middle of chapter 6,
has a nice discussion.
Just has the mathematical definitions
of what these things are.
Now, even in a steady state problem,
the problem can have a sort of directionality to it.
So if you have a problem that's very convective.
So you have a flow in a pipe.
And you have a nice, fast velocity flow down the pipe.
What's happening downwind is very much affected
by what happened upwind, but not vice versa.
So maybe you're doing some reaction upwind.
Say you're burning some hydrogen.
You're going to get water vapor made upwind
and it's going to flow downwind.
And you're going to have steam downwind.
If I squirt a little methane in downwind,
it doesn't actually do anything to what happened upwind.
I still have a hydrogen and oxygen flame up there
and it's been burning and making steam.
Do you see what I mean?
So it has a directionality, even though I
might do it as a steady state problem and not have time in it
at all.
And you've probably already done this, where you can
convert, like, a flow reactor.
You can write it, sort of, as a time thing or as a space thing.
Time and space are so related by the velocity of the flow.
We'll have similar kinds of phenomenon then.
It's almost like a parabolic time system
that, if the flow is fast enough,
the fusion backup, anything moving upstream is negligible.
It's all moving from upstream to downstream.
And so then, you'll have the same kind of issues.
In fact, you may even try to solve it
by just starting it at the upstream end and competing
stuff and then propagating what that output from that first one
does to the next one downstream, and then
do them, one at a time.
That would be one possible way to try
to solve those kinds of problems that are highly convective.
As you slow the velocity down, then the diffusion
will start to fight the velocity more and more.
If you make the velocity very slow,
the Peclet number very low, then you really
have to worry about things diffusing back upstream.
And then it's a different situation.
And that would be more like the elliptic problems we have.
Because things downstream would actually affect upstream,
if the flow is really small.
But typically for our problems, the velocity
is always really high.
Peclet numbers are really large and so not much goes
from downstream to upstream.
It's almost all upstream to downstream.
This leads to a famous situation.
If you look in the textbook at figures 6.7, 6.8,
it's kind of a very famous problem.
That if you just do this problem.
So this is just an ODE problem.
So I just have convection and diffusion, no reaction,
very simple.
This problem.
You'd think that I could just use my finite differences,
where I would put, I could use the center difference
for here and here and change this into this obvious looking
thing.
You think that that would be a reasonable thing to do, right?
And you can do the same one here.
You can write--
I'm terrible with minuses, sorry.
So these, the finite difference formulas.
And so, you can write these like this.
This stays equal to zero.
And now it's just an algebraic problem to solve.
A system of equations, you have equations
like this for every value of n.
And there's the mesh points.
You guys have done this before, yes?
Yes.
OK, so famous problem is, if you actually
try to solve it this way, unless you make delta z really tiny,
you get crazy stuff.
So on figure 6.7, they actually show
what happens for different values of delta z, which
is different values of what they call the local Peclet number.
So there's a thing called local Peclet number, which is--
OK?
If you make the local Peclet number too large, actually
anywhere bigger than 2, and you try
to solve this equation, what you get
is oscillations, unphysical oscillations.
They'll make the phi go negative.
It's crazy.
It looks like the simplest equation in the world, right?
It's a linear equation, so what's the problem with this?
But it doesn't work.
And so, then you could think, well, why doesn't it work?
And there's, like, multiple ways to explain this
about why this doesn't work.
I think that the best way to think about it
is about the information flow.
So if the local Peclet number is large,
that means that phi is just flowing downwind
and diffusion is not doing much.
Because the Peclet number is big.
And so, you really shouldn't use this formula
to calculate the flow.
Because what happens downwind, at point n plus one, I guess.
Here we have the flow this way.
So here's n minus one.
Here's my point n, n plus 1.
I really should not include, is this right,
n plus 1 in this formula.
Hope I have those signs right.
I apologize if I had them backwards.
Because I don't think that what happens
downstream should really affect this at all.
So really, I would do better to change from this center
difference formula to what's called the upwind difference
formula.
Where they would say, OK, instead,
let's use phi of n minus, the phi of n minus 1 over delta z.
Just use that instead of using this formula.
Now, you wouldn't think that that
would make much difference.
But it makes a gigantic difference.
And so if you look at figure 6.8 in the textbook,
you see you get perfectly stable solutions when you do that.
But where you do this one, you get
crazy, oscillatory, unphysical solutions,
unless you choose delta z really tiny.
Now, there's other ways to look at this.
So I think this is the correct way, the best
way to think about it is, it has to do with information flow.
If you try to make your solution depend on stuff
that it doesn't really depend on.
So upwind does not really depend on downwind.
If you force it to, by choosing to include
that information here, you end up with crazy stuff.
So that's one way to look at it.
Another way to look at it is that, if you really
want to compute these derivatives,
you've got to keep delta z small.
Because we're taking a differential.
We're turning it into a finite difference.
And if you take your delta z too big,
then it's a terrible approximation to do this.
And you can look at it.
You can see that what we're doing
is, when we have this problem where the local Peclet
number is getting too large, is we're actually choosing
delta z bigger than the, sort of,
the thickness of the solution.
So if you look at the solution, the analytical solution
of this problem, it's like this.
Where this has a very, very thin layer at the end.
It's like the flow has pushed all your stuff to one side.
And if you choose your delta z to be from here to here,
that's how big your delta z is, and then you're
computing the derivative of this point halfway along from here
to where this is.
You compute your derivative like that.
It's not a very good, not a very accurate representation of what
the real derivative is here.
And in fact, you really should have
a bunch of points here and here and here and here
and here to really resolve the sharp edge.
And so you did something really crazy
to use a big value of delta z to start with.
So that's another way to look at this problem.
Now, the fix, by doing this upwind differencing,
I think the best way to look at this is saying,
well, I'm just going to make sure I only depend,
I only make the equations depend on what's upwind,
with the convective term.
Because there's nothing convecting from downwind.
So just leave that out.
So that's a conceptual way to think about it.
Another way to think about it, which is in the text,
is that by doing this, you're introducing
another thing called numerical diffusion, where
you're actually, effectively increasing the diffusivity.
Because this is a very poor approximation
to the differential.
This is a asymmetrical finite difference formula.
It's not a very good value estimate of the derivative,
if delta z is big.
But you can write it out carefully and write this out
and say, oh, this actually is sort of
like this plus an effective diffusivity
chosen very carefully, so the terms cancel out just right.
And it turns out, if you look at it
that way, the effective diffusion you've
added, by using this instead of that,
is just enough to make the local Peclet number,
if you compute it with the effects of this d effective,
to stay less than 2.
Which is what you need as a condition to make
this numerically stable.
So I strongly suggest you read that.
It's only two pages long in the textbook.
It's very clear.
Just read, just look at it if you haven't seen this before.
There's a similar thing and it might be relevant to the COMSOL
problem that Kristyn showed you on Monday, is you can--
In that case, there was a region of the boundary here
that had some concentration c naught,
and then there was, over here, the concentration was zero.
The boundary is zero.
Here, it's some number, 17, whatever number it was.
And if you think of how to model what's going on,
one way to look at it is, I have a diffusive flux coming in
from the drug patch, diffusing the drug into the flow.
And I have this big, giant flow, flowing along here.
And the flow gets slower as I get closer to the wall because
of the friction with the wall.
And so if I looked somewhere right in here,
I could just draw a little control volume,
and look what's happening.
I have a diffusive flux coming in here.
I don't have any drug up here, so there's
no drug coming in here.
But I have drug leaving.
So basically, it's coming up and it's making a right-hand turn.
And so then I could think, well, if I
was going to try to figure out, what's
the steady state concentration of drug in there, I could say,
well, there's a certain amount of drug
entering from the diffusive flux, which
would be just this times dcdy.
This is the y direction.
And then this part over here would
be the velocity, actually, just times the concentration
I have in here, right?
So how much is flowing out.
And it's probably, my units are screwed up so there's--
Well, maybe not.
That's OK.
Is that all right?
Units are screwed up.
You guys will get it right.
It's probably an area.
Something to do with this size right here.
This one here, too.
Is that right?
So you could compute what the steady state
concentration would be if it's just coming in and flowing out.
And so, this is the finite volume view
of this kind of problem.
And so you can write the equations that way, instead.
And what you really care about are what the velocity flows
and the diffusive fluxes are on these boundaries
around the point.
So you don't try to compute things right at the point.
You compute around it.
And one way to look at this is, we just
did finite differencing.
We were looking at what was happening coming in.
And then, this is actually a center difference
around this interface, which is halfway between the two points.
And so actually, this is a good approximation
for the finite volume point of view.
This is all right?
Now, these are all different ways
to write down the equations.
All of them are correct in the limit
that delta z goes to zero.
If you have an infinite number of mesh points,
all these are exactly the same.
But they lead to really different numerical properties
when delta z is large.
And they're all inaccurate when delta z is large,
so it's all about what happens with the inaccuracies.
Now, we can go ahead and take any problem, even
really complicated 3-D problems or 4-D problems like this,
and we can write them with relaxation methods.
And what we're basically doing in a problem
like this is turning it into some function of y
is equal to zero, where y is the value of all the unknowns.
And there's a lot of equations.
So I might have 100 million equations
and 100 million unknowns that are the values,
or every state variable at every state point at all times.
And I can write it down, by finding differences,
by finding elements, by finding volumes.
Any way I want, I can write down an expression like this.
And I know in the limit, as I make
the number of elements of y go to infinity,
it will actually be the real solution.
And the problem is, my computer can't handle
infinite number of unknowns.
In fact, it's even going to have trouble
with a 100 million unknowns.
And so, I have to figure out what to do.
So how would I solve this normally,
is I would take the Jacobian of this and I would do,
I'd take an initial guess and I'd update it and have
the Jacobian times my change from the initial guess
to the equal negative of f, right?
You guys have done this a million times.
As Newton-Raphson is how you'd find
improved from an initial guess.
You have to provide some initial guess.
This is actually pretty hard if you have 100 million unknowns.
You have to provide the value of the initial guess for every one
of the 100 million.
But somehow you did it.
And now you want to refine that guess and this is
the formula for it.
And so, you just use backslash, right?
But the problem here is that now f is a 100 million long vector
and this is a 100 million squared matrix.
And a 100 million squared is pretty big,
10 to 16th elements in it.
So I'm not going to just write that matrix
down directly like this.
Now, we cleverly chose local basis functions,
so this is going to be sparse, this Jacobian.
So most of the numbers, are those 10 to the 16th numbers
in this matrix are zero.
So we don't have to represent them.
But there still might be quite a few.
Probably the diagonal events are non-zero,
so that's like 10 to the 8th of them right there.
And we've got, probably, a couple of other bands.
So I don't know, might be 10 to the 9th non-zero elements
inside this thing, which is pretty nutty.
And depending on how much memory they provided in the laptop
they gave you, you might be filling up the memory
already just to write down j.
And certainly, if you tried to do Gaussian elimination
to solve this, you're going to have a problem because you
get what's called fill in.
Do you guys remember this?
Even if you start from a very sparse Jacobian,
as you run the Gaussian elimination steps
the number of non-zero entries in the intermediate matrices
gets larger and larger.
Remember this?
And so, even if, initially, you can write down j and store it
in your computer, you could easily
overwhelm it with the third iteration or something.
So you're not going to just use Gaussian elimination.
So what do you have to do?
You want to find a method.
They're called direct.
Direct methods to solve this kind of problem,
in order to get your increment, which
is going to be your update, to get a better estimate of what
the solution is.
And the direct methods are ones that you don't actually
store the gigantic matrix.
So you're trading off CPU time for memory.
So you don't have enough memory to store the whole thing.
Which we, normally you would do this with Gaussian elimination
and huge steps.
It would solve it perfectly, right?
Right, Gaussian elimination's great, backslash.
You guys have used it a few times.
It's a really nice program.
It always works.
It's really nice.
But you're going to give up that because you can't afford it
because the memory, it's consuming too much memory
and your cheapo department didn't buy enough
RAM in your computer for you.
So you're going to instead try a direct method which
is, so it's trading off CPU time versus RAM.
So you're going to reduce how much memory you actually
actively have, but you're going to expect that it's
going to cost more CPU time.
Because otherwise, everybody would do something else besides
Gaussian elimination all the time, and they don't.
All right, so what's the method we're going to use?
It turns out that variance on the conjugate gradient method
turned out really well.
So we mentioned this earlier, briefly.
What you do is, instead of trying to solve this,
you try to solve this.
So you just try to minimize that, right?
You know at the solution this is going to be zero.
Right, jy plus [INAUDIBLE] zero when this is solved.
So you try and vary and find the delta y
that makes this go to zero by minimizing it.
And then this method, it turns out
that the conjugate gradient, if everything works great,
this is guaranteed in n iterations
to go directly to the solution of this kind of problem.
Now n is large, now, because we have 10 to the 8th unknowns.
So that's 10 to the 8th iterations.
You do anything in 10 to the 8th iterations,
you're not really sure it's really going to work.
So that's a warning.
But at least it should get, even after a few iterations,
you should get a smaller value of this than you had before.
And so this is the method.
And what's really nice about this method,
if you look into the details of how it works,
it never has to store any intermediate matrices.
So all it needs is the capability
to multiply your matrix times some vector.
And from that multiplication, it figures out a step direction.
That's the trick of this method.
And because it only has to do a forward multiplication,
you don't have to actually store j.
You can just compute elements of j, as needed, in order
to evaluate this top product.
This j times v. And you don't have to ever store
the whole j at once.
So it's really great for RAM to do this method.
Now, as I said, when you do 10 to the 8th iterations,
things don't always work out so well.
And so people have found that when you go more than maybe
10 or 20 iterations like this, you tend
to pick up numerical problems.
And so, they've worked out better methods.
And there's one that people use a lot for these problems
called bi cg stab.
Which is biconjugate gradient stabilized.
The applied mathematicians went crazy trying to figure out
better methods like this.
And this is the one that, I think currently,
people think is the best.
Though probably, there's a lot of research on this,
so maybe there's even better ones now.
But inside Matlab, I think this is the best on they have.
And this is a method for doing this.
Now, it's iterative.
Now, let's remember what we're doing.
We're trying to solve this problem.
We started with a guess.
We're going to break it up into an iterative problem like this,
where we're going to do some steps, delta y.
This is now doing iterative procedure
in order to compute delta y.
So we have an interim procedure and another iterative procedure
inside the first iterative procedure.
So this is going to be a lot more CPU time.
So just to warn you.
And it's also, it's not guaranteed
like Gaussian elimination, to beautifully come
to machine precision solution.
It's going to get you something.
But anyway, this is what people do.
So that's one good thing to try.
Now, this whole approach is sort of based on Newton-Raphson.
We know that doesn't have the greatest radius of convergence.
So you need to provide a pretty darn good initial guess.
So how are we going to get a good initial guess if you
have a 100 million unknowns?
So what methods do we know?
So one idea is you could do things
like homotopy, stuff like that.
If you could somehow take your PDE system,
get rid of the nonlinear terms.
You're really good at solving linear systems of PDEs.
You start with that and then you gradually
turn the non-linear term on.
Maybe you could coax it over.
So that's one possibility, if you're really good at PDEs.
Another possibility that I've run
into a lot is the problem you're trying to solve, for example,
might be a steady state problem like this.
Where this is zero and you try to figure out what
the steady state situation is in a flow reactor, for example.
And so, it's actually, your steady state problem
came from a time dependent problem.
And if you think that your time dependent problem really
converges to the steady state problem,
sort of without much concern about what
the initial guess is, then you could
put any random initial guess in and march it along
in time enough, and eventually, it should end up
at the solution you want.
So that's the idea.
So that's called a time marching idea.
So there's one idea, sort of the homotopy continuation
kind of idea.
That's one idea about what to do for initial guesses.
Another idea is the time marching idea.
And the key about it is, you have to really believe
that, physically, the system does march
to the steady state you want.
If it does, then this is a good idea.
If it doesn't, then you're not going to end up where you want,
because it's not going where you want.
And so, you have to watch out.
But in some situations, this is true, that things work.
So I do a lot of combustion problems.
A lot of these things, you light a spark on a stove burner.
No matter how you light it, it ends up
basically the same flame.
So that's a good one.
But I have some other ones where I light a candle
and I'm in a convective field, a turbulent field,
and it flickers.
And no matter what, it always flickers.
And I never get a steady state, so I'm never
going to find a steady state solution with that one.
And so I could time march until I retire.
My computer's still burning CPU time all that time
and it'll never get to a steady state solution.
So you really have to know what the real situation is.
But if you have one where it's going to converge,
then this is a reasonable idea.
Now let's talk about the time marching.
The time marching, you have, say, dy,
dt is equal to some f of y where this
is the discretized spatial vergence.
I took all the space dimensions and replaced them
with the finite difference expressions
for the derivatives.
Or I did colocation or something to convert this
into algebraic equations.
So the right-hand side is now algebraic equations, no
longer differentials.
And now, I'm going to time march this and I just--
This is very long.
It's a 100 million.
So how are we going to time march that?
Well there's kind of two schools of thought for this.
One school of thought is, if I can use an explicit method,
an explicit ODE solver, those are pretty nice.
Because all I have to do is evaluate f and then
I can march along and I never have to save any matrices.
Because remember, look at those formulas, y,
y new, by the explicit method, is
equal to some function g of y old.
And this is the update formula.
It can be forward Euler.
It could be RK 2.
It could be whatever.
But there's just some explicit formula.
I plug in the old vector.
I compute a new vector.
I never had to do any inversions of matrices or anything.
So that's good.
And this is, indeed, the way that the Navier-Stokes
equations are solved currently by the best possible solver.
So the best solver in the world discretizes in space
and then does ODE IVP using explicit method
for marching it.
They do problems that are pretty big.
So they'll have like 100 million mesh points and then 100
species.
So they have 10 to the 10th unknowns.
And at each time step, you're computing another vector
that's 10 to the 10th long.
So it's a pretty big problem.
And they use, like, 40% of one of the national supercomputer
centers will be running at one time for one job like this.
But it works.
You never have to store anything that's too big.
The biggest thing you have to store are the vectors.
Which, you're always going to have
to store some object the size of your solution,
right, if you're going to get your solution.
So that's one possibility.
Now, this is really good if an explicit solver can solve it.
So that means it's just not stiff, your problem.
And also, if you want a time accurate solution.
So in that case, they actually want the time accurate
solution.
But suppose we're trying to really solve the steady state
problem over here.
Then we really don't care about the time accuracy.
We're not going to report anything about our time steps.
In the end, we're only going to report our steady state
solution.
That's all we cared about.
So in that case, there's no reason
to try to have a really high accuracy, explicit formula
and try to make sure we all our time steps are exactly right.
We're just going to throw away all those time points anyway.
We'll just keep the final one as our initial guess
for a Newton-Raphson solve to find the steady state.
So in those cases, you might instead want to use--
So this is time accurate.
There's another method.
I don't know if I should call it time inaccurate.
If you don't care, if the solution is really,
what the solution yft is, all you
care about is where you get to, then you can do other methods.
And one of them is the backward Euler.
So you can say, well, y new is equal to y old
plus delta t times f of y new.
Now, the advantage of this is, I can make delta t pretty large.
This is guaranteed stable as long
as the ODE system is stable.
Remember?
And so, this will go to this, go to the solution pretty well,
to the steady state solution.
It won't go there in the right amount of time
because this is a very low order of formula.
It's not really accurate.
But it will end up to the real steady state solution.
This is what's actually used a lot.
However, the problem with this is, this is now
an implicit equation.
So we may have to solve it with some method
like Newton-Raphson.
But we got into this because we couldn't
solve our Newton-Raphson steps because we didn't
have a good initial guess.
So then the question is, can we find a good initial guess
for this problem?
And in this case, you can because you
can use the explicit formulas to again provide the initial guess
for this implicit formula.
So that's what they do, also.
And so now, if we're doing an iterative procedure in order
to get the initial guess for our other iterative procedure
ever there, which we're going to break up
into another iterative procedure inside it to solve every step.
But this is currently what people do.
So you guys are smart.
Maybe you can figure a better way to do it.
The world is waiting.
I think I'll stop there.
By the way, actually just names--
This kind of thing is called method of lines.
And if so, if anybody ever says I solved something
by method of lines, what they meant was,
they discretized in some of the coordinates
and they kept one of the coordinates as a differential.
And then they solved it with an ODE solver at the end.
And in certain systems like this, it's a smart thing to do.
You need to be kind of lucky that the set up is right,
so you can use it.
But if you can, it can be pretty good.
All right.
See you on Friday.