Last term, we talked about linear algebra and optimisation.
And this term, the other big half, if you like, of scientific computing is differential equations,
in particular ordinary and partial differential equations.
As always, we're going to go at a great rate, not very theoretically, but with a lot of hands on material.
There's a new URL for the course. That's how it works around here.
So all the lecture notes will be deposited there along with the assignments and so on.
The first assignment is due next Tuesday. As usual, the first one is smaller than the other three.
So the same pattern as last term, four assignments, the first one a little bit smaller, and then a mark at the end.
There is a pink handout, which is a variation on last term's blue handout.
Let's just quickly go through it to remind you of how things work.
So the course outline, as I say, is ODS and E's.
And in particular, I'm interested in easy methods for high accuracy and simple geometries called spectral methods.
The lectures me and then the two Tas are Adrienne Mont Anelli and Mikael Savitsky.
Are either of you here? I guess not. Okay. 12 lectures Tuesdays.
Fridays at 9:00. I've reminded you here under software tools of the collection of stuff we have online linking to many different things.
We won't talk about that again explicitly, but it's always there.
And remember, there are a couple of hundred pointers in that Web page, too, all sorts of interesting things.
So if you're new to the course, I definitely recommend you go to that Software Tools page and explore.
Regarding the assignments, the same pattern as last term as usual.
It's nice if you publish the things in MATLAB published, even though that's not the nicest system on earth.
But that's not required. What is required is that you turn in good solutions, including the programs you used,
of course, and that you don't just include programs and print outs.
It's crucial that you say a few words to explain what you've done.
One person left, turn by mistake, turned in assignment for it was 310 pages long.
He did this online somehow, so be sure to take a look at your assignments before you accidentally turned in 310 pages.
Okay. The last thing on the page is to remind you that the lectures are available online.
I realise most of you by definition don't need that, but they're there if you need them.
We haven't yet gotten around to putting last term's lectures in a publicly available place, but that's going to happen one of these days.
Okay. Any questions? Okay.
So let's begin with talking about. Ordinary differential equations.
So we call this part four of the course and I call that ody ease and nonlinear dynamics.
You know, those of you who've been here last term that I'm very interested in history always.
And I think historically speaking, not linearity was the most important development in mathematics in the last 50 years.
Essentially, computers got invented, and that made it possible to explore things that really couldn't be explored before.
So, for example, chaos is one of the many outgrowth of that new way of working with mathematical problems.
So with computers, we can explore all sorts of things, and the basic ideas we use are an ordinary differential equation means that you
have one independent variable as opposed to a partial differential equation,
which has more than one independent variable. And when I speak of an IVP, I mean an initial value problem as opposed to a boundary value problem.
So the starting point of ODI is a discussion of ODI vs.
And the basic thing to be said about those is that they can all be put in the same form.
We can always write. You Prime equals F of t u.
That's our universal standard form for an O.D. and we might as well suppose that it's posed on the interval to bigger than zero,
either the whole semi infinite interval or a finite part of it. And we have an initial condition for an initial value problem.
You not equal some prescribed data, you sub zero.
So there's the initial data for our IVP.
So that looks like a very specialised form. For one thing, it's a first order problem.
For another thing, it looks like a scalar, but it doesn't have to be a scalar.
The same notation works for systems of equations, so you could be either a scalar or a vector.
So this is scalar or vector, and I'll often call that an end vector, i.e. a system of equations.
And it's the usual story in scientific computing that sometimes you get systems of equations because you have many different things you want to track.
But other times you get systems because your disk criticising something continues.
So we saw last term that matrices arise of dimensions in the millions because we're doing finite difference or finite element approximations to eve.
Well, similarly, if you desc or ties at time dependent partial differential equation by sampling it in a lot of points,
you get a system of ordinary differential equations and that can easily have dimensions in the thousands or more.
So people really do care about systems of equations with very large dimension.
So it can be a scalar or it can be a system and it can be first order or higher order because.
Whenever you have a higher order problem, you can introduce additional variables to make it into a first order problem.
So, for example, if I had well, let's write one down.
Here's an example of the simple harmonic oscillator. The simplest ODI of all.
Let's write it with a w I could say w double prime equals minus w.
And I could give it some initial data. So for example, suppose I had w not equals one and w prime not equals zero.
So there is a simple example of a second order linear differential equation,
but I could turn it into a first order problem by defining a vector u with two components.
So I could define the first component u one to be W and the second component you to to BW Prime.
And then I got a system of two ordinary differential equations which I could even write in matrix form.
You one, you two, prime equals.
Well, because it's linear. I could write it as a two by two matrix.
It would be zero one minus one zero times you one, you two.
And there would be initial data. You need a vector of initial data so that the vector would be you one, you two.
At the time, zero is equal to one zero.
So any system of higher order odds can be written as a system of first order oddities.
That's amazingly convenient for the theory that most basic theorem of all about existence and uniqueness is always
proved in this first order context because it immediately translates into something for higher order problems.
At the software level, it's a subtler question. You can, in principle solve any O.D on a computer by translating it into a first order problem.
And that's certainly the way people do it 90% of the time, but it's not necessarily the best.
So if you really want optimal computations with a second order opening, maybe you don't want to reformulate it as a first order system.
Okay. So there's our simplest example. I emphasise that if it were a nonlinear problem, we couldn't write it in this matrix form.
We would need a function F, which would be a nonlinear function written in some more arbitrary way.
Let me also emphasise the word autonomous.
So when you say that it is autonomous, if it's independent of time now, what do I mean by independent of time?
Of course, the solution depends on time.
But I mean that the equation is independent of time, so there aren't any coefficients in the equation that depend on t.
But let me also make explicit what's implicit.
I guess that in applications, a lot of problems involve time and a lot of problems involve space.
Of course, mathematically it doesn't matter, but almost anything you do with odds, the independent variable is going to be one or the other of these.
So you could, for example, use a T here and an X here.
I won't be very systematic about making a distinction like that.
Another difference is that problems dependent on time usually not always tend to be IVP.
The initial value problems and problems dependent on space fairly often are boundary value problems.
But that distinction isn't absolute. Of course, there are boundary value problems in time, and there are initial value problems in space.
Here's another example. A nonlinear one. Consider the solar system.
The solar system, of course, has many bodies in it, but you could simplify it and say that it's got the sun.
Plus nine planets. Maybe it's eight planets. Maybe ten planets.
Nobody seems to be sure how many planets there are. And we could privilege our particular moon because that's a very nice moon.
So you could say that the solar system has 11 bodies. So 11 bodies.
But each one, you could say, has three coordinates each.
But each one also has three velocities because.
Newton's laws for the solar system are written in terms of force equals mass times acceleration, and acceleration is the second derivative.
So if we wanted to write the solar system as an oddity in the simplest fashion, we would have 11 bodies, each with six coordinates.
So there would be 66 dependent variables.
Assuming we're treating all of the bodies as simply points orbiting each other in the obvious fashion.
Or you could do it as a second order formulation with 33 independent variables.
Or you could pretend it's all in a plane, and then you'd have 22 independent variables.
But some of the planets aren't really in the plane. Okay.
Let's do an example. We will compute with everybody's favourite.
Simpler, nonlinear example of an everyday being is the so called Vanderpoel equation.
So the Vanderpoel equation is an oscillator like the simple harmonic oscillator,
except that it's damped and it's also nonlinear because the damping depends on the amplitude.
So whenever you have a physical process that depends on amplitude, then that probably means it's a nonlinear problem.
So it's a nonlinear oscillator. And even with without driving down the equation, I can describe it in words.
So it likes to oscillate. If the amplitude is large, then it is adapt system.
And so the amplitude decreases. If the amplitude is small, it's negatively damped.
So the amplitude increases. And so there's a preferred amplitude of the steady state solutions of the Vanderpoel equation.
So let's write it down w double prime plus a positive constant times w squared
minus one double w prime plus w equals zero and c is a parameter that is fixed.
So if we didn't have this term, it would simply be the harmonic oscillator.
But we do have this term. As long as this part of it is positive, you have damping.
Whenever this is negative, you have negative down. So if W is bigger than absolute one, an absolute value, it's damned if it's less than one.
It's negatively back. And we're going to compete with that in a moment.
What we find is that as time goes to infinity, the solutions converge.
To something of amplitude, approximately one.
And that's called our limit soccer. Now, of course, we can always write a higher order problem as a system of first order problems.
So let's do that one more time. Maybe the last time I'll do it.
If we write it in first order form, then the Vanderpoel equation becomes U one prime equals U of two, so w prime equals w prime and you two prime.
Is equal to minus u one minus the constant times u one squared.
Minus one. You too.
I hope I've written that correctly. Now, if I talk rather quickly and breezily about OBE as if it's a pretty easy subject.
The reason for that is it actually is a pretty easy subject. The odds are straightforward.
They're one of the most important things in mathematics, but they're not one of the difficult things in mathematics.
Let's play around a bit now on the computer.
Now there's one of the handouts, has some code on it, and the first bit of code at the top of that page is just a reminder of anonymous functions.
Such a convenient trick in MATLAB. So let me just remind you how anonymous functions work.
The sort of thing you can do is say F equals out of x x squared.
So that's an anonymous function. And I can now do things like F of three or F of three, four or five.
And MATLAB does the obvious thing. The standard way that the MATLAB guys work is to use anonymous functions to do various things.
So, for example, if you want to integrate an anonymous function or an M file, there's a code called Quad.
So Quad implements a numerical method which takes this anonymous function as
input and does the right thing to integrate it over the interval from 0 to 1.
As you know, I prefer the fun way of thinking where all functions are in the same setting of funds.
This more traditional way of thinking is that different operations like Quadrature Odds or whatever all have different interfaces,
but the net result is similar. So here's another example of a different interface.
Suppose I set a equals three. I'm just following the thing here.
Easy plot. I could now make another anonymous function.
So I'll say easy plot of A2 of x square root of abs, of f of x minus eight.
So you see, I'm piling up anonymous functions on anonymous functions.
Let's plot that. On the interval from -10 to 10.
What's happening? Yeah. So MATLAB has this command called plot or A-Z plot, which confuses people easy plot.
And it does the right thing to sample your function and draw a picture.
So that's an example of another interface to an anonymous function built into MATLAB.
Finally gets to emphasise the vector side of things. Suppose I made an anonymous function g.
Of several variables, so I could say g equals anonymous function of ABC and then I could say A two times B three times C.
So if I now say G of 111 with three arguments, I get one, two, three, or I could say F of g of one, one, one, and I get the squares of those.
So it's just a convenient way of doing a lot of coding in a single line.
Okay. So let's now assume we're good at anonymous functions and start playing with words.
The next code called Vendor Call. I'll just say a word or two about it and then we'll run it to see what the vendor policy later looks like.
So. M 25 Vanderpoel You see, first of all, we declare the initial condition, which is .01 for the function value and zero for the derivative.
So that's a vector of two components because we're going to convert it to first order form.
Then we give it a time span which goes from 0 to 50.
Then we input our damping parameter C and then we call one of the standard codes in MATLAB called ODI E four or five to solve this.
So it solves the equation defined by a function F which and you see there's an anonymous
function at the top that defines it on the given time span with the given initial condition.
And then it plots the result and it also makes a phase portrait. So let's run that to get a sense of things, I'll say.
And 25 Vanderpoel and it asks me for the damping constant.
So let's start in the middle and say one. So it's now chugging away.
And this is the result. Plenty of physics in that picture.
What you see is that. The initial condition is, of course, much smaller than a size 1.0.
So we have negative nine. So the negative damping makes the thing get much bigger.
And obviously, it's approaching a steady oscillation.
That's the limit factor whose average amplitude is about one that we see at maximum amplitude is about two.
If we press return, there's a pause. Now, if I press return, it will give me a face portrait,
which means we've now plotted this trajectory on the plane with you in that direction and you prime in that direction.
So as with most oscillations, you have an interplay between the kinetic and the potential energy, if you like.
Sometimes it's got a large amplitude and a small velocity.
Other times it's got a small amplitude and a large velocity, and it's moving around in that phase.
Let's change. See? To get a sense of dependence on that damn thing.
So if I change it to ten. So now we have much stronger damping parameter.
Then you can see the oscillation becomes much more visibly nonlinear there.
People have studied this a great deal. It's a fascinating case of slow motion and then a very sudden flip to the other state.
The face portrait you see also looks interesting.
Now at the other extreme. Suppose I run it with damping parameters zero.
Now I really do have the simple harmonic oscillator. The initial condition is point one and it just keeps going at .01 for at.
And the face portrait is a circle of radius point at one.
Finally, let's run it with a small but non-zero damping parameter.
So now we have C equals point three and what you see there is again, the negative damping makes it increase but more slowly.
Eventually it's reaching some kind of a steady oscillation and the phase portrait similarly starts
near the middle and it takes a lot of rotations around before it finally gets out to the limit cycle.
You see in the curve, it looks jagged. That's because MATLAB is by default plotting.
Piecewise linear is between data points. It doesn't mean that the data are that inaccurate.
They're probably quite accurate. But the plot doesn't happen to sample enough to show the accuracy.
Okay. So that's our simplest example of an oldie. And now let's say a bit about how these things are done.
That's a reasonable. Okay.
So the idea of solving Obis numerically is very, very old and it's one of the most concrete things in mathematics.
The basic idea, the starting point of an idea is that you're following some trajectory.
The body tells you the slope at a point.
So it's naturally enough to think of going a little bit in that distance, and then we get a new slope.
And it's the most obvious thing in the world to numerically do some kind of thing like that.
To use the Audi E as the starting point of a numerical approximation.
And this way of thinking was made famous by Euler, along with so much else that's going back to 1768, I guess.
I wouldn't say that Euler did a lot of really numerical mathematics, but a century later people began to do more.
And I want to mention the two big categories of numerical methods for obedience.
They're called wrong cutter methods. And multistep formulas.
These are the two big games in town. And they've been around for a long time.
So the multistep idea is really from the middle of the 19th century.
And the key guy there was John Adams.
In the 1850s. So about well, in the 1850s, that's when he did the work.
He didn't really publish it then. Adams was a great genius of numerical computing at Cambridge.
He was very young when he started this work, and he basically took Euler's low accuracy idea and looked at higher order accuracy analogues of that.
And these are the multi-step methods.
So this is the same guy who predicted the existence of the planet Neptune based on perturbations of the planet Earth.
In 1845, he was 26 years old at the time.
Sadly, the telescopes were not trained where he said they should be trained.
So a year later Neptune was actually discovered by La Vecchia in France.
Adams was an incredible genius. He was the senior wrangler, the best student at Cambridge in his year.
And in the examinations, the second best student.
This is back fourth got exam marks that were half those of Adams.
So he was really simply an extraordinary person who was much loved and certainly widely admired.
He was offered the post of Astronomer Royal and turned it down.
He was offered a knighthood and turned it down. Amazing guy. This guy best for his roommate.
I like to think they weren't actually roommates, I suppose, but they were pretty much college roommates.
Translated into the modern era, actually wrote down a lot of this stuff 30 years later in 1883.
So these this is the origin of the multistep methods.
And then it was a generation or so later that the wrong a cut of message came in.
So that was Carl Reiner, who is one of my heroes of numerical computation, did many, many things.
That was in 1895. He invented the fast Fourier transform, though nobody noticed, and many other things.
Another key name was Heine in 1900, and cut to himself was, I think, in 1901.
So basically that was the Brits and then the Germans.
And incidentally, regarding pronunciation in German and in English, the right thing to say is wrong to cut off or run a cut or something.
If you say runs, that's a complete error. However, in French they do say runes or something.
So if you're speaking French, you're allowed to say runes. Otherwise, I don't want to hear anyone say.
Right. Let me mention another interesting name.
Forest, Ray Moulton.
So there are formulas called atoms, best forth formulas, and there are formulas called atoms, molten formulas.
This guy, Molten, was an American, which was very unusual in those days.
America wasn't doing much academically in 1926.
And I think he may be related to Derek Moulton, who's a faculty member in this department.
So Derek isn't sure. But he says he had a great, great uncle who was supposedly very good at mathematics, and he thinks that may be the very good one.
One more name to mention, an early contributor is the very famous physicist Richard Farmhouses, who did many things.
And one of those was to get quite involved in these formulas in the 1930s.
Okay. So as I say, there are two big classes and.
The distinction between them is that the multi-step methods mar using formulas that coupled together many different steps.
Their only academy methods are, in some sense, one step method.
So now let me say what I mean by that. So what we always do in all of these methods is take the time, access and disparate it.
In the simplest setting. You imagine a uniform is for possession, though in practice, software will make that non-uniform.
So you imagine at time zero. You have t not.
And then t one we could call K. So K is our timestamp t two is two k and so so k is the timestamp.
Now you can do what Euler did and just can predict the future by little straight lines of length k.
Of course that works in some sense, but it's not an accurate method if the only converges linearly.
If you cut the time, step in half. The accuracy only doubles, which is not so great.
So the point is to get higher order second, third, fourth, or higher order methods.
So our goal in any of these formulas is to approximate the true solution.
You at the time, let's say t seven. Buy some number, which I could call this seven.
Which we're going to compute by some algebraic approximation to the differential equation.
So this will be computed somehow algebraically.
I tend to talk about numbers, but they don't have to be numbers in general.
You as a vector and V as a vector of the same size.
So we're approximating the vectors at these time values by.
Other factors that are approximations. T seven of course refers to ten times the times.
Let me just tell you the first four atoms formulas.
Or more precisely, the first for what are called atoms dash forth formulas.
So we're going back now 130 years. And the first one is Euler's formula.
So we write it like this. We say V and plus one equals the N.
Plus Hey Times fan and I better say what FM means.
F1 is an abbreviation for the function F evaluated at time t, n, and the current estimate of u.
V. So this formula is exactly this picture.
You're going along estimating your solution. You're at some point in order to do the next one, you take the current slope,
you apply F to your current estimate and march then a distant K in the time direction.
So that's Euler. And everybody calls that Euler's method.
But then Adams realised that if you play around with polynomials properly, you can get any order of accuracy you want.
So the next one is the end. Plus one equals VM plus k over two.
Times three F and minus F and minus one.
Now you're already seeing the multi step aspect of this to get the new value it's coupling together.
The current and the previous value and the order of accuracy has just improved by one order.
If you work with this one, you'll find its first order accurate. So the errors look like, okay, this one is second order accurate, case square,
and then I'll write down the next two, which are third order accurate and fourth order accurate.
So the third order Adams batch fourth formula is V and plus one equals V and plus k over 1223 times F,
n, -16 F and minus one plus five F and minus two.
And then the fourth order one is V and plus one equals v n.
Plus K over 24. You can see factorial hidden there times 55 F and -59 F and minus one plus 37.
And minus two. Minus nine and minus three.
And so on. It's an infinite family. You can take these forever.
So there is a 1,000th order out of the passport formula.
I'm sure nobody's ever written it down. You know, I like concreteness.
So the reason I write down these numbers is to make sure you have a sense of how concrete it all is.
Obviously, the numbers are not in any deeper sense. Very interesting.
But you can imagine there's a lot of very pretty mathematics designed to automate generation of such formulas.
And if you go back to the 1920s, 1930, 1940s, people were really good at that.
Once computers came along, the interest moved into somehow a little bit higher level questions.
Okay. So the good thing about add ons that's worth formulas is they're a very straightforward,
infinite set of formulas of whatever order accuracy you could ask for were a bit of a drawback is how to start them up.
Because in order to get the end plus one, I need not only v m but the previous three values with that form.
So if I want to actually use this in software, it's not enough to have an initial condition.
I need three more levels somehow,
so I need to use lower order formulas with very small time steps or something like that in order to generate the first few data values.
That's a little bit annoying on a uniform grade,
and it's more annoying if you're varying the time step because every time you vary it, you have to do something clever again.
So that's why I wrong and could have decided to construct a different set of formulas.
And let me show you some examples there. I'm going to show you two of.
So in the case of Adams formulas, I can say these are the first four and it's an infinite series.
Everything is conceptually very straightforward.
Wrong kind of formulas are not at all like that. There's no natural series.
There's more than one of any given order. They have all sorts of complexities.
You can't assume that the seventh order formula has seven stages in it, and you'll see what I mean by a stage.
Much more complicated and therefore maybe interesting mathematical theory.
The two I'm going to show you are a standard second order one and the standard fourth order.
So. Okay. Squared and that. I love K to the fourth.
As I say, there's more than one. But the one that I'll write down here is called Modified Euler.
And it looks like this. We take a look and define it to be the timestamp times F evaluated at the current point.
So that's F of key and the. And then we take B and define it to be the timestamp times f a value that t n plus a half of k.
And then plus a half of it.
So these are called stages. We've now done two stages.
We've used the function F twice, and then we combine the stages in modified order like this, we simply say the end plus one equals V and plus B.
So what's going on intuitively is that whereas Euler generates a slope from the left hand of the endpoint, thereby introducing I know of Keira,
this little bit of cleverness gives you some approximation to the average slope over the interval,
thereby giving you o of k squared and then you can crank up the order.
So here is the famous fourth order and another kind of scheme.
It's not the only fourth order running scheme, but it does have the property that it's often called the fourth order task.
And that one looks like this. And this one is very famous. I'm sure many of you have used it.
It has four stages. So you say AA equals K times F of ten then.
So we're starting as if we're doing Euler, but then we improve it.
We say B equals K times F of ten plus k over to V and plus over to.
And then amazingly, we do sort of the same thing. We now say see equals K times F of ten plus k over to V and plus B over to.
And then we do a fast one. We say D is equal to K times F of t, n plus k.
So now we're evaluating the function at fully the new time step V and plus C.
So we've done these four stages and then we put them together.
V And plus one equals v n +16 times a plus two.
B plus two C plus D.
So you can probably guess by looking at that that the analysis of these things is a little bit complicated.
It wasn't until the modern era that people were even able to generate a kind of formulas of order seven and eight and so on.
And there are all sorts of open conjectures about the kind of formulas.
The theory is actually quite beautiful. It involves trees in the sense of graphs.
And the person who is the great man who put the theory of wrong kind of methods into mathematics is called John Butcher.
He did this starting in the 1950s. He's still alive. He's a Kiwi from New Zealand.
Fascinating map. If you had one formula to take to a desert island, this is the one.
It's complicated to write down, but it's simple to use because there's no start problem.
You've got to do it. And it's got fourth order accuracy, which for many applications is very good accuracy.
So I'm curious. Raise your hand if you've used that form. That's amazing.
Okay, so here's a historical note.
Another great figure in numerical owning is Gerhard Vollmer from the University of Geneva.
And I heard him give a talk once about multi-asset versus market.
And he made the observation that over the years, Americans and British people tend to prefer multi-asset,
formal and Continentals tend to prefer kind of formulas.
And he actually argued convincingly that this difference ultimately goes back to the fight between lightness and Newton about calculus.
I forget the details of the argument, but somehow the Newtonian approach sort of led to this through Adams,
if you like, and this continental approach, the line, its approach through the Germans led to this.
Maybe true. Okay.
Let me say a word about the software available in MATLAB.
So this is officially Section 4.3.
Initial value problem codes in MATLAB and simulate.
So first of all, let me say all I'm going to say about Simula. You know, MathWorks is a very successful company with two or 3000 employees.
It all started from MATLAB, which was rooted in matrices.
But the real money-maker is differential equations and dynamics,
and the way most of the money is made is through this thing called simulant, which would be called a an environment for visual programming.
So in simulation you have you draw something on your screen and you connect various boxes and suitable ways to to construct the dynamics.
People like me, numerically, analysts never use simulated.
So I have actually never used it in my life, except once I tried the demo.
However, there's no denying that for many people it's the way to do it.
So I'm curious as you sit here. All right, have any used it a lot?
Yeah. So one day you'll give us a demo. Okay. But anyway, numerical analysts, mathematicians are a little more old fashioned.
We like to know more of what's going on. So we like to call codes to do things.
Let me just say 3 minutes about software. Oh, these are really the great subject in which numerical software took off 50 years ago,
because it has just that right level of complexity where nobody is going to solve an nobody on paper usually.
But to do it numerically is not that hard. So in the sixties, people started writing adaptive code.
So what does adaptive mean? For 50 years they've all been adaptive.
And there are two senses in which that word applies in the time steps or the space steps.
So that's always true. Basically, every ODP code in existence varies the time stamp depending on how the curve seems to be evolving.
We might say a little more about that later, but not today. It may also vary in the order, but not all of them do that.
So some of them will stick with a single formula like Ronnie Carter for something and vary the time step.
Others might go up and down a hierarchy of multistep formulas from order 1 to 13 or something, just to mention some classics in Fortran.
Going really back to the sixties. The two that you can find online of the oldest, best known codes are called Pack and Arcade Suite.
Easily found on the web, and they're also in our software tools directories.
A lot of these originate from the national labs in the US.
In MATLAB. There actually are eight codes that hardly any of us use, more than a couple of them.
And let me tell you the names.
The big three, really. r0de23034, five and O.D. 113.
These are both remarkable methods.
And the way they work, loosely speaking,
is that they run something with a fourth order method and a fifth order method and compare the two to get an error estimate to control the time step.
This is a multi-step one. Which actually has variable order going from order one to Order 13.
This tends to be my preference because I like high accuracy and the high order means that if you really want ten digits, this is usually best.
But if you only want five digits, this is fine. And if you only want two or three digits, this is hard.
Now, in the handout, if you look on the back of the code, I've copied the first couple of pages of the paper that introduced all of these codes.
It's by Larry Champagne and Mark Reichel Champagne.
I think at that time was at a national lab. No, he wasn't.
He was at Southern Methodist University, but he was connected with the labs.
Mark recalled was at MathWorks, I guess.
He's been at MathWorks for many years. He took some time off to found a company.
I think this is one of the two great papers published by the MathWorks team that made their software applicable to more people.
I've shown you last term. The other one about sparks matrices in MATLAB.
Well, this old paper described the suite of code that is still the basis of MATLAB computing 20 years later.
For your amusement. For my amusement. Notice the abstract.
And also the first sentence of the paper is a little curious.
The paper presents mathematical and software developments that are the basis for a suite of programs.
It doesn't say the paper presents a suite of programs.
Now, I don't know why, but I can guess it's the universal problem that publishing software papers is hard.
Nobody respects software.
So even when you have some software that a million people are going to use, you have to couch it in terms of mathematical development.
This is something we on the telephone team are very aware of.
We write a lot of papers and we're always careful to say that there's more going on in software.
And of course, there is more going on less software.
So let me finish guest have a couple of minutes by showing you how we can solve an initial volume problem in telephone.
The underlying methods are the same, basically.
In fact, our default is to call ody 113.
But the setting is all about the way of computing, you know.
Well. So let me just in the very little time available, just give you a sense of how this works.
Suppose I wanted to solve the Vanderpoel equation and said, fine,
I'm just going to do one version and I'll say C equals one, and then I'll say and equals what we call a chap.
So it's how about this? How we define a differential operator, we give it a domain.
So from 0 to 50 and the operator part of the differential operator is an anonymous function of time and the dependent variables.
So for example, for the Vanderpoel, it would be the second derivative of u comma two plus the constant times use.
Well, you don't need the dot u squared minus one.
Times differ view plus you so notice a key difference here from the standard MATLAB way of operating is
we're perfectly happy with higher order you don't have to reformulate it to have fun does that for you.
So there's a second order formulation, Chad. Fun ultimately is going to do it behind the scenes, convert that to first order.
Then, of course, you have to give it an initial condition and LBC stands for left boundary condition.
So if we wanted to set the function value in the derivative, we can do that.
And then to solve the problem we use backslash that we say and backslash zero.
So what we've just done is the first time you do anything, it takes a long time.
That's MATLAB. Let's do it again. Okay.
So as you know, in linear algebra, backslash is matlab notation for how you solve x equals B.
Well, in turbofan, we've overloaded that to solving a differential equation.
So what we've just done, we have just solved an of t you equals the right hand sorry times the unknown vector you equals the right hand side zero.
By methods I'm not talking about now. So if I played it, I get the same picture that we've seen before.
This is now a cheddar fan, so. It's actually a polynomial.
If I said the length of you, I'd see that it's a problem on your degree.
1070. Or if I said flat coefficients of you, I'd see that the Chevy shaft coefficients decrease down to machine precision like that.
And if I wanted to do the phase plane, I could say FIG. two and then plot you against the derivative here as I read my.
So the pictures end up very similar, but the setting is a little bit different.
Okay. See you on Friday.