1
00:00:00,270 --> 00:00:09,490
Okay. Good morning, everybody. This is a demo from Cleave Mueller's textbook, which I mentioned to you earlier, is freely available online.

2
00:00:09,510 --> 00:00:13,020
It's kind of an undergraduate textbook on numerical computing with MATLAB.

3
00:00:13,710 --> 00:00:17,040
This is the famous example of the double pendulum.

4
00:00:18,900 --> 00:00:25,440
There are obviously you have a fixed point in the middle and then a one pendulum attached to another pendulum.

5
00:00:25,620 --> 00:00:34,080
I think you can start the demo by clicking somewhere and if you click right there, it actually starts exactly vertical.

6
00:00:34,560 --> 00:00:37,920
So that's a fine example of an unstable equilibrium.

7
00:00:38,250 --> 00:00:41,670
And the code is set up so that it's exactly in equilibrium.

8
00:00:41,670 --> 00:00:52,880
Nothing will ever happen if I click near there. But then we're not quite exactly in the equilibrium and eventually it starts moving.

9
00:00:53,510 --> 00:00:56,960
The double pendulum is a great example of a chaotic system.

10
00:00:57,230 --> 00:01:02,240
And you remember our key feature of chaos is that it can't happen in the plane.

11
00:01:02,240 --> 00:01:09,410
It has to be in three D or more. If you have just two dependent variables, you can't have chaos.

12
00:01:09,710 --> 00:01:15,860
But here, depending how you count, there are at least three or four because the middle point is of course fixed.

13
00:01:15,860 --> 00:01:22,040
But then the end point has a position and a velocity, and then the other one also has a position and a velocity.

14
00:01:22,310 --> 00:01:29,719
So there are four parameters, at least in the simplest way you might formulate it.

15
00:01:29,720 --> 00:01:35,480
You could reduce that to three, but you couldn't reduce it to two. So it really is enough for chaos to be a possibility.

16
00:01:35,660 --> 00:01:40,460
Of course not. Every system with three variables is chaotic, but this one beautifully is.

17
00:01:41,810 --> 00:01:49,310
I think it goes a little too slowly. This demo one problem making demos is really difficult from the point of view of speed because you

18
00:01:49,310 --> 00:01:53,420
always make it on your own computer and then someone else's computer runs at a different speed.

19
00:01:53,690 --> 00:01:56,690
I'd like to be able to speed this up, but I haven't figured out how.

20
00:02:00,920 --> 00:02:14,840
Okay. Mueller is an old friend of mine, by the way, and know he invented MATLAB.

21
00:02:16,550 --> 00:02:23,270
He complains that they don't let him touch the code anymore. Okay.

22
00:02:24,470 --> 00:02:31,850
We're talking about ODIs. And today I wanted to say a word about stability regions which are a beautiful

23
00:02:31,850 --> 00:02:36,560
graphical tool to understand the behaviour of numerical methods for ODI use.

24
00:02:37,040 --> 00:02:40,970
So that's section 4.12.

25
00:02:59,110 --> 00:03:03,910
So let me remind you that our universal notation is that we have an O.D.

26
00:03:04,660 --> 00:03:15,130
These pens are long. Sorry. We have an odd initial value problem.

27
00:03:19,250 --> 00:03:25,530
Which we have been writing in the form U. Prine equals F of T and U.

28
00:03:26,660 --> 00:03:31,640
This demo here was an autonomous example where there was no explicit dependence on T,

29
00:03:32,450 --> 00:03:40,670
and we're approximating this by a finite, different formula, a run, a kind of formula or a Adam's formula.

30
00:03:42,020 --> 00:03:47,330
And let me remind you that sometimes they behave unstable.

31
00:03:47,720 --> 00:03:50,959
And a formula I gave you before was brand X.

32
00:03:50,960 --> 00:04:06,830
If you like the unstable formula V and plus one equals minus four V and plus five Z and minus one plus the time step.

33
00:04:08,640 --> 00:04:15,480
Times four F and plus two F and minus one.

34
00:04:16,560 --> 00:04:26,219
I remind you that the sub n refers to our approximate solution at time step and and f seven refers to f evaluated there.

35
00:04:26,220 --> 00:04:35,850
So f seven is shorthand for the function f evaluated that time and and at the approximate solution value.

36
00:04:35,850 --> 00:04:40,770
V7 So that was our example of an unstable formula.

37
00:04:43,280 --> 00:04:46,940
It's completely useless. Nothing you compute with that will ever be any good.

38
00:04:48,620 --> 00:04:55,070
However, if something interesting happens, even with good formulas, and that's sort of the theme of what we're going to talk about today.

39
00:04:55,340 --> 00:04:58,700
Even formulas that have names can have funny behaviour.

40
00:04:59,090 --> 00:05:03,410
So let me remind you of atoms back forth two,

41
00:05:04,310 --> 00:05:20,120
which is a stable formula that was V and plus one equals V and plus time step over two times three F and minus F and minus one.

42
00:05:23,070 --> 00:05:29,880
That's a very good formula. However, it too has some funny behaviour when the time step is too big.

43
00:05:30,210 --> 00:05:33,570
So this one will misbehave for all time, roughly speaking.

44
00:05:33,900 --> 00:05:39,420
This one will misbehave for large time sets and behave fine for small time steps.

45
00:05:39,720 --> 00:05:49,560
So I'm going to demonstrate that in the middle in a moment. We're going to just look at the simple model problem you prime equals minus you.

46
00:05:51,660 --> 00:06:01,110
Couldn't get simpler than that, and we'll do it on the interval from 0 to 10 with the initial condition, you not have one.

