1
00:00:00,150 --> 00:00:07,380
As you know, I find it beautiful how pervasive partial differential equations are throughout the scientific enterprise.

2
00:00:08,370 --> 00:00:15,210
And this is certainly the theme of the remainder of the course, exploring finite differences and other methods for solving pdes.

3
00:00:15,900 --> 00:00:21,990
Let me begin by passing around, actually for the second time, a pair of very nice books on this subject.

4
00:00:22,000 --> 00:00:29,430
So a few lectures ago I passed around a stack of books, two of which were also about PD ease.

5
00:00:29,430 --> 00:00:31,080
And I want to remind you of those two.

6
00:00:31,230 --> 00:00:38,340
One is by Aria Isserlis called a first course in the numerical analysis of differential equations, and the other is by Randy Leveque.

7
00:00:38,340 --> 00:00:41,730
Finite Difference Methods for Ordinary and Partial Differential Equations.

8
00:00:41,940 --> 00:00:49,140
These are both fantastic books by leaders in the field, which will tell you, of course, much more than we cover in these lectures.

9
00:00:52,320 --> 00:01:01,350
Now, whenever you talk about a numerical solution of pdes questions of stability and instability, numerical stability and instability come up.

10
00:01:01,590 --> 00:01:05,640
And that's what drives a lot of what we do. And it's pretty interesting.

11
00:01:05,880 --> 00:01:09,630
So the next section will call numerical instability.

12
00:01:10,530 --> 00:01:17,970
And today, you'll see how in a very simple way, it controls the kind of algorithms that work for Pdes.

13
00:01:19,110 --> 00:01:22,140
So let's see. Numerical instability.

14
00:01:25,110 --> 00:01:34,410
So to remind you, last lecture, we looked at the heat equation, which is you t equals you x a u x x.

15
00:01:36,060 --> 00:01:40,710
So the heat equation in one d would be u t equals u x x.

16
00:01:41,010 --> 00:01:45,960
And we disparities that and showed how. Sometimes you need a very small time step.

17
00:01:47,280 --> 00:01:54,300
So now let's consider the fourth order version of that u t equals.

18
00:01:56,590 --> 00:02:03,550
You x x x x. Except the way I've written, that would be an exponentially ill posed problem.

19
00:02:03,670 --> 00:02:10,730
You need a minus sign. The way it works with diffusion is four orders two, six, ten, 14, and so on.

20
00:02:10,750 --> 00:02:15,250
A plus sign is diffusive and four orders for eight, 12 and so on.

21
00:02:15,250 --> 00:02:19,480
A minus sign is diffusive. As you can see by plugging in a 48 mode.

22
00:02:19,870 --> 00:02:24,970
So this is the physically right approach to a fourth order diffusion problem.

23
00:02:26,500 --> 00:02:35,020
Now, let's suppose we want to discard ties that by a finite difference approximation, we do the obvious thing which will work.

24
00:02:37,770 --> 00:02:45,440
And well, we note that you x x x x is the same as u x x x x.

25
00:02:46,350 --> 00:02:51,630
So we can take the second difference of the second difference to get a pretty good approximation to that.

26
00:02:51,870 --> 00:02:56,910
So let's approximate you x x x x by.

27
00:03:01,870 --> 00:03:11,649
Something divided by eight squared and that's something we'll be three terms, one term minus twice another term plus a third term.

28
00:03:11,650 --> 00:03:24,460
And each of those will be a finite difference. So this will be the J plus two minus to the J plus one plus the J might be j over h squared.

29
00:03:25,870 --> 00:03:33,730
And this will be the j plus one minus 2vj plus the j minus one over h squared.

30
00:03:35,050 --> 00:03:44,200
And this will be the J minus to the j minus one plus v j minus two over h square.

31
00:03:46,390 --> 00:03:49,750
So we've taken the obvious. Second difference of the obvious.

32
00:03:49,750 --> 00:03:53,260
Second difference to get the obvious fourth difference.

33
00:03:54,550 --> 00:04:00,520
And that's a perfectly good way to describe the fourth derivative if you want to write it in matrix notation.

34
00:04:03,150 --> 00:04:09,209
You could think of that as if we're going to apply that to compute something.

35
00:04:09,210 --> 00:04:17,550
We're going to end up with a new time step, a vector of values being equal to a matrix times, the current time step.

36
00:04:18,300 --> 00:04:22,980
So this would be. The Matrix representation.

37
00:04:24,500 --> 00:04:28,070
Of an explicit as opposed to implicit.

38
00:04:29,270 --> 00:04:33,830
Finite difference approximation to our equation.

39
00:04:40,750 --> 00:04:43,990
So we've just arrived. A perfectly good method that will work.

40
00:04:44,200 --> 00:04:48,670
And now let's run it. You can see the code is called M 39.

41
00:04:49,510 --> 00:05:01,450
Fourth order diffusion. Just glancing at the code briefly.

42
00:05:02,350 --> 00:05:09,399
I use the usual trick in these demo codes that I actually do use a matrix which would be crazy if it were a dense matrix,

43
00:05:09,400 --> 00:05:13,750
but I use a sparse matrix, so it's not crazy. We set up a grid.

44
00:05:13,780 --> 00:05:23,049
The space size H is 1/40 and then you can see that we go through this time stepping with that time step of K and the code

45
00:05:23,050 --> 00:05:29,290
is set up so that the user inputs the time step so we can see the behaviour for different choices of the time step.

46
00:05:31,230 --> 00:05:34,860
And of course all the work is done just in the wire loop.

47
00:05:35,160 --> 00:05:40,590
You equals eight times you at every time step. So let's run that.

48
00:05:41,520 --> 00:05:46,080
I'll say m 39, fourth order diffusion.

49
00:05:47,730 --> 00:05:51,390
There's the initial distribution.

