Flash and JavaScript are required for this feature.
Download the video from Internet Archive.
Description: This video begins with questions from the students. These are followed by a discussion of stiffness, explicit and implicit schemes using an in-class example. Finally, the video introduces Runge-Kutta methods.
Instructor: Qiqi Wang
The recording quality of this video is the best available from the source.
Session 6: Numerical Integr...
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: OK. So let's start. So first of all, it looks like there are many questions regarding the reading and homework. So if possible, I like to answer some of your more [INAUDIBLE] questions. Anyone want to raise a question I can answer? Yes.
AUDIENCE: So something that we were running into was the definition of local ordered [INAUDIBLE]. Local and stationary and global [INAUDIBLE] and how to figure out [INAUDIBLE]. OK. So what we ended up doing was finding [INAUDIBLE] by taking dd plus 1 minus [INAUDIBLE].
PROFESSOR: Yeah.
AUDIENCE: And then saying that that was the [INAUDIBLE] location of it.
PROFESSOR: OK. Yeah.
AUDIENCE: And then that was equal to-- whatever order [INAUDIBLE] that was was equal to the order of delta t times [INAUDIBLE] plus 1 where t is the local order [INAUDIBLE].
PROFESSOR: That's exactly what you should be doing is finding the local order of axis. Did anybody else understand what [INAUDIBLE]? OK. [INAUDIBLE] saying. So that's exactly the right way of finding out the local order of axis.
Looking at the schema as [INAUDIBLE] n plus 1 minus the [INAUDIBLE] of the [INAUDIBLE]. That is you're tau. That is your [INAUDIBLE] area. And tau should be on the order delta t to the t plus 1 where p is a local order of [INAUDIBLE]. OK. That's basically repeating what [INAUDIBLE].
AUDIENCE: And how does that relate to the global order?
PROFESSOR: How does that relate to the global order of [INAUDIBLE]? Anybody want to attempt to answer that?
AUDIENCE: [INAUDIBLE]
PROFESSOR: So say that there is an if.
AUDIENCE: [INAUDIBLE]
PROFESSOR: If [INAUDIBLE] under what conditions does [INAUDIBLE] theorem hold?
AUDIENCE: [INAUDIBLE]
PROFESSOR: 0 stable, right? As long as t is great or equal to 1, which usually it is, you have a consistency. So consistency is usually is not that much of problem when you already have the local order of [INAUDIBLE]. 0 stable is something you have to separately. But as long as the schema is 0 stable, then the local order of accuracy is equal to the global order of accuracy as long as that order of accuracy is great of equal to 1. Is that clear?
AUDIENCE: Could you repeat that one more time?
PROFESSOR: If a [INAUDIBLE] is 0 stable, then the local order of accuracy is the same as the global order of accuracy.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Exactly. When a scheme is not 0 stable, it's [INAUDIBLE] equal to 0. And there is no global order of accuracy at all. Yes.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yes. If a schema is not 0 stable, then the scheme doesn't work. Remember when we had the supposedly most accurate two-step scheme, [INAUDIBLE] scheme? What happens when we rebuilt the delta t? The error, instead of decreasing, it increases.
So that's what happens with a not 0 stable scheme, where you refine your time step. Instead of the error decreasing, you have an increasing error. All right. Any other questions?
AUDIENCE: [INAUDIBLE]
PROFESSOR: OK. So the question is can we go over eigenvalues [INAUDIBLE] again? [INAUDIBLE] page-- yes. All right. So in doing eigenvalues, stability you plug in du dt equal to lambda u, right? So that is the question you use to analyze eigenvalue stability.
Then you write down an equation of Un plus 1. Again, if it's a midpoint rule, then the scheme I'm plugging is plus 2 delta t times f of Vn. And in this case, f of Vn is lambda times the Vn right? So that is the equation you get by plugging in the differential equation into a scheme.
Now, the next thing to do is to assume my Vn. And actually V n plus 1, Vn minus 1, and all other terms in your scheme satisfies an exponentially growing pattern. So when I whatever time stamp it is, it is equal to the 0 time stamp times a growth factor to whatever time stamp [INAUDIBLE].
So once you plug all of these into the scheme, you can find the quadratic equation for z if it's a two-step scheme. If it's a one-step scheme, you find just a linear equation for z. If it's a three-step scheme, you're going to find out with a cubic equation for z, et cetera.
And then the next is to see, does that equation allow any z, any solution, who's modulus is greater than 1? If that equation for particular combination of delta t and lambda-- delta t and lambda is always multiplied together. So for whatever delta t times lambda if this equation allows for z that has a modulus greater than 1, then it is not eigenvalue stable. The combination of the scheme and that lambda delta t is not eigenvalue stable. Yeah.
AUDIENCE: When you say the modulus greater than or equal to 1, I guess whenever I see modulus [INAUDIBLE].
PROFESSOR: OK. So the question is modulus is. Here the modulus is off a [INAUDIBLE] number.
AUDIENCE: Oh, so do you mean like [INAUDIBLE]?
PROFESSOR: Yeah. The [INAUDIBLE] to it is the real of z squared plus imaginary of d squared. What is the significance of that? The significance of that is that if that modulus is greater than 1, then z to the nth power is going to grow larger and larger as you make n larger. If the modulus is less than 1, then z to the nth power is going to become smaller as a take n to infinity.
What happens if the modulus of these equal to 1? Then the modulus of these always going to be equal to 1 whatever [INAUDIBLE]. OK. Any other question? Yeah.
AUDIENCE: [INAUDIBLE]
PROFESSOR: OK. So last week's [INAUDIBLE], why does du dt [INAUDIBLE] negative lambda times u. So in analyzing eigenvalue stability, you are looking at equation of this form. Right? So if I write down equal to lambda minus lambda times u, it is still in this form.
It's just the sign of the lambdas are kind of wrong. I mean, you can-- maybe I should use another symbol. Maybe I should say dudt equal to minus theta times u. In that case, then you just have to [INAUDIBLE] lambda equal to minus theta. And you can do the analysis the same way. Does that answer your question?
AUDIENCE: Yeah. [INAUDIBLE]
PROFESSOR: OK. Any other question? No? [INAUDIBLE] more. OK. OK. So for the reading recap, we have read about stiffness and the Newton-Raphson. Let me recap what they are. Stiffness are really orders of magnitude distance in timescales. OK.
So a stiff problem means there is one timescale that is a very, very small compared to the timescale of the main phenomenon we want to simulate. Maybe we want to simulate some phenomenon that is interesting. But in order to simulate that, there is something of a hidden timescale that is much, much faster than the one we are interested in.
And I'm going to show you two examples of what that very small timescale is. And so that is expressed in terms of physical phenomenon, two different timescales. we're presenting maybe two different physical phenomenon that you have to simulate at the same time. In terms of mathematics, where you write down OD, when you write an Ordinary Differential equation, you already kind of abstracted from the physics. So in that case, how do you interpret [INAUDIBLE]?
OK. We can say that when you transform a physical phenomenon into an ordinary differential equation, the timescale is really translated into eigenvalues. And we are also going to see examples. So if you have two timescales, one slow, one fast, we are going to see two eigenvalues where you write down the equation in its matrix form and [INAUDIBLE] do eigenvalue analysis of the matrix.
So we're going to see one large eigenvalue, one small eigenvalue. And the fast process responds to the larger one or the smaller one, what do you think?
AUDIENCE: Large?
PROFESSOR: The large, why?
AUDIENCE: Because it's [INAUDIBLE].
PROFESSOR: Exactly. So the large eigenvalue corresponds to the fast process. So if you think of dudt equal to lambda u, what unit does lambda have? That's another way to think about it. What unit does lambda have? In engineering, one of the most important things is to look at the units of everything, right? Huh?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yes. The unit of lambda is 1 over the unit [INAUDIBLE]. It's 1 over [INAUDIBLE]. Lambda really has a unit of frequency or rate, right? Actually, lambda is a-- imaginary lambda it is exactly the [INAUDIBLE] frequency of the solution, right? So lambda has the unit of frequency.
A very fast process is going to have a high frequency, therefore a large eigenvalue. And a very slow process is going to have a low frequency, and therefore a small eigenvalue. All right, does it make sense? So stiffness can either be interpreted in terms of order of magnitude difference in timescales or order of magnitude difference in eigenvalues.
OK. You also read about Newton-Raphson. And Newton-Raphson is the matter of solving nonlinear equations using increasing schemes. So the [INAUDIBLE] is, as you're going to see, is going to benefit a lot from using an implicit scheme instead of explicit scheme.
And when you have a nonlinear equation and you want to use an explicit scheme, as you dismay discover when you are doing the homework question, at every time stamp you need to solve a nonlinear equation, a nonlinear algebraic equation. In the homework, the equation happens to be a quadratic equation, and you can solve it using [INAUDIBLE]. But in general, that nonlinear equation you have solved every single time stamp is not going to happen to be an equation with an analytical solution. And a Newton-Raphson is a method people have invented centuries ago to solve such nonlinear equations numerically. All right, any questions before we dive into stiffness and using Raphson?
All right. Stiffness. The first [INAUDIBLE] when we illustrated stiffness is this guy. So we did in the first lecture, [INAUDIBLE]. Right? And this is a very, very similar set up. And to remind you, we have the same [INAUDIBLE] over here.
But instead of an aluminum cube, we have two metal plates sticked together like [INAUDIBLE]. So this blue [INAUDIBLE] is very [INAUDIBLE]. And one is copper, the other is aluminum. They are 2 inch by 2 inch big. And each only is a [INAUDIBLE], so one aluminum, one copper.
And I'm going to heat only this aluminum with a heat lamp. So I'm going to turn on the heat lamp. And there are two thermal [INAUDIBLE] that matches the [INAUDIBLE]. What differential equation would you use to model this problem?
When I'm heating only the aluminum, how would the [INAUDIBLE] of this and this going to behave? How would that be [INAUDIBLE]? You can discuss in pairs or in groups of three. And we'll see really answer the question of if we want to predict it, is it a stiff problem? If so, why? If not, why not? You can use a computer to look up any concepts or anything like that.
You can see the two colors, one is bronze-- sorry, sorry, one is copper. One is a copper temperature. One is a aluminum temperature. So we have some very good conclusions here. Why don't you just help us understand [INAUDIBLE]?
AUDIENCE: So this is a stuff problem, because there are two different things that are happening. We have the convection, which we [INAUDIBLE] the heat [INAUDIBLE] and the metal and then the metal in the air, and then the conduction between the two metals. The metals are both very conductive. So that process happens much faster than the rate of convection.
PROFESSOR: All right. Does that make sense? All right, thank you. So we have to physical phenomenon, convection and conduction. And if you'd think carefully, the two metals I exposed over here, within like the common metals, are the most heat conductive ones that you can reasonably [INAUDIBLE], aluminum and copper.
So I deliberately made the rate of conduction to be very fast. I deliberately made the conduction to be very fast. And I also made them very thick. So if you look at-- let me start writing some equations.
OK. So if you look at the conduction, rate of conduction, the rate of conduction, heat conduction, is equal to what? Can somebody tell me? There is definitely a k somewhere, the conductivity times the area times delta t over delta-- in this [INAUDIBLE]. Let's say the direction across the plates are x. So it is delta t over delta x.
OK. And changing the Q-dot is also equal to-- actually, I shouldn't say this. I should say the mass of each plate times the thermal capacitance, mc, times the partial T, partial T. So this is the rate of temperature change is also equal to Q-dot convection minus this is the conduction, Q conduction. So this is for the aluminum plate.
And if you just look at if, let's say, we isolate these two phenomenons if we only look at the conduction, what we get is partial T partial T equal to kA divided by mc delta T over delta x. So remember our eigenvalue analysis? This part would become our eigenvalue.
If we already have the conduction process, then this would be our lambda. OK. I think I get a minus sign here. So this will be our lambda. So we have a large k. OK We have our large A. But have a large [INAUDIBLE] compared to-- we have a large A compared to our delta x, two [INAUDIBLE].
So we have deliberately made these to be large. And if you plug in the numbers, this is significantly larger, about one to two orders of magnitude larger than the forced convection rate [INAUDIBLE] table. So this is a stiff process in terms of physics, because we have a very fast speed conduction and a relatively slow speed convection.
In terms of mathematics, let me say this is the aluminum. This is delta squared delta T is equal to the T aluminum minus T copper. So we have m aluminum c of aluminum equal to this. And we have aluminum equal to, let me say-- and let's set this k to be the average rate of conductivity in the aluminum and the copper.
So we replaced our delta T by TA minus Tc. All right. So this is conduction. And we also have convection. And the rate of convection, can somebody tell me what the rate of convection is?
We have the convection coefficient times A, area-- sorry-- times the difference between the T error. So let me just say to [INAUDIBLE] T error minus T aluminum and divided by aluminum mass and the heat capacitance. And our copper is going to satisfy a similar equation, where the heat conduction is the same. So what this is c instead. And we have TA minus TC here divided by delta x.
And here we have a removal of heat. So this is hot air. And we have a removal of heat. Cold air divided by Mc Cc. OK.
So as you're going to see, we are going to get a 2 by 2 matrix. We are going to get a 2 by 2 matrix. And the 2 by 2 matrix is going to have the 2 by 2 matrix. We write down the T of the aluminum temperature and copper temperature is going to be equal to if you regroup the terms.
And we know that this term is going to be large. So let me call this big A. And these terms are going to be small. Let me color these small terms by blue.
OK. So this term is going to be small. This term is going to be small. And let me color the last terms in red. So this term is going to be large.
So let me just write qualitatively here. We have a minus large term over here and a minus small term minus a small term in blue. So this is a multiplied by Ta and the Tc.
So on top of Ta, we have a negative large term and negative small term. On top of Tc here, we have a plus large term. And in the second equation, the evolution of the copper temperature, we have a large term in front of here and nothing in the small term. And in terms of the copper temperature, we have a minus large term and a minus small term.
So the matrix looks like this. That's how the matrix would look like when you write these two equations in terms of matrix form. Now, if you do the eigenvalue analysis of a matrix that looks like this, you can take the-- what happens to the eigenvalue if S infinitesimally small compared to L? Let's say what is the eigenvalue of minus L, minus L, L and L?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yes. Or it is L, because eigenvalues are linear. It's an L times the eigenvalue of this matrix. And what is the eigenvalue of that matrix? We can try to--
AUDIENCE: Well, at least one of them [INAUDIBLE].
PROFESSOR: OK. At least one of them are 0. Because this matrix is [INAUDIBLE], right? The first line is negative of the second line. What about the other one?
AUDIENCE: 0 minus 2.
PROFESSOR: 0 minus 2 exactly. Wait, is it minus 2? Yes, it is minus 2. OK. So it is 0 or minus 2, right? Minus 2L. So you can see the eigenvalues, one of them is 0. The other is minus 2.
And if you add these small terms the eigenvalue-- because the small items are just minus S on the diagonal, right? So how do they affect the eigenvalues? It's a test to your linear algebra.
AUDIENCE: [INAUDIBLE]
PROFESSOR: It is going to be put a minus S on both of the eigenvalues. It is the eigenvalue of an identity matrix times S times minus S. The eigenvalues of identity matrix is 1, right? So let me see-- [INAUDIBLE]. I'm not sure. Let me take it back. I'm not sure exactly what the eigenvalues are.
OK. So when I'm adding a [INAUDIBLE] matrix, I'm basically perturbing this by minus S. So this is going to be my eigenvalues. And because I have a small eigenvalue, I have a large eigenvalue, I would get a [INAUDIBLE].
So basically, we are looking at the example of a [INAUDIBLE] problem. We can analyze from the perspective of physics. We have a fast conduction of slow convection. And if we put these physically into equations and analyze the eigenvalues, we can get a slow eigenvalue and a fast eigenvalue.
The orders of magnitude difference in eigenvalues is also going to indicate the system is going to be stiff. Questions on this? Maybe you don't quite get through the mathematics, but it is going to be leading to the same conclusion.
AUDIENCE: Does the name derive from [INAUDIBLE]?
PROFESSOR: Good question. Does the name derive from stiff springs? So let me actually give you-- the answer is yes. The name does derive from stiff springs.
And it's actually going to be a homework question. But like I can kind of hint you a little bit right now. So what if you have a system like this? What if you have a system like this?
Instead of a thermal problem, let's think of a pendulum. So a pendulum is we have a mass over here, right? And a regular pendulum would have an inelastic string attached to the mass.
But what if you think of you have a spring here, so that the distance between the hinge and the mass can actually become longer or shorter. So if you have a system like this, in what condition-- so let's say the stiffness of the spring is K. Would you have a stiff system in terms of ODEs? Or you have a non-stiff system?
In what condition would you have a stiff system? In what condition would you have a non-stiff system.
AUDIENCE: [INAUDIBLE] stiffness [INAUDIBLE].
PROFESSOR: It depends on the stiffness of the spring, exactly. We have two physical processes over here. One is the motion of the pendulum driven by gravity. The other is the oscillation of the spring.
If the spring is very stiff, then the frequency of the oscillation is going to be very high, which means the timescale of the oscillation is going to be much faster than the motion of the pendulum due to gravity. And then we get two very different timescales. Therefore, we get a stiff system.
So that is an example of, like, how does the name stiff really comes into the picture of all this if you have, let's say, a large structure where you are interested in the motion of the whole. But inside the structure, there are some very stiff components, I mean, really stiff in terms of physically stiff components. When you derive the whole ODE, you discover two timescales.
One is the timescale you're interested in. But there is a much shorter timescale, much higher frequency timescale due to the stiff component inside the system. And therefore, you get a very difficult problem to solve. And [INAUDIBLE] this problem. Does that make sense? Yeah
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yeah. I guess my interpretation is that it's stiff, because of the stiffness of the spring. I mean, yeah, it's hard to solve. But it's because of the stiffness of the system it is hard to solve. It's like hard to solve is a result of it being stiff. Any other questions?
AUDIENCE: [INAUDIBLE]
PROFESSOR: You got it?
AUDIENCE: Yup.
PROFESSOR: OK. So before I give you another example of this system, let's talk about how to solve it. OK. As you read in your reading, in order to solve a stiff system, you really need implicit methods. Why is that the case?
Because in explicit method, it is only stable for a limited region of lambda delta t. And if there is some phenomenon in your system that is very fast, much, much faster than the phenomenons you want to resolve-- OK, for example, in this case, we may want to resolve the scale of how the temperature rise and how the temperature [INAUDIBLE]. That is the scale of the convection.
But in order for you to solve it, you have to resolve-- one of your lambdas is the lambda corresponding to a much, much faster process. What does this mean in terms of the plot? It means that if you have the stability region like this-- so this is the peanut shape of the active [INAUDIBLE], that's one of the best explicit schemes out there.
But if you have a lambda that is a thousand times the process, the frequency of the process you are interested in, that means you have to have a delta t that is a thousand times smaller than you would want to use. Because if you want to really resolve the process you want to resolve very well, you kind of only need a lambda times the smaller lambda times delta t to be in a reasonable range. But that wouldn't be stable, because of the existence of a very large eigenvalue.
So in order to make all the lambda delta t's to be within the stability region-- and remember, you have to make all the lambda delta t's within the stability region. Otherwise, you're scheme is going to be not eigenvalue stable. So in order for you to do that, you would have to have a delta t that is some times ridiculously small.
OK. So that is explicit. And for explicit schemes, all of the schemes only have a limited stability region. But if you have an implicit scheme, the story's much different.
You typically have a stability region that looks like this. It is stable for an infinitely large region. OK. You can easily make this scheme to be a stable for an infinitely region, so that even if some of the lambda are very large, you can still have a stable scheme. You can still use a delta t as large as the timescale you are interested in. Does it make sense?
Now, implicit method has a lot of benefits. But it also have some disadvantage, especially for [INAUDIBLE] your system. As you saw, [INAUDIBLE] system each step involves solving a nonlinear equation.
So for example, I'm using the most simple scheme backward [INAUDIBLE]. In backward [INAUDIBLE], the difference between V n plus 1 Vn minus over delta t is equal to the right-hand side of V n plus 1 instead Vn. If it's Vn, then it's forward [INAUDIBLE]. But if it's Vn plus 1, it's backward [INAUDIBLE]. It looks like a very innocent difference. But if f is nonlinear, this is hell a lot different.
Because the unknowns are also within here. So this is unknown. V n plus 1 is unknown. Only Vn is known.
If this is Vn, I can evaluate f of Vn. But if this is V n plus 1, I don't know what V n plus 1 is. So I cannot evaluate this. I have to solve the whole thing as a nonlinear equation.
OK. On the other the day, we said if f is a linear equation, if f is linear, we can express this in a matrix form and rearrange the equation and compute it by solving a linear equation. But if f is nonlinear, how do you do it? Of course the answer is Newton-Raphson.
But to understand what Newton-Raphson does, I think it benefits first to review what we would do in the linear case. In linear case, f of V n plus 1 can be represented using a matrix times the n plus 1. OK. Then what we can do is the following.
We can move all the known terms on the left-hand side, I over delta t minus A times V n plus 1 equals to the-- this is unknown equal to the knowns, which is Vn divided by delta t. So we are moving [INAUDIBLE] this is [INAUDIBLE] and moving this term, which is now here, onto the left-hand side. Once we arrive at a matrix form, we can use [INAUDIBLE] matrix and solve for the Vn plus 1 right?
Now, let's dive into the nonlinear case. OK. We don't know how to solve the nonlinear equation. So why don't we approximate it using a linear equation?
Why don't we approximate the nonlinear term f of V n plus 1. Try something linear. How do you do that?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Small perturbation Taylor series. This is another use of Taylor series completely different from how we used it in discretizing ordinary differential equations. We linearize f of V n plus 1 at some kind of initial guess.
So we first guess what V n plus 1 is. And usually, a good first initial guess is just Vn. Because if I know I'm taking small time steps, then I know the next time step I have a V n plus 1 that is really similar to what Vn is.
So I'm going to set f of Vn plus 1 equal to Vn equal to f of Vn plus delta V where delta V is V n plus 1 minus Vn times the derivative of V. So instead of writing it like this, I'm going to write a matrix, partial f partial V at Vn times delta V. So I'm going to write it like this.
Now, in this Taylor series expansion, I only can write a problem. Because I'm ignoring the [INAUDIBLE]. Now, in this Taylor series expansion, what is known? What is unknown? This is known, right?
And the unknown-- well, the unknown is only this. The matrix is also known. So somehow we reduced the nonlinear equation into a linear question.
Because all the terms that involves the unknown is linear, right? That's write this down and convince ourselves. Let's plug this into the original equation.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yes.
AUDIENCE: [INAUDIBLE]
PROFESSOR: So half of [INAUDIBLE] is unknown. But once we [INAUDIBLE] like this, we found the only unknown in terms of a linear term is [INAUDIBLE]. So let's do this. Let's plug into the scheme.
What we found is that so Vn plus 1 minus Vn is my delta v, right? My delta V is unknown. Delta V divide by delta t is equal to this. So the first term is known. And the Jacobi is known. The Jacobi is the derivative of f to V. And my delta V is over here.
OK. Now, I can solve this equation using exactly the same way which I solved the linear equation. I can say I, which is identity over delta t minus my Jacobi, partial f partial V at Vn times my delta V is equal to f of Vn. OK. Now, I can solve for the delta V. And I'm going to say my f Vn plus 1 [INAUDIBLE], which is my guess, is equal to Vn plus this quantity I just solved for, delta V.
This Vn plus 1 [INAUDIBLE] hopefully is a better approximation to V n plus 1 than my initial guess Vn. But how can I be sure? How can I be certain it is actually better?
AUDIENCE: [INAUDIBLE]
PROFESSOR: I calculate the residual, exactly. Now, I have this quantity. I can plug this quantity into the equation. I'm going to say is it equal to f of V n plus 1 tilde.
If this v tilde n plus 1, when I plug it in, is equal to the right-hand side, more approximately than if I plug in Vn as V n tilde plus 1. Then I know I have a better approximation. If this is almost exactly equal, then I know I already got a very good answer. What if I improved, but still didn't have an as accurate answer as I want?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yes. I keep iterating. Now, I'm going to define another delta n. So if not good enough, I'm going to say, OK, what I really want is equal to my guess plus my new delta V.
This is a new delta V. And do this again. And do this again. How do I do this again?
The unknown, which is f Vn plus 1, now is equal to f of v tilde n plus 1. Now, this is known, because I already computed the better approximation. Plus the derivative of f to V at this tilde V n plus 1 times my new delta V-- so this is iteration two.
And previously, this is iteration number one. And here is my iteration number two. In iteration number two, I constructed a new Taylor series expansion, I constructed a new approximation based on my updated guess.
And then the next thing is the same. I'm going to plug in this equation into the formula. And the formula is V n plus 1 minus Vn. Now, my V n plus 1 is equal to this, right?
So I actually have a V tilde n plus 1 minus Vn plus my delta V divided by delta t is equal to-- oh, let me actually color the terms that is known and unknown. These are known equal to f of v tilde n plus 1 plus the same thing times the unknown delta V. Sorry, should be blue. And again, I'm going to move all the known terms into one side and the unknown terms to the other side.
It's equal to f of V tilde n plus 1, which is known, minus delta t of [INAUDIBLE] minus Vn. I'm going to solve for that delta V again. Once I solve my delta V, I'm going to update my guess to be my first iteration guess plus the update again.
And then here goes my iteration number three. I'm going to keep iterating until the residual drops to into a tolerance level that I think is small enough. So the Newton-Raphson iteration is really another way of using Taylor series analysis by using the fact that we are able to solve these increases scheme equations.
So if we derive a algebraic equation using implicit scheme, we know how to solve that if it's a linear equation. We don't know how to solve that if it's a nonlinear equation. But we can approximate it using a linear equation and approximate it again and again. When I have a better estimate, then the Taylor series approximation is even better. And that should give me an even a better approximation at the next step.
All right is it clear? I mean, one of the Newton-Raphson iterations is one of the most useful numerical techniques that really goes well beyond this subject. It really accomplishes what you can-- if you go to MATLAB, you can use fsolve.
fsolve usually helps you solve nonlinear equations. But fsolve only works if you have a small amount of equations you want to solve. When we go to [INAUDIBLE] when we want to solve-- [INAUDIBLE] had [INAUDIBLE] you're going to find out maybe we have millions of equations we want to solve at the same time. And the fsolve it going to be helpful in the situation.
And in these cases, you really need Newton-Raphson equations [INAUDIBLE]. Any other questions on Newton-Raphson. Do you guys know how to implement a Newton-Raphson equation? Sure?
Because I just got some consensus right before class on some of the questions. When you try to solve-- let me say one word. When you try to solve this equation, you want to solve a matrix times delta v equal to a residual. So this is the residual, right?
It is really important to construct the A in the right way. Because if you don't be careful, you will get A transposed instead of A. OK. And when you solve delta V equal to A inverse times R, it is really important that A inverse is on the left-hand side of R instead of on the right-hand side. Because in linear algebra, the order really matters.
OK. And how the entries in a matrix A is placed is also important. Let me make a new page just to make my point. This I think can save a lot of headaches.
If delta V is arranged at V1, V2, et cetera to Vn and R is arranged as R1 et cetera to Rn, and the A is arranged as the derivative of R1 to V1 et cetera to derivative of R1 to Vn and the derivative of Rn to V1 derivative of Rn to Vn. It is important you [INAUDIBLE] rather than the transport of the matrix. Because it's correct Taylor series expansion is this matrix times this vector.
If you do the transpose of this matrix, then you would be modifying the derivative of Rn to [INAUDIBLE] Vn, which is going to be [INAUDIBLE]. And then you do delta V equal to A inverse times R. So just want to make sure you have the right denominator and the [INAUDIBLE] inside the matrix. Is this clear? When you do this in MATLAB, it's A backslash R. That gives you A inverse times R. Questions? No?
OK. So the rest of the 20 minutes, I want to have some fun. So I want to kind of have you try in your own computer a real example of a stiff system. And I have a joystick over here, so you can use.
And if some of you have Simulink installed on your computer, feel free to grab your computer, come to here, and see what is your own implementation of the system. So this is system is like this. So let' build a MATLAB flight simulator.
OK. [INAUDIBLE] build a flight simulation in MATLAB. In class, I won't be able to simulate the whole airplane. I'm only to be able to simulate the pitch motions, the longitudinal motion of the airplane, so no turning. All right.
So OK. So assume the [INAUDIBLE] coefficient is equal to 2 pi times alpha. So we're assuming we have [INAUDIBLE] air foil-- if you do think [INAUDIBLE] theory, that is exactly what I'm going to get, 2 pi times alpha angle of attack. And we have the lift and drag to be proportional to Cd and proportional to Cd and Cl.
And the dynamics is that dvdt, the derivative of my velocity, is equal to minus D, the drag times the gravity component in the backwards reaction. So if the my airplane has a drag, has a pitch of theta, I'm going to have deceleration mg cosine theta.
And the rate of change of theta, which is my pitch, is proportional to the lift minus the gravity component in the cosine direction. So that is my dynamic. So see [INAUDIBLE], I have two couple differential equations already where the terms d and l are expressing in terms of these.
Now, the rate of change of alpha, that is the angle of attack, I am going to model this as input minus alpha divided by Tau. So that is the model of the longitudinal stability. OK. So this is the case where I actually have a longitudinal stability.
Otherwise, the sign would be the other way. So alpha is my input. I'm going to use the [INAUDIBLE], add my input in. And alpha is the angle of attack.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Tau is a value I'm going to set. Tau is going to be-- so what do you think the value of Tau should be if I have a very stable system, if my Cd is way in front of my aerodynamic sample?
AUDIENCE: [INAUDIBLE]
PROFESSOR: So if I have a very stable system, Tau should be very large or very small?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Large? Are you sure? So what happens if I have a very stable system? If I pitch the airplane off, is it going to come back fast or come back very--
AUDIENCE: Fast.
PROFESSOR: Fast, right. So if I have a very stable system, Tau is going to be very small. If I have a marginally stable system, Tau is going to be very large. So Tau has the unit of pi. It's really kind of how long that it takes for the airplane to cut to static stability.
I have two files. One is the longitudinal [INAUDIBLE] n. So this function is really [INAUDIBLE] exactly the same differential equation I just showed on the slide. So I have four variables. They are the velocity, the pitch angle, the angle of attack, and my altitude.
I have all the [INAUDIBLE] defined here. I have my Cl equal to 2pi times alpha. I can do the drag in this. [INAUDIBLE]. And I have the dvdt, the theta dt, the alpha dt, and [INAUDIBLE]. So if you guys implement your own computer using [INAUDIBLE] style, [INAUDIBLE].
And here, I have Tau equal to 0.001. Does that mean I have a very stable system or a not very stable system? I got a very stable system, right. I get a very stable [INAUDIBLE].
AUDIENCE: [INAUDIBLE]
PROFESSOR: And if you succeeded [INAUDIBLE] you can [INAUDIBLE] input [INAUDIBLE] alpha [INAUDIBLE]. It's kind of like visualizing [INAUDIBLE]. For example, I can [INAUDIBLE]. It just kind of visualizes [INAUDIBLE]. This is a [INAUDIBLE]. I think I [INAUDIBLE].
AUDIENCE: [INAUDIBLE]
PROFESSOR: [INAUDIBLE] your own [INAUDIBLE] on the joystick, you can really control your aircraft. You can really control your [INAUDIBLE] simulated [INAUDIBLE]. So let's do this. Just for the rest of the five minutes, I'll just let you have fun for yourself. And one thing I want to show you is that-- so here, I [INAUDIBLE].
And the way I did is [INAUDIBLE] can save this file. And the way I [INAUDIBLE] the joystick input is by setting [INAUDIBLE] equal to [INAUDIBLE]. That actually gives you a joystick.
It's like a handle. It gives you a joystick handle. And in [INAUDIBLE] from the dropbox, what I did is the input equal to axis. So I did the input [INAUDIBLE] to axis [INAUDIBLE] 2. So this is the axis when I [INAUDIBLE].
And the [INAUDIBLE] is the third axis. And it happens to be this one. The first axis I think we [INAUDIBLE]. If you guys have time, that'd help me explain this to you something more like real simulator, that's be great.
OK. This is [INAUDIBLE]. [INAUDIBLE] my view which is the [INAUDIBLE] of the velocity, each angle, angle of attack, and altitude is equal to the last variable plus delta t times the function [INAUDIBLE]. And I've [INAUDIBLE] and did the display. The [INAUDIBLE] is asking [INAUDIBLE] the whole thing. But [INAUDIBLE] that [INAUDIBLE] very much depends on this thing, very much depend on my Tau.
So right now, my Tau is plus 1. It's stable, but not that stable. If I [INAUDIBLE] it to 0.001, it's going to be a more stable system. But what does that mean? That means it's [INAUDIBLE], right?
It has a smaller timescale that is much smaller than the one I'm interested in, which is the [INAUDIBLE] the longitudinal motion of aircraft. So if I set to be that value in MATLAB for [INAUDIBLE]. If I want to simulate [INAUDIBLE], what happens is my altitude is 10 to the 16.
And basically, the airplane goes unstable. Not unstable physically, but unstable numerically. It's not the airplane goes unstable. But it's my [INAUDIBLE] integrator goes unstable. It is particularly unstable in the pitch direction. Because in the pitch dimension, i have a timescale that is much, much smaller.
So if you guys are interested, take the [INAUDIBLE] I have and [INAUDIBLE] integrate on it using Newton-Raphson. That way if somebody is able to do that and bring the demo next class, that's fantastic. And that also saves you the homework this week. Because the homework this week is going to be about that.
So if you do that before next class, you can demo in class, and you already solved a big part of this week's homework. All right. OK. And what you [INAUDIBLE] is already in the [INAUDIBLE]. So I will see you all Wednesday.
