Okay. Good morning, everybody. This is a demo from Cleave Mueller's textbook, which I mentioned to you earlier, is freely available online.
It's kind of an undergraduate textbook on numerical computing with MATLAB.
This is the famous example of the double pendulum.
There are obviously you have a fixed point in the middle and then a one pendulum attached to another pendulum.
I think you can start the demo by clicking somewhere and if you click right there, it actually starts exactly vertical.
So that's a fine example of an unstable equilibrium.
And the code is set up so that it's exactly in equilibrium.
Nothing will ever happen if I click near there. But then we're not quite exactly in the equilibrium and eventually it starts moving.
The double pendulum is a great example of a chaotic system.
And you remember our key feature of chaos is that it can't happen in the plane.
It has to be in three D or more. If you have just two dependent variables, you can't have chaos.
But here, depending how you count, there are at least three or four because the middle point is of course fixed.
But then the end point has a position and a velocity, and then the other one also has a position and a velocity.
So there are four parameters, at least in the simplest way you might formulate it.
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.
Of course not. Every system with three variables is chaotic, but this one beautifully is.
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
always make it on your own computer and then someone else's computer runs at a different speed.
I'd like to be able to speed this up, but I haven't figured out how.
Okay. Mueller is an old friend of mine, by the way, and know he invented MATLAB.
He complains that they don't let him touch the code anymore. Okay.
We're talking about ODIs. And today I wanted to say a word about stability regions which are a beautiful
graphical tool to understand the behaviour of numerical methods for ODI use.
So that's section 4.12.
So let me remind you that our universal notation is that we have an O.D.
These pens are long. Sorry. We have an odd initial value problem.
Which we have been writing in the form U. Prine equals F of T and U.
This demo here was an autonomous example where there was no explicit dependence on T,
and we're approximating this by a finite, different formula, a run, a kind of formula or a Adam's formula.
And let me remind you that sometimes they behave unstable.
And a formula I gave you before was brand X.
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.
Times four F and plus two F and minus one.
I remind you that the sub n refers to our approximate solution at time step and and f seven refers to f evaluated there.
So f seven is shorthand for the function f evaluated that time and and at the approximate solution value.
V7 So that was our example of an unstable formula.
It's completely useless. Nothing you compute with that will ever be any good.
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.
Even formulas that have names can have funny behaviour.
So let me remind you of atoms back forth two,
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.
That's a very good formula. However, it too has some funny behaviour when the time step is too big.
So this one will misbehave for all time, roughly speaking.
This one will misbehave for large time sets and behave fine for small time steps.
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.
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.
So of course the solution is either the mining fee. And let me show you some instability.
So this code is called M 30 and 30 instability.
So if I say I'm 30 instability.
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.
So, for example, if I say two, then I'm just taking five steps from 0 to 10.
I'm sorry. I guess this code goes to 40, not to ten. Excuse me.
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.
So, of course, that's no good at all. Let's do it again.
Let's take time. Step 1.1. So there you see.
Well, something a little plausible for a little while, but it's clearly diverging.
So even stable only formulas can give garbage when the time steps are fairly large.
If you go to one as your time step, you see that's the neutrally stable time step.
It's obviously just on the edge of getting into the right territory with the time step of point nine.
There's a little bit of chatter, but it's sort of tracking the correct solution.
And then as the time step shrinks, you get nice convergence.
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.
Now, if you look at the code there, you'll see that I've included both the stable and the unstable formula.
So there's a comment character on the left.
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.
Let's confirm that I've got it preloaded as in 30 you it's called.
So this is now brand X if I run it with time. Step two same kind of thing.
1.11 it's always bad, no matter how small the time step is.
Notice it's telling you the maximum value. Ten to the 61.
We're getting exponential divergence to the 283.
So unstable formulas are bad for all time sets.
Stable formulas are good if the time step is small enough.
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
understand that and to decide how small a time step you need for a particular problem.
So now let's improve our model problem by putting in a parameter.
I'll say you prime equals few. So a is a constant.
So this is a model of the different frequencies in the system.
If A is positive, we have exponential growth, if it is negative with exponential decay.
If he's imaginary, we would have isolation. So by varying a, you can explore the different territories of the physics of a problem.
If you plug this in. Two A.D. formula, a discrete Audi E formula.
What you get is a recurrence relation. This is a linear problem.
The old formulas have a linear form, so you get a linear recurrence relation.
And by analysing that recurrence relation, you can learn a lot about whether things will blow up or not.
Loosely speaking, you're going to get blow up when this recurrence relation has some solutions that blow up.
So for example, let's do it.
For Adams Bass fourth to plug in you prime equals AEW and we get the end plus one equals van.
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.
So you see this is a three term recurrence relation involving some fixed coefficients,
three K over two and one, K over two and A relating V and plus one.
The N and the N minus one. So this is a three term recurrence.
And everything is known about the behaviour of recurrence relations.
Whether they have solutions that blow up or not is a matter of polynomial algebra.
So in order to analyse that, well, let me rewrite it first.
Rewrite it as van plus one equals one plus three k over to the N.
Minus K over to the end.
Minus one. And by the way, you always when you do this kind of analysis, get this product.
K l that's the parameter that appears. The timestamp.
Times the constant in your model. That's the scaling you always see.
So in order to find out the behaviour of a recurrence relation, what you do is you look for geometrically growing or decaying solutions.
So you make the answer that your solution is some number alpha or R to the N.
I grew up saying Anzacs, but my friends tell me that's politically incorrect because it assumes that people speak German.
So you're supposed to say Trial solution.
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.
Anyway, we tried to try a solution. The energy calls are to the end.
So to analyse the recurrence relation, you say you look for solutions that growth decay geometrically.
Some number race to the end. Now if you plug that in there, you very quickly get a conclusion.
Well, let's see there. You're going to get blow up.
If R is bigger than one in absolute value, you'll also get more gentle blow up if R is equal to one.
And it's a multiple solution in some sense. Which I'll explain in a second.
But our less than one is the K are bigger than one is blow up and then a fire is equal to one.
You need a more careful analysis. So in order to analyse the recurrence relation,
you plug in this unzips to the equation and you get a polynomial and you look at the roots of a polynomial.
So the values of r that matter.
Are the roots of what we would call the characteristic polynomial of the recurrence relation.
Whose aisle gets defined by showing you the example.
For this case the polynomial would be r squared minus one plus three r over 233k over two are plus k over two.
If there's a solution a root of that polynomial that's bigger than one, then there are exploding solutions of your recurrence relations.
If relation. If all of the roots are smaller than one, then all solutions decay.
And if there are some equal to one, then you might need a finer analysis.
In particular, our equals one is okay if it's a simple root.
Of the characteristic polynomial. But it's not okay if it's a multiple.
So let's do this for, let's say brand X for brand x.
We've kind of just wiped off the screen. There it is.
What you find is there is a root. With outside the unit.
His even as a ghost is zero.
So let's look at brand x with k equals zero.
Then it reduces to the recurrence relation r squared plus four, r minus five.
And we can factor that that's equal to R plus five times R minus one.
So what we see is that the brand X recurrence relation has a root that's far outside the unit circle.
And that is the mode of oscillation that you see in this picture.
At every step, roughly speaking, you get amplification by a factor of minus five.
So let's go back and. Look at that.
If, as the timestamp is sufficiently small, you're essentially getting amplification by a factor of minus five at each step here,
we've taken 20 steps, so that's something like five to the 20.
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.
But now let's do the more refined analysis for atoms.
That's fourth. So does the recurrence relation have a blowing up solution or not?
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.
Whenever K time Z is less than minus one.
On the other hand, if Kay is big or small enough that that's not less than minus one.
So if you have minus one less than K less than zero, then both routes are well behaved.
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.
Now the model we used that showed misbehaviour had a equal to minus one, didn't we?
Looked at the problem of you prime equals minus u.
So we looked at you. Prime equals minus you, i.e. a equals minus one.
So you would expect out of the best force to have trouble when.
The step size times minus one is less than minus one.
In other words, when the step size is bigger than minus one, so you have unstable, bigger than one, you have unstable behaviour.
Four steps size bigger than one unstable behaviour.
Forceps size less than one. I feel I'm doing this a little bit hastily.
But anyway, the point is that by analysing this model problem of just the linear equation,
you prime equal value, you can very quickly learn whether things are likely to blow up or not.
Now I mentioned a stability region. A stability region is a region in the complex plane of parameter space where you get stable behaviour.
So every formula like Adam's bash forth or Brand X has a stability region.
In the case of Brand X, it's an empty region. So given A.D. formula.
The formula. We define the stability region.
To be a region in the K plan.
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.
The characteristic polynomial.
Has all routes in the unit disc. Or more precisely.
All routes have to be less than or equal to one in modulus.
And if there have modulus one, they have to be simple.
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.
In a moment, I'm going to show you the album's best picture.
But in order to compute one, it's simpler to use an even simpler formula, the Euler formula.
So let's do that. As an example.
Then the first order Euler formula, which is items that fourth one, the first order items, that's fourth form.
So you remember that is the end plus one equals then plus k times s ten.
If we apply this to the model problem. Which is you.
Prime equals a you. Then we get the end.
Plus one equals the end. Plus k a the end.
In other words, vrn plus one equals one plus k times van.
So the characteristic polynomial. Is P of r equals R minus one plus k.
The root of the characteristic polynomial. There's only one root because it's a linear polynomial because it's a one step formula.
The root is one plus k a. When is that smaller than one in absolute value?
Well, one plus K is smaller than one of K is within one of the point minus one.
So there's our stability region. We're in the car.
Plane. The circle of radius one around the point minus one in the K plane is the stability region for the Euler formula.
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,
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.
B that's gets a linear theme, more complicated schemes have more beautiful pictures.
And let me show them to you just on the computer you can see here.
We're summarising 50 years of analysis of the schemes.
So I have a code here called and 31 regex which plots.
Some stability regions. I hope. So there's a lot of implication in these diagrams.
So let me talk about the three families illustrated there. We haven't discussed this one yet.
We have discussed these two.
So what you see here are pictures of stability regions for the first three atoms, vast fourth formulas and the first four rugby.
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.
Fourth formula has that stability region. We didn't arrive the red one, but we computed with it that atoms back forth two.
So you see that as the order of another is that the fourth formula increases.
The formula gets a little bit more finicky because the stability region is shrinking down a little bit,
which means you probably need somewhat smaller timestamps.
So for the red one, you probably need half the size of time step.
As for the blue, the green one half again, so slightly smaller timestamps in order to get good behaviour.
They're running a cut of formulas. Very loosely speaking, are similar, though the details are interestingly different.
The first one they had a formula is again the forward Euler method.
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,
and then the purple is the famous fourth order run the kind of formula.
So that purple curve is probably the most important curve in this picture.
That is the stability region for the fourth order run kind of form.
And one thing you note is that it has some nice behaviour on the imaginary axis.
It will do better with oscillatory problems whose physics is on the imaginary axis than the other.
Then at least the first order out of the batch. Fourth form.
As you increase the order of the rung, a cut of formulas notice the stability regions are getting a little bit bigger,
so they're slightly less finicky about step sizes than atoms vegetable.
Now these we're going to talk about in a minute. This picture is completely different.
It's a formula called backward differentiation. And I'm showing you the first six stability regions.
But what makes it completely different is that here the curves are the outer boundaries of the stability regions.
Here they're actually the inner ground. So the backward differentiation formulas have unbounded stability regions.
The blue one, which is the backward Euler formula.
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.
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.
Which is the end plus the end plus one equals V and plus k, f and plus one rather than n.
You find that the stability region. Is.
The exterior. Of the very same circle.
Sorry. I'm saying that wrong. It's not the very same circle. It's the circle on the right of the origin.
Apologies. So whereas these formula that is based forth and rugged cut need very small time steps to work this backward.
Euler method doesn't mean such a small time for that to work.
It's stable all the way out to infinity, and it will only blow up for tiny steps in that region.
I mentioned there's a lot of history, mathematical research that goes into these pictures, especially the ones on the right.
These are the boundaries of the curves for degree one, two, three, four, five and six.
Backward differentiation. It turns out the seventh one really looks different.
You can see the cyan colour. The light blue is bending towards the negative axis.
The seventh one crosses the negative axis, and that has big effects.
That means that people don't really use those formulas for orders higher than six.
Well, you're probably fairly confused at this point. So let me explain what we find from all this analysis.
We find that for some problems, these nice, simple methods are great, but that other problems need these fancier methods.
And those are the problems called stiff odds. So I want to talk now about the idea of a stiff old.
Talking. So if you go back to the 1950s, this man, Dahlquist, that I've mentioned,
established the basic theory of when these methods converge, when does a linear multistep method converge?
And it's pretty close to when does a wrong a kind of method converge? And the basic result of his stability theory.
Was. That all that matters is what happens as K goes to zero.
Remember I mentioned the theorem, that convergence.
Equals stability plus consistency. And any formula you derive intelligently will be consistent.
That just means it has order of accuracy, at least one.
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.
As all fruits.
Of the recurrence relation. Yeah our stable in the sense we've been talking about in the limit of zero timestamp.
So his great theory says that you don't have to look at these pictures.
All you have to do is look at the origin.
And so long as the origin is in the stability region, on the boundary of the stability region, your method will work.
Of course it's the theorem. So of course it's true.
If you have any stable, finite difference approximation to an O.D. and you use it in exact arithmetic without rounding errors,
then as the timestamp goes to zero, you will converge to the right answer.
But what happened is that very soon thereafter, people discovered that in many problems in practice they weren't getting convergence.
Or if they were getting convergence, the time step had to be really ridiculously small to get good answers.
And that's what the problem of stiffness was. So stiff problems.
Can be defined in many ways. One way to say it is there.
They tend to be problems with widely varying timescales.
And the effect of that. It turns out to be that you often need k very, very small to get reasonable solutions.
So it's a beautiful case of correct mathematics turning out to be overlooking part of what really mattered.
There was something going on in the problem that really mattered. And let me give an illustration of a step by step problem.
I say it's a problem with a widely varying time scale. Here's a simple example.
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.
And here's what I'll do. I'll say you prime equals -100 times, you minus cosine of t.
Minus sign up and you have zero.
Equals one. So that's a first order scaler.
ODI Of course, it's mathematically very simple and the solution can be written down.
It's cosine of tea. So solution.
U of t equals cosine of t.
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.
So there you are. So. The solution.
Is this cosine function? However, the equation has a much faster timescale in at this -100.
And that means that if you perturb the solution a little bit.
For example, if I started here instead of there, then there would be a very rapid transit.
Or if I started here, there would be a very rapid transit.
There's a physics in the problem that's 100 times faster than the timescale of the solution you're interested in.
So that's a classic example of stiffness.
If you're trying to compute this nice, smooth solution, you'd like to be able to take time steps on that timescale.
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.
We'd like to be able to take time steps that look something like this.
We don't want them to have to be 100 times smaller.
Let me illustrate that. So the illustration is the code called M 32.
And 32 stiffness. And what it does is solve this problem and I'm going to shrink the step size.
So that's the error when you solve it with a big step size of 0.25.
Now, of course, as you shrink the step size, you'd like the error to converge to zero.
So we'll cut the steps in half. And we're using here the atoms passport formula.
I cut the step size in half and the error gets worse. I cut it in half again and it gets worse.
So this is at its best fourth.
It's a convergent formula and yet you can see we're getting horrible behaviour and loosely speaking because of the one in there,
it's going to keep getting worse until we get down to a step size on the scale of 100.
Cut it in half again works. Cut it in half again.
Suddenly it's beginning to look plausible. It's still ten to the eighth, but it's heading the right way.
And now suddenly, look, the step size is now less than 100.
And suddenly I have good behaviour. So.
That's what often happens if you use firearms based on a problem like this.
With multiple time scales, you need a ridiculously small time step to make it work.
Here it's only 100, but you can easily imagine problems in the sciences with much bigger ratios.
And an example that people talk about a lot is in molecular dynamics.
Ideally, you have things happening on a femtosecond scale, ten to the minus 15 seconds that are in the physics,
and yet may be what your tracking is on a nanosecond timescale or a microsecond timescale.
You don't want to have to take one step outside ten to the -15.
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.
The physics, for the most part, moves at a relatively slow convective scale, and yet the equations allow for sound waves.
So that's a very different timescale. Of course, sound waves move very fast.
You don't want to have to take time steps that are on the scale, you know, small enough to resolve the soundwaves.
And yet, if you use the wrong numerical schemes, you might have to do that.
So these questions of stiffness are all over the place with ease and.
PD The basic solution I'm going to mention are these backwards differentiation formulas.
So let me show you that if you look at the code.
You see that the first block of code this is m 32 is using the second order Adams batch, fourth formula.
The next block of code is using the second order backward differentiation form for the big time step.
It's doing better. But then here's the amazing thing. Every time I cut the time, step in half, it just converts.
So you can see some formulas have totally different behaviour than others.
And the reason for that has to do with the stability routine that I plotted earlier.
So let's do that again. In fact, if I say FIG. and 31 regions.
So this is the second order backwards differentiation formula, which is the red curve in the right there.
The stability region is the exterior of that red curve.
So it's virtually the whole complex plain and in particular along the crucial negative, imaginary negative real axis.
It's it contains the whole negative, real axis.
So the second order of backward differentiation formula is, in some sense, it's going to behave strangely for all timestamps.
Whereas the second order at best fourth is different.
So anyway, this got discovered in the 1960s, largely.
So among the mathematicians, Dahlquist wrote a his second grade paper in 1963, in which he talked all about this stuff.
But it turns out that actually a crucial paper had been written a decade earlier by some eminent chemists, Curtis and Hirschfeld her.
Hershey Felder in 1952. And they wrote a paper which is really very special, in which they had actually pointed out this phenomenon.
There's no general convergence theory like Alquist, but they saw the point that things go wrong for certain formulas and not for others.
They tend to think it's always impressive when a paper in the sciences does more than one great thing.
So they identified the problem of stiffness.
They use that word, and then they introduce backward differentiation formulas to solve.
The word seems to have been coined by John Tukey, the great statistician who also coined the word [INAUDIBLE] and the word.
Well, maybe he coined the word software. At least he was almost the first person to use the word software.
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.
Do not include F and plus one and implicit formulas do include.
F and plus one. Now, why is that so important?
Well, if if this is in the formula, then you've got more work to do.
So, for example, the forward Euler formula is the end, plus one equals the M plus k F1.
The backward Euler formula is V and plus one equals v, n plus k and minus one plus one.
To implement that. To find the new value of V,
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.
So to implement an implicit formula, you have to solve a system of equations that every time step.
And in general, that will be a nonlinear system. So that's the downside.
More work is involved to implement them, but the upside is that they have potentially spectacularly better behaviour.
So for a fifth problem, it's worth it. For a non-native problem.
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,
you'll see that a little algebra has solved a scalar equation to get.
To solve this. But this was a linear problem, so that's easy for linear problems.
Of course, solving system of equations is no big deal. For nonlinear problems, you need some kind of iteration.
So that's the. Conclusion.
If you have a non sniff problem, probably explicit methods are best.
Because they're so easy to implement and work well. If you have a stiff problem, you probably want to use implicit methods.
E.g. backwards differentiation. When it comes for ODIs, this is a pretty good summary of the whole picture.
For PETA is the same issues are all there but everything is of course more complicated and they're in the land of pity.
You have all sorts of technologies that are sort of in between these two.
They're semi implicit, trying to find some compromise that enables you to take reasonable sized time steps without doing too much work
with pieces that the matrices could be so big that solving the fully nonlinear system would be impractical.
Let me show the next demo to illustrate this.
So Code three is illustrating the importance of this solvers by solving the Vanderpoel equation with two different sorts of methods.
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.
If I take ten then the solution looks like.
That. So you've seen things like this before.
Now, is that a stiff problem? Well, visually, it looks stiff because there's a very rapid transition compared with the other stuff.
So that is indeed a somewhat stiff problem. Different time scales.
And yet, if I solve it with MATLAB, stiff solver te23s,
which is a backward differentiation code, it takes a lot longer than using the non-state solver.
So this is not stiff enough to make much difference.
Suppose I do it again? With a constant of 100.
Now this is stiffer. We have the very rapid transition and then the slow timescale.
It's not really safe there.
You need small time steps there, but throughout this regime it is very safe because the equation has this big number 100 in it.
But the solution you care about is moving slowly.
Now you find a stiff solver is about three times quicker than the monster solver, which is not very exciting.
But there it is. However, let's go further.
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.
But then most of the time it's very stiff with a number 1000 in the equation that isn't reflected in the solution.
And now you see this big difference. It's 40 times faster using the backward differentiation on.
So if this matters hugely when you have big ratios, thousands or millions in your timescales.
The final thing I wanted to do was show you some more of the fun in action.
I think I won't do that because we're at time, but you'll see a little code at the bottom there.
Just to illustrate to have fun solving a boundary value problem for the homework assignment do in 11 days,
you're welcome to use MATLAB or have fun, whatever you like.
Personally, I find cheap fun much easier, but it's not appropriate for me to force you to do that.
By all means, use MATLAB.