47
00:06:03,670 --> 00:06:17,700
So of course the solution is either the mining fee. And let me show you some instability.

48
00:06:25,130 --> 00:06:29,000
So this code is called M 30 and 30 instability.

49
00:06:29,240 --> 00:06:31,490
So if I say I'm 30 instability.

50
00:06:33,950 --> 00:06:42,830
It asks me for a tiny step to solve this problem, and I'm going to start with bigger time steps and then smaller steps.

51
00:06:42,860 --> 00:06:47,990
So, for example, if I say two, then I'm just taking five steps from 0 to 10.

52
00:06:50,120 --> 00:06:53,210
I'm sorry. I guess this code goes to 40, not to ten. Excuse me.

53
00:06:54,170 --> 00:07:03,200
What you see here is the exact solution in red sampled at the appropriate times, and then the computed one with this atoms back forth formula.

54
00:07:03,530 --> 00:07:07,430
So, of course, that's no good at all. Let's do it again.

55
00:07:07,460 --> 00:07:12,730
Let's take time. Step 1.1. So there you see.

56
00:07:12,740 --> 00:07:16,850
Well, something a little plausible for a little while, but it's clearly diverging.

57
00:07:17,240 --> 00:07:23,090
So even stable only formulas can give garbage when the time steps are fairly large.

58
00:07:23,780 --> 00:07:29,989
If you go to one as your time step, you see that's the neutrally stable time step.

59
00:07:29,990 --> 00:07:37,850
It's obviously just on the edge of getting into the right territory with the time step of point nine.

60
00:07:38,480 --> 00:07:41,960
There's a little bit of chatter, but it's sort of tracking the correct solution.

61
00:07:42,200 --> 00:07:45,380
And then as the time step shrinks, you get nice convergence.

62
00:07:47,490 --> 00:07:56,790
So what I've shown you in that example is that even a stable formula like Adams Passport may require the time step to be small for good behaviour.

63
00:07:59,310 --> 00:08:05,250
Now, if you look at the code there, you'll see that I've included both the stable and the unstable formula.

64
00:08:05,820 --> 00:08:10,410
So there's a comment character on the left.

65
00:08:10,500 --> 00:08:18,510
If I get rid of that character, then it will be running Brand X instead of having the best fourth two, which should be unstable for all values of K.

66
00:08:18,690 --> 00:08:23,220
Let's confirm that I've got it preloaded as in 30 you it's called.

67
00:08:24,570 --> 00:08:29,400
So this is now brand X if I run it with time. Step two same kind of thing.

68
00:08:30,060 --> 00:08:35,190
1.11 it's always bad, no matter how small the time step is.

69
00:08:36,480 --> 00:08:40,110
Notice it's telling you the maximum value. Ten to the 61.

70
00:08:40,560 --> 00:08:46,680
We're getting exponential divergence to the 283.

71
00:08:48,530 --> 00:08:52,250
So unstable formulas are bad for all time sets.

72
00:08:52,460 --> 00:08:55,910
Stable formulas are good if the time step is small enough.

73
00:08:56,510 --> 00:09:10,959
And that's the theme of today's lecture. And I want to describe this technique of stability regions that people use in order to

74
00:09:10,960 --> 00:09:15,400
understand that and to decide how small a time step you need for a particular problem.

75
00:09:17,020 --> 00:09:21,309
So now let's improve our model problem by putting in a parameter.

76
00:09:21,310 --> 00:09:27,280
I'll say you prime equals few. So a is a constant.

77
00:09:31,870 --> 00:09:35,410
So this is a model of the different frequencies in the system.

78
00:09:35,420 --> 00:09:39,790
If A is positive, we have exponential growth, if it is negative with exponential decay.

79
00:09:40,000 --> 00:09:47,200
If he's imaginary, we would have isolation. So by varying a, you can explore the different territories of the physics of a problem.

80
00:09:48,430 --> 00:09:58,320
If you plug this in. Two A.D. formula, a discrete Audi E formula.

81
00:09:58,620 --> 00:10:03,690
What you get is a recurrence relation. This is a linear problem.

82
00:10:03,700 --> 00:10:08,730
The old formulas have a linear form, so you get a linear recurrence relation.

83
00:10:10,730 --> 00:10:18,200
And by analysing that recurrence relation, you can learn a lot about whether things will blow up or not.

84
00:10:18,560 --> 00:10:23,450
Loosely speaking, you're going to get blow up when this recurrence relation has some solutions that blow up.

85
00:10:23,660 --> 00:10:26,780
So for example, let's do it.

86
00:10:27,970 --> 00:10:38,800
For Adams Bass fourth to plug in you prime equals AEW and we get the end plus one equals van.

87
00:10:40,570 --> 00:10:55,090
Plus K over two and then our three F, minus F and minus one term becomes three times a VF, minus a V and minus one.

88
00:10:57,330 --> 00:11:02,459
So you see this is a three term recurrence relation involving some fixed coefficients,

89
00:11:02,460 --> 00:11:08,700
three K over two and one, K over two and A relating V and plus one.

90
00:11:08,820 --> 00:11:13,230
The N and the N minus one. So this is a three term recurrence.

91
00:11:18,450 --> 00:11:21,660
And everything is known about the behaviour of recurrence relations.

92
00:11:22,110 --> 00:11:28,140
Whether they have solutions that blow up or not is a matter of polynomial algebra.

93
00:11:30,620 --> 00:11:34,730
So in order to analyse that, well, let me rewrite it first.

94
00:11:35,330 --> 00:11:45,620
Rewrite it as van plus one equals one plus three k over to the N.