50
00:05:51,660 --> 00:05:55,950
And then it asks me what time step I want. So I'm going to put in a very small time step.

51
00:05:56,160 --> 00:06:00,750
I'll put in one E minus eight.

52
00:06:01,560 --> 00:06:12,800
And let me make that visible. Okay.

53
00:06:15,650 --> 00:06:19,760
So we put in a timestamp of ten to the minus eight. And here's what happens.

54
00:06:22,020 --> 00:06:27,240
The usual principle of diffusion that it's interesting a little bit and then it gets boring pretty fast.

55
00:06:27,510 --> 00:06:28,979
So there's fourth order diffusion.

56
00:06:28,980 --> 00:06:34,170
You can see if you're interested in the physics of these things that it is a little different from the second order case.

57
00:06:34,170 --> 00:06:42,120
The magnet tenacity properties are rather different. But anyway, there we have a perfectly reasonable solution of sorts.

58
00:06:42,270 --> 00:06:46,590
Let's increase the time step. So I'll change that to.

59
00:06:52,890 --> 00:07:00,150
Four times ten to the minus eight. Same thing.

60
00:07:00,330 --> 00:07:06,840
Time step four times bigger. Now I'll change it to 4.8, eight times ten to the minus eight.

61
00:07:08,250 --> 00:07:12,720
4.88 D minus eight. Little wiggly.

62
00:07:13,680 --> 00:07:21,990
This is actually just barely stable, maybe. But of course, you can see that there are oscillations that aren't going away very quickly.

63
00:07:22,560 --> 00:07:33,780
So that's very marginally stable. And of course, I'm now going to do it with 4.90 times ten to the minus eight, and the Wiggles will grow.

64
00:07:35,720 --> 00:07:40,010
So that's about as marginal an instability as you could ask for in a linear problem.

65
00:07:40,520 --> 00:07:44,540
And instability will grow at a good, clean exponential rate.

66
00:07:44,570 --> 00:07:46,550
There you have a good clean exponential growth.

67
00:07:49,040 --> 00:07:55,490
And of course, if I take a bigger time step, it'll only get worse if I took ten to the minus seven, for example.

68
00:07:55,850 --> 00:08:03,570
Then it blows up. So we see there, again, the kind of instability that one sees in all sorts of schemes.

69
00:08:03,960 --> 00:08:10,670
What's new is that because it's a fourth order problem. The critical time step is even smaller than before.

70
00:08:10,680 --> 00:08:17,520
It's very, very small. Ten to the minus eight. And what's going on there is that this is a stiff problem.

71
00:08:18,270 --> 00:08:21,809
The fourth derivative couples together, things at all frequencies.

72
00:08:21,810 --> 00:08:25,710
And there's the high frequencies cannot be forgotten.

73
00:08:25,740 --> 00:08:30,600
They're moving so quickly that it's causing the problems we associate with stiffness.

74
00:08:31,650 --> 00:08:35,860
So something has to be done in order to improve that before we do anything.

75
00:08:35,880 --> 00:08:52,480
Let's analyse it and understand what's going on. So let's do the standard for a trick that people use to analyse instability.

76
00:08:54,810 --> 00:08:59,190
So this idea basically was introduced by von Neumann in the 1940s.

77
00:08:59,460 --> 00:09:05,880
Computers were invented. Los Alamos was discovering all sorts of things, sort of during and then after the war.

78
00:09:06,090 --> 00:09:08,940
And von Neumann was a key figure.

79
00:09:11,910 --> 00:09:20,330
He was really one of the great geniuses of the 20th century, and he knew everything about pure and applied mathematics.

80
00:09:20,340 --> 00:09:23,560
It's amazing how he did the whole range, from pure to applied.

81
00:09:23,790 --> 00:09:31,290
He knew Fourier analysis very well, so when he saw instabilities like that, he of course knew that the thing to do was for analysis.

82
00:09:31,530 --> 00:09:36,990
Here's how we do it. We're going to plug in a Fourier mode and see what happens to it as time passes.

83
00:09:38,370 --> 00:09:45,830
So. Let's do our trial solution and we'll say, what if?

84
00:09:47,250 --> 00:09:50,460
At one time? Steppe We have a pure forage mode.

85
00:09:51,840 --> 00:09:58,050
So suppose that step end the solution looks like this.

86
00:09:58,410 --> 00:10:10,160
The J at step n. I'll call it a sine wave, but I'll actually write down a complex exponential suppose at time at step and it's e to the i.

87
00:10:10,850 --> 00:10:17,150
C. Times. SJ Which is another way of saying E to the.

88
00:10:18,240 --> 00:10:27,140
I see. JH So we're doing very standard for analysis in space here.

89
00:10:27,440 --> 00:10:32,960
C is the wave number. Which tells you it's the spatial frequency.

90
00:10:37,000 --> 00:10:43,660
And we want to see how different wave numbers evolve as we take steps with the finite, different scheme.

91
00:10:43,930 --> 00:10:48,670
So we plug this into the finite difference scheme and we see what the and plus one looks like.

92
00:10:54,050 --> 00:10:57,650
So at step end plus one.

93
00:11:00,080 --> 00:11:06,290
Well, the miracle of Fourier analysis is that this sine wave will be multiplied by a constant.

94
00:11:06,560 --> 00:11:13,950
What we'll find is that the end plus one will be simply equal to some number, depending on the wave number.

95
00:11:13,970 --> 00:11:18,230
So we call it JFC Times the and.

96
00:11:20,400 --> 00:11:23,860
And G of C is known as the amplification factor.

97
00:11:25,290 --> 00:11:30,720
So stability analysis in this business is about studying amplification factors.

98
00:11:31,890 --> 00:11:43,170
And I need a new red pen. Yeah.

99
00:11:48,400 --> 00:11:51,610
Okay. Well, so let's figure out what that amplification factor is.

