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.
WILLIAM GREEN: So welcome to our special class with half the
people gone for tonight's hearing.
You guys are the chosen few, so I'm pleased you're here.
AUDIENCE: [INTERPOSING VOICES]
WILLIAM GREEN: So I'm just going to say a little bit
about models versus data and then
I'm going to go back to undervalue problems and talk
a little bit about how you handle really big ones.
Yeah?
AUDIENCE: Are these going to be [INAUDIBLE]..
WILLIAM GREEN: Yeah.
Yes, sir.
AUDIENCE: All right.
WILLIAM GREEN: And in fact, these
are even going to be on the video,
because the guy just this fixed the video,
so it will be recorded.
All right.
Yeah, I'm sorry.
You said the video was broken last time
so that recording didn't work.
All right.
So when we're talking about models versus data
a really critical question is whether the model
is really consistent with the data or not.
And how do we know this is really true?
So we talked about how we can figure out
some sort of chi-squared maximum that we would tolerate.
So the maximum that we tolerate, we
call it Q, capital Q in a lecture before.
And so if we can find a chi-squared of theta that's
less than this, then it's at least plausible
that our model is true.
Right?
Or, alternatively, if all chi-squared thetas
for any value of thetas are greater than chi-squared max,
then we would say the model is false,
or something else is really wrong with our experiment.
Right?
That's what we say when we're saying
we have a maximum chi-squared that we would tolerate.
So a critical question with us is
whether we really found the very best fit
value of the parameters.
In some simple problems you might have linear dependence
on the parameters, and then there's only one
best fit usually, as Long as you have
a non-singular set of equations and everything
we find there's no trouble.
However, more commonly you have equations
which are singular or close to singular,
and also you normally have non-linear dependence
between the model predictions and the parameters.
So, for example, I do kinetics.
We always care about the ray coefficients, the K's and they
always come in inside a differential equation.
Right?
You always have like d concentrations
dt is equal to some terms, and some of them
are, you know, k concentration I times concentration j.
whatever one it is, ij.
Right?
You have a lot of like terms like this by molecular reaction
terms.
And these K's, some of them I don't know.
And so these would be one of my thetas.
So this would be, you know, theta 17
would be one of these K's that I want to determine.
And when I have a differential equation like this,
when I integrate it, the K is always
going to come in non-linearly in the solution.
Right?
So when I integrate this solution,
even in the simplest case you get things like c of t
is equal to ae to the negative kt.
That might be a really super simple case
if you have a linear kinetics.
The K is still coming in non-linearly
to the predictions.
Right?
And if you have anything complicated
you can't even write it down in a closed form expression,
and you have to integrate with oed45 or ode15s
and so the relationship between K
and the concentrations that you measure
can be really complicated.
So it's always non-linear.
In my world all the models are non-linear.
And I think there's many other situations like this.
So if you try to measure a diffusivity,
it might be possible to set up an experiment
where you'd expect a linear dependency to something
you'd measure in a diffusivity, but most likely not.
So most likely it will come in non-linearly somehow
and you'd have some model for it,
and you'd have a non-linear dependence
between that number and the numbers you would measure.
Right?
And almost all the properties are like that.
So because you have a non-linear dependence,
then you have to worry about how many minima
do you have in your surface?
So chi-squared of theta is a non-linear function
of a lot of variables.
And typically it won't just have one minimum,
you'll have a whole lot of local minima.
Now one of those is the global minimum.
And that's the one you want, that's
the very best possible fit.
But you normally won't know how many minima it has
and so no matter how many minima you
find you don't know if you're done,
and you don't know if you really found the best fit.
OK?
So just because all the ones you found so far bigger
than chi-squared max you might not
be ready to write that article in Nature
to claim that the laws of physics
are incorrect because you're not sure you really
found the very best fit.
OK?
So that's one of the many problems
of trying to figure out if a model of the data
are really consistent.
Now that particular one, you can overcome
if you can find a guaranteed way to find the global optimum,
to find the very best possible that there is.
And some people spend their life working on this,
and one of those people is professor Barton
in this department.
And so he teaches a class on how to do this,
and there are actually methods to do it.
So if you want to see the application of this
to kinetic models there's a paper
by me and professor Barton.
And the lead author is one of his former students,
Adam Singer and J. Phys.
Chem. and it just shows how you can use global optimization
to guarantee you have the very best possible
via the parameters, and then they'll show you what happens.
So this is a case where the experimental data was measured
by one of my students, James Taylor, and that's
the thing with the error bars, and he did a lot of hard work
to make sure he had pretty good error bars
and he really knew what they were.
And so he had those points.
So he measured a lot of points.
This measuring a lot of points is a key thing.
So it means his degrees of freedom is very large.
Because he uses many, many different time
points along the way there and he
has error bars on all of them.
And so because the degrees of freedom are very large
it actually makes it a very good constraint on the model.
It tested very thoroughly.
Where if you only measured one point,
then that's not a very good constraint.
Now what James had done, and we published a paper in 2004 where
he did it, where he said, you know, when I do my fits
this is what I get.
I get these kind of purple and green things,
and therefore they're all outside the error bars,
and therefore I publish in the Journal of Physical Chemistry,
I claim this model is wrong based on my experimental data.
OK?
So we published that in 2004, it was accepted by the referees.
It's in the literature, you can read it.
Don't believe it though, because Adam singer went back and did
global optimisation using the same model
and his red line goes right through all the points really
well.
Do you see that?
You can see the red line there.
See this red curve?
That's the true best fit.
And so all the fits that James had found
were all local minima in the chi-squared squared surface.
And when he actually found the best fit
with Adam's fancy method, then it
turns out that actually the model
matched the data perfectly.
So that's definitely true.
So what we published the first time
was completely wrong, because we said
the model disproved the data.
But in fact the model proved the data.
I mean the data proved the model,
or was completely consistent with the model
in every way in great detail.
So I was a little embarrassed.
So then James, also in the same paper, had found--
and they measured another condition
in a slightly different case, and actually
was looking at a different solvent.
And in this case he found this green curve as his best fit,
and it looked pretty good by eye,
though you can see it sort of skittering
on the tops of the error bar bands.
And so this is one where he said, you know,
we do the statistical test, there's like a 5% chance
that you would measure data like he
measured if the green model was the truth.
And so we said some wishy-washy words in the paper,
basically because we weren't really
if it was consistent or not consistent with the data.
But in that case Adam used his fancy global optimization
techniques and found that James' fit was actually not
that far from the global optimum,
and he got the red curve there, marked by the red arrow,
and that one looks pretty good, it actually
goes through the data.
But it doesn't really do it right,
and found that the confidence intervals are only 16%.
So it's only a 16% chance, one chance out of six,
that James would have measured the data he measured
if that model was true.
And now we're sure that we have the global best fit,
so now we're, like, thinking well,
you know, maybe we should say that the model is not
true for this case.
Because it's an 84% chance it's not true anyway.
So actually we got it wrong both ways.
So in the first case we said that our data
disproved the model.
And in the second case, we said that our model was pretty good
with the data, but actually they were both wrong.
So one was wrong that actually--
the data did prove the models correct for the first case
and it actually looks a little fishy for the second case.
So anyway this is just a concrete example
to make you be a little bit more worried about what you do when
you do these fitting things.
Now there's a variety of ways to try to do the least course
fitting to try to find optima.
There was a guy named Carr from University of Minnesota.
And he published this [INAUDIBLE] way
to just try a whole lot of different initial guesses,
and that usually works to find the global minimum.
But if you really want to have a guaranteed global minimum
that you're absolutely sure is definitely
is the global minimum, then you really
need to use global optimization methods
and you really should talk to professor Barton.
And for these kinetic models he actually has
a distributed code called GDOC that will just
find the global optimization as long as you don't have
too many kinetic parameters.
I think if you go to about 6--
around 6 adjustable parameters it can handle.
If you have more than that then usually typically not.
Because the problem with global optimisation methods
is that they scale exponentially with the number
of adjustable parameters, or the number of things
you're trying to optimize.
Now on the other hand, if you have
a model that you need to use more
than six adjustable parameters to fit your data,
then you might want to try to think
of a different experiment.
Because you know the famous saying is like, you know,
give me N parameters I can fit an elephant.
So if you have more than six parameters
and you're trying to fit some data set,
probably if you do it right you'll find a fit,
but it may not really prove anything.
So from that point of view professor
Barton's code, GDOC, is pretty convenient for doing
these non-linear optimizations.
So that's a code for non-linear optimization with differential
equations embedded.
And they can be stiff and they can be large,
it doesn't matter.
What really matters it's just the number
of adjustable parameters, that's the cost.
All right now there was a really good question,
that you asked I think, in class last time, maybe.
Maybe the time before.
And it was what happens--
so let's look at different case you get.
So one case you get, you know, if chi-squared square,
theta best, your best fit parameters is a lot bigger,
think chi-squared max, then you're ready to say the model
disproves it.
So that one you're OK.
If chi-squared of your best fit is less than,
or significantly less than chi-squared max--
it has be much less than it, really less, not equal.
--then you're really happy to say
that the models are consistent.
The model and data are consistent to each other.
And, in this case, what you would normally do next
is determine the confidence intervals on the thetas
that you adjusted.
Right?
And if you read the chapter in Numerical Recipes,
for example, some of the notes online,
it could tell you how to do this.
Right?
And so that was OK?
It's OK.
So the case that's really problematic
is you find chi-squared theta best,
and it's less than chi-squared max,
so you think the model is consistent with data,
but it's pretty darn close to theta-squared max.
And so then you have to figure out, well, the weird aspect
of this is I now have--
here's my famous plot of theta one and theta two.
So I have two parameters I'm adjusting.
If they want theta two, here's my best fit.
My best fit.
If chi-squared max-- if the curve
of constant chi-squared max is something like this,
this is chi-squared of theta equals chi-squared max.
If it's like that then I have some big region
of a lot of different theta values
that give pretty good fits, and I'm
happy to draw confidence intervals,
like draw lines here and here, for example, and say
that's the confidence interval on theta one,
and I'm all right.
OK?
And it's sloped here so I might want
to tell people about the co-variance between theta
one and theta two.
All right?
But I know what to do, it's fine.
Now suppose this is the case that's like this.
Where chi-squared theta best is significantly less
than chi-squared max.
Now suppose instead chi-squared best
was very close chi-squared max.
So then here is chi-squared max.
This is chi-squared of theta equals chi-squared max,
this is the good case.
Here's the bad case.
So in the bad case, there is some region of thetas
that the model is consistent with the data,
but it's very tiny.
There's a little tiny range of theta values
that would make the models consistent with the data.
Now what's weird about this is that my best fit, it's not
very good.
The chi-squared is almost up to the maximum
chi-squared I would tolerate.
So it means my deviations between my model and data
are pretty bad.
Like this plot actually.
OK?
And so I'm, like, you know, it could-- it's possible.
I was a little bit unlucky and my data set
is pretty far off but still within the possible range.
And so maybe the model data are consistent.
But then when I try to compute the confidence intervals
on my parameters, there's only a really small range
of parameters that are going to give good fits here.
Because the best possible fit is this red line.
There's no choice of the parameters
that's going to go through all the circles with this model.
So now I'm like, what's going on here?
So if I went ahead and do the same way
I did before, I draw my confidence intervals like this
and say, OK.
Say, wow, I've done an amazingly great job
of determining theta one.
I know theta one very, very, very precisely.
And the same for theta two.
Here's theta two, it's got to be in that range
and I've really done a great job.
So I'm awesome.
But I'm only awesome because my original best
fit wasn't very good.
So this might make you feel a little sick in your stomach
if you're getting ready to publish in nature.
Right?
And then if the reviewers are on their game,
you might get something really scathing
back from the reviewers.
So there's kind of two ways to play it at this point.
One way is, you are absolutely sure that your equations are
correct.
You are Sir Isaac Newton, you have just
written the law of gravitation.
You are measuring the orbits of the moons of Jupiter,
and you are sure that your equations are right.
And so, you know, my problem is my antiquated telescopes
of 1600 are not that great, and so I
think I have maybe-- maybe I misestimated my error
bars on my measurements.
Maybe I have some aberration or something,
the angles are off, right?
My twisted telescope looks as if the star moves or something,
right?
OK?
So that's one possible way to go,
say I really believe the model, and so I
believe this is the truth, and I really
believe my chi-squared max to be what it is,
and I double checked the error bars and I think it's true.
So I think have done a great job and I should publish it
in Nature because I have determined
these parameters with 17 more decimal places
than anybody ever could before.
OK?
So that's one way to look at it.
The other way to look at it is what's
recommended in the Numerical Recipes book.
And they say don't use--
instead of writing the equation to be chi-squared of theta
has to be less than chi-squared max,
they say when you're computing the confidence intervals,
instead write chi-squared of theta
has to be less than chi-squared best fit plus chi-squared max.
So they say don't use a stinky little circle here,
draw a nice big one around there.
And what is this saying?
I guess what this is saying is I don't believe
that my poor agreement between the model and data
justifies me declaring very narrow parameter values.
So
I just use the kinds of parameter value confidence
intervals that it would have actually computed a priori when
I do my experimental design.
I think that's really how well I can determine stuff.
And I think something is wrong that my chi-squared best
fir is so bad.
Maybe I didn't do my global optimization right
and so I have the wrong parameters.
Maybe I mis-estimated my error bars
and so my chi-squared are all off.
Maybe-- I don't know what else happened.
Something happened, and I don't believe I'm that good,
so therefore I'm going to quote much broader confidence
intervals, which means much less precision in my determination
of parameters.
And that's what the Numerial Recipe
guys recommend you should do.
But you know, it's up to you.
And I see it both ways.
So the people who won the Nobel Prize for the dark energy,
do you guys know about this?
So if you ever saw the data there.
So that's one where they said--
they were using Einstein's theory
of relativity, General relativity.
They say we believe Einstein, we believe that equations true.
And we're not going to, you know, the fact
that the data only matches when you choose
a certain value of the cosmological constant,
we're just going to go with it.
And we're reporting that's the value
of the cosmological constant.
Now some other people might say that's kind of goofy theory.
Einstein changed his mind six times about
whether to put the cosmological constant
in the equation at all, so I'm not so confident
and I would report something different.
But they went with their gut.
They said we believe Einstein.
They went with it.
The Nobel Committee agreed, and they got the Nobel prize.
So you know, you're taking your chances.
The Numerical Recipe guys are more like engineers,
they say, oh, we don't believe models that much anyway.
Let's just use more conservative confidence intervals.
Any questions about this?
OK.
All right.
So now let's change topics.
And we'll change to [INAUDIBLE] splitting.
So if you remember back when we were talking about boundary
value problems, we're often trying
to solve problems that are like this.
All right.
So this is the conservation equation
for a species in a reacting flow situation.
So you have some convection, you have some diffusion,
you have some reactions.
Now in my life I have a lot of models.
I have about 200 species in them.
So this K is 200 different variables.
And then you have these equations, and on top of them
you have [INAUDIBLE] equations for that momentum conservation
and continuity.
And so this is actually a pretty bad set of equations.
And typically the chemistry is stiff.
Because I have chemistry time scales
it'll be, like, in picoseconds, sub-nanosecond time scales.
And I'm caring about simulating-- well in my case
maybe I have an advanced new engine burning some fuel.
And so we're running over a simulation of a piston cycle,
maybe 10 milliseconds, but the chemistry is sub-nanosecond,
so then you multiply it out.
So that's like I'm going to 10 to minus two seconds down to 10
to minus 10 seconds.
That's eight orders of magnitude difference in time scales.
So I'm in a bad way trying to solve this equation.
And I have 200 of these guys coupled to each other.
And then I have to think about what my mesh looks like.
And I'm trying to simulate inside a piston,
actually I have a moving mesh because the pistons
compressing, so the equations are kind of complicated even
with what the mesh is doing.
And then I need a bunch of mesh points
there because the-- do you guys hear the Kolmogorov scale?
You guys learned this yet?
It's a scale like the minimum eddy size
if you have a turbulent flow.
And that's pretty darn small also.
So it might be sub-micron.
And so I have a natural, lens scale
that's like microns maybe.
So 10 to minus six meters, but the piston diameter
is couple centimeters, so again I
have like four orders of magnitude between the mesh
size I need to resolve the physics
and the size of the thing.
So I need about 10,000 points across the cylinder,
but it's a 3-D problem, so I need 10 to the fourth, times 10
to the fourth, times ten to the fourth,
so I need 10 to the 12th mesh points.
And then I have 200 variables, say variables at each mesh
point, so I have 10 to the 14th stayed variables.
So how much memory you have in your computer?
How much ram?
You have probably like 16 gigabytes.
So you have ten to the ninth- tenth to the ninth bytes,
is that right?
No-- what's giga?
Ninth?
OK, so you have ten to the tenth.
Sorry.
Ten to the tenth byes, but I have 10
to the 14th stayed variables.
So I have a problem.
So then you need to figure out how you solve it.
So there's two schools of thought of how
to solve these equations.
The school thought which has actually been the most
successful so far, but is not readily available,
is to get the government to build you
the biggest possible computer in the world with as much memory
as possible.
And then you can actually keep track of your state vector,
all 10 to the 14th elements, and you
reserve this whole computer for yourself for, your know,
a month.
And just run it.
And you figure out how to parallelize
all the calculations, and you just use an explicit solver.
And so that way it's just straightforward.
All you're doing is you're computing,
you know y of T plus delta T is equal to something,
and it's just an explicit formula
ode45, fancy version of that.
OK?
And it's just algebra, and you just add them up
and it's no problem.
But this is a limited feasibility thing.
This only became possible about three years ago
when they built a computer that actually had enough memory
to store the state factor.
It's also-- if the equations are too stiff it still
doesn't work, because the step size and time step you need
to use is so small--
so say it's 100 picoseconds and you're
trying to get it out to a millisecond,
so you need 10 to the eighth something
time steps, which means you need to double precision numbers all
the way through.
But if the government builds you a 64-bit 128-bit machine
then you can do it.
But anyway this is like a very special kind of approach.
So now most people, including me I
don't have access to that kind of computer power.
Actually if you want to read about that
there's a woman who does that, who is Jacqueline Chen,
and you might want to look on the internet and see about her.
And somehow she became buddies with the people
in charge of the super computers,
and so they're happy to give her a month of, you know,
the best computer in the world every year,
and she runs some really awesome calculation
that's exactly right.
Solves everything, has all the state variables
ten to the 14th, and that's fantastic.
And currently that's the benchmark
for how people test approximate methods for how to solve this.
Because she actually has the full solution numerically
converged.
But this again is a limited applicability.
So a lot of people are working on approximation methods
to try to solve this equation if I don't have access
to such awesome computer resources.
And the general idea of them--
well the two ideas.
So one idea is somehow get rid of the turbulence.
OK?
And so what they do is they use a mesh coarser than what
you what you would need to resolve the Kolmogorov scale,
and then they use models for what happens
for itty bitty tiny eddies.
And those are called sub-grade scale models,
and this whole approach is called large eddy simulation.
And the idea is you keep track of the large eddies,
the big recirculation zones that's modelled explicitly,
but you don't keep track of all the itty bitty tiny
swirling off little eddies that lasts
for a little second, a little bit of a period of time
before they dissipate because of diffusion.
OK?
So that's the concept.
And there's, I don't know, half of all
of the mechanical engineers who work in large eddy simulations
are trying to figure out how to do this.
Maybe Jim as well, I don't know.
No.
OK, but at least he knows some of them probably.
So that's one whole branch.
It's to try to get rid of the turbulence.
And basically that's trying to reduce
the density of mesh I need.
So trying to get rid of this--
use the spatial mesh.
But the other issue is the time mesh.
And the time is even worse, right?
Because I said I have ten to the eighth time points
and I only need ten to the fourth spatial points
in the dimension.
So in some ways, if I can do something with the time,
that really can buy me a lot.
I have eight orders of magnitude I'd
like to try to reduce to three orders of magnitude if I could,
that would be really nice.
That idea is to say, well let's separate
the fast timescales from the slow time scales.
And in fact all the fast time scales
are mostly all chemistry, particularly after I changed
the eddy simulation where I get rid of the really small spatial
scales, then the time scale of diffusion over this big mesh
is kind of slow.
And so I can go to a time scale that's completely
controlled by chemistry.
So we already told you, how do you solve stiff differential
equations?
Use a stiff solver, right?
And there are special methods, there are implicit methods,
and there's specialized solvers.
And what they do is they allow you to take larger time steps
than you would be allowed to use with explicit methods
where you would run into numerical instabilities.
If you try to use the explicit method with a big time stamp,
it's much bigger than the real physical time scale,
you'll get some crazy oscillations.
But if we use the stiff solvers you
can get away with time steps that
are much, much larger than the real physical scale.
So that allows you to go from 100 picoseconds up
to maybe 100 nanoseconds.
Get rid of three orders of magnitude of time step
by using a stiff solver.
So that's pretty good, three orders magnitude reduction?
Well go for it.
It's not only reducing the CPU time potentially,
but it's even more importantly reducing
the number of points you have to save,
because when you run this calculation what Jacqueline
Chen ends up with is a list, y at every x, xn, yn, el, times
e, something like this.
Alpha, I don't know.
That's the output of her program.
OK, so how big is this object?
This is 200 times 10 to the fourth, times 10 to the fourth,
times 10 to the fourth, times 10 to the eighth.
Right?
That's how many numbers she's going
to get out of her simulation.
So, actually, her biggest problem
is not actually solving the simulation.
As long as she can get the one month of CPU time she's good.
She's got to solve.
Her Problem is how to analyze her data.
Because now she has this much data
that she has to figure out to say
something rational about it.
And so she has like a whole army of post-docs
whose whole job is trying to invent new day interrogation
methods.
And also this data set is so large it cannot be stored
on a hard disk.
In fact it cannot be stored in the biggest farm of hard disks.
So she has to, as the calculation is running,
send the data to different farms of hard disks
around the world and filling them completely up,
and then sending them to the next one and filling them up.
So she has fragments of this data
set stored all over the place.
And then if she wants to make a graph
she has to have a special graphs software that
looks inside this data farm, get some numbers,
look inside this data farm and puts it
all together in some way.
So this is a pretty hard project.
OK?
Now if I can use a stiff solver I
don't need to use so many time stamps.
So I can get this down, say 10 to the fifth at least
I got rid of three orders of magnitude in my data.
Right?
So now I just need one data farm maybe,
instead of, you know, using all the ones at Amazon and Google
own or something.
OK?
So this is a pretty big advantage to do it this way.
Now the disadvantage, as you know,
is solving implicit differential equations
is way more expensive, right?
Because you have to solve an implicit non-linear algebraic
equation at every time step.
And in this case it can be a pretty complicated equation.
Because if I have this many state variables, so 10
to the fourth, times 10 to the fourth, times 10 to the fourth,
so 10 to 12th times 200, so 10 to the 14th at each time step,
that's my current state.
And I'm trying to solve a non-linear algebraic equation,
10 to the 14th unknowns, you need a pretty good solver.
I mean, I know backslash is good and F [? solved ?] and stuff,
but I really don't think it's going to work.
So that's a really serious problem.
So the next trick is to say, OK, well this stiffness is only
happening because of the chemistry.
So I don't--
I could solve all of the finite volumes separately.
So the idea is to split the problem up.
I'm going to solve the chemistry using a stiff solver only
at one point a volume at a time.
And then I'm going to do the transport some other way.
So what I'm going to do is I'm going
to rewrite this equation as My unknown is dy to t,
is equal to some kind of transport term
and some kind of reactant term.
And this term is local.
All this in this equation, this is really
same dy at position x depends on this at position x.
Whereas if you compute these guys
I need to do finite differences or something
from adjoining finite volumes.
Right?
I have to do fluxes from the adjoining volumes.
So these guys are sensitive not just
to my local value of y at this finite volume,
but also to the neighboring values of y.
But if I don't worry about them, if I break it up
I have this local thing and this non-local thing.
This one is local, it's stiff.
This one is non-local, but it's not so stiff.
And often it's linear or nearly linear.
Because the transport terms are mostly linear.
And so I can use--
I have specialized solvers that people
have worked on for many years for each of these things.
So for the stiff thing I have things like ode15s,
I have DAEPACK and [? DASPAK ?] and [? DASL ?] and sundials
and all those programs that people have written and spent
a lot of effort to solve this problem.
And similarly there's a whole bunch
of people who have spent their whole life trying
to solve the convection diffusion problems,
the non-reacting flow problem.
And so I can use specialized solvers that they have.
So I want to split this up.
And what I really would like to solve
is I'd prefer to solve dy/dt equals r and dy/dt equals t
completely separately if I could.
Because I can use specialized methods
that are the best for each.
So not to figure out, well this is not really
the real equation because it has both of them added together.
How can I handle that?
So this operator splitting recipe
is solve dy/dt equals to t for a while,
and at the end of this process I'll have a [? yat ?] final.
And then I'll solve this equation with y equals y0.
That's the equation system of solving, of transport.
I do some transport for a while, I get some values for the y,
then I take these guys and I solve
the reaction equation starting from y equals y final.
How to write this?
One or two 0.
OK?
And I solve this one for a while and I'll get a different,
you know, a somewhat different [? yrt ?] final.
All right?
And the simplest way of putting things
is just keep going back and forth.
Boom, boom, boom, boom, boom, boom, boom, boom.
And in the limit where I make this t final minus
t0 to be very small, then maybe I
might expect that they might converge to the same solution
as the couple of questions.
All right?
So that's the concept.
Now, if you do it this way it doesn't work very well.
You need to make this delta t incredibly small
to get any accuracy.
However you can do a slightly more clever version, which
is I saw dy/dt equal t starting from yt0, so y0,
and I solve this to me y at t plus delta-t over 2.
Then I solve dy/dt is equal to r, starting from y of t0
equal to, call this y star.
y star and this I integrate out to y
of t plus the full delta t.
I'll call that y double star.
And then I solve this one again, and I
solve starting from y of t0 equal to y double star,
and I integrate that out for another delta t over 2.
And this one I called my y final.
And that's the one I use.
So I take a half of a transport step,
sort of like I'm starting to do the transport,
then I suddenly-- oh, hold on.
I got you my reaction.
I do my reaction, then I do the other half of my transport
after that.
And if you can do the analysis of this,
it turns out this thing converges
to the true solution of the couple equations
with second order accuracy.
So you might make this delta t, if I cut it
in half the precision increased by factor of four.
If I cut if by a factor of 10 it increased by a factor of 100.
So this is looking more promising as a way to do it.
This is the most popular way to do it.
So this recipe here for operator splitting
is called Strang splitting.
Named after Gill Strang, he was a professor
in applied mathematics department here.
Former head of CIAM, the international body
for applied mathematics.
And also author of a fantastically great
linear algebra textbook which maybe some of you guys
might have seen some videos of his.
So he invented this 1968, I was five years old.
He didn't know that his invention
was of any practical use, being a mathematician.
So he just published that paper.
Now in the meantime, this paper is like a citation classic.
It's cited like, I don't know, 20,000 times.
He had no idea.
He never looks at citation list, he's a mathematician.
So I was working on this problem and--
let's see-- let's back up.
So this is how he usually solves it.
Right?
We usually solve the problems, the--
splitting, sorry.
Just be able to see it.
And we like to solve it this way because I
can use specialized solvers.
I use a transport solver, I use a stiff solver here.
This one I can parallelize because I
can do a different ode solution at each finite volume.
So if I have a big computer, a big parallel computer,
say I have 10,000 nodes, I can solve 10,000 finite volumes
simultaneously in this step, so the CPU time
is reduced a lot because it parallelizes easily.
So this is really popular.
This is the main way that people do things.
And then in the ideal case you do it
with a certain guess of delta t, then
you try delta t divided by five.
And if you get the same result then you publish it.
OK?
And you've converged your solution.
However, what I found doing this is I was trying to find flames.
I was trying to calculate the positions of flames.
And what I found was that with the delta t's that I would try,
sometimes I would get a qualitatively different
solution depending on the delta t I chose.
And so in particular what I found
was that I was solving burner flames.
So I have, like, a burner here, and it had a flame on top of it
like this.
And there's two importantly different cases.
One is where the flame is floating above the burner.
And one is where the flame kind of goes all the way
and touches the metal.
And you've probably seen both those cases,
like a bunsen burner or something like that.
Sometimes if you get the flow rates high enough,
the whole thing lift off and be floating there.
And sometimes if you slow the flow rates down
they say that they will anchor to the metal.
Some piece the flame will actually touch the metal.
Usually it's like if you have a--
a lot of times you have a case like this.
You have like a metal piece, metal piece,
here's your mixture coming in, out here is the air,
and a lot of times the flame will actually
touch this metal piece here.
You know, the flame sort of looks like that.
Have you ever seen a flame like this?
Anyways, this is a common thing, anchored-flame.
If you just have a slow speed burner that's what you get.
If you go high speed on the burner
you can make the whole thing lift up and sort of
dance around if you look at the flame.
All right?
You could probably see the same if you look
at a piece of wood burning.
If it's burning well the flame will be lifted off the wood,
you'll see a gap between the flame and the wood.
If it's not burning that well you can see bits of the flame
actually touch the pieces of wood.
OK?
So that's a very important thing for a commission guy.
So was really-- to me this is a big deal.
So when I solved it with one choice of delta t
I get the solution.
When I solved with another choice of delta t
I get the solution, with a nice big gap.
Now they both can't be right.
And so I was a little bit alarmed
because it's like, it's really different, right?
So I thought, you know, I didn't really
care that much about the dynamics of my problem.
I just wanted to make sure that the steady state flame I ended
up with was the right one.
So if you're interested in steady state problems
this is a really important thing.
You want to make sure that your solution really converges
to the true steady state.
So what I demonstrated primarily is the Strang splitting method
is not reliable to converge steady state.
Now I believe I made the delta t tiny enough,
probably I would converge the [INAUDIBLE] steady state.
But I don't know how tiny I have to go,
and I can't really judge my convergence criteria
very well, because at some point of delta t
is it suddenly jumps from one solution to another solution
that's qualitatively and completely different.
So the convergence is not smooth and continuous,
it's like it's converging, converging, converging,
jumps to some other solution.
And I don't know where the jump happens.
And I don't know if I keep going further,
will it jump to something else?
So this makes me very loathe to publish the paper.
Right?
So I don't know.
So I was agitated about that.
So I start working on trying to fix this splitting.
The problem was that Strang is really smart.
So his splitting method is really good.
It's stable, it's pretty accurate,
it's really hard to beat it actually.
And so I published the papers and they were that good
because they really didn't really improve it that much.
But then I had a really smart post-doc
named Ray Speth, now works in the AeroAstro department.
And he came up with a different way to look at this.
And he said, you know, if I had a problem that's like this,
I can always add some constant to this
and subtract the same constant from this,
and I'm really getting the same problem.
And so now the thing is, can I cleverly
choose this constant so that, for example, my steady state
solution is going to be good?
So how to figure out what constants to choose?
And so the methods where you add these constants
are called balancing methods.
So it's trying to balance the operator splitting scheme.
And what Ray's method--
pretty simple balancing method.
So this is called simple balancing.
Let's see if I have a slide that shows it.
There we go.
Here we go.
Simple balancing.
You choose a c to be the average of your right hand
side at the previous time step.
OK, now idea here is that as you're
getting close to a steady state, the previous time step
is not going be that much different than the current time
step.
And so this average number is going
to stay more or less study as you as you get close.
And so what you're really doing is adding the, sort of,
the average of these two guys to each of them
from the previous step, and that has
the consequence of really changing the solution.
So let me just show you.
So the left hand side is what Strang splitting does.
The blue chips are the transport,
and the green steps are the reaction.
And what's happening is I did a half a transport step,
then I do a reaction.
The reaction sort of jumps me off the trajectory,
then brings me back to the trajectory,
then overshoots because I do a whole delta t step of reaction.
And then I have to do a half step of transport back.
And so I get from the black-- one black dot
is the first point, and the other black dot
is what I call y final.
And that's showing what happens.
And it just repeats like that over and over again.
And it makes sense.
If you're close to a steady state,
then if I'd just take one of these terms,
without the balancing theta transport,
it's going to head away from steady state,
because at steady state the transport and the reaction
are balancing against each other.
Typically you have something transporting in
and is being consumed by the reaction.
Or something is being created by reaction
and is being transported away.
And they have opposite signs.
If I only put one in that has a certain sign,
then I'm going to move away from the steady state.
And then I need the other one to bring it back.
And the reason why Strang splitting has trouble
is it's doing these huge excursions away
from the steady state at each time step,
with the hope that it ends up back at the steady state.
Right?
That's the way it's supposed to work out.
And to second order it does, that's
because it's pretty clever about how it works.
But it does it's crazy excursions.
Now if you have a non-linear model,
the further you do excursions from the real trajectory
the bigger the area you get.
So if you do that balance method,
you can see what happens is that when you get to steady state
all the steps are just moving you along the steady state
trajectory.
They're basically maintaining the steady state,
and that's because I added--
it's like transport plus half of transport
plus r or something like that.
And what is it?
There's probably a minus sign here.
And so I'm getting rid of half of my transport
and I'm adding in half of this.
So I changed-- before I had dy/dt is equal to t,
but now I have dy/dt is equal to 1/2 of t.
I feel I need a sign to stick here.
So half of t minus f of r, something like that.
And the-- this can't be right.
I must have a sign.
Maybe this is not the average it's a half.
It should be minus.
It's like the-- these things should be added,
and these guys are you cancelled out,
so basically [INAUDIBLE] is going to be zero,
which it should be as a steady state.
So I'm sorry, that previous slide had the sign wrong.
This is what happens as the solution.
If you do Strang splitting on the left
you get more accuracy as you cut delta t, but it's only--
it's not getting any better.
It doesn't matter if it's near steady state
or away from steady state, it treats everything equally,
the Strang splitting.
In the balance splittings, as you get to steady state,
all of a sudden the errors drive exponentially.
You see this litte logarithm of the error?
Logarithmic plot.
So you get tremendously great convergence
to the real steady state once you get anywhere close
to the steady state.
Now the cost of this is that because the formula, which
has a sign mistake I think, it depends on yn,
it makes the whole thing more numerically unstable.
Because it's like an explicit thing.
We're using what happened at a previous time
to correct our equation for future times.
And so you are more likely to get numerical instabilities.
So then Ray invented an implicit version
of that, and also a higher order method,
and that's called re-balance splitting.
The equation is way too complicated for me
to write down, since I can't even write down the minus signs
correctly in this one.
But you're happy to read the paper, it's in CIAM.
And he has the equations for that.
And when you put the second order one in, that's implicit,
then the thing is just as stable as Strang splitting,
but it also has this great property
that as you get close to steady state
it exponentially converges to the solution, which
is sort of like how Newton's method is really so good.
It's sort of like that.
You get a lot more decimal places really fast
once you get close to steady state.
All right, that's enough.
Enjoy Thanksgiving.
I posted a bunch of readings and things
about metropolis Monte Carlo and stochastic stuff
that we're talking about next week, and also a little bit
more about models versus data.
So do the reading when you get a chance.
Say hi to your families, say hi to your friends.
Have a nice time.