95
00:11:48,920 --> 00:11:52,999
Minus K over to the end.

96
00:11:53,000 --> 00:12:00,459
Minus one. And by the way, you always when you do this kind of analysis, get this product.

97
00:12:00,460 --> 00:12:04,330
K l that's the parameter that appears. The timestamp.

98
00:12:06,190 --> 00:12:11,290
Times the constant in your model. That's the scaling you always see.

99
00:12:13,390 --> 00:12:20,390
So in order to find out the behaviour of a recurrence relation, what you do is you look for geometrically growing or decaying solutions.

100
00:12:20,440 --> 00:12:26,800
So you make the answer that your solution is some number alpha or R to the N.

101
00:12:27,580 --> 00:12:34,720
I grew up saying Anzacs, but my friends tell me that's politically incorrect because it assumes that people speak German.

102
00:12:35,230 --> 00:12:37,570
So you're supposed to say Trial solution.

103
00:12:40,380 --> 00:12:47,610
Actually, when I took physics as a freshman, we had a wonderful Austrian lecturer and he said with a good thick accent, and that was wonderful.

104
00:12:48,330 --> 00:12:53,790
Anyway, we tried to try a solution. The energy calls are to the end.

105
00:12:55,200 --> 00:13:00,779
So to analyse the recurrence relation, you say you look for solutions that growth decay geometrically.

106
00:13:00,780 --> 00:13:07,860
Some number race to the end. Now if you plug that in there, you very quickly get a conclusion.

107
00:13:08,850 --> 00:13:10,890
Well, let's see there. You're going to get blow up.

108
00:13:12,940 --> 00:13:27,010
If R is bigger than one in absolute value, you'll also get more gentle blow up if R is equal to one.

109
00:13:28,180 --> 00:13:38,230
And it's a multiple solution in some sense. Which I'll explain in a second.

110
00:13:40,210 --> 00:13:44,530
But our less than one is the K are bigger than one is blow up and then a fire is equal to one.

111
00:13:44,530 --> 00:13:51,040
You need a more careful analysis. So in order to analyse the recurrence relation,

112
00:13:51,040 --> 00:13:56,680
you plug in this unzips to the equation and you get a polynomial and you look at the roots of a polynomial.

113
00:13:56,920 --> 00:14:00,700
So the values of r that matter.

114
00:14:06,890 --> 00:14:14,020
Are the roots of what we would call the characteristic polynomial of the recurrence relation.

115
00:14:20,480 --> 00:14:22,910
Whose aisle gets defined by showing you the example.

116
00:14:23,150 --> 00:14:45,200
For this case the polynomial would be r squared minus one plus three r over 233k over two are plus k over two.

117
00:14:47,540 --> 00:14:55,730
If there's a solution a root of that polynomial that's bigger than one, then there are exploding solutions of your recurrence relations.

118
00:14:56,000 --> 00:15:00,380
If relation. If all of the roots are smaller than one, then all solutions decay.

119
00:15:00,650 --> 00:15:05,420
And if there are some equal to one, then you might need a finer analysis.

120
00:15:05,420 --> 00:15:12,170
In particular, our equals one is okay if it's a simple root.

121
00:15:15,650 --> 00:15:23,030
Of the characteristic polynomial. But it's not okay if it's a multiple.

122
00:15:31,250 --> 00:15:37,610
So let's do this for, let's say brand X for brand x.

123
00:15:42,840 --> 00:15:49,530
We've kind of just wiped off the screen. There it is.

124
00:15:55,580 --> 00:16:01,719
What you find is there is a root. With outside the unit.

125
00:16:01,720 --> 00:16:06,870
His even as a ghost is zero.

126
00:16:06,880 --> 00:16:10,120
So let's look at brand x with k equals zero.

127
00:16:12,280 --> 00:16:25,450
Then it reduces to the recurrence relation r squared plus four, r minus five.

128
00:16:26,410 --> 00:16:32,440
And we can factor that that's equal to R plus five times R minus one.

129
00:16:34,610 --> 00:16:42,410
So what we see is that the brand X recurrence relation has a root that's far outside the unit circle.

130
00:16:47,400 --> 00:16:51,660
And that is the mode of oscillation that you see in this picture.

131
00:16:54,440 --> 00:16:59,450
At every step, roughly speaking, you get amplification by a factor of minus five.

132
00:16:59,720 --> 00:17:03,560
So let's go back and. Look at that.

133
00:17:04,670 --> 00:17:11,510
If, as the timestamp is sufficiently small, you're essentially getting amplification by a factor of minus five at each step here,

134
00:17:11,510 --> 00:17:19,010
we've taken 20 steps, so that's something like five to the 20.

135
00:17:19,010 --> 00:17:27,950
If I'm not doing a careful analysis here, that's where the blow up comes from, from these unstable solutions of the recurrence relation.

136
00:17:29,270 --> 00:17:33,190
But now let's do the more refined analysis for atoms.

137
00:17:33,200 --> 00:17:44,450
That's fourth. So does the recurrence relation have a blowing up solution or not?

138
00:17:44,490 --> 00:17:58,680
Now there's the formula. And what you find is. I'm not going to write down the algebra, but you find that there is a root outside the unit circle.

139
00:18:01,600 --> 00:18:05,140
Whenever K time Z is less than minus one.

140
00:18:11,440 --> 00:18:16,540
On the other hand, if Kay is big or small enough that that's not less than minus one.

141
00:18:16,780 --> 00:18:29,680
So if you have minus one less than K less than zero, then both routes are well behaved.