100
00:11:51,730 --> 00:12:00,190
We just plug the sine wave into the finite difference formula and here's what we get g of T is equal to.

101
00:12:00,430 --> 00:12:03,610
And I'll just write it down once here. Practised with this stuff.

102
00:12:03,610 --> 00:12:07,270
It's amazingly quick. We're going to plug this into there.

103
00:12:08,260 --> 00:12:15,940
And what we find is that our galaxy is equal to one minus two times step over the fourth power of the space step.

104
00:12:15,940 --> 00:12:21,340
And then we get various terms from those various things.

105
00:12:21,520 --> 00:12:39,580
The v j plus two term gives you an e to the 2ij2i see h the v j plus one terms give you a minus four e to the i c h the vee-jay term in the middle.

106
00:12:40,250 --> 00:12:44,740
It gives you a plus six e to the zero.

107
00:12:46,630 --> 00:12:52,930
The the j minus one term gives you a minus four e to the minus.

108
00:12:52,930 --> 00:12:59,290
I see h and finally, the leftmost term gives you a plus one.

109
00:13:00,670 --> 00:13:04,120
E to the minus two i. C.

110
00:13:04,120 --> 00:13:13,660
H. I'm a little embarrassed that I never wrote down up above.

111
00:13:15,820 --> 00:13:24,310
So let me quickly do that. I never simplified that formula, which would have given me.

112
00:13:26,780 --> 00:13:41,630
The J plus two minus for the J plus one plus six the J minus for v j minus one plus v j minus two.

113
00:13:41,990 --> 00:13:51,700
So all of these collapse into that. And it's when I plug the sine wave into there.

114
00:13:52,270 --> 00:13:55,390
Here are these numbers. One minus four, six minus four, one.

115
00:14:01,970 --> 00:14:11,450
So there's the amplification factor as a function of the wave number C four for the fourth order disk realisation.

116
00:14:12,560 --> 00:14:20,600
And again, of course we could well what we do is we look at all possible wave numbers and see how they behave.

117
00:14:20,780 --> 00:14:24,710
Are they amplified, are they decayed? So let's simplify this a bit.

118
00:14:28,140 --> 00:14:43,170
These things turn into one minus K over H to the fourth times, the fourth power of E to the i c h over two minus e to the minus i.

119
00:14:43,710 --> 00:14:49,560
C h over to. That's just algebra.

120
00:14:50,640 --> 00:14:56,550
And if you look at these things, you see a sine wave, right? It's an E to the eye, something minus E to the minus eye something.

121
00:14:56,760 --> 00:15:03,120
So this in turn is one minus K over H to the fourth times.

122
00:15:03,120 --> 00:15:10,110
The sign of c h over two squared to the fourth.

123
00:15:15,250 --> 00:15:19,570
So we have the fourth power of a sine wave dependence on the wave number.

124
00:15:20,020 --> 00:15:25,600
And of course, the fourth power is always a positive number. So we have one minus a positive number.

125
00:15:26,440 --> 00:15:33,850
And the question for stability is, is that thing less than one in which case will get decay or is it bigger than one?

126
00:15:34,000 --> 00:15:37,030
In which case we'll get instability and growth.

127
00:15:38,930 --> 00:15:44,270
So that's the key question of von Neumann analysis is this.

128
00:15:46,200 --> 00:15:53,530
Less than one in absolute value. For all wave numbers.

129
00:15:58,390 --> 00:16:05,480
If he's good, if no bad. Well, that's pretty easy to figure out.

130
00:16:05,840 --> 00:16:10,180
The sign as ranges over numbers between minus one and one.

131
00:16:10,190 --> 00:16:13,430
So the fourth Power Rangers, our numbers between zero and one.

132
00:16:13,820 --> 00:16:17,660
So this will range over numbers to the left of one.

133
00:16:20,000 --> 00:16:28,100
On the real axis, there's one. As sea varies, that will always take values one or two, the left.

134
00:16:28,100 --> 00:16:35,210
And how far to the left will it go? It will go to one minus K over H to the fourth.

135
00:16:36,660 --> 00:16:45,970
So there is this interval of the real axis. That will be the amplification factor for different wave numbers.

136
00:16:46,780 --> 00:16:50,080
We want to know, is that interval always less than one an absolute value?

137
00:16:50,320 --> 00:16:58,960
So the question becomes is, is it the case that one minus K over to the fourth is bigger than minus one?

138
00:17:00,730 --> 00:17:04,690
That is our stability condition based on von Neumann analysis.

139
00:17:07,250 --> 00:17:22,640
So in other words gets to simplify that is to bigger than 16 K over to the fourth i.e. is K less than or equal to H to the fourth over eight.

140
00:17:24,620 --> 00:17:29,790
That is our von Neumann stability condition. But.

141
00:17:36,490 --> 00:17:39,600
16. Where did the 16 come from?

142
00:17:39,960 --> 00:17:43,560
There should have been a two to the fourth. I forgot the two to the fourth, didn't I?

143
00:17:43,590 --> 00:17:48,010
Thank you. So this is to high times the sign.

144
00:17:48,030 --> 00:17:55,380
It's not just the sign. It's two high times the sign. Thank you. Two are well spotted.

145
00:17:56,880 --> 00:18:05,660
So hopefully when you work that out. So of course I've done this very quickly.

146
00:18:05,990 --> 00:18:13,370
But in outline, what you can see is a very simple procedure for analysing these irregular grids and linear problems.

147
00:18:13,730 --> 00:18:16,730
Plug in for your mode to see if they can grow or not.

148
00:18:17,090 --> 00:18:20,550
You reach a conclusion, and that's a pretty pessimistic conclusion.

149
00:18:20,570 --> 00:18:25,520
It says that if you have the space step, you have to shrink the time step by a factor of 16.

150
00:18:25,910 --> 00:18:31,820
So you're going to need very small time steps, which, of course, is exactly what we found with the experiment.

