1
00:00:01,740 --> 00:00:05,460
As you can see, there are a bunch of handouts from the coffee table book.

2
00:00:05,820 --> 00:00:13,180
Let me first of all mention the one with the beautiful spots and scrolled picture of the great dot equation.

3
00:00:14,250 --> 00:00:18,240
A beautiful set of equations which you're all working with on Assignment three.

4
00:00:18,690 --> 00:00:21,060
So this will give you a little bit of the background.

5
00:00:21,690 --> 00:00:31,829
These equations have become one of people's favourite illustrations of how tuning a few parameters can change the qualitative behaviour of solutions.

6
00:00:31,830 --> 00:00:35,490
So you see, this picture shows different regions of parameter space.

7
00:00:35,940 --> 00:00:38,519
In many regions solutions are relatively boring,

8
00:00:38,520 --> 00:00:45,839
but then there are other regions where they have interesting structure and sometimes the structure looks like long rolls,

9
00:00:45,840 --> 00:00:51,660
and in other places it's more like spots. And in the assignment you'll be seeing some pictures like that.

10
00:00:53,100 --> 00:00:55,890
Let me also alert you to a few things here.

11
00:00:56,190 --> 00:01:02,670
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.

12
00:01:03,450 --> 00:01:07,649
We have an assignment due next Tuesday and then I will be handing back assignment

13
00:01:07,650 --> 00:01:10,889
to at the end of this lecture and I'll probably forget unless you remind me.

14
00:01:10,890 --> 00:01:12,480
So I trust you'll remind me.

15
00:01:14,340 --> 00:01:23,250
So today we're going to talk about higher order difference approximations and in particular the limit when the order goes to infinity,

16
00:01:23,430 --> 00:01:33,270
which is called spectral district ization. So 5.9 Fourier Spectral Description zation.

17
00:01:35,270 --> 00:01:39,680
And to talk about that, we're going to begin with finite differences to motivating.

18
00:01:49,430 --> 00:01:52,669
And I'm going to do everything today with the wave equation.

19
00:01:52,670 --> 00:02:04,309
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,

20
00:02:04,310 --> 00:02:09,020
one to the left and the right. And two of the handouts actually are about the wave equation.

21
00:02:09,020 --> 00:02:18,860
So let's glance at them. The one D Wave Equation handout describes how you can represent a general

22
00:02:18,860 --> 00:02:23,240
solution to the wave equation in terms of left going and right going solutions.

23
00:02:24,710 --> 00:02:29,660
I hope you're reading these coffee table handouts because they're very carefully put together.

24
00:02:30,620 --> 00:02:36,560
You're getting an exposure to all sorts of PD that is very hard to get elsewhere.

25
00:02:38,240 --> 00:02:45,710
The other one about the wave equation, the just with the title wave equation talks about how things change in higher dimensions.

26
00:02:46,010 --> 00:02:49,610
So of course you see a picture of a propagating wave there.

27
00:02:49,610 --> 00:02:54,769
You could imagine it's a gravitational wave, if you like, which would be correct, by the way,

28
00:02:54,770 --> 00:03:01,820
because the gravitational wave, by the time it is a few light years away from those black holes, might as well be linear.

29
00:03:02,090 --> 00:03:04,400
So it really is essentially the wave equation.

30
00:03:06,200 --> 00:03:14,270
One interesting mathematical fact about wave equations is that dimensions and dimensions are different from even dimensions,

31
00:03:14,720 --> 00:03:21,800
but the case end equals one is special. So in odd dimensions all information travels at speed exactly one.

32
00:03:22,190 --> 00:03:27,590
And it even dimensions. All information travels at speed less than or equal to one.

33
00:03:27,920 --> 00:03:34,880
Except that dimension one is special. It's like an even dimension where the information travels at speed less than or equal to one.

34
00:03:35,150 --> 00:03:40,670
That's called Huggins principal. So having fun reading those, please.

35
00:03:49,620 --> 00:03:56,700
Now what I want to do is immediately show you some simulations by finite difference methods of this wave equation.

36
00:03:57,360 --> 00:04:00,630
And we're going to do the second order.

37
00:04:00,660 --> 00:04:04,050
Fourth order and sixth order disorganisation.

38
00:04:04,620 --> 00:04:09,420
So the second order then. So let's get the leapfrog sort of representation.

39
00:04:09,690 --> 00:04:16,530
So that would be a second order formula. And this would be a fourth order formula.

40
00:04:18,060 --> 00:04:21,150
And this would be a sixth order formula.

41
00:04:24,270 --> 00:04:30,840
I hope those pictures are enough to tell you everything you need to figure out the equation.

42
00:04:31,170 --> 00:04:37,709
Let's get right down one of them. So in this case, what I mean is that V and plus one,

43
00:04:37,710 --> 00:04:49,440
minus two the and plus V and minus one all evaluated at the middle space point divided by k squared equals

44
00:04:49,440 --> 00:05:00,870
v j plus one minus to v j plus v j minus one all evaluated at the middle time step divided by h square.

45
00:05:02,720 --> 00:05:07,580
So that's the second order leapfrog style, discrete ization of the way equation.

46
00:05:08,000 --> 00:05:12,650
And in these cases, we do the same thing in time, but something different in space.