142
00:18:38,430 --> 00:18:47,320
So the recurrence relation tells you that you would expect the other the best fourth formula to misbehave when K is less than minus one.

143
00:18:47,340 --> 00:18:53,670
Now the model we used that showed misbehaviour had a equal to minus one, didn't we?

144
00:18:54,030 --> 00:18:58,560
Looked at the problem of you prime equals minus u.

145
00:19:01,350 --> 00:19:06,150
So we looked at you. Prime equals minus you, i.e. a equals minus one.

146
00:19:07,560 --> 00:19:11,190
So you would expect out of the best force to have trouble when.

147
00:19:13,250 --> 00:19:17,330
The step size times minus one is less than minus one.

148
00:19:17,360 --> 00:19:24,830
In other words, when the step size is bigger than minus one, so you have unstable, bigger than one, you have unstable behaviour.

149
00:19:29,900 --> 00:19:35,360
Four steps size bigger than one unstable behaviour.

150
00:19:37,310 --> 00:19:44,530
Forceps size less than one. I feel I'm doing this a little bit hastily.

151
00:19:44,530 --> 00:19:51,519
But anyway, the point is that by analysing this model problem of just the linear equation,

152
00:19:51,520 --> 00:19:57,040
you prime equal value, you can very quickly learn whether things are likely to blow up or not.

153
00:20:05,450 --> 00:20:13,880
Now I mentioned a stability region. A stability region is a region in the complex plane of parameter space where you get stable behaviour.

154
00:20:14,450 --> 00:20:19,650
So every formula like Adam's bash forth or Brand X has a stability region.

155
00:20:21,290 --> 00:20:27,380
In the case of Brand X, it's an empty region. So given A.D. formula.

156
00:20:27,390 --> 00:20:35,620
The formula. We define the stability region.

157
00:20:39,460 --> 00:20:42,820
To be a region in the K plan.

158
00:20:48,030 --> 00:21:05,210
So it's the seven points. Or the set of values, in this case, a complex plane such that the recurrence is stable, which means.

159
00:21:07,190 --> 00:21:10,400
The characteristic polynomial.

160
00:21:13,460 --> 00:21:17,120
Has all routes in the unit disc. Or more precisely.

161
00:21:20,790 --> 00:21:23,909
All routes have to be less than or equal to one in modulus.

162
00:21:23,910 --> 00:21:27,360
And if there have modulus one, they have to be simple.

163
00:21:43,870 --> 00:21:52,660
So that gives you a picture which you can draw in the complex plane to find out whether a loaded formula will seem to blow up or not.

164
00:21:55,430 --> 00:21:58,310
In a moment, I'm going to show you the album's best picture.

165
00:21:58,610 --> 00:22:03,500
But in order to compute one, it's simpler to use an even simpler formula, the Euler formula.

166
00:22:04,520 --> 00:22:08,879
So let's do that. As an example.

167
00:22:08,880 --> 00:22:16,310
Then the first order Euler formula, which is items that fourth one, the first order items, that's fourth form.

168
00:22:17,640 --> 00:22:24,150
So you remember that is the end plus one equals then plus k times s ten.

169
00:22:26,310 --> 00:22:31,710
If we apply this to the model problem. Which is you.

170
00:22:31,880 --> 00:22:37,430
Prime equals a you. Then we get the end.

171
00:22:37,430 --> 00:22:44,330
Plus one equals the end. Plus k a the end.

172
00:22:47,810 --> 00:22:54,800
In other words, vrn plus one equals one plus k times van.

173
00:22:59,310 --> 00:23:13,520
So the characteristic polynomial. Is P of r equals R minus one plus k.

174
00:23:18,160 --> 00:23:24,400
The root of the characteristic polynomial. There's only one root because it's a linear polynomial because it's a one step formula.

175
00:23:24,730 --> 00:23:35,550
The root is one plus k a. When is that smaller than one in absolute value?

176
00:23:35,730 --> 00:23:42,060
Well, one plus K is smaller than one of K is within one of the point minus one.

177
00:23:42,420 --> 00:23:48,320
So there's our stability region. We're in the car.

178
00:23:48,390 --> 00:24:07,720
Plane. The circle of radius one around the point minus one in the K plane is the stability region for the Euler formula.

179
00:24:16,190 --> 00:24:25,069
So to remind you what that means, it means if you solve a nobody with this formula and your time constants are such that K sort of lives in there,

180
00:24:25,070 --> 00:24:33,469
then you'll probably get good behaviour if your time constants are bigger so that K lives outside there, you'll probably get explosively bad.

181
00:24:33,470 --> 00:24:39,830
B that's gets a linear theme, more complicated schemes have more beautiful pictures.

182
00:24:39,830 --> 00:24:43,520
And let me show them to you just on the computer you can see here.

183
00:24:45,250 --> 00:24:48,730
We're summarising 50 years of analysis of the schemes.

184
00:24:50,550 --> 00:24:56,700
So I have a code here called and 31 regex which plots.

185
00:24:58,580 --> 00:25:16,860
Some stability regions. I hope. So there's a lot of implication in these diagrams.

186
00:25:16,860 --> 00:25:21,650
So let me talk about the three families illustrated there. We haven't discussed this one yet.

187
00:25:22,070 --> 00:25:23,690
We have discussed these two.

188
00:25:24,650 --> 00:25:33,420
So what you see here are pictures of stability regions for the first three atoms, vast fourth formulas and the first four rugby.

189
00:25:33,420 --> 00:25:40,570
A kind of formula in the case of out of the first fourth, we just defined we get to derive the blue one, the first order out of the past.