151
00:18:32,150 --> 00:18:40,970
In fact. If H equals 1/40, which is what we had in the code.

152
00:18:44,060 --> 00:18:56,240
Then what you find is that the stability restriction is less than 4.88, three times ten to the minus eight.

153
00:18:58,800 --> 00:19:02,520
So 4.88 is stable. 4.90 is unstable.

154
00:19:06,860 --> 00:19:12,080
So we have a stiff problem here. The high frequencies are giving us trouble in some sense.

155
00:19:12,230 --> 00:19:18,050
And actually, next week, we're going to talk in a general way about how one addresses the stiffness of pdes.

156
00:19:18,890 --> 00:19:27,920
For the moment, let's talk about the basic tool that we have for that which is implicit as opposed to explicit, different thing.

157
00:19:32,640 --> 00:19:42,540
Just like he's. We want to use implicit ness to allow ourselves to be less constrained in the time step.

158
00:19:43,170 --> 00:19:47,280
So implicit one the finite differences.

159
00:19:52,510 --> 00:19:56,950
And the basic thing that you start by doing is exactly what we did with ODS.

160
00:19:57,220 --> 00:20:00,610
You use a backward rather than a forward formula.

161
00:20:01,000 --> 00:20:10,900
Now, I didn't draw the stencil, in fact, but I could have done that when we wrote down our finite difference formula.

162
00:20:11,680 --> 00:20:13,150
Where did we write down our finite?

163
00:20:13,210 --> 00:20:26,140
I guess we never quite wrote it down, but I hope you understood that we were coupling together five points at step RN and computing from them.

164
00:20:27,290 --> 00:20:32,090
A point at step and plus one. So that's the stencil of this explicit formula.

165
00:20:33,800 --> 00:20:38,810
Well, to get around the stability restraint,

166
00:20:38,810 --> 00:20:50,560
we want to use an implicit formula that will change that picture to find a different thing in the nude in the future, coupled with the present.

167
00:20:50,570 --> 00:21:06,190
That will be our new stencil. And we can immediately write down the formula.

168
00:21:06,220 --> 00:21:20,320
Let's do that. So our implicit backward Euler style disorganisation of the fourth order problem will be v, j plus v and plus one equals v.

169
00:21:20,320 --> 00:21:29,190
N. And then we have minus because it's minus the fourth derivative K over H to the fourth.

170
00:21:30,510 --> 00:21:46,290
And then we have these five terms, the J plus two minus for the J plus one plus six, the J minus for the J, minus one plus the J minus two.

171
00:21:47,910 --> 00:21:56,220
The explicit formula had nn, but the implicit formula will have N plus one and plus one and so on.

172
00:22:01,780 --> 00:22:09,860
So that's the. Formula for the obvious implicit descriptors in the backward Euler style disorganisation.

173
00:22:10,040 --> 00:22:11,660
To analyse it, we do the same thing.

174
00:22:11,900 --> 00:22:22,460
We plug in a R and that's of a four way mode, a sine wave, and we see by what constant is it multiplied as you go from one step to the next.

175
00:22:22,730 --> 00:22:29,750
And I won't spell that out, but you get the very same result, except that the sign terms are in the denominator rather than the numerator.

176
00:22:29,990 --> 00:22:36,300
You find that the amplification factor. Is now equal to one over.

177
00:22:38,330 --> 00:22:43,970
One plus 16 K over H to the fourth.

178
00:22:45,520 --> 00:22:50,320
Times. The sign of c h over to.

179
00:22:51,630 --> 00:22:56,620
To the fourth. Hope I've got that right. You've got me nervous now.

180
00:22:57,400 --> 00:23:04,960
All right. What did we get before? Yes.

181
00:23:04,980 --> 00:23:12,180
So although I neglected the 16, as you pointed out, this is we before we had that in the numerator and now we've got it in the denominator.

182
00:23:13,800 --> 00:23:17,610
But the good news is this is less than one no matter what.

183
00:23:19,050 --> 00:23:22,620
For all see, regardless of K.

184
00:23:32,670 --> 00:23:39,330
In fact, the amplification factor, as you can see, are always one or bigger sorry, one, one or smaller.

185
00:23:39,540 --> 00:23:43,540
And so nothing can ever blow up.

186
00:23:43,560 --> 00:23:49,320
This will be a stable formula for any timestamp, and the phrase people use is unconditionally stable.

187
00:23:56,120 --> 00:24:01,849
No conditions on the time stamp. So of course, that doesn't mean you can use any time stamp,

188
00:24:01,850 --> 00:24:06,290
because if your time step is too big, it'll still be inaccurate, but at least it won't blow up.

189
00:24:06,800 --> 00:24:13,430
And we know the sort of general theme of numerical computing with differential equations is that to get convergence,

190
00:24:13,430 --> 00:24:21,980
you need accuracy and you need stability. Well, this means that the only thing we have to check is the accuracy to get a good solution.

191
00:24:24,910 --> 00:24:32,260
So of course we have a code for that. If you look at M40, you can see it's just like M 39.

192
00:24:33,400 --> 00:24:40,990
And the only change is that we've changed a times U in the main loop to be backslash U in the main loop.

193
00:24:41,710 --> 00:24:46,570
Now, with that high level of abstraction, it's a change of about two characters.

194
00:24:47,380 --> 00:24:51,700
Of course, B is marginally different, but not very different. The signs are just a little bit different.

195
00:24:53,050 --> 00:24:57,310
So we've turned matrix multiplication into solving a system of equations.

196
00:24:57,490 --> 00:25:04,330
And that's doable because in this one dimensional problem, it's a penta diagonal system of equation.

197
00:25:05,020 --> 00:25:15,850
So this thing here. We have V and plus one times B equals VM.

198
00:25:16,780 --> 00:25:20,860
So we're solving a system of equations to get the new value.