47
00:05:14,030 --> 00:05:23,090
And in this case, the coefficients turn out to be minus a half for thirds, minus five has.

48
00:05:24,760 --> 00:05:32,890
4/3. Minus a half. And in this case, those coefficients turn out to be 1/90.

49
00:05:34,580 --> 00:05:38,480
Minus 3/20 three has.

50
00:05:40,020 --> 00:05:43,340
-40 9/18.

51
00:05:44,460 --> 00:05:48,050
And three have. Uh.

52
00:05:48,340 --> 00:05:54,400
I forgot. My signs are on. Three has minus three 20th and 1/90.

53
00:05:57,650 --> 00:06:01,940
And you know where to find such numbers in the tables that I handed out from Feinberg.

54
00:06:02,270 --> 00:06:11,030
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.

55
00:06:12,620 --> 00:06:18,829
So today's lecture is really about what happens as you continue to higher and higher order,

56
00:06:18,830 --> 00:06:22,070
and in particular, making sense of the limit of infinite order.

57
00:06:22,460 --> 00:06:31,220
Just to give you an indication of what happens, the number two here and then five has and then 49,

58
00:06:31,220 --> 00:06:38,810
18th, those numbers in the limit of infinite order of accuracy become pi squared over three.

59
00:06:40,860 --> 00:06:51,590
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.

60
00:06:55,790 --> 00:07:03,980
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.

61
00:07:04,910 --> 00:07:13,600
Let's run it. So I'll say m 49 wave equation and it asks me how many points I want in my grid.

62
00:07:13,620 --> 00:07:19,040
So I'll say 64. And then it asks, Do I want the second, fourth or sixth order formula?

63
00:07:19,430 --> 00:07:26,510
So I'll say second. And then it computes a wave propagating to the right at speed one.

64
00:07:29,410 --> 00:07:35,360
That's pretty bad, isn't it? There's a lot of spurious affiliation there that shouldn't be there.

65
00:07:35,380 --> 00:07:39,400
Let's improve it by taking a higher order approximation.

66
00:07:40,900 --> 00:07:45,430
I'll take four. Four. That's better.

67
00:07:45,940 --> 00:07:49,750
There's still a little bit of dispersion being introduced by the numerical scheme.

68
00:07:51,220 --> 00:08:06,030
Let's try a sixth order. That's better, but there's still a bit of dispersion being introduced by the numerical scheme.

69
00:08:06,240 --> 00:08:11,640
So that little tail of oscillation should not be there. It ought to be just a cleanly propagating way.

70
00:08:11,970 --> 00:08:15,150
And I need to say a word about why it's moving to the right.

71
00:08:15,690 --> 00:08:19,800
The equation is completely symmetrical. Things can move left or right.

72
00:08:20,310 --> 00:08:25,050
But in this code, I've given initial data corresponding to moving to the right.

73
00:08:25,290 --> 00:08:32,310
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.

74
00:08:33,030 --> 00:08:36,059
So I set a value, you old, and then I set a value.

75
00:08:36,060 --> 00:08:42,870
You add the fact that that value you involved X minus K is what makes it move to the right.

76
00:08:43,140 --> 00:08:54,070
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.

77
00:08:54,310 --> 00:09:01,900
You have the spurious dispersion effect which, as you can imagine, could confuse you a great deal in an application.

78
00:09:02,440 --> 00:09:11,860
There's a surprise that comes up in wave simulations that linear problems are kind of harder than nonlinear problems because nonlinear waves,

79
00:09:12,340 --> 00:09:15,680
thanks to the non linearity, tend to be self-organising.

80
00:09:15,700 --> 00:09:25,870
A solid time, for example, keeps it shape. And numerically, that sometimes means that nonlinear simulations are less prone to this kind of problem.

81
00:09:26,200 --> 00:09:29,530
But a linear simulation, this sort of thing happens all the time.

82
00:09:30,970 --> 00:09:40,000
Let's refine the grid a little bit. So instead of 128 points, let's take that sorry, 64 and let's take 128.

83
00:09:40,510 --> 00:09:45,430
So with the second order. You see, things are better.

84
00:09:47,470 --> 00:09:50,890
But there are oscillations. Let's try the sixth order.

85
00:09:53,620 --> 00:10:02,460
128/6 order. That's looking very good now, isn't it?

86
00:10:03,150 --> 00:10:07,680
There's really hardly any visible evidence of the dispersion. There's a little bit there, but not too much.

87
00:10:08,010 --> 00:10:13,110
On the other hand, 128 is a lot of points to resolve a wave as boring as that.

88
00:10:13,680 --> 00:10:20,250
We have at least ten or 15 points in the way, which seems a lot of points to resolve a wave accurately.

89
00:10:22,020 --> 00:10:29,610
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.

90
00:10:32,600 --> 00:10:36,490
And indeed, let's go back and do a plot that shows it more clearly.

91
00:10:36,500 --> 00:10:41,120
So I'll do I'm 49 again, 128 points.

92
00:10:41,120 --> 00:10:57,060
But second one. So I want to say a word about that word and where it comes from.

93
00:11:02,370 --> 00:11:12,420
You remember that diffusion is when things decay, basically, and dispersion is when things are spread into a pulse of different wave numbers.

