Hand handouts. And as you can see, one of them is a questionnaire. I hope you'll feel that out.
And in particular, we're curious to know which bits of software that have appeared in this course you might have used in your own work.
There's also the assignment number three, which is do 11 days from today, Tuesday in lecture.
And it's basically MATLAB, though. You'll see that on the third question, the biggest one,
you're very welcome to use check find if you prefer, but you certainly aren't required to do that.
So today's a pretty amazing day in the news because of the gravitational waves.
I couldn't resist looking at the paper and having looked at the paper.
I couldn't resist making copies of it for you. It's truly amazing.
I'm sure you've all been hearing about this on the news. So one thing I note is that the actual paper is much more exciting than the news coverage.
On the news coverage, you see things like, you know, stills from the Big Bang Theory and stuff, which is very cute.
But to actually look at the paper is extraordinary.
So let's just glance at that for a moment. This is a historic occasion, the discovery of gravitational waves.
It's a beautifully written paper. You should really sit down and read this.
It's fantastic. This will probably turn out to be the most sensational scientific discovery of your time in graduate school.
Some of you are physicists. Some of you aren't. But this belongs to everybody.
This is huge. It involves certainly all of physics and all of mathematics and all of computational mathematics.
It's just a great achievement of international science.
If you just leaf through, I haven't yet read it all carefully.
But you can see, first of all, it says so much about how science, big science is done these days.
There are hundreds of authors and even hundreds of institutions.
It's all part of the Lego collaboration, which is international, although I guess the two main places are in the USA.
The beginning just describes perfectly the history of this problem.
And then if you turn to the second page, you get these sensational plots.
So well, actually, before you turn to the second page, go back to the abstract.
With its sensational beginning on September 14th, 2015, at 950 and 45 seconds.
The two detectors simultaneously detected pretty much the same signal.
So a billion years ago, two black holes merged in the course of a few seconds.
And that merger is energetic enough to shake the Einstein and gravitational space,
creating waves big enough to be detected a billion years later here on Earth.
So they travel for a billion years at the speed of light, and then they go right through the United States in these detectors.
So in these sensational plots, you see a curve from Hanford, Washington, and another curve from Livingston, Louisiana, and they match beautifully.
Incredible. I don't know any of the story behind this,
but I can sort of imagine the excitement in these two labs when they discovered they had each detected the same signal and then the two curves below.
I don't know any of the details,
but clearly the relativists have numerical models that predict what sort of a signal colliding black holes should send out.
So these predictions match the experiments magnificently.
It's quite incredible. I wouldn't say most of the credit here should go to numerical computation.
There's much more going on with instrumentation and real physics.
But obviously a part of the problem is numerics. And I don't know any of the details of what sort of equations they solve to get these amazing curves.
But it is clear from reading the paper that once they get the signal, they know what they're looking for.
And they were quickly able to look at the papers that generated such a curve and back out from that estimates of the masses of the black holes,
for example. So it's truly an extraordinary bit of science.
And let me mention one other aspect of it that relates to this course and to what we're doing right now.
Of course, we're talking about partial differential equations.
I mentioned that the first great example of a scientific advance that really only made sense in the context of a PD was Maxwell's equations.
So this is going back 150 years. Of course, light waves and radio waves already were there.
You didn't have to discover light wave.
But Maxwell, as a the crowning achievement of his equations, was able to show that the equations had solutions in the form of waves.
And of course, we know and he knew that those are light and radio.
Waves. And very quickly, all those pieces came together. Well, this is a much more complicated, much longer story involving Einstein.
So, you know, Einstein did special relativity in 1905, and that doesn't particularly involve these deformations of the universe.
But then general relativity involving gravity came along a decade later.
And in 1916, he realised that the solutions of his pdes should have, should include waves as possibilities.
Now his PD eves are infinitely pdes are infinitely more complicated than Maxwell's.
They're nonlinear, and they're ten of them instead of two or four of them.
But they are these and they have these wave solutions. And he predicted the existence of these solutions 199 years ago.
And. In the news, you'll see happy stories about what a great man Einstein is and what a great man Hawking is.
But I think the even deeper story is how extraordinary it is that physical and mathematical reasoning can lead to this kind of a prediction.
The equations suggested certain things, and a hundred years later, we discover those things are true.
It's just extraordinary that mathematics can have such magical ability to predict how things work.
So it's a very exciting day for all of us, and the rest of the lecture will be infinitely less exciting than what has just happened.
I'm curious. Raise your hand. First of all, if you're in physics, raise your hand.
Wow, that's impressive. And second of all, if you have any particular connection to this kind of physics, raise your hand.
Okay. Well, of course, this event is going to affect people's career choices.
That's complicated. I view with horror the thought of being on a team this big.
What's it like for an individual to be, you know, author 146 on this paper?
I don't know, but it's quite something. Okay.
So you'll note I included the whole paper, including the authors and the institutions, because somehow that's part of the drama.
I think now we're talking about a solution of relatively easy equations.
And today the first thing I want to talk about is to say a word about what you do when your grid is not a uniform grid.
So let's talk about finite different thing in general grids for.
Now by general, I mean, it's not just equally face points.
There are all sorts of reasons one might wish to have more complicated grid distribution.
One is that the geometry of the problem you're working with might be complicated.
You might want to do special things near boundaries. Another is that the solutions you're working with might have singularities or
something so that you want to adaptively refine the grids in certain areas.
So there are all sorts of reasons to use more complicated grids. So let's say our question is how should we derive a finite difference formula?
A finite difference approximation to, let's say, the case derivative.
On an arbitrary grid. And I'm thinking in one dimension, but of course, in multiple dimensions, the same question comes up.
So we have some grid and maybe it's gradually getting coarser, for example.
So there would be the kind of grid that could easily come up where maybe you have
a boundary here and you're resolving things more closely near the boundary.
And maybe you want the third derivative approximation at some point on this grid.
What's the right way to get that?
Well, there's a basic principle that nine times out of ten is the one people use, and that's polynomial interpolation.
So let's suppose, for example, that we want at that point an approximation to the case derivative.
Well, basically the standard thing to do amounts to picking a set of adjacent points.
You know, maybe these interpolating those data by a polynomial and then differentiating that interpolate.
So let's write that down. Here's the principle we're going to use.
So at each point, each grade point, we decide which data we're going to use.
So the data will go presumably to the left and to the right.
So I could say which data the sub j minus R dot dot dot up to v j plus s to u.
So in this case, this is the j grid point.
And this is the case where our equals two in this picture and s equals two.
Because I've gone two points to the left and two points to the right.
Then the principle is we now imagine interpolating that data by a polynomial.
So we interpolate the data by a polynomial and it would be the unique interpolate of degree r plus s.
Because if you have our plus X plus one data points.
There's a unique polynomial of that degree through that data.
And then we differentiate the internal and k times and evaluate the result at the point we care about.
In other words, our finite difference approximation. Is going to be the case.
Derivative of p evaluated at the point we're interested in.
So it's a very simple principle polynomial interpolation differentiation.
But this really is the basis of most of what people do when they describe ease,
at least when they disparities them by finite differences as opposed to finite elements.
So let's do an example.
Well, I should say that every finite, different formula I've shown you does fit this pattern, except I've shown you just simple ones on regular grids.
So, for example, if you have a uniform grid and you want to get a second derivative approximation based on three points,
you do exactly that if you interpolate by a parabola and differentiate the interpolate twice.
So let's work out a case like that. I'll do a different one.
Suppose I want the one sided first derivative?
So I want to approximate the first derivative. And let's suppose that I have equally spaced points.
But what's different is that it's a one sided derivative I want.
So it's an x not equals zero and x one equals some grid size H and x two equals two h.
And the motivation here would be that we have a boundary and a grid coming in from the boundary.
And at the boundary point we only have data to one side.
So we want a one sided formula. Well, we're going to interpolate by a parabola and differentiate the interpolate once.
So here's the interpolate. I'll write it in what's called the Newton form.
In turbulent through three data points is a parabola. And it turns out that equal to P of x equals the value at the left grid point plus.
V one minus v not x over h where v one is the value at the second grade point plus v to minus two v one plus the not.
Times X. Times X minus H over two x square.
So you can easily check that that polynomial fits the data.
It takes the value of not at zero and the value of one H and the value of two at h2h.
It's a quadratic polynomial degree one there and degree two there.
So we differentiate that and we find that P prime is equal to the one minus the
nought divided by H plus v two minus two v one plus v nought divided by three times.
To explain minus H. Over to a square.
The interpreter is a parabola, so its derivative is linear.
And then we just evaluate that derivative at zero if we want the one sided formula.
So we compute that this was P prime of X, we compute that p time of zero is equal to the one minus V,
not over h, and then when x is equal to zero, we get that thing with a minus.
So that's minus. V two minus two.
The one plus d nought divided by two h.
And that in turn is equal to minus a half.
The two. Plus 2v1 minus three halves the not.
Divided by age. So that's it.
We've just done the algebra and figured out the coefficients of the 3.1 sided first derivative.
And in fact, we've seen these coefficient before. Minus two have two.
Minus three has those are the coefficients of the second order Adams batch fourth formula because the Adam's back fourth formulas have this
one sided nature because they're explicit formulas and it's a second order one which relates to parabola interpolating in three points.
So that's why we see this. Now, what happens in computation, of course, is that you don't grind all this stuff out.
It's all done in advance a priori or on the fly using clever algorithms.
You don't write down anything on paper anymore. That's much too 19th century.
Let me point you to two remarkable papers in this business by Bengt Fineberg.
So. Surprisingly recently, that is to say, only a couple of decades ago.
He has two key papers. One of them in mathematics of computation in 1988.
And another one in Siam review.
In 1998. And the details in the names are in our lecture notes.
The first one is called Generation of Finite Difference Formulas on equally on arbitrarily spaced grids.
And the second is called calculation of weights and finite difference formulas.
These two surprisingly late find elegant formulations and stable algorithms to do a very general computations like this very rapidly and cleanly.
Of course, people would have known what to do 50 years ago, but using somewhat less clean and robust algorithms.
And I've got a handout from his first paper here.
So if you take a look at the handout,
you see on what I've got from the 1988 paper called Generation of Finite Difference Formulas on arbitrarily spaced grids.
And let's just glance at his four tables there.
So in table one. You can see it's centred approximations at a grid point.
So the notion there is that you have an equally spaced grid and you want to approximate a derivative at the middle.
So this is the most basic thing one would do all the time and you should recognise at least the first coefficient.
So suppose you want the zeroth derivative.
That's the first line of the table. Well, the zeroth derivative just means the interpolate.
So you interpolate these five data points by a polynomial of degree four and you evaluate that polynomial at this point.
Of course, you get the value at this point. So that's why the weight is just the number one.
The next four lines of the table are the first derivative.
So suppose you wanted the first derivative here based on three points or five points or seven points.
Well, here are the coefficients.
So, for example, here we have five points and that would be minus the 12th, minus two thirds, plus two thirds, and so on.
So he's offering these.
Partly because the numbers are interesting, but also partly to illustrate how systematically one can regard these problems of generating coefficients.
The Next Table. Table two is centred approximations at a halfway point.
So here the idea is that you have an even number of points and you want the value there.
Now, why would you want that? Well, the answer is that very often, going back many years, people have described using so-called staggered grids.
So that's the phrase staggered grids.
What you find is that very often in fluid mechanics, for example, you have several variables you're trying to solve.
For that, each depends on the other. And it turns out that, for example, suppose we had a grid like this, let's draw it in two dimensions.
So here's a two dimensional grid, and maybe at each of these points we have an estimate of the pressure.
Well, it turns out that the velocity depends on a derivative of the pressure, so it makes sense to compute the velocity at in points.
So if we were in one dimension, I would put a velocity grid right there.
In two dimensions. I might even put it here. Let's do one. That's simpler.
So a staggered grid might evaluate the velocity at these points.
So that the pressure is that some kind of a derivative of the velocity and the velocity of the derivative of the pressure.
That kind of staggering of a grid is an old trick that people use.
So of course, table two gives you coefficients for that and the top entry in the table we instantly
recognise as if you want to approximate the value here based on these two numbers,
of course you would take the mean of those two numbers. So we see a half and a half for the zero for the derivative.
On the other hand, suppose you wanted to approximate the first derivative based on these two numbers.
Then in the next four lines you see minus one and one and so on.
Now let's look at the other two tables. The other two tables are analogous behaviour added boundary.
So in table three, you see one sided approximations at a grid point.
And in table four, you see one sided approximations valuated there.
So of course, the motivation is what to do at a boundary in a simulation, and the motivation here is what to do at a boundary.
When you have some kind of a staggered grip. So for example, these formulas here, we can recognise that's in table three.
If you want, if you only have two points at which to estimate a first derivative, well, all you can do is take their difference.
But if you have three points, you've got the minus three has two minus a half.
And that's the second order. Adam's back fourth. And then I've also written down the next few lines because that's the third and fourth order.
Adam's back fourth. So -11, sixth, three minus three has one third.
Those are the Adam's best fourth coefficients for the third order form and so on.
Now these are all with regular grids.
But Feinberg's algorithms, especially in the second paper, give you very slick, robust ways to do it for arbitrary grids.
Now we're going to come back to these issues of how you construct derivatives.
Shortly when we turn to spectral methods. But for the moment, that's all I'm going to say about them.
And then the rest of today, I want to switch to multidimensional problems.
So this is 5.8. Multiple space dimensions.
Well, of course, computational science happens in all sorts of dimensions in terms of physical space.
We tend to live in a three dimensional world.
Sometimes one might regard that as a four dimensional space time continuum if gravity waves are at issue.
But even higher dimensions come up all the time.
First of all, in physics, for example, you could have a state space if you have a bunch of particles in many dimensions or in finance.
People are always working with multivariate databases involving different financial instruments.
There are all sorts of abstract context in which one would even like to work in a hundred or a thousand dimensions.
That stuff is actually pretty hot these days.
There's a lot going on in high dimensions, but more classically, 2D and 3D are the kind of simulations we focus on.
So a typical standard big problem of computational science might be a three dimensional,
fluid mechanical simulation of the high Reynolds number, where you might have several variables involved.
There would be three components of velocity, maybe under pressure.
All sorts of things might come in at density.
So you would very typically have a three dimensional grid with a number of different variables at each grid point.
And that's the kind of motivation I'm thinking of here. So there are sort of two facts, a good fact and a bad fact, if you like.
So the good news is. Little changes conceptually.
So 2D and 3D are based on the same mathematics as one D.
And what do I mean by that? Well, you can approximate derivatives in ways that are obvious extensions from one day to two or three.
You set up a discrete system of equations. You need to solve them.
All the issues we've talked about arise.
For example, if you have a hyperbolic problem with a finite wave speed, then you can probably use an explicit method.
No problem. If you have a diffusive problem with an infinite wave speed, you probably need an implicit method.
If you have a stiff PD, you might be able to mix together an implicit method for the linear part and an explicit method for the nonlinear parts.
All of that is still true in two or three dimensions.
The bad news is. Linear algebra.
And what I mean by that is all of these methods, as soon as you have a system of equations to solve, it's no longer going to be abandoned system.
So of course, in one dimension, we're used to matrices that look like this,
and those are so easy to solve that I haven't even discussed how one solves them.
But in two dimensions the matrices will get more complicated and in three dimensions much more complicated.
And so nobody would ever have predicted it. But linear algebra really is the heart of a lot of computational science.
Because of this, if all equations were hyperbolic, linear algebra would be less important.
But because we care about diffusive processes which require implicit miss, linear algebra is at the heart of what's done.
So this constrains the kind of simulations that are feasible in all sorts of ways, and it leads.
We've talked about this in the first term.
Of course, it leads to tremendously advanced methods, both in the direct category and in the iterative category.
And when you're doing iterative methods, you're always talking about preconditions.
So these are just huge areas of study.
Very roughly speaking, you might say that direct methods often win in 2D, but in 3D, the Sparsity structure is too complicated to take advantage of.
So maybe iterative methods win, but it's not really so simple.
Both are used in all sorts of way. So linear algebra changes the picture very much.
We're not going to say too much about that for the moment. We're going to focus on the part that isn't so different in multiple dimensions.
And let's take our example, the 2D heat equation.
So suppose we have you t equals the le plus c interview.
So that's u x x plus u y y and let's put it on a square.
So I'll say zero less than X and y less than one.
And since it's a very simple problem, I'll simply say that you equal zero on the boundary.
So we have some initial heat distribution which eventually is going to flow out through the boundary,
and we'd like to know how that evolves as a function of time.
Well, let's do the obvious discrimination.
So we'll take our grid, we'll set the values to zero on the boundary, but then in the interior, we'll set up a grid of points.
And we'll do finite differences. To approximate the equation.
So let me write down the obvious finite difference approximation of it there.
Well, if I write down the obvious backward Euler approximation, it could look like this.
We're going to set the iPhone X coordinate to be I times a space step and we'll set the j y coordinate
to be j times the same space step and we'll set the time position to be in times of time step.
And the notation I'll use is that h is one over some integer j.
So we have a j by j grid basically. So the backward Euler disparate ization.
I hope you can all write that down without my telling you. But I'll tell you, the stencil could look like this at the new time step.
We sort of have X and Y there. So this is timestamp and plus one.
We'll compute our derivatives there and then because it's backward Euler.
We'll do the time derivative along there. Back to step. Yeah.
So the formula will be that the time derivative is v and plus 1ij minus z n.
I. J divided by k.
So I hope the notations obvious the x coordinate, the j coordinate so their time approximation to the time derivative.
And that will be the sum of two terms, the double derivative in x and the double derivative in y.
So these will be the I plus one j minus two the i j plus the I minus one j at time step and plus one.
Divided by X squared. And then similarly the i j plus one minus 2vij plus the i j minus one at time step and plus one.
Divided by square. So there we have our approximation to u t equals u x x plus u.
Y y. It's a paradox of mathematics that the policy is, of course, isotropic.
It doesn't treat X and Y differently from 45 degrees, and yet it can be written in this form, which is completely nano.
So our discrete possession of it is of course not isotropic.
And yet we're discounting an isotropic operator.
So that's the backward Euler approximation and I won't work it out, but it's unconditionally stable.
By which we mean there's no time step restriction.
Of course, if the time stamp is too big, you'll get inaccuracy, but you won't get explosions.
And by the way, if we do forward Euler. Which would be everything the same except these and plus ones would be ends.
Then you do get a time step restriction which turns out to be K has to be less than H squared over four.
So backward. Euler, as usual, allows big time steps forward.
Euler requires small time steps, the accuracy of both formulas.
His O of Kate plus eight squared.
So if K is proportional to x squared, then time and space errors would be of comparable magnitude.
If K is bigger than the time, errors would be the biggest.
So if you think about what's going on, when something like this happens,
we're using matrices, except they're sort of four dimensional matrices, right?
If you have a1d problem, then you say vector equals matrix times vector, or maybe you solve a system of equations matrix times vector equals vector.
In a 2D problem. I need to draw a four dimensional picture, don't I?
It's a tensor. So what we'd like to say is that a matrix is equal to a some kind of a 4D tensor times another matrix, or to solve a problem.
We'd like to say that for D tensor. Time's matrix equals matrix.
There's no doubt that that's what's logically going on when we in the game like this, and yet nobody thinks that way.
This has always bothered me because it must be fruitful to think more like that than we do.
Yet when we design our computers, we don't think that way very much.
When we write down our algorithms, we don't think that way. Of course, people do talk about TensorFlow.
They're not unaware of this interpretation, but most of computational science turns everything into a matrix.
You take this four dimensional object and you string it out so it becomes a two dimensional object.
Of course, in 3D we would want to have a three dimensional cube would be equal to a six dimensional tensor, times a three dimensional cube.
So ignoring cancers and just thinking like linear algebra, it's.
Abstractly, we can write down that the backward Euler formula will have our usual format of matrix times.
Vector at the new time step equals vector at the current time step.
So this is so abstract that formula could apply to the tensor picture as well as to a matrix interpretation.
But certainly in practice, people think in terms of matrices. And B is equal to the identity.
Minus. The time stamp times our discrete LA Plus.
The street. LE Plus, is this part of the picture?
So the whole linear algebra turns out to involve solving a system which basically is a perturbation of the identity by the discrete.
And as we know from the first term in MATLAB, you can write the discrete in the plus in the form of a chronic product I comedy.
Plus the chronic product in the other direction d com.
I where d is equal to h to the minus two times the tri diagonal matrix whose nonzero entries are one minus two and one.
So I think actually these formulas are a complete description of the linear algebra of backward Euler.
You set up this J by J tri diagonal matrix.
Then you make this j squared by j squared matrix from the chronic r product or the sum of the two chronic or product is called a chronic or sun.
As we discussed last term, you then have this big matrix of dimension J squared.
You subtracted off from the big identity.
So let's say this is a J by J identity, but this is a j squared by j squared identity, the tensor identity, if you like.
That's the whole linear algebra. So this is iconic.
Her son. And chronicler sums are all about taking TEDsters and squashing them down into two dimensions.
Now the linear algebra side is all about the ordering of the unknowns and how you can exploit that ordering at least four direct methods.
So let me look at the ordering by focusing on the case j equals four where we have three unknowns in each direction.
So we now have our unknowns in a three by three grid, which is interior to the boundary, the boundaries out here with taking the value zero.
But the unknowns are in this three by three grid. So if I think of this one as the first and this is the second and this is the third unknown.
Five, six, seven, eight.
Nine. You can number them however you want mathematically,
but the way you number them will affect how you put them on a computer if you string them out in order.
So if I number them like this telegraph typewriter ordering one, two, three,
four, five, six, seven, eight, nine and then put them on the computer that way.
That implies something about the structure of the matrices that we'll be working with.
So, in fact, if I use that ordering.
Then the discrete le plus in. Is a mapping from a vector of length nine to a vector of length nine.
And that mapping. So let's see. I want a nine by nine matrix.
Right. One, two, three, four, five, six, seven, eight, nine. That nine by nine matrix will have block three by three structure.
And it will contain numbers corresponding to couplings between adjacent?
No. So for example, this point is coupled to those because of the second derivative in X and to those because of the second derivative in Y.
But it's not coupled to the because we're not doing any different thing along the diagonal.
So you end up with a block matrix and each block has a special structure.
And in fact, if you look at the numbers, it's minus four, minus four, minus four, one, one, one, one.
That's on the diagonal blocks.
So it's a what we would call a block tri diagonal matrix whose diagonal blocks are try diagonal on the super diagonal blocks.
We have little identities. And then of these three things, you have zero.
So the big structure is block diagonal, and then within each block further structure.
Of course, in 3D you get analogous things with more complicated structure and in end dimensions you'd have a sort of recursive structure of some kind.
So the wording is that we say L is black, tri diagonal.
With tri diagonal blocks. And for the last few minutes, I'd like to use that technology to solve a few problems.
So now let's look at the codes I handed out there. There are three of them 46, seven and eight.
The first two are for the heat equation, and the third is for the wave equation in two dimensions.
So looking first at em 46 backward. Euler It's a two dimensional backward Euler simulation for the heat equation and it will do exactly this.
And you see the key line is where I say B you equals B backslash u.
So in that line you equals B backslash u. MATLAB is being told to solve a system with this structure.
Now, this structure isn't that bad. It's not as narrow a band as you might like, but it's still a banded matrix with a bandwidth of on the order of J.
So. This is the band.
The band did this is the Matrix. Outside their everything is zero and MATLAB exploits that quite easily.
So this two dimensional simulation goes quite smoothly, even though linear algebra is involved at each step.
So let's run that. So if I type m 46 backward Euler it asks me what size grid I want and I'll begin by saying I'd like a ten by ten grid.
So it shows you The Matrix and you can see there are ten blocks each.
There's a block. This is the lowest, right try diagonal block with the ten by ten matrix.
This is the lowest right of diagonal block, which is a ten by ten identity.
So there's the structure of the whole matrix. It's bounded enough to be easy to solve.
And so every time I take a step. Solve the problem.
That's not too exciting. Let's make the grid finer. Let's say I'm 46.
I'll do it on a 90 by 90 grid. So there's again the Matrix.
It's of dimension 90 by 90.
So 8100. But plenty of structure.
So that's easy for your computer to solve. There's the initial condition and every time I take a time step.
The heat flow. There are a couple of things I always like to do to fiddle around with codes.
So in this one, you'll notice there's a line that says periodic boundary conditions with a percent sign.
If I remove that percent sign, then it will change the boundary conditions to be periodic.
So let's do that. It's called and 46 PR.
So if I now say let's try a 30 by 30 grid, there's the grid and notice same tri diagonal structure.
But now these things have appeared. Our matrix is now. It wraps around.
It's a circular matrix, but MATLAB still is able to solve that without much difficulty.
So the initial condition looks like that. And now, as time goes by, the physics is a bit different because it's periodic boundary conditions.
It's not exciting to watch, but you can see that at the boundary the temperature is rising a little bit, so we're not losing heat anymore.
It's merely equalising throughout the domain. But if you tried to take a big time step with this, you'd lose accuracy pretty quickly.
So the Crank Nicholson approach would be to increase the order of accuracy in time from K to k squared.
I won't write it down because you all know what it means, but the picture tells you all.
So the second code m 47, which incidentally is the year that the crack Nicholson paper was published in 1947, I think shows you the picture.
We're going to disparities, both of the new time step and of the old time step and treat those equally.
So the symmetry now gives us k squared rather than k accuracy.
It won't particularly appear in the pictures, but let's just run it for fun.
If I say m 47. There's very little difference in the code.
I've chosen a different initial condition just for variety, so you can see it works fine.
Let's do a periodic version of that. You can see a couple of lines change their.
Not much to say, but that works.
Maybe something more interesting is to look at the penultimate line of code.
Where it says removing the parents is catastrophic.
So we have you equals B backslash a time to use.
Why do we have that? Well, with Jack Nicholson, we're just criticising at the old time step.
That's the A. And then we're also disguising at the new time step.
That's the B. So we're solving a system of equations involving B, but the right hand side of those equations depends on it.
Now, mathematically, you could do things either way.
It's perfectly legitimate to talk about the matrix.
B Inverse time. Z That's a matrix. And if you asked MATLAB to compute that matrix, it will do it.
The trouble is that to actually compute that whole matrix can be very time consuming as opposed to solving a system of equations with B,
which is much quicker. So let's try that. I think I've called it m 47 slow.
So I've now removed the parentheses.
So MATLAB is just doing that line from left to right and trying to compute the matrix B inverse.
And if I say ten. No problem.
If I say 30, it's a bit slower.
I don't know if you notice the difference now let's say 60.
So now it's trying to compute. It's doing it. That's pretty impressive.
It's computed this matrix or actually at every time step it's computing this matrix of size 3600.
I won't try 90. That would really grind to a halt. Let's in our last 30 seconds, let's quickly do the leapfrog for the wave equation.
So there you have the two dimensional wave equation and we're solving it by a leapfrog formulation.
Maybe bigger grid is nicer. And here the picture says it all about how we're discontinuing Leapfrog is suitable for an explicit
it's an explicit scheme suitable for a hyperbolic problem with finite speed of propagation.
So we're evaluating the second derivative in the middle, and then the time derivative goes from timestamp and minus one to time,
step and plus one so that you get symmetry for the higher order, but also explicitness to avoid having to do any linear algebra.
So that's to finish the lecture. We can imagine that that's a gravitational wave.
And the timescale, the space scale is 1 billion light years from the left to the right.
Okay. See you on Friday.