190
00:25:40,640 --> 00:25:49,060
Fourth formula has that stability region. We didn't arrive the red one, but we computed with it that atoms back forth two.

191
00:25:49,610 --> 00:25:53,780
So you see that as the order of another is that the fourth formula increases.

192
00:25:54,470 --> 00:26:00,050
The formula gets a little bit more finicky because the stability region is shrinking down a little bit,

193
00:26:00,410 --> 00:26:03,500
which means you probably need somewhat smaller timestamps.

194
00:26:03,830 --> 00:26:07,580
So for the red one, you probably need half the size of time step.

195
00:26:07,580 --> 00:26:14,480
As for the blue, the green one half again, so slightly smaller timestamps in order to get good behaviour.

196
00:26:16,510 --> 00:26:24,350
They're running a cut of formulas. Very loosely speaking, are similar, though the details are interestingly different.

197
00:26:24,560 --> 00:26:27,920
The first one they had a formula is again the forward Euler method.

198
00:26:27,920 --> 00:26:35,540
So that blue circle is exactly the same as that. But then a second order among a kind of formula has that red behaviour and third order,

199
00:26:35,540 --> 00:26:39,200
and then the purple is the famous fourth order run the kind of formula.

200
00:26:39,590 --> 00:26:44,930
So that purple curve is probably the most important curve in this picture.

201
00:26:45,740 --> 00:26:49,480
That is the stability region for the fourth order run kind of form.

202
00:26:50,180 --> 00:26:54,590
And one thing you note is that it has some nice behaviour on the imaginary axis.

203
00:26:55,070 --> 00:27:01,970
It will do better with oscillatory problems whose physics is on the imaginary axis than the other.

204
00:27:02,690 --> 00:27:05,000
Then at least the first order out of the batch. Fourth form.

205
00:27:07,170 --> 00:27:12,300
As you increase the order of the rung, a cut of formulas notice the stability regions are getting a little bit bigger,

206
00:27:12,630 --> 00:27:17,850
so they're slightly less finicky about step sizes than atoms vegetable.

207
00:27:20,450 --> 00:27:26,030
Now these we're going to talk about in a minute. This picture is completely different.

208
00:27:28,690 --> 00:27:34,030
It's a formula called backward differentiation. And I'm showing you the first six stability regions.

209
00:27:34,240 --> 00:27:41,050
But what makes it completely different is that here the curves are the outer boundaries of the stability regions.

210
00:27:41,650 --> 00:27:48,460
Here they're actually the inner ground. So the backward differentiation formulas have unbounded stability regions.

211
00:27:48,910 --> 00:27:51,970
The blue one, which is the backward Euler formula.

212
00:27:54,130 --> 00:28:11,110
Let me draw that. This is the picture of the stability region of the forward Euler formula, which is the blue curve on the left or in the middle.

213
00:28:13,990 --> 00:28:25,930
But if you look at the first order of backward differentiation formula, which might be called B, D, this is the same as backward Euler.

214
00:28:30,060 --> 00:28:42,690
Which is the end plus the end plus one equals V and plus k, f and plus one rather than n.

215
00:28:44,640 --> 00:28:48,830
You find that the stability region. Is.

216
00:28:49,720 --> 00:28:53,430
The exterior. Of the very same circle.

217
00:28:55,910 --> 00:29:00,140
Sorry. I'm saying that wrong. It's not the very same circle. It's the circle on the right of the origin.

218
00:29:00,660 --> 00:29:27,209
Apologies. So whereas these formula that is based forth and rugged cut need very small time steps to work this backward.

219
00:29:27,210 --> 00:29:30,660
Euler method doesn't mean such a small time for that to work.

220
00:29:30,660 --> 00:29:36,630
It's stable all the way out to infinity, and it will only blow up for tiny steps in that region.

221
00:29:40,560 --> 00:29:47,010
I mentioned there's a lot of history, mathematical research that goes into these pictures, especially the ones on the right.

222
00:29:49,530 --> 00:29:54,929
These are the boundaries of the curves for degree one, two, three, four, five and six.

223
00:29:54,930 --> 00:29:59,220
Backward differentiation. It turns out the seventh one really looks different.

224
00:29:59,250 --> 00:30:03,570
You can see the cyan colour. The light blue is bending towards the negative axis.

225
00:30:03,780 --> 00:30:07,970
The seventh one crosses the negative axis, and that has big effects.

226
00:30:07,980 --> 00:30:12,960
That means that people don't really use those formulas for orders higher than six.

227
00:30:22,690 --> 00:30:31,360
Well, you're probably fairly confused at this point. So let me explain what we find from all this analysis.

228
00:30:31,960 --> 00:30:39,790
We find that for some problems, these nice, simple methods are great, but that other problems need these fancier methods.

229
00:30:39,790 --> 00:30:45,730
And those are the problems called stiff odds. So I want to talk now about the idea of a stiff old.

230
00:31:06,070 --> 00:32:01,209
Talking. So if you go back to the 1950s, this man, Dahlquist, that I've mentioned,

231
00:32:01,210 --> 00:32:09,070
established the basic theory of when these methods converge, when does a linear multistep method converge?

232
00:32:09,070 --> 00:32:15,730
And it's pretty close to when does a wrong a kind of method converge? And the basic result of his stability theory.

233
00:32:19,800 --> 00:32:26,210
Was. That all that matters is what happens as K goes to zero.

234
00:32:27,380 --> 00:32:30,440
Remember I mentioned the theorem, that convergence.

235
00:32:34,290 --> 00:32:44,819
Equals stability plus consistency. And any formula you derive intelligently will be consistent.