199
00:25:21,160 --> 00:25:24,790
But it's a painted diagonal system of equations.

200
00:25:26,670 --> 00:25:35,550
Five non-zero diagonals, all the others zero and MATLAB is of course clever enough to do a good algorithm for that.

201
00:25:35,820 --> 00:25:39,780
It's o of n operations for an n by n point a diagonal matrix.

202
00:25:40,080 --> 00:25:46,590
So this very simple sparsity pattern is easily taken advantage of to get a fast solution.

203
00:25:49,200 --> 00:25:52,440
So we run the code. And of course, I hope you know what we'll see.

204
00:25:52,650 --> 00:26:09,570
Nothing will ever blow up. So this one is called M 40 implicit.

205
00:26:12,530 --> 00:26:20,210
And again, I can choose the time step. So if I say ten to the minus eight, the solution will be almost identical to what we saw before.

206
00:26:21,170 --> 00:26:26,210
And the speed the speed of computation isn't much different. It's not that hard to solve the system of equations.

207
00:26:28,770 --> 00:26:33,659
You can sense as that evolves that it's an absurdly small time step.

208
00:26:33,660 --> 00:26:37,680
It doesn't make sense to be doing all this work to get from one time step to the next.

209
00:26:40,810 --> 00:26:44,920
If we increase that, we can go as far as we want.

210
00:26:45,400 --> 00:26:50,670
So, you know, pick something ten to the minus five. So now we have a thousand times bigger time step.

211
00:26:51,820 --> 00:26:57,580
Now you can see how the thing really evolves for a longer time, which we effectively couldn't see at all before.

212
00:27:00,110 --> 00:27:04,040
So in terms of stiffness, remember the general definition of stiffness?

213
00:27:07,010 --> 00:27:12,440
Is that a problem is if it contains time scales much shorter than the ones you are tracking.

214
00:27:12,830 --> 00:27:19,610
So is this problem stiff? Well, to be precise, at the very beginning, it's not when you begin the calculation.

215
00:27:23,350 --> 00:27:30,400
Because that is not a smooth curve. For the first few milliseconds, there's something going on, on a very small timescale step.

216
00:27:30,640 --> 00:27:38,200
So for the first few billionth of a time unit, you would say that this is not a stiff PD.

217
00:27:38,440 --> 00:27:45,700
But then after that, things are smooth. And so the equation has these fast time scales in it that are irrelevant to the thing you tracking.

218
00:27:45,940 --> 00:27:49,820
So it quickly becomes a stiff problem. So let's just do that again.

219
00:27:49,840 --> 00:27:56,920
If I say ten to the minus five again. You can see very quickly, it's gotten to the point.

220
00:27:58,860 --> 00:28:03,390
Where the fast timescales have nothing to do with the evolution of the curve.

221
00:28:04,140 --> 00:28:07,290
And of course, it's unconditionally stable. So I could take point one.

222
00:28:07,980 --> 00:28:15,420
Yeah, it works. So we have this just absolutely huge difference between explicit and implicit.

223
00:28:26,170 --> 00:28:29,620
Now everything becomes much more fun when you have nonlinearity there.

224
00:28:29,950 --> 00:28:35,770
And it turns out the same ideas are at the heart of what is a practical way to cope with nonlinearity.

225
00:28:37,930 --> 00:28:41,620
And we're going to talk about that a bit now with an example and then more next week.

226
00:28:42,640 --> 00:28:47,800
So I want to tell you about one of my favourite PD e's, the motto Sufficiency Equation.

227
00:28:59,200 --> 00:29:15,069
This is a prototype of a non-linear equation with stiff behaviour and the equation is u t equals minus u x x minus you.

228
00:29:15,070 --> 00:29:22,610
X x x x. Minus you squared over two x.

229
00:29:25,610 --> 00:29:33,500
So of course it's a one dimensional nonlinear PD and there's a page from the coffee table book there, which I hope you'll read when you come to it.

230
00:29:34,550 --> 00:29:41,330
To remind you, the coffee table book has about 40 of these nicely finished pages and you can find them online at my website.

231
00:29:41,330 --> 00:29:51,320
Just click on books. So this equation has three terms on the right and everyone has a different physics.

232
00:29:51,560 --> 00:29:59,480
And the interplay of those three bits of physics is what makes it so beautiful, just to say in words what those three things are.

233
00:30:01,600 --> 00:30:06,370
This is negative linear diffusion because of the minus sign.

234
00:30:09,550 --> 00:30:15,310
So if you had just that equation, you two equals minus x x.

235
00:30:15,310 --> 00:30:18,520
That's called the backward heat equation. It's completely opposed.

236
00:30:18,910 --> 00:30:22,870
Things just blow up. And the higher the wave numbers, the faster they blow up.

237
00:30:24,380 --> 00:30:27,620
This one. Has positive diffuse.

238
00:30:35,600 --> 00:30:43,160
And now let's ignoring the nonlinear term. Let's think about what would happen to a Fourier mode plugged into that equation.

239
00:30:43,370 --> 00:30:49,220
Well, this act to increase it exponentially and this acts to decrease it exponentially.

240
00:30:49,460 --> 00:30:56,370
So let's try the trial solution for the nonlinear part.

241
00:30:57,370 --> 00:31:11,830
All right. Sorry for the linear part of the equation. Suppose our trial solution is that you is equal to a four year mode, which I'll say is e to the.

242
00:31:12,080 --> 00:31:19,130
I see x. And I'll give it some sigma time t.

243
00:31:19,970 --> 00:31:27,020
I'm not putting an eye there because we're going to get not wavelike behaviour but explosively increasing or decreasing behaviour.

244
00:31:27,290 --> 00:31:37,969
If you plug that into these two terms, well the minus u double x term we have Ike C squared,

245
00:31:37,970 --> 00:31:41,480
so we get minus x squared and then with another minus sign.