94
00:11:13,890 --> 00:11:21,480
And the way people analyse dispersion is with dispersion relations, which are all about tracking and how frequency depends on wave number.

95
00:11:21,780 --> 00:11:29,850
So the procedure throughout physics and mathematics, what you do is you look at solutions to the wave equation.

96
00:11:32,820 --> 00:11:39,750
Of the farm U of x and t equals an exponential or a sine wave.

97
00:11:39,750 --> 00:11:42,870
But mathematically, it's always simpler to write an exponential.

98
00:11:42,870 --> 00:11:48,540
So I'll say frequency times. Time plus wave number times space.

99
00:11:49,700 --> 00:11:53,120
So i omega t plus c x.

100
00:11:57,710 --> 00:12:06,970
And when people analyse dispersion, what they do is plug this into the equation and see how Omega relates to see.

101
00:12:08,900 --> 00:12:13,040
So let's make sure you know the words their frequency refers to time.

102
00:12:13,850 --> 00:12:19,070
It's the dual variable of time. And wave number refers to space.

103
00:12:19,850 --> 00:12:29,360
It's the dual variable of space. So if you plug that into an equation, you get a dispersion relation.

104
00:12:31,740 --> 00:12:36,390
If you plug it in to u t t equals u x x.

105
00:12:37,260 --> 00:12:52,260
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.

106
00:12:53,280 --> 00:12:57,690
So that's the dispersion relation for the wave equation.

107
00:13:04,780 --> 00:13:09,670
And from that dispersion relation, you can infer how different ways will propagate.

108
00:13:10,150 --> 00:13:18,700
Let's plot the dispersion relation. If this is wave number and this is frequency, then it has two branches.

109
00:13:18,700 --> 00:13:23,590
One of them looks like that and one of them looks like that.

110
00:13:24,790 --> 00:13:29,530
And from the slopes of those curves, you can figure out what the wave speeds are.

111
00:13:29,560 --> 00:13:32,590
So this has slope plus one.

112
00:13:35,290 --> 00:13:46,070
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.

113
00:13:49,790 --> 00:13:54,710
So they're in the dispersion relation. You can see the right going and the left going solutions to the wave equation.

114
00:13:55,790 --> 00:14:00,650
So there's no dispersion there when the dispersion relation is a straight line.

115
00:14:00,920 --> 00:14:05,930
Everything goes at the same speed and that's all there is to say. The wave equation doesn't do very much.

116
00:14:08,420 --> 00:14:16,850
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.

117
00:14:16,850 --> 00:14:20,690
It's only the transients in some sense that go at a different speed.

118
00:14:22,930 --> 00:14:27,970
Now here's what happens for a finite difference approximation.

119
00:14:30,830 --> 00:14:38,690
If you plug in this hypothesis to a finite difference approximation, you get a different dispersion relation.

120
00:14:42,290 --> 00:14:50,300
And let's just do the second order one. If you take good old second order leapfrog.

121
00:14:54,620 --> 00:14:57,890
In other words, the one with this stencil.

122
00:15:01,120 --> 00:15:06,219
You plug that trial solution in to the formula for Leapfrog?

123
00:15:06,220 --> 00:15:08,560
I won't go through the algebra. Here's what you find.

124
00:15:09,310 --> 00:15:26,770
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.

125
00:15:36,080 --> 00:15:41,150
So on the one hand, you have this formula for the actual PD.

126
00:15:41,660 --> 00:15:45,860
On the other hand, you have this formula for the finite difference approximation.

127
00:15:46,130 --> 00:15:53,060
And let me plot what that looks like. Now, the curves are not straight lines.

128
00:15:54,710 --> 00:15:59,120
And in fact, they're periodic. So I'll show you the fundamental picture.

129
00:16:01,410 --> 00:16:04,890
In wave number and frequency.

130
00:16:07,470 --> 00:16:12,660
Here you have zero zero. Here you have wave number pi over H.

131
00:16:13,470 --> 00:16:22,860
Here you have wave number minus pi over H. Here you have frequency minus pi over K and frequency pi over k.

132
00:16:25,480 --> 00:16:32,380
And the dispersion curve. Looks like.

133
00:16:33,670 --> 00:16:39,570
This. And then it keeps going periodically, so it keeps going forever.

134
00:16:50,890 --> 00:16:56,670
So what's happening here is that in the middle of the picture, you get what you should have.

135
00:16:56,680 --> 00:17:00,430
I'm sorry, I've only drawn one of them. I've just drawn this line.

136
00:17:00,430 --> 00:17:07,180
There's also another set for this line in the middle of the picture that looks like that.

137
00:17:07,690 --> 00:17:14,560
So it's doing the right thing when Omega and C are small, but when Omega and see are big, it's not correct.

138
00:17:14,800 --> 00:17:23,320
So if Omega and see are small, that means many points per wavelength.

139
00:17:27,740 --> 00:17:35,770
Part. So the wave is well resolved and everything behaves approximately as it should.

140
00:17:37,420 --> 00:17:42,850
But when Omega and C are not small, the behaviour is completely different.

141
00:17:43,210 --> 00:17:48,690
This limit out here would correspond to the smallest number of meaningful points you could possibly have.

142
00:17:48,700 --> 00:18:01,000
If C is pi over H, then you're in a sawtooth situation, a completely minimally resolved wave if you like.

