As you can see, there are a bunch of handouts from the coffee table book.
Let me first of all mention the one with the beautiful spots and scrolled picture of the great dot equation.
A beautiful set of equations which you're all working with on Assignment three.
So this will give you a little bit of the background.
These equations have become one of people's favourite illustrations of how tuning a few parameters can change the qualitative behaviour of solutions.
So you see, this picture shows different regions of parameter space.
In many regions solutions are relatively boring,
but then there are other regions where they have interesting structure and sometimes the structure looks like long rolls,
and in other places it's more like spots. And in the assignment you'll be seeing some pictures like that.
Let me also alert you to a few things here.
So Friday's lecture will not be in this room, but in L two, and then we come back here for the final week next week.
We have an assignment due next Tuesday and then I will be handing back assignment
to at the end of this lecture and I'll probably forget unless you remind me.
So I trust you'll remind me.
So today we're going to talk about higher order difference approximations and in particular the limit when the order goes to infinity,
which is called spectral district ization. So 5.9 Fourier Spectral Description zation.
And to talk about that, we're going to begin with finite differences to motivating.
And I'm going to do everything today with the wave equation.
So in fact, it'll be just the one DX wave equation will be today's illustration u t t equals u x x which has solutions that travel at speed,
one to the left and the right. And two of the handouts actually are about the wave equation.
So let's glance at them. The one D Wave Equation handout describes how you can represent a general
solution to the wave equation in terms of left going and right going solutions.
I hope you're reading these coffee table handouts because they're very carefully put together.
You're getting an exposure to all sorts of PD that is very hard to get elsewhere.
The other one about the wave equation, the just with the title wave equation talks about how things change in higher dimensions.
So of course you see a picture of a propagating wave there.
You could imagine it's a gravitational wave, if you like, which would be correct, by the way,
because the gravitational wave, by the time it is a few light years away from those black holes, might as well be linear.
So it really is essentially the wave equation.
One interesting mathematical fact about wave equations is that dimensions and dimensions are different from even dimensions,
but the case end equals one is special. So in odd dimensions all information travels at speed exactly one.
And it even dimensions. All information travels at speed less than or equal to one.
Except that dimension one is special. It's like an even dimension where the information travels at speed less than or equal to one.
That's called Huggins principal. So having fun reading those, please.
Now what I want to do is immediately show you some simulations by finite difference methods of this wave equation.
And we're going to do the second order.
Fourth order and sixth order disorganisation.
So the second order then. So let's get the leapfrog sort of representation.
So that would be a second order formula. And this would be a fourth order formula.
And this would be a sixth order formula.
I hope those pictures are enough to tell you everything you need to figure out the equation.
Let's get right down one of them. So in this case, what I mean is that V and plus one,
minus two the and plus V and minus one all evaluated at the middle space point divided by k squared equals
v j plus one minus to v j plus v j minus one all evaluated at the middle time step divided by h square.
So that's the second order leapfrog style, discrete ization of the way equation.
And in these cases, we do the same thing in time, but something different in space.
And in this case, the coefficients turn out to be minus a half for thirds, minus five has.
4/3. Minus a half. And in this case, those coefficients turn out to be 1/90.
Minus 3/20 three has.
-40 9/18.
And three have. Uh.
I forgot. My signs are on. Three has minus three 20th and 1/90.
And you know where to find such numbers in the tables that I handed out from Feinberg.
One of his tables gave you all of these numbers. And of course, they can be generated in a systematic way by various bits of mathematics.
So today's lecture is really about what happens as you continue to higher and higher order,
and in particular, making sense of the limit of infinite order.
Just to give you an indication of what happens, the number two here and then five has and then 49,
18th, those numbers in the limit of infinite order of accuracy become pi squared over three.
And you'll see what I mean by that in a moment. So let's look at the computer to get a sense of how these things work.
I'm looking at code M 49 now called the wave equation and not much point in looking at it because it's just like all our other code.
Let's run it. So I'll say m 49 wave equation and it asks me how many points I want in my grid.
So I'll say 64. And then it asks, Do I want the second, fourth or sixth order formula?
So I'll say second. And then it computes a wave propagating to the right at speed one.
That's pretty bad, isn't it? There's a lot of spurious affiliation there that shouldn't be there.
Let's improve it by taking a higher order approximation.
I'll take four. Four. That's better.
There's still a little bit of dispersion being introduced by the numerical scheme.
Let's try a sixth order. That's better, but there's still a bit of dispersion being introduced by the numerical scheme.
So that little tail of oscillation should not be there. It ought to be just a cleanly propagating way.
And I need to say a word about why it's moving to the right.
The equation is completely symmetrical. Things can move left or right.
But in this code, I've given initial data corresponding to moving to the right.
And in particular, if you see at the beginning or the middle, there's a block of code that says time steps, zero and one.
So I set a value, you old, and then I set a value.
You add the fact that that value you involved X minus K is what makes it move to the right.
If I had said x plus K, it would have moved to the left. So this is the problem that you see in numerical simulations, especially of waves.
You have the spurious dispersion effect which, as you can imagine, could confuse you a great deal in an application.
There's a surprise that comes up in wave simulations that linear problems are kind of harder than nonlinear problems because nonlinear waves,
thanks to the non linearity, tend to be self-organising.
A solid time, for example, keeps it shape. And numerically, that sometimes means that nonlinear simulations are less prone to this kind of problem.
But a linear simulation, this sort of thing happens all the time.
Let's refine the grid a little bit. So instead of 128 points, let's take that sorry, 64 and let's take 128.
So with the second order. You see, things are better.
But there are oscillations. Let's try the sixth order.
128/6 order. That's looking very good now, isn't it?
There's really hardly any visible evidence of the dispersion. There's a little bit there, but not too much.
On the other hand, 128 is a lot of points to resolve a wave as boring as that.
We have at least ten or 15 points in the way, which seems a lot of points to resolve a wave accurately.
So spectral methods are all about getting beyond that. But before we get beyond that, I want to tell you a little more about this dispersion effect.
And indeed, let's go back and do a plot that shows it more clearly.
So I'll do I'm 49 again, 128 points.
But second one. So I want to say a word about that word and where it comes from.
You remember that diffusion is when things decay, basically, and dispersion is when things are spread into a pulse of different wave numbers.
And the way people analyse dispersion is with dispersion relations, which are all about tracking and how frequency depends on wave number.
So the procedure throughout physics and mathematics, what you do is you look at solutions to the wave equation.
Of the farm U of x and t equals an exponential or a sine wave.
But mathematically, it's always simpler to write an exponential.
So I'll say frequency times. Time plus wave number times space.
So i omega t plus c x.
And when people analyse dispersion, what they do is plug this into the equation and see how Omega relates to see.
So let's make sure you know the words their frequency refers to time.
It's the dual variable of time. And wave number refers to space.
It's the dual variable of space. So if you plug that into an equation, you get a dispersion relation.
If you plug it in to u t t equals u x x.
What you get is minus omega squared equals minus c squared minus omega squared equals minus c squared, i.e. omega equals plus or minus C.
So that's the dispersion relation for the wave equation.
And from that dispersion relation, you can infer how different ways will propagate.
Let's plot the dispersion relation. If this is wave number and this is frequency, then it has two branches.
One of them looks like that and one of them looks like that.
And from the slopes of those curves, you can figure out what the wave speeds are.
So this has slope plus one.
And this, of course, has a slope minus one. And that turns out to correspond to a wave speed of minus one and a wave speed of plus one.
So they're in the dispersion relation. You can see the right going and the left going solutions to the wave equation.
So there's no dispersion there when the dispersion relation is a straight line.
Everything goes at the same speed and that's all there is to say. The wave equation doesn't do very much.
I mentioned in the handout the Higgins phenomenon, and I said that not everything goes at the same speed, but every sine wave goes at the same speed.
It's only the transients in some sense that go at a different speed.
Now here's what happens for a finite difference approximation.
If you plug in this hypothesis to a finite difference approximation, you get a different dispersion relation.
And let's just do the second order one. If you take good old second order leapfrog.
In other words, the one with this stencil.
You plug that trial solution in to the formula for Leapfrog?
I won't go through the algebra. Here's what you find.
You get the sine squared of omega k over two equals k squared over x squared times the sine squared of c h over two.
So on the one hand, you have this formula for the actual PD.
On the other hand, you have this formula for the finite difference approximation.
And let me plot what that looks like. Now, the curves are not straight lines.
And in fact, they're periodic. So I'll show you the fundamental picture.
In wave number and frequency.
Here you have zero zero. Here you have wave number pi over H.
Here you have wave number minus pi over H. Here you have frequency minus pi over K and frequency pi over k.
And the dispersion curve. Looks like.
This. And then it keeps going periodically, so it keeps going forever.
So what's happening here is that in the middle of the picture, you get what you should have.
I'm sorry, I've only drawn one of them. I've just drawn this line.
There's also another set for this line in the middle of the picture that looks like that.
So it's doing the right thing when Omega and C are small, but when Omega and see are big, it's not correct.
So if Omega and see are small, that means many points per wavelength.
Part. So the wave is well resolved and everything behaves approximately as it should.
But when Omega and C are not small, the behaviour is completely different.
This limit out here would correspond to the smallest number of meaningful points you could possibly have.
If C is pi over H, then you're in a sawtooth situation, a completely minimally resolved wave if you like.
Whereas when C is small, you have a well resolved way.
Now I know some of you have studied dispersion relations and others haven't.
But all sorts of neat things turn up as soon as you get a curve like this.
Then it's called a dispersive equation and different wave numbers travel at different speeds.
And those speeds depend on this. And in fact, suppose we have a wave number here, for example, or suppose we have that wave number,
frequency pair, then two different numbers turn out to be important.
The. The slope of that line is equal to the negative of the group velocity.
So so-called. And the slope of the line that passes through the origin.
Is the negative of the phase velocity. Now some of you have studied those things and others haven't been true.
But the idea of a group velocity and a phase velocity is that if you have some kind of a wave package,
then the individual wave crests travel at the phase of velocity determined by this red slope, whereas the packet as a whole.
Travels at the group velocity, which is the slope of the green curve for a non dispersive equation, they're identical.
You don't get the distinction between phase and group velocity, but for a dispersion equation they're different and that's what's going on here.
So we have a wave which isn't as well resolved as we'd like it to be.
Roughly speaking, it has energy in some region like this, so a lot of the energy's near the origin and it does just the right thing.
But a bit of the energy is at a higher wave number and that energy you can see has a slope of the green line that's less than what it should be.
Instead of the slope one, it has slope point nine or something.
That's why the group velocity is less than one.
You see that these waves are going too slowly and the higher the wave number, which is to say the shorter the wavelength, the more too slowly they go.
So as you get fewer and fewer points per wavelength, shorter wavelength, you get greater and greater lags behind.
So, in fact, let's look at that one more time in motion.
If I run that again. Let's try a new configuration for fun.
I'll take 256 points and second order. So there's a bit of energy lagging behind.
You see, it's going slightly slower than it should. That wasn't dispersive enough to be interesting, was it?
Let's take fewer points. 128 Second order.
So you see that the higher wave numbers are going a little bit too slowly.
That's the dispersive effect in finite different schemes. So spectral methods are all about avoiding that by going to formula infinite order.
Oh, I should tell you my favourite reference on this stuff. It's by me written when I was your age.
So when I was a graduate student, I wrote a paper on this published in 1982.
About group velocity and finite differences. So spectral methods are the limit of infinite order.
And what does that mean? Well. I've told you that to devise finite difference formula formulas, you interpolate by a polynomial or something,
and then differentiate the interpolate in spectral methods, you use a global interpolation rather than a local one.
That's what it means. So whatever points we have, we're going to interpolate globally and then differentiate that interpolate.
And it will turn out that that method gives you something which has no particular finite order.
It's in some sense infinite order if you have a periodic domain.
Then your global in turbulence will be trigonometric.
And the methods that result are, of course, called Fourier spectral methods.
If you have a non periodic domain. Then the interpreters you use are just polynomials.
And that's called Chubby Chef Spectrum.
And broadly speaking, what you'll find is that the best method for accurate solutions of all sorts of problems,
if there is not very significant geometry, are the spectral methods because they're global.
They're not good at dealing with a complicated shape of a region.
They like boxes. They like intervals or boxes or cubes that they're good at.
If you're in an interval or a box or a cube,
then these disparate associations can be really magnificent and they give you the most efficient ways to solve all sorts of problems.
So let me tell you what we get in the limit of.
Our Leapfrog formula, if you like. But where?
We go to more and more points approaching infinity and that limit does make sense, can be derived in various ways.
So as the number of points goes to infinity, here are the coefficients we get.
We get a differentiation stencil if you like, which is one over h times.
Well, I already gave you the number in the middle. It's minus pi squared over two.
Three. And then the others are just rational numbers.
So here you get two and then minus two fourth and then 2/9 and minus to 16th and so on.
And on the left it's the same to minus two, fourth to ninth, minus 2/16.
And so. So this is the limit.
These are the numbers that the coefficients converge to.
And then as you take more and more points in your grid. And this is for second order differentiation.
Second order. Differentiation.
And there's no need for hand-waving and explaining what these numbers mean.
So let me say it precisely. If you go through this prescription of interpolating by a polynomial on a grid centred here a polynomial of degree two,
four, six, eight, 16 and so on, three, three, four, six, eight, ten, 12 and so on.
And in each case you differentiate the polynomial twice and evaluate the derivatives on the grid.
Then you'll get a sequence of numbers, the ones we wrote down here.
And those numbers, as the grid goes to infinity, will approach these numbers.
The best way to derive that fact is by Fourier analysis, which I might do later.
So this is the formula, infinite version of spectral differential of differentiation.
This picture would apply on an infinite grid. Now, of course, nobody really uses an infinite grid.
You might use a periodic grid. So let's write down what happens on a finite periodic grid.
So this was an infinite grid. Makes sense mathematically, but not computationally on a periodic grid.
The principle is the same, except now because of the periodic you interpolate by trigonometric polynomials.
So you have your domain and you sample it in 3.5.7 points,
always periodically interpolate by a trigonometric polynomial involving signs and cosine differentiate.
You get something analogous to this, but not quite the same. And here's what you get.
Turns out you get one half. Then in the middle, the number in the middle instead of two becomes the square of the cosecant of.
Sorry, I'm not writing down the right one. Apologies.
The number of the middle corresponding to pi squared over three turns out to be minus two.
Pi squared over three h squared minus a third.
And I hope I've done that right. Hope I haven't made an error.
And then this becomes cosecant squared of over two and this becomes minus cosecant
squared of of two over two and then cosecant squared of three over two and so.
And similarly on the left, cosecant squared over two minus cosecant squared of two, h over two and so on.
So this is the one that's elegant mathematically. This is the one that you would use in practice.
You would use a matrix with these coefficients to differentiate your periodic signal.
Let me pass around a book I have on this called Spectral Methods in MATLAB.
There are about a 15 books on spectral methods. This is the easiest.
It's not the most serious academically, but it's good fun.
So. In order to try this, we can set up a matrix which applies these coefficients to each entry of a vector.
And let's let's in fact look at the code, which is m 50.
So if you look at M 50 wave equation 48, what you see is that the vector C contains the numbers I've just written down.
And then the vector, the command D equals Toeplitz is where we construct the matrix.
So a toeplitz matrix. Is a matrix that constant along diagonal.
So the Toeplitz command in MATLAB will take the vector C and make that your first row and then it will just continue constant along these diagonals.
Because that's what you need for a periodic process.
It's translation invariant where disparate ising on a periodic grid.
So there you have a four point periodic grid. Here's an eight point periodic grid.
Here's the 16 point periodic grid.
Whether you disparate rise on that grid because everything is periodic, the process should be translation invariant.
And that means that the numbers in your matrix should be the same in every row, but just shifted.
Moreover, they wrap around. So the numbers here wrap around too, down there.
So it's in fact, a circular matrix.
Which means a toeplitz matrix where things wrap around.
So I'm being I'm just using words there. But what does it mean?
Mathematically, a hopeless matrix is one in which a i j only depends on the difference between the coefficients.
And you can write down a similar algebraic condition for that.
So the code makes a toeplitz matrix.
The reason it's actually more than that and circle is the reason the coefficients wrap around.
Has to do with the entries in that vector. The cosecant have that property that they have the appropriate periodicity.
So this code makes that Toeplitz matrix and then at every step it does leapfrog multiplying by that Toeplitz matrix.
So notice that as we increase the order.
The Sparsity changes in the order to case we have a tri diagonal matrix, and in the order for case we would have a penta diagonal matrix.
But as we go to order infinity, as it were, it's a dense matrix.
And that's a general property of spectral methods.
The matrices are dense, but they tend to be smaller in dimension than with finite differences or finite elements.
So this code just solves the same problem using this spectral technology.
Let's see what happens. It's called I'm 50.
64 points. And there's no dispersion at all.
It's doing exactly the right thing for every wave number.
So there's not much to say. You can see a little jagged miss in the curve.
That's only because 64 points are being correct, connected by straight lines.
If you interpolated them properly by a trigonometric interpolate, you wouldn't even see that jagged Ms.
So spectral methods can get away with sort of two points per wavelength instead of ten points per wavelength.
If I do it with more points, of course it'll look good again.
Same thing. Bit slower than before I notice because I'm multiplying by a dense matrix.
Okay. So that's the basis of spectral methods.
But in practice, they're often not always implemented by taking a 48 transform.
I've shown it to you as a matrix process, but the other way to do it is in 48 space where you act on the waves wave number by wave number.
So I'm going to show you how that works. And this is where the word spectral comes from.
So spectral analysis sounds like breaking something into a different frequency, different wave numbers.
This word was coined back around 1970 by Steve Orzag, as it happens.
So let's talk then about Fourier. Spectral different possession.
Via FFP. So by means of the foreign aid transfer.
I'll just remind you, I'm sure you all know the Fast Furious transform was made famous in 1965 by Julian Tookey.
It had actually been invented long before, about five different times.
The first time was by Carl Gauss in 1805.
So it's an old idea. The 51st Invented in 1805 and last invented in 1965.
There's a law called Stigler Law that says nothing is ever named after the actual inventor.
And certainly this is a good example of that. It's actually in Tukey in 65 who made it famous.
So the idea is to act on each individual wave number in the fashion that your PD wants to act on those wave numbers.
So we simply note that if we had, let's say, a pure sine wave.
So suppose we had e to the i j x for sound j to eyes the square root of minus one, and j is a number.
Then you prime of x. Well, it's a times j times you.
So when you look at different wave numbers, differentiation becomes a trivial operation.
Differentiating a sine wave is just multiplying it by a number.
So this complicated operation of differentiation become just a scalar multiplication in 48 space.
And that's the idea. Take a signal, convert it to space.
And then differentiation just becomes a lot of multiplications.
So if you have X is a superposition of exponentials.
Then let's write it as a sum. So suppose we write you of x equals a sum over indices j of sum coefficients capital use of j times e to the i j x.
Then. We can take the derivative thanks to linearity by differentiating each piece in the sun.
So the derivative becomes the sum of use of j.
And then i. J e to the i.
J x. And similarly you double prime would be the sum of you j times minus j square e to the i j x etc.
So that gives you a very slick way to perform differentiation via the 48 transform.
And in fact, Capital U is the discrete Fourier transform.
So you is the discrete Fourier transform.
Discrete for the transform of the.
Sampled values of little new. So I really should say.
Have you sampled on a grill? So this is the idea of FFP implementations of spectrum methods.
You do whatever you want to do in foyer space and then convert back to space space to find your solution.
So let's write it down and call it an algorithm if you like.
The algorithm for spectral differentiation is.
Given a vector you.
Compute its discrete eight transform. So it's discrete for transform capital u equals in standard notation 50 of you.
That's how MATLAB and many other codes would write it. So FFP is an algorithm.
DFT is a mathematical object computed by that algorithm.
So we say that the fast 48 transform is a fast algorithm for computing the discrete Fourier transform.
So you compute the discrete for a transform in this fast fashion.
And then. If I want to take two derivatives.
To be analogous to what I've done in the code so far. I would then multiply.
By the vector minus J Square.
So in other words, I'm constructing a new vector.
Call it capital W whose jth entry is minus j squared times the jth entry of capital U and then we just transform back.
And in MATLAB or many other languages that would be the command i f f t so little w equals i f f t of capital w the inverse discrete for a
transform and the forward for a discrete discrete for a transform are mathematically essentially the same with some constant factors modified.
So this is your algorithm. The amount of work is and log in.
As compared with N squared for multiplying by the matrix.
So that's the point of using the FFE you convert n squared to and log n so mathematically.
We have w equals a matrix behind you.
If we do this directly via the circulate matrix.
Without in any way exploiting its circle and structure then its o of and squared operations.
If we do it indirectly via the FFP its analogue in.
One of those historical oddities that always interests me.
This is why spectral methods were invented in 1965.
People invented the 52,
three discovered the 50 spectral methods came along within the next five years because people were dazzled by this great improvement.
Suddenly, methods involving Fourier transform became tractable on the computer.
But what's ironic historically is it turns out that for the sort of small grids that people were mostly using back then anyway,
the direct method is often as fast as the indirect one. If you're on a 64 by 64 grid, it may not make much difference.
Once you're up to 128 to 56, then the F 50 is a big win.
But in practice then people often use the four AA transform for spectral methods, but by no means always.
So let's take a look at the next code, which is 51.
So in M51 you can see as always, I'm doing things in time by a leapfrog democratisation.
But in space it's a spectral method. And now we're using the Fast 48 transform.
So I'm 51 and I'm 50. Are ISO morphic in terms of what they do, but the method they use is different and in particular,
you can see the time stepping at every step for loop at the bottom.
We take the five little you to get big you.
We multiply by a suitable vector to get capital w then we take little w equals the inverse transform.
You notice it says reel of 50. That's just to get rid of rounding errors.
So mathematically, the thing you're taking the real part of is already real.
But in practice it will have ten to the -16 times imaginary components.
So you strip them away. So let's see how that does.
Let's do em 50 again and time it. So if I take a greater say to 56.
One, two, three, four, five, six, seven, eight, nine, ten, 11, 12, 13, 14, 15, 16, 17, 18, 19.
Okay, and let's do em. 51. Which is the first version, same grid.
One, two, three, four, five, six, seven, eight, nine, ten, one, two.
So it wasn't even twice as fast. There are all sorts of reasons why that comparison isn't fair.
On grids of this size, actually, the FFP really is a win, but you can see that it may not be as spectacular as you first expect.
One reason the comparison isn't fair is that a lot of the time here is spent on the plot.
While we're talking about that, I want to draw to your attention perhaps the mathematically most interesting of the coffee table book handouts here.
So if you look at handout 19, the one way, way the question.
You'll see this is an idea that's of enormous importance in various applications.
The idea is to take a two dimensional dispersion relation or a three dimensional dispersion relation and chop it in half.
So on the left, in this diagram, you see circles which correspond to the dispersion relation of the two dimensional second order wave equation.
And then in some applications you'd only like semi-circle. So you'd like the waves to go to the right, but not to the left for various reasons.
Achieving that is mathematically, in some sense, impossible.
At least it can't be done with a partial differential equation.
You need a pseudo differential equation, but you can come very close to it with suitable approximations.
And this particular page is describing how you can derive those approximations to get waves that propagate just like ordinary waves to the east,
while not propagating at all to the west. Now I want to show you one more thing.
I wanted to mention to you that in Champ Fun, these 48 discrete positions are available.
And let me illustrate what that means. So suppose I take a periodic function, I'll say F equals R, let's define the function.
I'll say f, f equals. The E to the sign of pi times x.
Now, if I say F equals seven of F, then what's happening there is that the function is being represented as a chubby chef series,
which we haven't really talked about much. If I plot F, it will look periodic, of course, because it is periodic.
And the representation there that's being used is a Chevy Chef series.
But Chevy Fun also enables you to represent periodic functions by Fourier series, trigonometric series.
So suppose I said f trig equals.
And then I'd say have fun of the same thing and give it the flag trig.
Well, that's the same function. But now it really is a four year series.
So for example, if I plot the coefficients of F three.
That's hard to see. Let me do it like this. These are four year expansion coefficients for this function.
The middle number. Wave number. Wave number zero corresponds to the constant term.
And then to the right and the left. You have E, to the I, X and to the minus and so on.
The mathematically, the 48 series would go forever, but instead fun things always go to 16 digits.
So you can see that to 16 digits. It takes about 30 coefficients to represent this function.
And you can solve differential equations and so on. So Turbofan has all this Fourier stuff built into it.
Now, we should stop there because I need to hand back to the assignments.
I don't think they're in alphabetical order. So what I'll just do is spread them out on the floor up here and please come up and collect.
Remember, the lecture on Friday is an L two. Not in this round.