246
00:31:41,480 --> 00:31:48,300
So we get plus c squared. So we find that. You get?

247
00:31:50,830 --> 00:31:56,890
You at s t. Well, you get that.

248
00:31:58,520 --> 00:32:02,840
Sigma is related to see by this equation.

249
00:32:06,190 --> 00:32:11,530
And this is a kind of a dispersion relation or a diffusion relation would be a better term,

250
00:32:11,530 --> 00:32:16,210
but nobody ever speaks of diffusion relations at the dispersion relation.

251
00:32:18,430 --> 00:32:22,300
For the behaviour of the linear part of the promoter of efficiency equation.

252
00:32:22,720 --> 00:32:31,390
It tells you that small wave numbers where this term dominates have a positive sigma so they'll grow.

253
00:32:31,810 --> 00:32:36,310
And large wave numbers where this term dominates have a negative sigma, so they'll decay.

254
00:32:36,550 --> 00:32:43,600
And in fact, if I plot this as a function of C, you'll get this kind of behaviour.

255
00:32:43,600 --> 00:32:46,990
When C equals one, you get a crossover.

256
00:32:49,270 --> 00:32:58,360
Between a positive sigma and a negative sigma. So this is the regime of negative diffusion and this is the regime of positive diffusion.

257
00:32:58,630 --> 00:33:07,000
The high wave numbers diffused. So if we just had that linear equation without the nonlinear term, it would be opposed.

258
00:33:07,690 --> 00:33:10,750
So solutions, technically speaking, wouldn't really exist.

259
00:33:10,870 --> 00:33:18,190
But anyway, the behaviour in question would be that the low wave numbers would simply explode and the high wave numbers would decay.

260
00:33:18,430 --> 00:33:23,170
So you would see bigger and bigger, fairly smooth waves.

261
00:33:25,470 --> 00:33:36,930
However. This term has the beautiful effect of moving energy from the low wave numbers to the high wave numbers.

262
00:33:38,130 --> 00:33:40,170
So this is the nonlinear term.

263
00:33:46,000 --> 00:33:54,370
What it does is it takes these relatively smooth components and stiffens them up, like the Burger's equation, which we mentioned earlier.

264
00:33:54,820 --> 00:34:02,430
And I'm not going to try to write that down mathematically, but in effect, it takes low, smooth energy and turns it into less smooth energy.

265
00:34:02,440 --> 00:34:06,100
So it's moving energy from this part of the picture to that part of the picture.

266
00:34:06,340 --> 00:34:12,490
And because of that, you end up with a well posed PD with wonderfully complicated behaviour.

267
00:34:12,760 --> 00:34:19,690
The you get the structures that are fairly smooth because this is where the energy is amplified.

268
00:34:19,900 --> 00:34:27,550
But those structures don't blow up because as they get bigger, energy is siphoned off into the high wave numbers where it decays.

269
00:34:27,760 --> 00:34:33,100
So it's a kind of an engine moving energy around, growing it in one region, decaying it in another.

270
00:34:34,950 --> 00:34:44,000
Now to solve that numerically. This is an equation you don't particularly care about, although it's pretty,

271
00:34:44,010 --> 00:34:50,340
but it really is a prototype of a lot of computational science with pdes to solve it numerically.

272
00:34:51,310 --> 00:34:57,170
You need to treat the linear terms implicitly and the nonlinear terms explicitly.

273
00:34:57,190 --> 00:34:59,500
You don't need to, but that's the best way to do it.

274
00:34:59,980 --> 00:35:07,660
And that's almost the central, most important practical fact about how people solve pdes like this.

275
00:35:08,950 --> 00:35:16,210
And let me summarise why it's important to do the linear terms implicitly and the nonlinear terms explicitly.

276
00:35:20,490 --> 00:35:24,930
So the good strategy for this and many other problems.

277
00:35:37,340 --> 00:35:44,300
Is that somehow or other? You want to do the linear terms implicitly and the nonlinear terms explicitly.

278
00:35:50,050 --> 00:35:51,850
And why is that? Well.

279
00:35:58,980 --> 00:36:07,440
Because it's the linear terms that are causing the stability constraint, this very high frequency and this fairly high frequency.

280
00:36:07,470 --> 00:36:13,800
Those are the ones that would force you to have an absurdly small timestamp if you used an explicit formula.

281
00:36:14,070 --> 00:36:24,060
So the reason the linear terms should be implicit is to avoid a very tight stability restriction.

282
00:36:28,880 --> 00:36:35,480
A very tight limit on the timestamp. Why are the nonlinear terms explicit?

283
00:36:35,840 --> 00:36:41,330
Because whenever you have an implicit formula, you have to solve some equations.

284
00:36:41,510 --> 00:36:47,230
If they're linear equations, that's just linear algebra. But if they're nonlinear equations, that's more than linear algebra.

285
00:36:47,240 --> 00:36:51,180
You have to do a Newton iteration or something crazy. That can be hard.

286
00:36:51,530 --> 00:37:00,890
So explicit for the nonlinear terms to avoid having to solve a nonlinear problem at each step.

287
00:37:15,600 --> 00:37:22,350
Now, the reason this strategy is so fruitful is that for all sorts of problems, the linear terms are the stiff ones.

288
00:37:22,380 --> 00:37:27,300
They're the ones with the high derivatives, with the fast timescales.

289
00:37:27,840 --> 00:37:31,140
The nonlinear terms are operate on a slower timescale.

290
00:37:31,380 --> 00:37:37,350
In principle, it could be the other way around. If you had the nonlinear terms with the high derivatives, then you'd be in big trouble.

291
00:37:37,350 --> 00:37:42,990
Then this strategy would get you now where you basically have to be implicit and do a lot of work at each time step.

292
00:37:43,290 --> 00:37:54,390
But luckily a lot of problems are in this category where the stiff terms are linear and it's kind of amazing how easy it is,

