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.
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.