143
00:18:02,170 --> 00:18:07,330
Whereas when C is small, you have a well resolved way.

144
00:18:10,790 --> 00:18:15,710
Now I know some of you have studied dispersion relations and others haven't.

145
00:18:16,160 --> 00:18:20,090
But all sorts of neat things turn up as soon as you get a curve like this.

146
00:18:20,090 --> 00:18:25,790
Then it's called a dispersive equation and different wave numbers travel at different speeds.

147
00:18:26,090 --> 00:18:34,610
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,

148
00:18:34,610 --> 00:18:40,490
frequency pair, then two different numbers turn out to be important.

149
00:18:42,470 --> 00:18:53,420
The. The slope of that line is equal to the negative of the group velocity.

150
00:18:53,420 --> 00:19:02,990
So so-called. And the slope of the line that passes through the origin.

151
00:19:06,740 --> 00:19:18,750
Is the negative of the phase velocity. Now some of you have studied those things and others haven't been true.

152
00:19:18,770 --> 00:19:24,290
But the idea of a group velocity and a phase velocity is that if you have some kind of a wave package,

153
00:19:25,070 --> 00:19:34,070
then the individual wave crests travel at the phase of velocity determined by this red slope, whereas the packet as a whole.

154
00:19:35,830 --> 00:19:43,000
Travels at the group velocity, which is the slope of the green curve for a non dispersive equation, they're identical.

155
00:19:43,000 --> 00:19:50,560
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.

156
00:19:51,430 --> 00:19:57,580
So we have a wave which isn't as well resolved as we'd like it to be.

157
00:19:58,240 --> 00:20:07,840
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.

158
00:20:07,840 --> 00:20:16,720
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.

159
00:20:16,960 --> 00:20:20,530
Instead of the slope one, it has slope point nine or something.

160
00:20:21,550 --> 00:20:24,610
That's why the group velocity is less than one.

161
00:20:24,970 --> 00:20:35,920
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.

162
00:20:36,640 --> 00:20:44,380
So as you get fewer and fewer points per wavelength, shorter wavelength, you get greater and greater lags behind.

163
00:20:45,970 --> 00:20:49,420
So, in fact, let's look at that one more time in motion.

164
00:20:53,880 --> 00:20:59,710
If I run that again. Let's try a new configuration for fun.

165
00:20:59,890 --> 00:21:09,200
I'll take 256 points and second order. So there's a bit of energy lagging behind.

166
00:21:10,580 --> 00:21:15,260
You see, it's going slightly slower than it should. That wasn't dispersive enough to be interesting, was it?

167
00:21:15,590 --> 00:21:21,050
Let's take fewer points. 128 Second order.

168
00:21:24,270 --> 00:21:28,410
So you see that the higher wave numbers are going a little bit too slowly.

169
00:21:28,680 --> 00:21:43,440
That's the dispersive effect in finite different schemes. So spectral methods are all about avoiding that by going to formula infinite order.

170
00:21:47,290 --> 00:21:52,450
Oh, I should tell you my favourite reference on this stuff. It's by me written when I was your age.

171
00:21:54,550 --> 00:21:59,050
So when I was a graduate student, I wrote a paper on this published in 1982.

172
00:22:01,730 --> 00:22:11,980
About group velocity and finite differences. So spectral methods are the limit of infinite order.

173
00:22:18,260 --> 00:22:29,860
And what does that mean? Well. I've told you that to devise finite difference formula formulas, you interpolate by a polynomial or something,

174
00:22:29,860 --> 00:22:36,309
and then differentiate the interpolate in spectral methods, you use a global interpolation rather than a local one.

175
00:22:36,310 --> 00:22:45,730
That's what it means. So whatever points we have, we're going to interpolate globally and then differentiate that interpolate.

176
00:22:45,970 --> 00:22:50,950
And it will turn out that that method gives you something which has no particular finite order.

177
00:22:50,950 --> 00:22:56,470
It's in some sense infinite order if you have a periodic domain.

178
00:22:59,530 --> 00:23:02,530
Then your global in turbulence will be trigonometric.

179
00:23:05,490 --> 00:23:09,600
And the methods that result are, of course, called Fourier spectral methods.

180
00:23:14,100 --> 00:23:27,030
If you have a non periodic domain. Then the interpreters you use are just polynomials.

181
00:23:32,140 --> 00:23:34,210
And that's called Chubby Chef Spectrum.

182
00:23:44,650 --> 00:23:51,970
And broadly speaking, what you'll find is that the best method for accurate solutions of all sorts of problems,

183
00:23:52,570 --> 00:23:58,570
if there is not very significant geometry, are the spectral methods because they're global.

184
00:23:59,140 --> 00:24:02,590
They're not good at dealing with a complicated shape of a region.

185
00:24:03,040 --> 00:24:07,090
They like boxes. They like intervals or boxes or cubes that they're good at.

186
00:24:07,360 --> 00:24:09,760
If you're in an interval or a box or a cube,

187
00:24:10,000 --> 00:24:17,020
then these disparate associations can be really magnificent and they give you the most efficient ways to solve all sorts of problems.

188
00:24:18,730 --> 00:24:25,150
So let me tell you what we get in the limit of.

189
00:24:27,570 --> 00:24:30,810
Our Leapfrog formula, if you like. But where?

