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: Ready to go.
So first of all, to say thank you to several people
in the class who sent codes to find these wavefronts when
a shock, a step function is moving along with velocity c.
I had roughly drawn what these three methods produced,
but it was much better to see it in Mr. Sekora's movie that's
on 18.086 site.
And sure enough, behind the moving front is an overshoot,
that reminds us of the Gibbs phenomenon, from Lax-Wendroff,
and I think it's sort of unavoidable if we have
a negative -- because one of the three coefficients
at the previous time is negative for Lax-Wendroff,
that's how it gets its second order accuracy,
and we see that better accuracy in a much steeper,
much steeper profile.
Then Lax-Friedrichs, if you've noticed,
has a sort of jagged profile.
Maybe you've realized why that would be.
Somehow Lax-Friedrichs kind of operates on a staggered grid.
We'll see, actually with -- I've drawn a staggered grid here,
so I can say what I mean.
Lax-Friedrichs computes this value
from the two earlier values there and then
at the next level again and so forth,
so that actually, Lax-Friedrichs is using two staggered grids
separately.
I think -- normally, with Lax-Friedrichs,
we call that n plus 1 and n plus 2 and n plus 3,
but the picture is still right, so I think that's probably
responsible for what you see, that there's two waves,
two fronts, right, one delta x away from each other.
And finally, upwind is also smeared out.
I guess I was so pleased to see these actually showing
on the website and that makes me think of more questions.
Maybe you too.
First of all, as time goes forward,
we see this bump beginning to appear.
So am I right in thinking that it reaches a certain height
and stays there?
So I guess I'm thinking about sort of further questions,
because every time you do a numerical experiment,
good additional question show up.
So I'm wondering, as this wavefront forms,
this characteristic profile for these three methods,
first of all, does that profile settle down
to a steady profile?
Ultimately it'd be great to be able to predict
what it would be.
In particular, does that bump reach a certain height?
Because you can see it just barely forming and then
you see it growing to a certain point.
And then, when the movie ends, fortunately it
plots this picture and you really see the profile,
and you can actually see, if you look closely,
that it goes above 1 again and probably more.
That would be, does it approach a steady profile?
Does that height settle down?
It'd be fantastic to be able to predict that height.
For example, you always want to know, how do things
depend on the parameter?
The one parameter we have in this problem is r --
c delta t over delta x.
If I change r, I'm changing all the coefficients.
If I changed r to be exactly 1, then the profile
would be perfect, right?
That's the golden ratio, r equal 1.
I would be computing that value exactly from this value,
at exactly the right slope, on the characteristic,
and perfection.
Now as r goes away from 1, I presume --
it goes down below 1.
Of course, if it goes above 1, catastrophe.
If it goes below 1, I presume the bump begins to appear,
goes higher.
I don't know.
Could you see how the dependence on r appears?
And then I can even suggest one more thing, which is --
let's put a fourth guy onto this picture,
which would be leapfrog.
Leapfrog for first-order equations,
so let me just put down below what that would be.
So it would be on that kind of a staggered grid.
A new value at time n plus 1 would
come from these at time n, really, and this one, at time n
minus 1.
So leapfrog -- you'd think of it right away.
I use a centered difference in time, U_(j, n+1) minus U_(j,
n-1) over 2 delta t equals c times a centered difference
in space, U_(j+1, n) minus U_(j-1, n) over 2 delta x.
So you see right away, the 2's cancel.
The delta x times c -- the delta t times c over delta x is
our old ratio r, so this is easily written in terms of r
again.
It does have this two time step issue
and of course, we have to think about stability.
I'm writing down a fourth candidate for the equation u_t
equals c*u_x, of course, before I come back to the second-order
equations, or systems.
I just wanted to put in these comments,
partly because I'm thinking about projects that are --
and I hope you are -- so some people will have theses
or problems from other courses which involve solving PDEs
and that'd be perfectly natural to base some project
on, but if you don't have some specific application,
then here are a whole lot of interesting questions,
to just go further with that.
So the dependence on r would be a --
the profile dependence on r, c delta t over delta x,
would be a question.
Let me just emphasize, what's the goal?
When I say that this upwind one is smeared out,
I mean that it doesn't satisfy what you really hope for.
You want to capture that shock, that discontinuity --
not perfectly, you can't expect that --
but you want to capture it, maybe in, let's say, 2 delta x.
A good method captures the shock within 2 delta x
and does it without much undesirable oscillation,
physical oscillation like that.
So these are methods where we're trading off a good shock
capture versus a smeared one, which is not satisfactory,
really, in a lot of applications.
So capturing the shock within 2 delta
x, roughly, is highly desirable, and you
might have to go add a artificial viscosity that we'll
discuss to get there.
So that's past things, but still very much present,
if I can say.
Here I've introduced a new leapfrog,
a one-way equation leapfrog, I'll
call this, where the last lecture was
about two-way leapfrog.
Leapfrog for u_tt equals c squared
u_xx, a two-way wave equation.
I might just mention about this one-way leapfrog,
once you put in e to the i*k*x as we always do,
or e to the i*k*j delta x in the discrete case --
so you would put in an e to the i*j delta x times k,
and then find the G, the growth factor in --
that'll depend on k and will tell us the time dependence.
It'll be a G to the n-th.
I think you'll find that for this guy,
the absolute value of G is exactly 1 while it's stable.
Of course, if the Courant condition is violated,
then there's no way it could be stable
and G's will grow, be bigger than 1.
But, I think in the stable range, you will have --
this G will have absolute value exactly 1.
So it'll be on the unit circle.
You don't have much leeway, right?
We found upwind was, like, well inside;
Lax-Friedrichs was well inside; Lax-Wendroff,
because it had higher accuracy for low frequencies,
was closer; and I think this guy is right on.
So this is the first-order methods.
This would be Lax-Wendroff and this would
be this second-order method.
So you're really playing with fire there
and maybe this is not a favorite, just
because it's too close to the edge, in fact, right
on the edge.
So those are some comments really generated
by your numerical experiments and I think that's the way this
course should -- homeworks and projects should develop out
of numerical experiment.
Now I'll come back to complete my lecture
on the second-order wave equation.
So we studied leapfrog for that.
I won't repeat that, but then with second-order wave
equation, we have another option,
which is to reduce a second-order equation with two
time derivatives to a system of two equations that are
first order in time and space.
I've chosen to use the letters H and E as a kind of hint
that we're touching a really big application, which
is electromagnetism.
So E for the electric field, H for the magnetic field but here
in 1-D -- so it allows me just to mention that in 3-D,
when E has three components, H has three components,
this is really the curl and this is really minus the curl.
And the notes will show the beautiful Yee mesh,
staggered mesh that copies for finite difference methods,
the properties of the curl, our interpretation
of and the properties of -- that Maxwell discovered the physics
of, you know, current, magnetic fields going around a wire that
along which one an electric field is going.
But let me stay here with one dimension.
So first of all, can we recover the wave equation from that?
Let's just see that we haven't left behind the wave equation.
So I just take the second derivative --
so the second derivative -- I just want to recover the wave
equation so that I haven't lost it here -- so this will be,
if I -- I'm taking the time derivative of the first
equation, so I get d second h dt dx, right?
But by the useful, extremely useful fact that this is
the same as taking the t derivative first and the x
derivative second, now I have the x derivative of this
equation, so this is c -- now I have the x derivative,
or c times it, so now I'm up to c squared,
the x derivative of the x derivative, right?
So we eliminated H, got back to one equation, second order,
and it's the wave equation.
I want to show the same thing now for finite differences.
It's just -- it's such a simple and useful device and our point
will be that it works for -- just as well for differences.
That allows us to do this staggered grid idea.
So the staggered grid idea is, I'll do E on this grid,
at these grid points, H at the half points,
E again at the integer points, H again at the half points.
Now, so my difference method, leapfrog,
on this staggered mesh, will copy --
so the first equation will give E at the new time from
the difference in H at this time.
So, you see it.
dE/dt is going to be replaced by that time difference.
dH/dx is going to be replaced by this space difference,
and that will determine E at the new time.
So I need time 0 and time 1/2 to get started.
Maybe I put those down here, somewhere.
I have to use the initial condition,
the initial velocity to get started on E and H
and then I go onwards.
I've got to this point and this point.
Then this guy comes from this one, this one and this one.
Of course, that similarly comes from the same equation
shifted over.
Now I know these.
Now I'm ready for this one.
Coming from -- what am I approximating now?
I'm on the staggered part of the grid, the half steps.
I'm approximating this on the half steps, where H is defined.
So that minus that is the time difference.
This minus this -- which, I know these now --
is the space difference in this equation,
and that's what allows me to compute the new one.
You see, it's a simple idea.
It's close but different from my scalar
leapfrog, one-way leapfrog here.
Here it was the same U at all the points.
Here it is E at the integer time levels and H at the half
integer time levels.
It's natural to say -- I'm calling this leapfrog.
Could I do this for the difference equation?
I just want to see that if I eliminate H --
you saw what happened there.
By using the identity that the second derivative in t and x
is always the same as the second derivative in x and t,
I eliminated H and got an equation that involved E only.
What's going to be the corresponding identity?
Trivial, but just -- trivial things are not bad things
to notice -- for the different method.
This is what I wanted to be sure of, for the H part.
This H -- so the H were -- do you mind if I draw yet another
symbol?
I want to use this idea to eliminate H. So
let me put two H's there.
Is that right?
Yeah.
I think what I'd like to know -- let me just say what I believe.
I believe that the time difference of the space
difference of H at mesh points is
equal to the space difference of the time difference of H
at the same mesh points.
I meant to put these as subscripts.
It's not delta x.
It's delta sub x.
Difference in the time direction, difference
in the space direction.
Let's just see what this means.
Let me redraw these four points.
These are H. H at all those points.
So this says, figure out the space differences -- so this,
this minus this and this minus this --
and then take the time difference of those.
Can I see what I'll get if I do that?
So this means I take the space differences,
so I have 1 minus 1 and that would be 1 minus 1,
but now I'm going to take the time difference of those.
So that'll be this one minus this one
and this one plus this one.
Do you see that?
What we have here is H at this point minus H
at this point minus H at this point
plus H at this point plus 1.
I claim that we have the same thing here,
that if I take the time differences first --
can I do that?
Instead of taking the space differences first,
I'll take this one, this time difference,
this minus this and I'll subtract --
I'm taking now the space difference -- I'll subtract --
are you with me?
You see, it's just the same four values --
two 1's and two minus 1's, whether we go first -- in fact,
this actually is a -- I think maybe throws a little light
on the calculus identity, that we can exchange the order
of the derivatives, because we know we can exchange the order
of the differences, because we just got numbers like 1,
minus 1, minus 1, and plus 1.
If we let the spacing, the mesh, go to 0, this --
and divided by the delta x and the delta t,
I suppose we would prove the continuous case.
So in a way, this throws light on something that's probably
something that we take for granted,
but we wouldn't want to be pressed on too closely.
So having that identity allows me
to do exactly what I did with the continuous case there
and show that this, when I eliminate the half levels,
I've got ordinary leapfrog.
So when I eliminate half levels, I will be --
when I eliminate this half level, let's say,
then I'll be connecting -- when I eliminate half levels,
I'm going to be connecting these at three levels -- one, two,
three.
It'll be leapfrog.
So since it's been a few days since we wrote down
two-way leapfrog or second-order equation leapfrog, let me just
write it again, and that's what we would get.
The time difference squared of E will
be equal to c squared times the space difference squared of E.
And of course, here we have a delta x squared
and here a delta t squared.
So altogether, an r squared.
That's the formula we know.
So if I simplify this to make it look good
and put the delta t squared up there,
you see that it's just r squared.
Maybe it's worth realizing that, again, the magic r equal 1 --
The magic ratio r equal 1 gives the exact U,
agreeing exactly with u.
This equation, if r is 1, would be satisfied exactly
by the continuous solution and so it's the same as this one.
So that's leapfrog.
So can I summarize leapfrog?
Leapfrog was, first of all, that led us to these --
it worked for second-order equations.
It led us to the same stability condition
that r had to be less or equal to 1,
but it took a quadratic equation for G. Remember,
we had G squared minus 2G times some quantity
a that involved r's and cosines plus 1 equals 0.
The reason it was a quadratic was there
were two time steps involved.
We handled that.
We got stability, if r was less than or equal to 1.
That led us to a, in magnitude, less or equal 1.
And that led us to G -- I'll just remember those steps --
a less or equal 1, and that let us to G less than or equal
to 1.
Yes, the notes, including more about Maxwell's equation
and Yee's beautiful method, with a figure that shows his mesh,
will go up -- I mean, every section is getting upgraded
and this one is the next one to be upgraded.
Probably by tomorrow, you'll have a good section
on second-order equations.
For me, that's what I wanted to say about the wave
equation alone, with no diffusion
and now I'm ready to move to the heat equation.
So what do we have coming up?
Maybe important to say what's coming up.
So we've done one-way wave.
Now we've done two-way waves.
Next comes the heat equation, heat/diffusion,
and then will come a mixture of the two,
convection-diffusion -- and by the way,
so I'll go straight to these -- in a couple of weeks,
we'll have a guest lecturer who will do the applications
to financial mathematics, mathematical finance.
Specifically, the Black-Scholes equation.
So you probably know that a lot of effort and a lot of money
are going into the valuation of options
and other financial derivatives and that this Black-Scholes
equation beautifully produced the heat equation,
and additional correction terms will
produce a convection-diffusion equation.
So this is an application that we
wouldn't have made some years ago,
but now it's worth focusing on that.
I think it's interesting in itself.
So this is what's coming and then
you could say one-way wave equations, but nonlinear.
So nonlinear wave equations and those
are called conservation laws.
That's after these linear cases.
All right.
So I'm ready to tackle the heat equation, starting
then on to the next section of the notes.
So let me tackle it here.
So the heat equation.
Let me take the constant to be 1; u_t equals u_xx.
Because we already see the big difference.
Now, there's a second x derivative, but only
a first derivative in time.
The dimensions of t and x are no longer the same.
We won't have delta t and delta x comparable.
Now it will be delta t comparable with delta x
squared, and that's a very significant difference,
because that's -- if delta x is small,
delta x squared is extremely small.
So when delta t is constrained by delta x squared,
we're looking at small time steps.
We're looking at stiff systems, absolutely.
This is going to be a stiff problem,
because the low frequencies and the high frequencies
differ enormously in the decay rate.
In fact, why don't we see it right away,
starting as always with exponentials.
So if I plug in u of x and t is some G of t e to the i*k*x,
separating variables as always.
Look for a pure exponential solution,
plug that in and we get dG/dt times e to the i*k*x
on the left, and now I have to take two x derivatives,
so that brings i*k down twice, so that's minus k squared --
i squared giving the minus 1 -- times G, times e to the i*k*x.
And now, of course, I can cancel that, which is never 0,
and I've separated out the time part,
the G part, dG/dt is minus k squared G. So G
is e to the minus k squared t.
Notice that it's k squared, and notice
that it's real and negative.
See, that's a very big difference.
Compare e to the i*k*c*t, which it was for the one-way wave
equation.
The difference is enormous here.
That's magnitude 1, energy conserved.
This is magnitude smaller than 1, energy dissipated.
We had a ratio delta t over delta x.
Now we're going to see delta t over delta x squared
because somehow x and t are -- now have --
it's x squared that matches t now.
So now, when I put that into here -- let me just do it --
G of t is e to the minus k squared t.
So there is the pure exponential.
Simple as always, because we have constant coefficients,
but highly informative.
I guess another thing that I read off of that is,
what frequencies are decaying fast?
High frequencies?
If k is large, k squared is very large
and e to the minus k squared t is decaying very quickly.
So high-frequency noise in this system,
roughness, discontinuities are going to get smoothed out,
strongly, by moving forward in time.
Physically, what's happening is, we
have our wavefront, which stayed a perfect step
function in the wave equation, will instantly smear.
So diffusion smears discontinuity.
Heat travels.
If we were back over here, if the temperature starts out as 1
on this side and 0 on this side, then in an instant,
some heat flows from right to left.
In fact, it flows with infinite speed.
In a short delta t, there's a little heat all the way out
there -- very little, of course, because it's --
I guess we want to see what that behavior is.
So I now want to write down -- staying with the differential
equation, what do we want to do?
We did the exponential case and learned a good bit
from that, the fast decay of high frequencies.
The next step would be to put different frequencies,
different k's together.
By linearity, we can add these solutions,
we can multiply that by any numbers we want,
depending on k; we still have solutions.
And we can integrate over k, so we will, actually.
That that will be the next step.
Let me write down what that does.
I'm going to jump, take a little jump here
to put on the board what I just said in words.
I can multiply this, e to the minus k squared t e
to the i*k*x -- that's what's happening to the pure frequency
k.
If I multiply that by the amount of that frequency which is
in the initial condition -- so here is the initial condition,
u_0, and I've taken its Fourier transform.
This tells me how much of frequency k is in the problem,
is in the initial problem.
This tells me how much of frequency k
is in the initial problem.
That amount multiplies this, with this rapid decay as time
goes forward, but at any later time, I'm just adding up --
and since k -- we're on the whole line --
k is a continuous -- so we add up from k equal minus infinity
to infinity and I think we have to remember the 2*pi that
depends on our particular definition of that transform,
but I'm going to stay with this one.
That's the answer.
That's u of x and t.
So I've said that in words, but actually, we
could just check it.
We would like to know that this left-hand side solves the heat
equation, and we would like to know
that it has the correct start, the correct u of x and 0.
Let's check that.
At t equals 0, at the start, this is 1.
So at the start, this formula, without that thing
there, is just the inverse Fourier transform
that takes the Fourier transform, multiplies by this,
integrates to recover u_0.
So it does start out correctly.
Should I say that again?
It starts correctly because at t equal 0,
this is just the inverse Fourier transform of u_0,
so it produces u_0.
It inverts this transform that we took.
Secondly, does it solve the heat equation?
Sure, because for each k, we've just
found that that solves the heat equation
and now we're just putting different k's together.
So that's the answer.
So you might say, OK, perfection.
Do we need finite differences?
I guess we do.
This is typical of formulas that give the answer exactly,
but in terms that are not numerically convenient.
Because, to use his formula, first of all,
this this requires an integral, an infinite integral to find --
integral to find -- for u_0, u hat.
To use this formula, we'd have to find this transform,
and then we have to do this integral
for the inverse transform, this integral, for the answer, u.
Two infinite integrals.
And, not to mention the fact that on a finite interval,
when x doesn't go all the way from minus infinity
to infinity, we have to think all over again.
So we have a nice formula for the answer, but --
it looks nice, but it's not numerically terrific.
We can get a lot of information out of this
and there's one other special solution that's
also extremely informative.
What if -- so I'm going to do a special case,
but it's the big case.
What if the initial conditions is a delta function?
What if the initial condition is a delta function?
An impulse at time 0, x equal 0, an instant source of heat
at a point -- a point source with finite strength,
you could say.
An impulse at times 0.
For the wave equation, what happened?
That was like turning on a flash of light
and it went along the characteristics.
For the heat equation, it'll be quite different.
There aren't characteristics here,
but this is an extremely important solution,
and I think you could call it the fundamental solution
to the heat equation.
You may know what it is.
You may have seen this formula.
It's fantastic that we get a formula for this one.
You could say, we've got a formula.
What's the Fourier transform, what's u_0 hat of k
for the delta function?
1, right?
The Fourier transform of the delta function is a constant 1.
All frequencies are there in equal amounts.
This is a constant, so I can remove it.
It's 1.
I'm left with a specific integral to do.
And the question is, can I do it?
And the answer turns out to be yes.
You can do that integral when this is a constant
and that will give you the answer
with that starting value.
Let me just say what that answer is.
u of x and t -- this is the fundamental solution,
u fundamental, maybe I should say,
because it's such an important case -- is, turns out to be --
it's not obvious what that integral is.
By the way, normally -- let's see -- normally it's --
to get the answer that I'm going to write down,
we take a kind of end run on the integral,
and I'm going to take a serious end run and write down
the answer exactly.
Here's the key part of the answer.
It's a solution which starts from a delta function, spreads
of course, and the shape as it spreads is a bell.
It's the famous bell-shaped Gaussian distribution of e
to the minus -- so is it e to the minus x squared, I guess?
And then, always here comes -- that scalar tells us the width
of the bell.
How far has it spread at time t?
And what goes there is 4t.
That's the key observation.
And notice how x squared and t are coming in.
It's that ratio, x squared to t, that's crucial.
That's the parameter again.
Of course, we need to multiply by a constant.
How do I know I need a constant?
Because the total heat is conserved.
In other words, the integral of the delta function,
the original source of heat, was 1.
The integral from minus infinity to infinity
of the delta function, which is all totally concentrated
at one point, is 1.
So I need to put on whatever constant
it takes so that the integral of this thing
shall stay 1 for later times.
This is one integral from minus infinity to infinity that we
can do and it turns out to need 1 over the square root
of 4*pi*t.
That solves the heat equation.
It's not a lot of fun to plug that into the heat equation,
but it's possible, right?
The x derivative, the time derivative
is going to be messy.
It'll have two terms.
The x derivative will have two terms and they'll agree.
The notes will give one way to do the integral
and get that answer.
You see that that integral, it's --
integrals with e to the minus k squared in them are impossible
on a finite interval, but here we're going all the way from
minus infinity to infinity, which you might think makes
the problem harder, but actually it makes it a great deal easier
and we can get an explicit answer, a beautiful answer.
So we learn that the solution to the heat equation,
coming from a point source, is a bell-shaped Gaussian curve that
gets wider and wider as t increases,
but I guess you would want to look to see,
what happens if t comes down to 0?
Does that really approximate -- I meant to say,
does it really converge to -- this should, as t goes to 0,
this should go to the delta function, and in some way,
it does.
In some way, it does.
It's wonderful.
The delta function, of course, is 0 away from the origin.
So why does this approach 0 away from the origin?
Suppose x is 1.
Then do you see this thing going to 0 as t goes to 0?
Look what's happening.
As t goes to 0, this part is blowing up, but kind of weakly.
That's not a disastrous blowup.
Compared to this e to the minus -- let's say 1 over 4t.
As t goes to 0, that's e to a very negative exponent.
That's going to 0 very fast, much
faster than this is going to 0.
If you like to think of l'Hopital's rule or ratio
or something, this quantity is going to 0 way
faster than this one is, and the result
is 0, except at the one point x equals 0, where this is 1.
The fact that the total amount of heat stays at 1
tells us that the delta function is the limit there.
Today was the heat equation, continuous case,
with these explicit answers.
Then next time is -- tomorrow, I guess -- is the heat equation,
finite differences.
Again, we'll have several difference methods
and to compare them would be terrific.
Please think about that, comparison with the heat
equation and/or what I was saying
at the beginning of class, a deeper look
at the comparison for the wave equation.
See you tomorrow for difference equations for this problem.
Thanks.