236
00:32:44,820 --> 00:32:47,460
That just means it has order of accuracy, at least one.

237
00:32:47,880 --> 00:33:00,930
So Dahlquist theory from 1956 said that basically to find out if a formula works, you have to find out if it's stable and stability was defined.

238
00:33:03,050 --> 00:33:06,350
As all fruits.

239
00:33:08,490 --> 00:33:16,860
Of the recurrence relation. Yeah our stable in the sense we've been talking about in the limit of zero timestamp.

240
00:33:22,970 --> 00:33:27,290
So his great theory says that you don't have to look at these pictures.

241
00:33:27,290 --> 00:33:28,910
All you have to do is look at the origin.

242
00:33:29,120 --> 00:33:35,210
And so long as the origin is in the stability region, on the boundary of the stability region, your method will work.

243
00:33:37,920 --> 00:33:40,200
Of course it's the theorem. So of course it's true.

244
00:33:41,820 --> 00:33:50,010
If you have any stable, finite difference approximation to an O.D. and you use it in exact arithmetic without rounding errors,

245
00:33:50,370 --> 00:33:54,180
then as the timestamp goes to zero, you will converge to the right answer.

246
00:33:56,010 --> 00:34:03,780
But what happened is that very soon thereafter, people discovered that in many problems in practice they weren't getting convergence.

247
00:34:04,200 --> 00:34:09,960
Or if they were getting convergence, the time step had to be really ridiculously small to get good answers.

248
00:34:10,380 --> 00:34:16,320
And that's what the problem of stiffness was. So stiff problems.

249
00:34:20,460 --> 00:34:23,880
Can be defined in many ways. One way to say it is there.

250
00:34:24,780 --> 00:34:28,170
They tend to be problems with widely varying timescales.

251
00:34:37,580 --> 00:34:54,570
And the effect of that. It turns out to be that you often need k very, very small to get reasonable solutions.

252
00:35:04,130 --> 00:35:10,300
So it's a beautiful case of correct mathematics turning out to be overlooking part of what really mattered.

253
00:35:10,310 --> 00:35:16,640
There was something going on in the problem that really mattered. And let me give an illustration of a step by step problem.

254
00:35:17,000 --> 00:35:21,380
I say it's a problem with a widely varying time scale. Here's a simple example.

255
00:35:26,490 --> 00:35:34,920
So what I'm going to do is I'm going to have an example with a smooth solution that nevertheless has some physics in it on a faster timescale.

256
00:35:35,130 --> 00:35:44,820
And here's what I'll do. I'll say you prime equals -100 times, you minus cosine of t.

257
00:35:47,610 --> 00:35:51,930
Minus sign up and you have zero.

258
00:35:53,630 --> 00:35:58,620
Equals one. So that's a first order scaler.

259
00:35:59,370 --> 00:36:06,079
ODI Of course, it's mathematically very simple and the solution can be written down.

260
00:36:06,080 --> 00:36:09,410
It's cosine of tea. So solution.

261
00:36:11,880 --> 00:36:15,660
U of t equals cosine of t.

262
00:36:16,260 --> 00:36:26,310
That's obvious. You plug it in. This thing cancels for the whole term with 100 and it isn't there, and the derivative of the cosine is minus the sign.

263
00:36:26,580 --> 00:36:32,440
So there you are. So. The solution.

264
00:36:34,180 --> 00:36:42,430
Is this cosine function? However, the equation has a much faster timescale in at this -100.

265
00:36:42,760 --> 00:36:45,850
And that means that if you perturb the solution a little bit.

266
00:36:48,220 --> 00:36:54,190
For example, if I started here instead of there, then there would be a very rapid transit.

267
00:36:55,390 --> 00:36:58,180
Or if I started here, there would be a very rapid transit.

268
00:36:58,450 --> 00:37:06,670
There's a physics in the problem that's 100 times faster than the timescale of the solution you're interested in.

269
00:37:08,450 --> 00:37:10,820
So that's a classic example of stiffness.

270
00:37:12,230 --> 00:37:18,470
If you're trying to compute this nice, smooth solution, you'd like to be able to take time steps on that timescale.

271
00:37:19,100 --> 00:37:29,090
But because of this, if you use the wrong finite difference formulas, you may discover you have to take the time steps that are really microscopic.

272
00:37:29,420 --> 00:37:32,930
We'd like to be able to take time steps that look something like this.

273
00:37:33,740 --> 00:37:36,770
We don't want them to have to be 100 times smaller.

274
00:37:39,700 --> 00:37:55,700
Let me illustrate that. So the illustration is the code called M 32.

275
00:37:58,280 --> 00:38:09,770
And 32 stiffness. And what it does is solve this problem and I'm going to shrink the step size.

276
00:38:10,100 --> 00:38:15,710
So that's the error when you solve it with a big step size of 0.25.

277
00:38:16,880 --> 00:38:20,990
Now, of course, as you shrink the step size, you'd like the error to converge to zero.

278
00:38:24,080 --> 00:38:29,270
So we'll cut the steps in half. And we're using here the atoms passport formula.

279
00:38:32,490 --> 00:38:40,770
I cut the step size in half and the error gets worse. I cut it in half again and it gets worse.

280
00:38:41,310 --> 00:38:43,590
So this is at its best fourth.

281
00:38:43,590 --> 00:38:51,810
It's a convergent formula and yet you can see we're getting horrible behaviour and loosely speaking because of the one in there,

282
00:38:52,620 --> 00:38:57,180
it's going to keep getting worse until we get down to a step size on the scale of 100.

283
00:38:59,280 --> 00:39:03,720
Cut it in half again works. Cut it in half again.