190
00:24:32,750 --> 00:24:39,380
We go to more and more points approaching infinity and that limit does make sense, can be derived in various ways.

191
00:24:41,120 --> 00:24:47,720
So as the number of points goes to infinity, here are the coefficients we get.

192
00:24:48,830 --> 00:24:55,340
We get a differentiation stencil if you like, which is one over h times.

193
00:24:56,180 --> 00:25:00,590
Well, I already gave you the number in the middle. It's minus pi squared over two.

194
00:25:02,800 --> 00:25:08,580
Three. And then the others are just rational numbers.

195
00:25:08,910 --> 00:25:19,260
So here you get two and then minus two fourth and then 2/9 and minus to 16th and so on.

196
00:25:20,280 --> 00:25:28,020
And on the left it's the same to minus two, fourth to ninth, minus 2/16.

197
00:25:28,740 --> 00:25:34,480
And so. So this is the limit.

198
00:25:35,170 --> 00:25:38,200
These are the numbers that the coefficients converge to.

199
00:25:38,200 --> 00:25:43,480
And then as you take more and more points in your grid. And this is for second order differentiation.

200
00:25:45,040 --> 00:25:49,050
Second order. Differentiation.

201
00:25:55,650 --> 00:25:59,930
And there's no need for hand-waving and explaining what these numbers mean.

202
00:25:59,940 --> 00:26:11,760
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,

203
00:26:11,820 --> 00:26:17,790
four, six, eight, 16 and so on, three, three, four, six, eight, ten, 12 and so on.

204
00:26:18,510 --> 00:26:24,960
And in each case you differentiate the polynomial twice and evaluate the derivatives on the grid.

205
00:26:26,550 --> 00:26:31,440
Then you'll get a sequence of numbers, the ones we wrote down here.

206
00:26:33,180 --> 00:26:39,810
And those numbers, as the grid goes to infinity, will approach these numbers.

207
00:26:42,030 --> 00:26:47,310
The best way to derive that fact is by Fourier analysis, which I might do later.

208
00:26:48,990 --> 00:26:54,810
So this is the formula, infinite version of spectral differential of differentiation.

209
00:26:56,130 --> 00:27:00,690
This picture would apply on an infinite grid. Now, of course, nobody really uses an infinite grid.

210
00:27:00,870 --> 00:27:07,860
You might use a periodic grid. So let's write down what happens on a finite periodic grid.

211
00:27:09,270 --> 00:27:17,140
So this was an infinite grid. Makes sense mathematically, but not computationally on a periodic grid.

212
00:27:22,120 --> 00:27:28,780
The principle is the same, except now because of the periodic you interpolate by trigonometric polynomials.

213
00:27:29,140 --> 00:27:34,990
So you have your domain and you sample it in 3.5.7 points,

214
00:27:34,990 --> 00:27:41,920
always periodically interpolate by a trigonometric polynomial involving signs and cosine differentiate.

215
00:27:42,340 --> 00:27:46,300
You get something analogous to this, but not quite the same. And here's what you get.

216
00:27:47,920 --> 00:28:01,100
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.

217
00:28:02,330 --> 00:28:10,650
Sorry, I'm not writing down the right one. Apologies.

218
00:28:11,230 --> 00:28:18,630
The number of the middle corresponding to pi squared over three turns out to be minus two.

219
00:28:20,070 --> 00:28:25,380
Pi squared over three h squared minus a third.

220
00:28:27,650 --> 00:28:36,050
And I hope I've done that right. Hope I haven't made an error.

221
00:28:37,730 --> 00:28:45,590
And then this becomes cosecant squared of over two and this becomes minus cosecant

222
00:28:45,590 --> 00:28:54,680
squared of of two over two and then cosecant squared of three over two and so.

223
00:28:57,840 --> 00:29:06,000
And similarly on the left, cosecant squared over two minus cosecant squared of two, h over two and so on.

224
00:29:11,530 --> 00:29:16,240
So this is the one that's elegant mathematically. This is the one that you would use in practice.

225
00:29:16,660 --> 00:29:23,740
You would use a matrix with these coefficients to differentiate your periodic signal.

226
00:29:27,090 --> 00:29:31,050
Let me pass around a book I have on this called Spectral Methods in MATLAB.

227
00:29:31,260 --> 00:29:34,920
There are about a 15 books on spectral methods. This is the easiest.

228
00:29:35,970 --> 00:29:39,090
It's not the most serious academically, but it's good fun.

229
00:29:47,740 --> 00:30:01,450
So. In order to try this, we can set up a matrix which applies these coefficients to each entry of a vector.

230
00:30:04,530 --> 00:30:08,100
And let's let's in fact look at the code, which is m 50.

231
00:30:09,120 --> 00:30:19,800
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.

232
00:30:23,360 --> 00:30:29,899
And then the vector, the command D equals Toeplitz is where we construct the matrix.

233
00:30:29,900 --> 00:30:38,490
So a toeplitz matrix. Is a matrix that constant along diagonal.

234
00:30:39,800 --> 00:30:49,940
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.

235
00:30:52,460 --> 00:30:55,610
Because that's what you need for a periodic process.

236
00:30:55,610 --> 00:30:59,840
It's translation invariant where disparate ising on a periodic grid.