293
00:37:54,780 --> 00:37:57,990
at least for this one D problem to take advantage of that.

294
00:37:58,350 --> 00:38:02,220
So what we're simply going to do is use a stencil.

295
00:38:05,020 --> 00:38:13,500
That's implicit. For the fourth order term and the second order term.

296
00:38:20,860 --> 00:38:24,460
But the nonlinear thing will not be at that level at all.

297
00:38:26,460 --> 00:38:32,350
And we in feed in, if you like. The nonlinear term.

298
00:38:34,940 --> 00:38:47,450
At times step n. So this, of course, as always, is the end plus one part of the stencil, and that's the time stamp and part of the stencil.

299
00:38:47,750 --> 00:38:52,639
So you blend these together and let's just look at the code.

300
00:38:52,640 --> 00:38:59,240
It's so simple what you need to do. So the code is called M 41 Courier Moto Z MUSZYNSKI Take a look at that.

301
00:39:00,530 --> 00:39:05,870
You'll see that it's got the same flavour as all these last few codes I've been showing you.

302
00:39:05,870 --> 00:39:11,570
It sets up a grid. It sets up some matrices. And then the key line in the code.

303
00:39:14,100 --> 00:39:23,870
In the wire loop. Is this one that says you equals B backslash.

304
00:39:25,130 --> 00:39:29,390
You minus a times you squared over two.

305
00:39:37,030 --> 00:39:39,280
Now that is doing just what I mentioned.

306
00:39:39,970 --> 00:39:48,610
The solution of a system of equations involving that point, a diagonal matrix is doing the second and fourth order linear turns.

307
00:39:51,080 --> 00:39:58,970
And then the A matrix which is not being solved, it gets multiplying by something is of course the nonlinear term.

308
00:40:03,720 --> 00:40:11,670
So given a vector you of values at time step in, we can evaluate that and then we can do a backslash.

309
00:40:11,670 --> 00:40:18,360
And we've got our discrete approximation that the next time step and this works beautifully.

310
00:40:18,360 --> 00:40:36,010
And let me show you the picture. The reason it's I like this equation so much is that the solutions are chaotic.

311
00:40:37,420 --> 00:40:41,350
So I'll say M. 41.

312
00:40:43,780 --> 00:40:49,890
Oops. Komodo said Wyszynski. So that's the initial condition.

313
00:40:50,100 --> 00:40:54,360
And it asks me, what time do I want to go to? Well, let's make it a nice big time.

314
00:40:54,360 --> 00:41:00,590
I'll say a thousand. So as I go out type return to see solution.

315
00:41:00,600 --> 00:41:06,450
Okay, so it's evolving and the initial condition is doing complicated stuff.

316
00:41:08,200 --> 00:41:11,860
And then you see this beautiful, seemingly random motion.

317
00:41:16,370 --> 00:41:21,810
And notice. The interplay of randomness with tremendous structure.

318
00:41:22,080 --> 00:41:27,090
So, of course, that looks random, but the structure is that it's clearly got a well-defined wavelength,

319
00:41:27,090 --> 00:41:34,020
it has a preferred wavelength, and that will be determined by this curve.

320
00:41:34,050 --> 00:41:42,570
The preferred wavelength is the one that's amplified the most. So energy gets gets increased at that fairly low wave number, fairly smooth wave.

321
00:41:43,860 --> 00:41:50,880
But it can't get amplified too much because the nonlinear term then sends it off to the higher wave numbers where it gets diffused away.

322
00:41:52,200 --> 00:41:56,669
And. People have proved that these solutions truly are chaotic.

323
00:41:56,670 --> 00:42:01,350
There's a strange attractor. All the requisite features of chaos are there.

324
00:42:02,890 --> 00:42:06,640
Including the exponential dependence on initial conditions.

325
00:42:06,850 --> 00:42:11,620
So that picture is not in any sense. Well, it's not in the literal sense, correct?

326
00:42:11,890 --> 00:42:18,459
At times 626, this is not the right curve because little errors along the way will have been

327
00:42:18,460 --> 00:42:22,900
amplified exponentially according to the up and off exponent for this chaotic problem.

328
00:42:23,230 --> 00:42:25,809
So none of these curves are literally correct.

329
00:42:25,810 --> 00:42:32,740
After the first time, a few time units, we have that chaotic situation where the overall behaviour is correct,

330
00:42:33,070 --> 00:42:36,760
but the precise detail at a precise time is not correct.

331
00:42:39,400 --> 00:42:43,420
I can just stare at that stuff all day, but luckily we're only going to 1000.

332
00:42:53,370 --> 00:42:57,959
Let's, however, run it a couple more times just for a brief time interval.

333
00:42:57,960 --> 00:43:02,190
So suppose I run that again. If I guess, go to time ten, for example.

334
00:43:02,850 --> 00:43:08,130
What you see is that the hump kind of increases and then it begins to wiggle.

335
00:43:08,730 --> 00:43:14,850
So the increase is this linear effect of energy at low wave numbers being amplified.

336
00:43:15,060 --> 00:43:21,090
But then the beginning of the Wiggles is this non-linear effect of it being transferred to higher wave numbers.

337
00:43:21,330 --> 00:43:25,680
So wiggling that's that's the Burger's equation, the steepening up of a front.

338
00:43:26,370 --> 00:43:34,219
Let's go a little further. Suppose I took it to term 20. So it grows and then wiggles.

339
00:43:34,220 --> 00:43:39,060
Come in and then. You're beginning to sense complexity.

340
00:43:39,090 --> 00:43:44,530
And of course, then it slowly gets chaotic. Now I want to show you.

341
00:43:46,180 --> 00:43:53,440
The latest thing the club fund team has been doing, or rather our Adrian Martinelli has been doing, which is related to problems like this.

342
00:43:53,440 --> 00:43:57,040
I'm going to show you how you can solve this problem. And Jeb, fun.