284
00:39:04,560 --> 00:39:10,200
Suddenly it's beginning to look plausible. It's still ten to the eighth, but it's heading the right way.

285
00:39:11,310 --> 00:39:14,370
And now suddenly, look, the step size is now less than 100.

286
00:39:15,120 --> 00:39:21,370
And suddenly I have good behaviour. So.

287
00:39:23,830 --> 00:39:29,290
That's what often happens if you use firearms based on a problem like this.

288
00:39:31,330 --> 00:39:36,400
With multiple time scales, you need a ridiculously small time step to make it work.

289
00:39:37,060 --> 00:39:42,010
Here it's only 100, but you can easily imagine problems in the sciences with much bigger ratios.

290
00:39:42,010 --> 00:39:45,520
And an example that people talk about a lot is in molecular dynamics.

291
00:39:45,760 --> 00:39:52,300
Ideally, you have things happening on a femtosecond scale, ten to the minus 15 seconds that are in the physics,

292
00:39:52,750 --> 00:39:58,450
and yet may be what your tracking is on a nanosecond timescale or a microsecond timescale.

293
00:39:59,140 --> 00:40:02,710
You don't want to have to take one step outside ten to the -15.

294
00:40:03,980 --> 00:40:16,300
So you need some solution to this problem. Another example is in atmospheric calculations when people are computing the state of the air around us.

295
00:40:16,870 --> 00:40:26,890
The physics, for the most part, moves at a relatively slow convective scale, and yet the equations allow for sound waves.

296
00:40:27,370 --> 00:40:31,120
So that's a very different timescale. Of course, sound waves move very fast.

297
00:40:31,600 --> 00:40:39,130
You don't want to have to take time steps that are on the scale, you know, small enough to resolve the soundwaves.

298
00:40:39,370 --> 00:40:42,730
And yet, if you use the wrong numerical schemes, you might have to do that.

299
00:40:43,210 --> 00:40:46,640
So these questions of stiffness are all over the place with ease and.

300
00:40:46,640 --> 00:40:53,350
PD The basic solution I'm going to mention are these backwards differentiation formulas.

301
00:40:53,650 --> 00:40:57,220
So let me show you that if you look at the code.

302
00:41:00,960 --> 00:41:07,910
You see that the first block of code this is m 32 is using the second order Adams batch, fourth formula.

303
00:41:10,370 --> 00:41:17,260
The next block of code is using the second order backward differentiation form for the big time step.

304
00:41:17,270 --> 00:41:23,860
It's doing better. But then here's the amazing thing. Every time I cut the time, step in half, it just converts.

305
00:41:28,850 --> 00:41:33,920
So you can see some formulas have totally different behaviour than others.

306
00:41:34,310 --> 00:41:39,320
And the reason for that has to do with the stability routine that I plotted earlier.

307
00:41:39,980 --> 00:41:46,100
So let's do that again. In fact, if I say FIG. and 31 regions.

308
00:41:54,340 --> 00:42:00,850
So this is the second order backwards differentiation formula, which is the red curve in the right there.

309
00:42:01,120 --> 00:42:04,210
The stability region is the exterior of that red curve.

310
00:42:04,240 --> 00:42:11,500
So it's virtually the whole complex plain and in particular along the crucial negative, imaginary negative real axis.

311
00:42:11,920 --> 00:42:15,490
It's it contains the whole negative, real axis.

312
00:42:15,820 --> 00:42:22,980
So the second order of backward differentiation formula is, in some sense, it's going to behave strangely for all timestamps.

313
00:42:24,800 --> 00:42:28,190
Whereas the second order at best fourth is different.

314
00:42:31,730 --> 00:42:36,350
So anyway, this got discovered in the 1960s, largely.

315
00:42:42,670 --> 00:42:52,600
So among the mathematicians, Dahlquist wrote a his second grade paper in 1963, in which he talked all about this stuff.

316
00:42:52,810 --> 00:43:00,550
But it turns out that actually a crucial paper had been written a decade earlier by some eminent chemists, Curtis and Hirschfeld her.

317
00:43:02,600 --> 00:43:14,650
Hershey Felder in 1952. And they wrote a paper which is really very special, in which they had actually pointed out this phenomenon.

318
00:43:14,690 --> 00:43:22,850
There's no general convergence theory like Alquist, but they saw the point that things go wrong for certain formulas and not for others.

319
00:43:23,300 --> 00:43:30,140
They tend to think it's always impressive when a paper in the sciences does more than one great thing.

320
00:43:30,440 --> 00:43:33,680
So they identified the problem of stiffness.

321
00:43:34,730 --> 00:43:43,400
They use that word, and then they introduce backward differentiation formulas to solve.

322
00:43:45,750 --> 00:43:53,040
The word seems to have been coined by John Tukey, the great statistician who also coined the word [INAUDIBLE] and the word.

323
00:43:53,220 --> 00:43:57,570
Well, maybe he coined the word software. At least he was almost the first person to use the word software.

324
00:43:59,150 --> 00:44:10,010
So the idea of backward differentiation formulas is to they're highly implicit and I haven't really used these words yet, but explicit formulas.

325
00:44:14,900 --> 00:44:23,270
Do not include F and plus one and implicit formulas do include.

326
00:44:26,770 --> 00:44:29,800
F and plus one. Now, why is that so important?

327
00:44:31,420 --> 00:44:34,540
Well, if if this is in the formula, then you've got more work to do.

328
00:44:34,810 --> 00:44:45,360
So, for example, the forward Euler formula is the end, plus one equals the M plus k F1.