237
00:31:02,000 --> 00:31:06,110
So there you have a four point periodic grid. Here's an eight point periodic grid.

238
00:31:06,620 --> 00:31:08,600
Here's the 16 point periodic grid.

239
00:31:11,320 --> 00:31:18,400
Whether you disparate rise on that grid because everything is periodic, the process should be translation invariant.

240
00:31:18,610 --> 00:31:23,470
And that means that the numbers in your matrix should be the same in every row, but just shifted.

241
00:31:23,680 --> 00:31:29,980
Moreover, they wrap around. So the numbers here wrap around too, down there.

242
00:31:30,460 --> 00:31:33,640
So it's in fact, a circular matrix.

243
00:31:37,020 --> 00:31:41,040
Which means a toeplitz matrix where things wrap around.

244
00:31:44,750 --> 00:31:49,700
So I'm being I'm just using words there. But what does it mean?

245
00:31:49,730 --> 00:31:59,720
Mathematically, a hopeless matrix is one in which a i j only depends on the difference between the coefficients.

246
00:32:01,970 --> 00:32:05,000
And you can write down a similar algebraic condition for that.

247
00:32:06,380 --> 00:32:09,470
So the code makes a toeplitz matrix.

248
00:32:10,010 --> 00:32:14,060
The reason it's actually more than that and circle is the reason the coefficients wrap around.

249
00:32:14,360 --> 00:32:21,590
Has to do with the entries in that vector. The cosecant have that property that they have the appropriate periodicity.

250
00:32:21,950 --> 00:32:29,660
So this code makes that Toeplitz matrix and then at every step it does leapfrog multiplying by that Toeplitz matrix.

251
00:32:30,470 --> 00:32:36,110
So notice that as we increase the order.

252
00:32:37,470 --> 00:32:46,350
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.

253
00:32:46,590 --> 00:32:50,910
But as we go to order infinity, as it were, it's a dense matrix.

254
00:32:53,940 --> 00:32:56,220
And that's a general property of spectral methods.

255
00:32:56,550 --> 00:33:02,760
The matrices are dense, but they tend to be smaller in dimension than with finite differences or finite elements.

256
00:33:03,870 --> 00:33:08,580
So this code just solves the same problem using this spectral technology.

257
00:33:08,820 --> 00:33:13,930
Let's see what happens. It's called I'm 50.

258
00:33:17,180 --> 00:33:23,889
64 points. And there's no dispersion at all.

259
00:33:23,890 --> 00:33:26,920
It's doing exactly the right thing for every wave number.

260
00:33:32,210 --> 00:33:37,190
So there's not much to say. You can see a little jagged miss in the curve.

261
00:33:37,340 --> 00:33:41,150
That's only because 64 points are being correct, connected by straight lines.

262
00:33:41,420 --> 00:33:46,460
If you interpolated them properly by a trigonometric interpolate, you wouldn't even see that jagged Ms.

263
00:33:46,820 --> 00:33:52,550
So spectral methods can get away with sort of two points per wavelength instead of ten points per wavelength.

264
00:33:53,050 --> 00:33:56,330
If I do it with more points, of course it'll look good again.

265
00:33:58,110 --> 00:34:04,260
Same thing. Bit slower than before I notice because I'm multiplying by a dense matrix.

266
00:34:13,260 --> 00:34:16,500
Okay. So that's the basis of spectral methods.

267
00:34:18,330 --> 00:34:23,610
But in practice, they're often not always implemented by taking a 48 transform.

268
00:34:23,940 --> 00:34:32,640
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.

269
00:34:33,600 --> 00:35:32,890
So I'm going to show you how that works. And this is where the word spectral comes from.

270
00:35:33,400 --> 00:35:38,710
So spectral analysis sounds like breaking something into a different frequency, different wave numbers.

271
00:35:39,160 --> 00:35:43,180
This word was coined back around 1970 by Steve Orzag, as it happens.

272
00:35:44,590 --> 00:35:50,590
So let's talk then about Fourier. Spectral different possession.

273
00:35:55,320 --> 00:36:04,340
Via FFP. So by means of the foreign aid transfer.

274
00:36:09,190 --> 00:36:16,840
I'll just remind you, I'm sure you all know the Fast Furious transform was made famous in 1965 by Julian Tookey.

275
00:36:17,230 --> 00:36:20,770
It had actually been invented long before, about five different times.

276
00:36:21,100 --> 00:36:24,280
The first time was by Carl Gauss in 1805.

277
00:36:24,700 --> 00:36:34,420
So it's an old idea. The 51st Invented in 1805 and last invented in 1965.

278
00:36:35,230 --> 00:36:40,809
There's a law called Stigler Law that says nothing is ever named after the actual inventor.

279
00:36:40,810 --> 00:36:46,330
And certainly this is a good example of that. It's actually in Tukey in 65 who made it famous.

280
00:36:49,470 --> 00:36:58,650
So the idea is to act on each individual wave number in the fashion that your PD wants to act on those wave numbers.

281
00:36:58,950 --> 00:37:06,000
So we simply note that if we had, let's say, a pure sine wave.

282
00:37:06,360 --> 00:37:15,180
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.

283
00:37:18,950 --> 00:37:26,230
Then you prime of x. Well, it's a times j times you.

