Flash and JavaScript are required for this feature.
Download the video from iTunes U or the Internet Archive.
Description: Students learned to solve partial differential equations in this lecture.
Instructor: William Green
Session 22: Partial Differe...
The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high quality educational resources for free. To make a donation, or to view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.
PROFESSOR: So today we're going to talk about partial differential equations. Let me just give a quick review where we are. So we're going to have a PDE lecture today. Then Professor Swan is going to do a review on Friday. There's no lecture on Monday, but you have the quiz, too, Monday night. Wednesday, next week is Veterans Day so it's a holiday. There's no classes at MIT. So the next time I'll be lecturing to you will be a week from Friday. The homework following the quiz will be posted and so you guys can get started on that.
That homework involves a COMSOL problem which, for those of you who had trouble with COMSOL, might take some extra time, since it's the first time you're doing it. So I would suggest you at least give it a try early. I'm sure the TAs would be very happy to help you to force COMSOL to cooperate. Bend COMSOL to your will.
All right, so we're going to talk about partial differential equations. And, as it sounds, what's different about a partial differential equation is, it has partial derivatives in it. And so we have a lot of equations like this. So for example, I work a lot on this equation. It's a poor choice, isn't it?
OK, so this is the conservation equation for a chemical species in a reacting flow. So the change in the mass fraction of the species z, at a point, xyz at time, t. That differential is given by a diffusion term, a convection term, and a reaction term. So this is a pretty famous equation. This is part of the Navier-Stokes equations for the whole reactive flow system, which I hope you guys have seen already, and for sure you will see, if you haven't seen them yet.
And there's a lot of other partial differential equations. I work a lot, also, on the Schrodinger equation. That's another partial differential equation. And what's special about the partial differential equations is that, in this case, this partial derivative is respect to time, holding all the spatial coordinates fixed. And these partial derivatives are respect to space, holding time fixed.
A very important thing about partial derivatives, which you probably have encountered in 1040 already, is, it really depends what you hold fixed. And you get different results if you hold different things fixed. Now the implication, when you write an equation like this, is that whatever partials you see-- When I wrote this, if I have this and it has a term that's like d squared, dx squared, is one of the terms in this thing. This says t. This says x. The convention is, oh boy, I better hold x fixed. And here must be, I must hold t fixed. Because there's another one in the same equation.
Now, you can change coordinates however you want. Just really watch out that you carefully follow the rules about what you hold fixed, when you change things. Otherwise you can cause all kinds of craziness. Now you guys probably took multivariable calculus at some point in your distant past. And they told you very carefully what the rules were. So follow them. And I suppose, when you're doing thermo right now? Have you started to do all this partial derivative stuff? Yeah, so you've seen how completely confusing it can be with negative signs showing up all over the place. And other things like this. I don't know if you've encountered this yet? At least, I thought it was confusing when I took it. So just be aware. It's the same thing.
Now, often you do want to change the equations. So you can write down an equation like this. But, for example, if you have a special symmetry of the problem, like it's cylindrical, for example, then you might want to change the cylindrical coordinates. And then you have to just be careful that you write that correct Laplacian and cylindrical coordinates. Make sure it doesn't mess up anything about what you're holding fixed.
Sometimes you might want to be crafty and use some weirdo coordinates. For example, if you wanted to solve the Schrodinger equation for h2 plus, it turns out there's a special coordinate system called the elliptic coordinate system. And in that coordinate system, the partial differential equation is separable. And so, if you do it in that coordinate system, you can actually solve analytical solution, which is pretty unusual for the Schroedinger equation.
But it's a pretty goofball coordinate system and you really have to watch out, when you convert, to make sure you do it just right. But it's just the rules, chain rule and the rules about how to keep track of what's being held constant. And when you change, will tell them, how, being held constant, how to change it.
Now in principle, this problem is like no different than the ODE BBP problems we were just doing. You just have a few more terms. And so, you can use finite element. You can use finite differences. You can use basis set methods. And it's really no different. So all those methods basically look the same. You just get equations that have more unknowns in them. Because, for example, if you're using finite element methods or basis set methods with a local basis set, you need to put points both in the x direction and the y direction, and in this case, the z direction, and then maybe also in the time direction. So now you have a lot of mesh points in all different directions. And so you get a lot of mesh points. And so fundamentally, there's no problem, but in practice, there's a big problem if the number of unknowns gets to be so large. That's the main problem.
So think about it. You have some kind of equation like this. It's in three spatial dimensions and in time. How many points do you think you need to discretize your solution in the spatial coordinate? There might be somewhere, might be 100, I don't know, points might be enough. And so now I have x and y and z. So I'll need, well, just in COMSOL, you saw this on Monday. Suppose I just have two dimensions. Suppose I just have x and z. And I start putting points in. Well, I want to have some points along this way. And I want to have some points along this way. And so, I actually have an unknown here and an unknown here, and one there and one there, and one there and one there. There's quite a few of them.
I think of how many points I have. Each of those is a point where I have a little basis function. And each one of those has a magnitude of how much weight I have on that basis function, things they called the d's before in the basis functions. And I have a lot of them. Now, it's going to be, if I have 100 points this way and 100 points this way, I have 10,000 of these little points. That means I have 10,000 unknowns I have to solve for. 10,000 unknowns are a lot harder to solve for than 100 unknowns. and if you saw, on some of the problems we did with the ODE BBP problems, even with 100 unknowns, we were having a lot of trouble, sometimes, to get the real solution. All right? So you can just see how it just gets to be, like, a big numerical problem.
So that's two dimensions. The third dimension, you'll get another 100 times as many points. You might have a million unknowns. And once you get up to a million unknowns, we're starting to talk serious math, now. OK? Even for your good computers that you have. In addition, it's not just you have the number of mesh points, but it's the number of mesh points times the number of variables at each mesh point. So in the Navier-Stokes equations, what variables do you have?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Louder.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Velocity and pressure, OK. So you have pressure. We have the different components of velocity. What else we got?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Sorry, what?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yeah, though, yeah. Density or temperature. Yeah, temperature maybe, so. So maybe temperature would be what you might use. And then you'd have a state function maybe with these guys.
What else would you have? So then you'd have all these species, right,? The mesh fractions are all the species. So you have, like, y, species one. Mesh fraction, species two. However many you got. So that's a pretty big set of variables and you need to know all these numbers at each point. So this might be-- how many have we've got? One, two, three, four, five, plus however many species we have? So this is n species plus 5 times the mesh. And I told you the mesh was what, a million? So maybe I have 10 million unknowns.
So the issue, to a large extent, has to do with just, can I solve it? Can I solve a system of 10 million unknowns? 10 million equations, 10 million unknowns. All right, so that's going to be one big issue. Now the other issues, we still have the same problems we had in ordinary differential equations, that we have to worry about numerical stability. We have to worry about actually, is the whole problem even stable to begin with? Physically? Because if it's not stable physically, then probably we're not going to get a stable solution.
So you have all those kinds of problems that we had before in the ODE case, and now it's just sort of amplified by the fact, now we have boundary conditions in more than one dimension. So we have multiple chances to screw it up. And we're trying to figure out how to set the boundary conditions correctly. And we have to worry about stability in all the different directions. So those are the different kinds of things we have to worry about. But it's all the same problems you had before, just amplified.
Now, because of this problem of the very large number of unknowns, currently, if you have a problem that has two dimensions, you can kind of very routinely solve it. And things like COMSOL on your laptop, you can usually do those. When you get up to three dimensions, you might be able to solve it or might not be able to solve it. And it could be, like, a really big challenge sometimes to solve it. And then in reality, we have problems like the Navier-Stokes is four dimensions, right? There's time, as well. So then that gets to be really hard. And then the Schrodinger equation has three times the number of electrons dimensions in it. So that gets impossibly hard.
And so, specialized methods are developed to handle those systems with a lot of dimensionality. So, like, for the Schrodinger equation, people almost always do it with basis set methods. And there's been a huge effort, over decades, to figure out really clever basis sets that are really close to the real solutions, so that you don't need to use too many basis functions to get the good solution. So that's the way to keep the number of variables down.
In the Navier-Stokes world, I'll show you what people do to try to deal with that. There's many other problems. I don't know all the problems in the world. But in general, when you have three or more dimensions in your PDE system, you're looking up special tricks. You're looking up, how do people in this field deal with this particular PDE? Special ways that are good for that kind of PDE.
All right, so I said, one problem was that you could have, that the PDE system, itself, can be intrinsically unstable. So one example I know like that is, detonation problems. So if I have a system that has some explosive in it, or a mixture of hydrogen and oxygen even, and I have it in a tube or something, if I change my initial conditions, I can go from the case where it's stable, that it just sits there, a hydrogen and oxygen. Or it could turn into a regular burning. Or if I change the conditions a little bit differently, it could change into a detonation, where it actually sends a shockwave down the tube faster than the speed of sound, totally different situation.
So those kind of systems can be really sensitive to the physical initial conditions. A very small spark can make a really big difference in a system like that. And also, in your numerical solution, if the thing is physically sensitive, it'll typically be sensitive also to numerical noise. If you get numerical noise at some point, that might be enough to kick it over from, say, not igniting at all to having a detonation, which is a completely different solution.
But another case, like in biology, if you're trying to study the growth of a cell culture or the growth of a tumor, as you know, if somebody has a tumor and one of the cancer cells mutates and does something different, it can have a really different outcome of what happens to the patient, what happens to the tumor. The whole morphology can change completely. Everything about it can change. So that's the kind of system, also, that can be very sensitive to small details. A small change in one cell, for example, could be a really big difference. So this is not just in combustion problems you have these kinds of instabilities. All kinds of problems have this kind of thing. You can have just regular chemical reactors, like in multiple steady states, little tiny fluctuations can make it jump from one steady state to another one. That's another very famous kind of problem.
There's another class of PDEs where it's stable in some sense, but it causes a problem. We call those hyperbolic. And those are like wave equations. The solutions are propagating waves. So the equations of acoustics, the equations of electromagnetism, are wave equations where, basically, a wave propagates more or less without any damping. So what happens, then, is you have some numerical noise that will make a little wavelet that will propagate, more or less, without any damping, too. And it'll keep bouncing back and forth inside your domain numerically. And if you keep introducing numerical noise, eventually the whole thing is just full of noise. And it doesn't have much relation to the real physical solution.
So hyperbolic systems of differential equations need special solution methods. That in the numerical method, it's carefully trying to damp out that kind of propagating wave that's coming from the noise. And they set them up a special way. I'm not going to talk about how you solve them in this class, but there's a special, whole group of solution methods that people use who are solving acoustics equations and solving electromagnetic equations. If you get into solving shockwave equations-- I actually work for Professor [? Besant. ?] He has, like, electrochemical shocks in some of his systems. You have the same kind of thing. You have to use special mathematics to correctly handle that, special numerical tricks.
But we're not going to get into that in this class. If you get into those things, you should take the PDE class for hyperbolic equations and they'll tell you just all about the solvers for those kinds of things and special tricks. And there's journal papers all about it.
Here, what we're going to focus primarily on two kinds of problems: elliptic problems and parabolic problems. And those problems are ones where there's some dissipation that's very significant. And so, in an elliptic problem, you introduce some numerical noise and as you solve it, the numerical noise kind of goes away. Sort of like intrinsically stable all throughout. Now the real definition of these things, parabolic, hyperbolic, elliptic, has to do with flow of information. So in a hyperbolic system, information-- there's regions of the domain that you can make a change here and it doesn't make any effect on some domain part over there. And so that causes-- because that is the real situation, that also propagates into the numerical solution method.
So for example, if I'm modeling a shockwave that's moving faster than the speed of sound, I can introduce any pressure fluctuation I want behind the shockwave, and it has no effect on the shockwave. Because the pressure waves from that pressure fluctuation are traveling at the speed of sound. They won't catch up to the shock. So they will have no effect whatsoever. So if I was judging what was happening by looking at what's happening in the shock, I can introduce any kind of random noise I want back over, downstream of the shock, I guess. Yeah? And it won't make any difference at all. And so, that's part of the reason why I need a special kind of solver.
If I have elliptic equations, every point affects every other point. A famous case for that is, like, this steady state heat equation. You have some temperature source here, a heat source here, and you have some coolness source over here. And there's some heat flowing from the hot to the cold. And the heat flow is kind of slow enough that everything sort of affects everything else, since you actually have, the heat is trying to flow by, sort of, every path. And there's some insulation that's resisting its flow. And it has been in the connectivity of flow a certain way. And those kinds of equations are called elliptic. And the methods that we use for the relaxation mesh that we showed are really good for those kinds of methods, for those kind of problems.
And then another kind of problem that we get into a lot is like this one. This is called parabolic. And typically, in the cases we see that are parabolic, we have time as part of the problem. Always what we think should happen is that what happens in the future depends on what happened in the past, but not vice versa. So you would expect that when you're computing what's happening at the next time point, it's going to depend on what happened at earlier time points. But if you try to do it the other way around, it would be kind of weird, right? You wouldn't expect that what in the future really affected stuff in the past.
So it naturally leads us to sort of pose it as an ODE IVP, or IVP kind of problem. And we know initially, some time zero something, and we're trying to calculate what happens at some later time. And so those kinds of problems we call parabolic. And the methods that we're going to talk about mostly are ones for solving parabolic and elliptic problems. And you can read, in the middle of chapter 6, has a nice discussion. Just has the mathematical definitions of what these things are.
Now, even in a steady state problem, the problem can have a sort of directionality to it. So if you have a problem that's very convective. So you have a flow in a pipe. And you have a nice, fast velocity flow down the pipe. What's happening downwind is very much affected by what happened upwind, but not vice versa. So maybe you're doing some reaction upwind. Say you're burning some hydrogen. You're going to get water vapor made upwind and it's going to flow downwind. And you're going to have steam downwind.
If I squirt a little methane in downwind, it doesn't actually do anything to what happened upwind. I still have a hydrogen and oxygen flame up there and it's been burning and making steam. Do you see what I mean? So it has a directionality, even though I might do it as a steady state problem and not have time in it at all. And you've probably already done this, where you can convert, like, a flow reactor. You can write it, sort of, as a time thing or as a space thing. Time and space are so related by the velocity of the flow. We'll have similar kinds of phenomenon then. It's almost like a parabolic time system that, if the flow is fast enough, the fusion backup, anything moving upstream is negligible. It's all moving from upstream to downstream.
And so then, you'll have the same kind of issues. In fact, you may even try to solve it by just starting it at the upstream end and competing stuff and then propagating what that output from that first one does to the next one downstream, and then do them, one at a time. That would be one possible way to try to solve those kinds of problems that are highly convective. As you slow the velocity down, then the diffusion will start to fight the velocity more and more. If you make the velocity very slow, the Peclet number very low, then you really have to worry about things diffusing back upstream. And then it's a different situation. And that would be more like the elliptic problems we have. Because things downstream would actually affect upstream, if the flow is really small.
But typically for our problems, the velocity is always really high. Peclet numbers are really large and so not much goes from downstream to upstream. It's almost all upstream to downstream.
This leads to a famous situation. If you look in the textbook at figures 6.7, 6.8, it's kind of a very famous problem. That if you just do this problem. So this is just an ODE problem. So I just have convection and diffusion, no reaction, very simple. This problem. You'd think that I could just use my finite differences, where I would put, I could use the center difference for here and here and change this into this obvious looking thing.
You think that that would be a reasonable thing to do, right? And you can do the same one here. You can write-- I'm terrible with minuses, sorry. So these, the finite difference formulas. And so, you can write these like this. This stays equal to zero. And now it's just an algebraic problem to solve. A system of equations, you have equations like this for every value of n. And there's the mesh points. You guys have done this before, yes? Yes.
OK, so famous problem is, if you actually try to solve it this way, unless you make delta z really tiny, you get crazy stuff. So on figure 6.7, they actually show what happens for different values of delta z, which is different values of what they call the local Peclet number. So there's a thing called local Peclet number, which is-- OK? If you make the local Peclet number too large, actually anywhere bigger than 2, and you try to solve this equation, what you get is oscillations, unphysical oscillations. They'll make the phi go negative. It's crazy. It looks like the simplest equation in the world, right? It's a linear equation, so what's the problem with this? But it doesn't work.
And so, then you could think, well, why doesn't it work? And there's, like, multiple ways to explain this about why this doesn't work. I think that the best way to think about it is about the information flow. So if the local Peclet number is large, that means that phi is just flowing downwind and diffusion is not doing much. Because the Peclet number is big. And so, you really shouldn't use this formula to calculate the flow. Because what happens downwind, at point n plus one, I guess. Here we have the flow this way. So here's n minus one. Here's my point n, n plus 1.
I really should not include, is this right, n plus 1 in this formula. Hope I have those signs right. I apologize if I had them backwards. Because I don't think that what happens downstream should really affect this at all. So really, I would do better to change from this center difference formula to what's called the upwind difference formula. Where they would say, OK, instead, let's use phi of n minus, the phi of n minus 1 over delta z. Just use that instead of using this formula.
Now, you wouldn't think that that would make much difference. But it makes a gigantic difference. And so if you look at figure 6.8 in the textbook, you see you get perfectly stable solutions when you do that. But where you do this one, you get crazy, oscillatory, unphysical solutions, unless you choose delta z really tiny.
Now, there's other ways to look at this. So I think this is the correct way, the best way to think about it is, it has to do with information flow. If you try to make your solution depend on stuff that it doesn't really depend on. So upwind does not really depend on downwind. If you force it to, by choosing to include that information here, you end up with crazy stuff. So that's one way to look at it.
Another way to look at it is that, if you really want to compute these derivatives, you've got to keep delta z small. Because we're taking a differential. We're turning it into a finite difference. And if you take your delta z too big, then it's a terrible approximation to do this. And you can look at it. You can see that what we're doing is, when we have this problem where the local Peclet number is getting too large, is we're actually choosing delta z bigger than the, sort of, the thickness of the solution. So if you look at the solution, the analytical solution of this problem, it's like this. Where this has a very, very thin layer at the end. It's like the flow has pushed all your stuff to one side.
And if you choose your delta z to be from here to here, that's how big your delta z is, and then you're computing the derivative of this point halfway along from here to where this is. You compute your derivative like that. It's not a very good, not a very accurate representation of what the real derivative is here. And in fact, you really should have a bunch of points here and here and here and here and here to really resolve the sharp edge. And so you did something really crazy to use a big value of delta z to start with. So that's another way to look at this problem.
Now, the fix, by doing this upwind differencing, I think the best way to look at this is saying, well, I'm just going to make sure I only depend, I only make the equations depend on what's upwind, with the convective term. Because there's nothing convecting from downwind. So just leave that out. So that's a conceptual way to think about it.
Another way to think about it, which is in the text, is that by doing this, you're introducing another thing called numerical diffusion, where you're actually, effectively increasing the diffusivity. Because this is a very poor approximation to the differential. This is a asymmetrical finite difference formula. It's not a very good value estimate of the derivative, if delta z is big. But you can write it out carefully and write this out and say, oh, this actually is sort of like this plus an effective diffusivity chosen very carefully, so the terms cancel out just right. And it turns out, if you look at it that way, the effective diffusion you've added, by using this instead of that, is just enough to make the local Peclet number, if you compute it with the effects of this d effective, to stay less than 2. Which is what you need as a condition to make this numerically stable.
So I strongly suggest you read that. It's only two pages long in the textbook. It's very clear. Just read, just look at it if you haven't seen this before. There's a similar thing and it might be relevant to the COMSOL problem that Kristyn showed you on Monday, is you can-- In that case, there was a region of the boundary here that had some concentration c naught, and then there was, over here, the concentration was zero. The boundary is zero. Here, it's some number, 17, whatever number it was.
And if you think of how to model what's going on, one way to look at it is, I have a diffusive flux coming in from the drug patch, diffusing the drug into the flow. And I have this big, giant flow, flowing along here. And the flow gets slower as I get closer to the wall because of the friction with the wall. And so if I looked somewhere right in here, I could just draw a little control volume, and look what's happening. I have a diffusive flux coming in here. I don't have any drug up here, so there's no drug coming in here. But I have drug leaving. So basically, it's coming up and it's making a right-hand turn.
And so then I could think, well, if I was going to try to figure out, what's the steady state concentration of drug in there, I could say, well, there's a certain amount of drug entering from the diffusive flux, which would be just this times dcdy. This is the y direction. And then this part over here would be the velocity, actually, just times the concentration I have in here, right? So how much is flowing out. And it's probably, my units are screwed up so there's-- Well, maybe not. That's OK. Is that all right? Units are screwed up. You guys will get it right. It's probably an area. Something to do with this size right here. This one here, too. Is that right?
So you could compute what the steady state concentration would be if it's just coming in and flowing out. And so, this is the finite volume view of this kind of problem. And so you can write the equations that way, instead. And what you really care about are what the velocity flows and the diffusive fluxes are on these boundaries around the point. So you don't try to compute things right at the point. You compute around it.
And one way to look at this is, we just did finite differencing. We were looking at what was happening coming in. And then, this is actually a center difference around this interface, which is halfway between the two points. And so actually, this is a good approximation for the finite volume point of view. This is all right?
Now, these are all different ways to write down the equations. All of them are correct in the limit that delta z goes to zero. If you have an infinite number of mesh points, all these are exactly the same. But they lead to really different numerical properties when delta z is large. And they're all inaccurate when delta z is large, so it's all about what happens with the inaccuracies.
Now, we can go ahead and take any problem, even really complicated 3-D problems or 4-D problems like this, and we can write them with relaxation methods. And what we're basically doing in a problem like this is turning it into some function of y is equal to zero, where y is the value of all the unknowns. And there's a lot of equations. So I might have 100 million equations and 100 million unknowns that are the values, or every state variable at every state point at all times. And I can write it down, by finding differences, by finding elements, by finding volumes. Any way I want, I can write down an expression like this. And I know in the limit, as I make the number of elements of y go to infinity, it will actually be the real solution.
And the problem is, my computer can't handle infinite number of unknowns. In fact, it's even going to have trouble with a 100 million unknowns. And so, I have to figure out what to do. So how would I solve this normally, is I would take the Jacobian of this and I would do, I'd take an initial guess and I'd update it and have the Jacobian times my change from the initial guess to the equal negative of f, right? You guys have done this a million times.
As Newton-Raphson is how you'd find improved from an initial guess. You have to provide some initial guess. This is actually pretty hard if you have 100 million unknowns. You have to provide the value of the initial guess for every one of the 100 million. But somehow you did it. And now you want to refine that guess and this is the formula for it. And so, you just use backslash, right? But the problem here is that now f is a 100 million long vector and this is a 100 million squared matrix. And a 100 million squared is pretty big, 10 to 16th elements in it. So I'm not going to just write that matrix down directly like this.
Now, we cleverly chose local basis functions, so this is going to be sparse, this Jacobian. So most of the numbers, are those 10 to the 16th numbers in this matrix are zero. So we don't have to represent them. But there still might be quite a few. Probably the diagonal events are non-zero, so that's like 10 to the 8th of them right there. And we've got, probably, a couple of other bands. So I don't know, might be 10 to the 9th non-zero elements inside this thing, which is pretty nutty. And depending on how much memory they provided in the laptop they gave you, you might be filling up the memory already just to write down j.
And certainly, if you tried to do Gaussian elimination to solve this, you're going to have a problem because you get what's called fill in. Do you guys remember this? Even if you start from a very sparse Jacobian, as you run the Gaussian elimination steps the number of non-zero entries in the intermediate matrices gets larger and larger. Remember this? And so, even if, initially, you can write down j and store it in your computer, you could easily overwhelm it with the third iteration or something.
So you're not going to just use Gaussian elimination. So what do you have to do? You want to find a method. They're called direct. Direct methods to solve this kind of problem, in order to get your increment, which is going to be your update, to get a better estimate of what the solution is. And the direct methods are ones that you don't actually store the gigantic matrix. So you're trading off CPU time for memory. So you don't have enough memory to store the whole thing. Which we, normally you would do this with Gaussian elimination and huge steps. It would solve it perfectly, right? Right, Gaussian elimination's great, backslash. You guys have used it a few times. It's a really nice program. It always works. It's really nice.
But you're going to give up that because you can't afford it because the memory, it's consuming too much memory and your cheapo department didn't buy enough RAM in your computer for you. So you're going to instead try a direct method which is, so it's trading off CPU time versus RAM. So you're going to reduce how much memory you actually actively have, but you're going to expect that it's going to cost more CPU time. Because otherwise, everybody would do something else besides Gaussian elimination all the time, and they don't.
All right, so what's the method we're going to use? It turns out that variance on the conjugate gradient method turned out really well. So we mentioned this earlier, briefly. What you do is, instead of trying to solve this, you try to solve this.
So you just try to minimize that, right? You know at the solution this is going to be zero. Right, jy plus [INAUDIBLE] zero when this is solved. So you try and vary and find the delta y that makes this go to zero by minimizing it. And then this method, it turns out that the conjugate gradient, if everything works great, this is guaranteed in n iterations to go directly to the solution of this kind of problem. Now n is large, now, because we have 10 to the 8th unknowns. So that's 10 to the 8th iterations. You do anything in 10 to the 8th iterations, you're not really sure it's really going to work. So that's a warning. But at least it should get, even after a few iterations, you should get a smaller value of this than you had before.
And so this is the method. And what's really nice about this method, if you look into the details of how it works, it never has to store any intermediate matrices. So all it needs is the capability to multiply your matrix times some vector. And from that multiplication, it figures out a step direction. That's the trick of this method. And because it only has to do a forward multiplication, you don't have to actually store j. You can just compute elements of j, as needed, in order to evaluate this top product. This j times v. And you don't have to ever store the whole j at once. So it's really great for RAM to do this method.
Now, as I said, when you do 10 to the 8th iterations, things don't always work out so well. And so people have found that when you go more than maybe 10 or 20 iterations like this, you tend to pick up numerical problems. And so, they've worked out better methods. And there's one that people use a lot for these problems called bi cg stab. Which is biconjugate gradient stabilized. The applied mathematicians went crazy trying to figure out better methods like this. And this is the one that, I think currently, people think is the best. Though probably, there's a lot of research on this, so maybe there's even better ones now. But inside Matlab, I think this is the best on they have.
And this is a method for doing this. Now, it's iterative. Now, let's remember what we're doing. We're trying to solve this problem. We started with a guess. We're going to break it up into an iterative problem like this, where we're going to do some steps, delta y. This is now doing iterative procedure in order to compute delta y. So we have an interim procedure and another iterative procedure inside the first iterative procedure. So this is going to be a lot more CPU time. So just to warn you. And it's also, it's not guaranteed like Gaussian elimination, to beautifully come to machine precision solution. It's going to get you something. But anyway, this is what people do. So that's one good thing to try.
Now, this whole approach is sort of based on Newton-Raphson. We know that doesn't have the greatest radius of convergence. So you need to provide a pretty darn good initial guess. So how are we going to get a good initial guess if you have a 100 million unknowns? So what methods do we know? So one idea is you could do things like homotopy, stuff like that. If you could somehow take your PDE system, get rid of the nonlinear terms. You're really good at solving linear systems of PDEs. You start with that and then you gradually turn the non-linear term on. Maybe you could coax it over. So that's one possibility, if you're really good at PDEs.
Another possibility that I've run into a lot is the problem you're trying to solve, for example, might be a steady state problem like this. Where this is zero and you try to figure out what the steady state situation is in a flow reactor, for example. And so, it's actually, your steady state problem came from a time dependent problem. And if you think that your time dependent problem really converges to the steady state problem, sort of without much concern about what the initial guess is, then you could put any random initial guess in and march it along in time enough, and eventually, it should end up at the solution you want.
So that's the idea. So that's called a time marching idea. So there's one idea, sort of the homotopy continuation kind of idea. That's one idea about what to do for initial guesses. Another idea is the time marching idea. And the key about it is, you have to really believe that, physically, the system does march to the steady state you want. If it does, then this is a good idea. If it doesn't, then you're not going to end up where you want, because it's not going where you want. And so, you have to watch out.
But in some situations, this is true, that things work. So I do a lot of combustion problems. A lot of these things, you light a spark on a stove burner. No matter how you light it, it ends up basically the same flame. So that's a good one. But I have some other ones where I light a candle and I'm in a convective field, a turbulent field, and it flickers. And no matter what, it always flickers. And I never get a steady state, so I'm never going to find a steady state solution with that one. And so I could time march until I retire. My computer's still burning CPU time all that time and it'll never get to a steady state solution. So you really have to know what the real situation is. But if you have one where it's going to converge, then this is a reasonable idea.
Now let's talk about the time marching. The time marching, you have, say, dy, dt is equal to some f of y where this is the discretized spatial vergence. I took all the space dimensions and replaced them with the finite difference expressions for the derivatives. Or I did colocation or something to convert this into algebraic equations. So the right-hand side is now algebraic equations, no longer differentials. And now, I'm going to time march this and I just-- This is very long. It's a 100 million.
So how are we going to time march that? Well there's kind of two schools of thought for this. One school of thought is, if I can use an explicit method, an explicit ODE solver, those are pretty nice. Because all I have to do is evaluate f and then I can march along and I never have to save any matrices. Because remember, look at those formulas, y, y new, by the explicit method, is equal to some function g of y old.
And this is the update formula. It can be forward Euler. It could be RK 2. It could be whatever. But there's just some explicit formula. I plug in the old vector. I compute a new vector. I never had to do any inversions of matrices or anything. So that's good. And this is, indeed, the way that the Navier-Stokes equations are solved currently by the best possible solver. So the best solver in the world discretizes in space and then does ODE IVP using explicit method for marching it.
They do problems that are pretty big. So they'll have like 100 million mesh points and then 100 species. So they have 10 to the 10th unknowns. And at each time step, you're computing another vector that's 10 to the 10th long. So it's a pretty big problem. And they use, like, 40% of one of the national supercomputer centers will be running at one time for one job like this. But it works. You never have to store anything that's too big. The biggest thing you have to store are the vectors. Which, you're always going to have to store some object the size of your solution, right, if you're going to get your solution. So that's one possibility.
Now, this is really good if an explicit solver can solve it. So that means it's just not stiff, your problem. And also, if you want a time accurate solution. So in that case, they actually want the time accurate solution. But suppose we're trying to really solve the steady state problem over here. Then we really don't care about the time accuracy. We're not going to report anything about our time steps. In the end, we're only going to report our steady state solution. That's all we cared about. So in that case, there's no reason to try to have a really high accuracy, explicit formula and try to make sure we all our time steps are exactly right. We're just going to throw away all those time points anyway. We'll just keep the final one as our initial guess for a Newton-Raphson solve to find the steady state.
So in those cases, you might instead want to use-- So this is time accurate. There's another method. I don't know if I should call it time inaccurate. If you don't care, if the solution is really, what the solution yft is, all you care about is where you get to, then you can do other methods. And one of them is the backward Euler. So you can say, well, y new is equal to y old plus delta t times f of y new. Now, the advantage of this is, I can make delta t pretty large. This is guaranteed stable as long as the ODE system is stable. Remember? And so, this will go to this, go to the solution pretty well, to the steady state solution.
It won't go there in the right amount of time because this is a very low order of formula. It's not really accurate. But it will end up to the real steady state solution. This is what's actually used a lot. However, the problem with this is, this is now an implicit equation. So we may have to solve it with some method like Newton-Raphson. But we got into this because we couldn't solve our Newton-Raphson steps because we didn't have a good initial guess. So then the question is, can we find a good initial guess for this problem? And in this case, you can because you can use the explicit formulas to again provide the initial guess for this implicit formula. So that's what they do, also.
And so now, if we're doing an iterative procedure in order to get the initial guess for our other iterative procedure ever there, which we're going to break up into another iterative procedure inside it to solve every step. But this is currently what people do. So you guys are smart. Maybe you can figure a better way to do it. The world is waiting. I think I'll stop there.
By the way, actually just names-- This kind of thing is called method of lines. And if so, if anybody ever says I solved something by method of lines, what they meant was, they discretized in some of the coordinates and they kept one of the coordinates as a differential. And then they solved it with an ODE solver at the end. And in certain systems like this, it's a smart thing to do. You need to be kind of lucky that the set up is right, so you can use it. But if you can, it can be pretty good. All right. See you on Friday.