329
00:44:45,880 --> 00:44:54,340
The backward Euler formula is V and plus one equals v, n plus k and minus one plus one.

330
00:44:56,880 --> 00:44:59,870
To implement that. To find the new value of V,

331
00:45:00,240 --> 00:45:10,830
you have to solve a system of equations which in principle in general will be nonlinear because F and plus one is an abbreviation for a function.

332
00:45:23,700 --> 00:45:29,880
So to implement an implicit formula, you have to solve a system of equations that every time step.

333
00:45:32,400 --> 00:45:41,420
And in general, that will be a nonlinear system. So that's the downside.

334
00:45:41,900 --> 00:45:52,280
More work is involved to implement them, but the upside is that they have potentially spectacularly better behaviour.

335
00:45:52,880 --> 00:45:56,420
So for a fifth problem, it's worth it. For a non-native problem.

336
00:45:56,420 --> 00:46:05,809
It isn't worth it. And if you look at the code we just ran M32, if you look in detail at what's going on there,

337
00:46:05,810 --> 00:46:12,200
you'll see that a little algebra has solved a scalar equation to get.

338
00:46:13,240 --> 00:46:17,740
To solve this. But this was a linear problem, so that's easy for linear problems.

339
00:46:18,430 --> 00:46:23,500
Of course, solving system of equations is no big deal. For nonlinear problems, you need some kind of iteration.

340
00:46:25,000 --> 00:46:29,110
So that's the. Conclusion.

341
00:46:32,750 --> 00:46:37,640
If you have a non sniff problem, probably explicit methods are best.

342
00:46:41,920 --> 00:46:50,380
Because they're so easy to implement and work well. If you have a stiff problem, you probably want to use implicit methods.

343
00:46:52,830 --> 00:47:12,630
E.g. backwards differentiation. When it comes for ODIs, this is a pretty good summary of the whole picture.

344
00:47:12,840 --> 00:47:19,100
For PETA is the same issues are all there but everything is of course more complicated and they're in the land of pity.

345
00:47:19,180 --> 00:47:22,589
You have all sorts of technologies that are sort of in between these two.

346
00:47:22,590 --> 00:47:31,560
They're semi implicit, trying to find some compromise that enables you to take reasonable sized time steps without doing too much work

347
00:47:31,860 --> 00:47:39,810
with pieces that the matrices could be so big that solving the fully nonlinear system would be impractical.

348
00:47:44,270 --> 00:47:47,780
Let me show the next demo to illustrate this.

349
00:47:48,980 --> 00:48:04,640
So Code three is illustrating the importance of this solvers by solving the Vanderpoel equation with two different sorts of methods.

350
00:48:10,100 --> 00:48:22,700
So if I run in 33, it's the Vanderpoel equation which we know and love and I get to adjust the parameter that controls the non-linear down.

351
00:48:23,210 --> 00:48:28,430
If I take ten then the solution looks like.

352
00:48:30,230 --> 00:48:33,410
That. So you've seen things like this before.

353
00:48:35,000 --> 00:48:41,420
Now, is that a stiff problem? Well, visually, it looks stiff because there's a very rapid transition compared with the other stuff.

354
00:48:41,720 --> 00:48:45,410
So that is indeed a somewhat stiff problem. Different time scales.

355
00:48:46,940 --> 00:48:52,370
And yet, if I solve it with MATLAB, stiff solver te23s,

356
00:48:52,790 --> 00:48:58,160
which is a backward differentiation code, it takes a lot longer than using the non-state solver.

357
00:48:58,880 --> 00:49:02,090
So this is not stiff enough to make much difference.

358
00:49:04,790 --> 00:49:12,580
Suppose I do it again? With a constant of 100.

359
00:49:15,050 --> 00:49:21,110
Now this is stiffer. We have the very rapid transition and then the slow timescale.

360
00:49:21,110 --> 00:49:22,270
It's not really safe there.

361
00:49:22,280 --> 00:49:29,150
You need small time steps there, but throughout this regime it is very safe because the equation has this big number 100 in it.

362
00:49:29,450 --> 00:49:32,120
But the solution you care about is moving slowly.

363
00:49:32,600 --> 00:49:38,239
Now you find a stiff solver is about three times quicker than the monster solver, which is not very exciting.

364
00:49:38,240 --> 00:49:41,810
But there it is. However, let's go further.

365
00:49:43,190 --> 00:49:57,880
Suppose the constant is a thousandth. So now there's a very rapid transit that's not stiff because in that regime you really need small time steps.

366
00:49:58,210 --> 00:50:05,590
But then most of the time it's very stiff with a number 1000 in the equation that isn't reflected in the solution.

367
00:50:05,800 --> 00:50:11,920
And now you see this big difference. It's 40 times faster using the backward differentiation on.

368
00:50:14,050 --> 00:50:20,710
So if this matters hugely when you have big ratios, thousands or millions in your timescales.

369
00:50:21,280 --> 00:50:26,160
The final thing I wanted to do was show you some more of the fun in action.

370
00:50:26,170 --> 00:50:30,190
I think I won't do that because we're at time, but you'll see a little code at the bottom there.

371
00:50:30,400 --> 00:50:38,830
Just to illustrate to have fun solving a boundary value problem for the homework assignment do in 11 days,

372
00:50:39,610 --> 00:50:42,700
you're welcome to use MATLAB or have fun, whatever you like.

373
00:50:43,090 --> 00:50:48,819
Personally, I find cheap fun much easier, but it's not appropriate for me to force you to do that.

374
00:50:48,820 --> 00:50:50,230
By all means, use MATLAB.