284
00:37:34,210 --> 00:37:40,120
So when you look at different wave numbers, differentiation becomes a trivial operation.

285
00:37:40,360 --> 00:37:44,200
Differentiating a sine wave is just multiplying it by a number.

286
00:37:44,560 --> 00:37:51,790
So this complicated operation of differentiation become just a scalar multiplication in 48 space.

287
00:37:52,000 --> 00:37:56,130
And that's the idea. Take a signal, convert it to space.

288
00:37:56,470 --> 00:37:59,620
And then differentiation just becomes a lot of multiplications.

289
00:37:59,950 --> 00:38:05,560
So if you have X is a superposition of exponentials.

290
00:38:09,480 --> 00:38:28,080
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.

291
00:38:31,270 --> 00:38:41,110
Then. We can take the derivative thanks to linearity by differentiating each piece in the sun.

292
00:38:41,810 --> 00:38:45,709
So the derivative becomes the sum of use of j.

293
00:38:45,710 --> 00:38:49,310
And then i. J e to the i.

294
00:38:49,610 --> 00:39:06,800
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.

295
00:39:10,420 --> 00:39:16,630
So that gives you a very slick way to perform differentiation via the 48 transform.

296
00:39:16,840 --> 00:39:21,360
And in fact, Capital U is the discrete Fourier transform.

297
00:39:21,370 --> 00:39:27,190
So you is the discrete Fourier transform.

298
00:39:27,520 --> 00:39:30,730
Discrete for the transform of the.

299
00:39:34,030 --> 00:39:41,890
Sampled values of little new. So I really should say.

300
00:39:41,900 --> 00:39:54,510
Have you sampled on a grill? So this is the idea of FFP implementations of spectrum methods.

301
00:39:55,500 --> 00:40:02,490
You do whatever you want to do in foyer space and then convert back to space space to find your solution.

302
00:40:02,670 --> 00:40:07,920
So let's write it down and call it an algorithm if you like.

303
00:40:09,240 --> 00:40:12,780
The algorithm for spectral differentiation is.

304
00:40:15,390 --> 00:40:18,450
Given a vector you.

305
00:40:21,570 --> 00:40:36,430
Compute its discrete eight transform. So it's discrete for transform capital u equals in standard notation 50 of you.

306
00:40:36,550 --> 00:40:42,700
That's how MATLAB and many other codes would write it. So FFP is an algorithm.

307
00:40:43,510 --> 00:40:47,200
DFT is a mathematical object computed by that algorithm.

308
00:40:47,530 --> 00:40:54,610
So we say that the fast 48 transform is a fast algorithm for computing the discrete Fourier transform.

309
00:40:56,290 --> 00:41:00,010
So you compute the discrete for a transform in this fast fashion.

310
00:41:00,460 --> 00:41:06,960
And then. If I want to take two derivatives.

311
00:41:08,950 --> 00:41:13,570
To be analogous to what I've done in the code so far. I would then multiply.

312
00:41:17,080 --> 00:41:20,170
By the vector minus J Square.

313
00:41:22,530 --> 00:41:26,399
So in other words, I'm constructing a new vector.

314
00:41:26,400 --> 00:41:38,100
Call it capital W whose jth entry is minus j squared times the jth entry of capital U and then we just transform back.

315
00:41:46,970 --> 00:41:59,710
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

316
00:41:59,720 --> 00:42:08,420
transform and the forward for a discrete discrete for a transform are mathematically essentially the same with some constant factors modified.

317
00:42:09,410 --> 00:42:13,400
So this is your algorithm. The amount of work is and log in.

318
00:42:18,510 --> 00:42:22,230
As compared with N squared for multiplying by the matrix.

319
00:42:22,680 --> 00:42:30,030
So that's the point of using the FFE you convert n squared to and log n so mathematically.

320
00:42:33,400 --> 00:42:38,000
We have w equals a matrix behind you.

321
00:42:40,150 --> 00:42:44,740
If we do this directly via the circulate matrix.

322
00:42:47,500 --> 00:42:53,680
Without in any way exploiting its circle and structure then its o of and squared operations.

323
00:42:55,700 --> 00:42:59,300
If we do it indirectly via the FFP its analogue in.

324
00:43:10,990 --> 00:43:14,890
One of those historical oddities that always interests me.

325
00:43:17,230 --> 00:43:20,610
This is why spectral methods were invented in 1965.

326
00:43:20,620 --> 00:43:22,590
People invented the 52,

327
00:43:22,630 --> 00:43:29,920
three discovered the 50 spectral methods came along within the next five years because people were dazzled by this great improvement.

328
00:43:30,160 --> 00:43:34,780
Suddenly, methods involving Fourier transform became tractable on the computer.

329
00:43:35,230 --> 00:43:42,400
But what's ironic historically is it turns out that for the sort of small grids that people were mostly using back then anyway,

330
00:43:43,150 --> 00:43:49,990
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.

331
00:43:50,230 --> 00:43:54,160
Once you're up to 128 to 56, then the F 50 is a big win.

332
00:43:55,060 --> 00:44:01,690
But in practice then people often use the four AA transform for spectral methods, but by no means always.

333
00:44:04,450 --> 00:44:09,490
So let's take a look at the next code, which is 51.

334
00:44:14,190 --> 00:44:20,230
So in M51 you can see as always, I'm doing things in time by a leapfrog democratisation.