343
00:43:58,470 --> 00:44:05,760
So let's. Well, the code there is m 41 seven, but I'm actually going to type it in so I can talk about it.

344
00:44:07,980 --> 00:44:12,090
So the code we're speaking of is called spin.

345
00:44:14,170 --> 00:44:18,730
And the reason it's called spin is that that stands for a stiff PD integrator.

346
00:44:18,940 --> 00:44:24,370
It's not currently in the standard telephone, but you can get it by going to the branch feature spin up.

347
00:44:24,370 --> 00:44:27,310
So the it tells you here what to do if you want to play with it.

348
00:44:29,080 --> 00:44:39,070
So what you do is you declare what interval you're working on and the interval for this problem is 16 pi times minus one one.

349
00:44:39,700 --> 00:44:46,829
So that's our domain. And then in fun style, you make a suitable operator.

350
00:44:46,830 --> 00:44:51,060
We call it a spin up. This PDP operator defined on that domain.

351
00:44:52,730 --> 00:45:00,380
And it has a linear part and a nonlinear part. So the linear part would be defined.

352
00:45:03,240 --> 00:45:12,360
As an anonymous function in MATLAB, the usual fashion, and then in the standard seven fun style, we have these overloads of MATLAB commands.

353
00:45:12,780 --> 00:45:22,640
So. The linear part is an anonymous function that has the negative of the second derivative and the negative of the fourth derivative.

354
00:45:25,120 --> 00:45:33,489
And then the nonlinear part. Has to do that you squared effect.

355
00:45:33,490 --> 00:45:40,370
So that will be an anonymous function. If anyone sees me type something wrong, please tell me.

356
00:45:42,520 --> 00:45:45,950
What did they do in part? Oh yes.

357
00:45:45,950 --> 00:45:49,070
And annoyingly, we have fights on the Champs team about this.

358
00:45:49,280 --> 00:45:53,870
Currently, I'm losing. So you do have to get the case properly.

359
00:45:55,050 --> 00:46:00,540
S dark non-linear part is equal to minus five times the.

360
00:46:01,710 --> 00:46:05,310
Derivative of u squared and I. Yeah. Okay.

361
00:46:08,260 --> 00:46:12,579
So in fact, there's now a thing called a spin up that we have now set up.

362
00:46:12,580 --> 00:46:18,219
And to solve the problem, we, of course, needed the initial condition, which of course we give in the form of a telephone.

363
00:46:18,220 --> 00:46:23,890
So I could say to have fun and then I could say cosine of x over 16.

364
00:46:26,140 --> 00:46:30,310
Times one plus sign of the x over 16.

365
00:46:32,250 --> 00:46:34,379
And then that's on the domain we specified.

366
00:46:34,380 --> 00:46:40,170
And because it's periodic, I'll tell everyone that it's periodic by saying Trig, though, you don't have to do that.

367
00:46:40,380 --> 00:46:43,950
And then we're going to evolve over some hour interval.

368
00:46:43,950 --> 00:46:49,170
So that doesn't matter. Let's get hard. Line it. I'll say you.

369
00:46:50,410 --> 00:46:58,120
Solution equals spin. So the spin up and then a time interval so we could say zero up to ten.

370
00:47:01,330 --> 00:47:11,739
And then the initial condition you not so assuming I haven't mistyped anything more that will now call a certain algorithm,

371
00:47:11,740 --> 00:47:17,560
which I'm not telling you about just now. And we will get much as what we got before.

372
00:47:17,590 --> 00:47:21,850
Let's see, there's 1020 already on the screen. So let's change this to 20.

373
00:47:23,980 --> 00:47:28,510
So we computed that with the finite differences, and now we're doing it with the spin code.

374
00:47:28,520 --> 00:47:33,430
So please take an image in your brain of that curve and we'll see whether we get the same curve.

375
00:47:37,400 --> 00:47:41,000
Okay. Type space when ready? I'm ready. Is that the same?

376
00:47:41,600 --> 00:47:45,469
I hope so. It won't be the same at time. 500 because of the chaos.

377
00:47:45,470 --> 00:47:51,650
But of course, for the initial times it should be the same. Now, this is actually a much more powerful system.

378
00:47:53,460 --> 00:48:03,720
It's faster and more accurate. So although you can't really see just by looking at the picture, this is a much better solution than we had before.

379
00:48:04,290 --> 00:48:08,220
I think it's stable enough that we could go to large times and nothing would go wrong.

380
00:48:09,300 --> 00:48:12,030
Let's go to a slightly large time and see if that's true.

381
00:48:17,440 --> 00:48:22,690
So next week I'm going to tell you about the algorithms that are used for this a little bit anyway.

382
00:48:22,870 --> 00:48:25,989
And it also works in 2D and 3D.

383
00:48:25,990 --> 00:48:29,410
So you can do various fun things with spin.

384
00:48:29,680 --> 00:48:35,320
If you're interested in reaction to fusion equations in more dimensions, you can say things like spin.

385
00:48:35,830 --> 00:48:43,300
Well, I could say Spin of Chaos, and I'd get the hardcore coated, comatose efficiency demo.

386
00:48:43,360 --> 00:48:47,990
There it is. But then various other things are in there.

387
00:48:48,000 --> 00:48:55,180
I could say KTV. Or I could say spin too.

388
00:48:55,180 --> 00:49:00,790
And do something in multiple dimensions, say. I think that's the Ginsburg Landau equation.

389
00:49:00,790 --> 00:49:04,430
Does that work? Did I make a mistake there?

390
00:49:05,090 --> 00:49:13,520
Guess we have a two dimensional version of this evolving against Burton Lando equation that again has this linear,

391
00:49:13,760 --> 00:49:17,360
non-linear structure and we're going to stop there for today.

392
00:49:17,360 --> 00:49:18,020
So thanks.