335
00:44:20,490 --> 00:44:24,330
But in space it's a spectral method. And now we're using the Fast 48 transform.

336
00:44:24,600 --> 00:44:32,669
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,

337
00:44:32,670 --> 00:44:37,350
you can see the time stepping at every step for loop at the bottom.

338
00:44:37,710 --> 00:44:40,800
We take the five little you to get big you.

339
00:44:41,400 --> 00:44:49,020
We multiply by a suitable vector to get capital w then we take little w equals the inverse transform.

340
00:44:49,230 --> 00:44:53,760
You notice it says reel of 50. That's just to get rid of rounding errors.

341
00:44:54,000 --> 00:44:58,620
So mathematically, the thing you're taking the real part of is already real.

342
00:44:58,980 --> 00:45:03,000
But in practice it will have ten to the -16 times imaginary components.

343
00:45:03,270 --> 00:45:09,160
So you strip them away. So let's see how that does.

344
00:45:14,820 --> 00:45:21,740
Let's do em 50 again and time it. So if I take a greater say to 56.

345
00:45:24,060 --> 00:45:38,060
One, two, three, four, five, six, seven, eight, nine, ten, 11, 12, 13, 14, 15, 16, 17, 18, 19.

346
00:45:38,250 --> 00:45:45,620
Okay, and let's do em. 51. Which is the first version, same grid.

347
00:45:48,210 --> 00:45:56,640
One, two, three, four, five, six, seven, eight, nine, ten, one, two.

348
00:45:57,150 --> 00:46:02,460
So it wasn't even twice as fast. There are all sorts of reasons why that comparison isn't fair.

349
00:46:03,300 --> 00:46:11,310
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.

350
00:46:11,880 --> 00:46:15,840
One reason the comparison isn't fair is that a lot of the time here is spent on the plot.

351
00:46:18,450 --> 00:46:26,550
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.

352
00:46:26,910 --> 00:46:31,170
So if you look at handout 19, the one way, way the question.

353
00:46:35,000 --> 00:46:39,320
You'll see this is an idea that's of enormous importance in various applications.

354
00:46:39,680 --> 00:46:47,480
The idea is to take a two dimensional dispersion relation or a three dimensional dispersion relation and chop it in half.

355
00:46:47,810 --> 00:46:56,090
So on the left, in this diagram, you see circles which correspond to the dispersion relation of the two dimensional second order wave equation.

356
00:46:56,720 --> 00:47:03,920
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.

357
00:47:04,880 --> 00:47:08,030
Achieving that is mathematically, in some sense, impossible.

358
00:47:08,210 --> 00:47:10,880
At least it can't be done with a partial differential equation.

359
00:47:10,880 --> 00:47:17,270
You need a pseudo differential equation, but you can come very close to it with suitable approximations.

360
00:47:17,510 --> 00:47:26,930
And this particular page is describing how you can derive those approximations to get waves that propagate just like ordinary waves to the east,

361
00:47:27,410 --> 00:47:33,850
while not propagating at all to the west. Now I want to show you one more thing.

362
00:47:38,410 --> 00:47:44,860
I wanted to mention to you that in Champ Fun, these 48 discrete positions are available.

363
00:47:44,860 --> 00:47:54,519
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.

364
00:47:54,520 --> 00:48:04,040
I'll say f, f equals. The E to the sign of pi times x.

365
00:48:05,600 --> 00:48:16,819
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,

366
00:48:16,820 --> 00:48:24,530
which we haven't really talked about much. If I plot F, it will look periodic, of course, because it is periodic.

367
00:48:27,090 --> 00:48:30,680
And the representation there that's being used is a Chevy Chef series.

368
00:48:31,230 --> 00:48:37,380
But Chevy Fun also enables you to represent periodic functions by Fourier series, trigonometric series.

369
00:48:37,680 --> 00:48:41,909
So suppose I said f trig equals.

370
00:48:41,910 --> 00:48:46,200
And then I'd say have fun of the same thing and give it the flag trig.

371
00:48:47,670 --> 00:48:53,770
Well, that's the same function. But now it really is a four year series.

372
00:48:53,770 --> 00:48:57,310
So for example, if I plot the coefficients of F three.

373
00:48:59,850 --> 00:49:15,470
That's hard to see. Let me do it like this. These are four year expansion coefficients for this function.

374
00:49:15,650 --> 00:49:19,970
The middle number. Wave number. Wave number zero corresponds to the constant term.

375
00:49:20,240 --> 00:49:25,700
And then to the right and the left. You have E, to the I, X and to the minus and so on.

376
00:49:26,540 --> 00:49:32,480
The mathematically, the 48 series would go forever, but instead fun things always go to 16 digits.

377
00:49:32,750 --> 00:49:39,350
So you can see that to 16 digits. It takes about 30 coefficients to represent this function.

378
00:49:40,310 --> 00:49:45,650
And you can solve differential equations and so on. So Turbofan has all this Fourier stuff built into it.

379
00:49:46,670 --> 00:49:49,880
Now, we should stop there because I need to hand back to the assignments.

380
00:49:50,660 --> 00:49:57,110
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.

381
00:49:57,560 --> 00:50:01,460
Remember, the lecture on Friday is an L two. Not in this round.