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.
WILLIAM GREEN, JR.: Shall we get started?
Yes.
Apologies for the notes being posted so late.
I had trouble with the new two factor authentication system,
which I need to get into Stellar and put things up there.
My phone wasn't working so I had to find my iPad someplace
and use that to get credentials to log in.
So if you don't have them don't worry about it.
The full notes are posted online, no missing pieces.
So you can utilize those afterwards
if something is unclear.
All right?
So your quizzes are graded.
The TAs have them.
I think overall, we were really pleased
with the outcome of this quiz.
It is very conceptual in nature and we
thought people really understood a lot of the material
quite well.
And you can see that from the distribution of scores.
So the mean on the quiz was about 77,
a standard deviation of 12.
I think everyone did really quite well.
And on a number of the problems, perfect solutions were posted.
So I'll take a second here.
If there are any questions about the material that
was on the quiz, or the format of it
that you'd like to have addressed,
I'm happy to answer them.
No?
Does it seem like it took too much time
to complete all three questions in the two hours allotted?
Or was it about the right length of material for the time slot?
So we vetted it with the two TAs.
And it took each of them less than an hour to complete it.
So that was a good sign that we thought
you guys should be able to get through most
of the material in two hours.
So let's recap.
We're going to move on today to a topic called differential
algebraic equations.
But before doing that, there are a couple of concepts
I want to review from numerical integration,
or actually concepts that weren't covered when we talked
about numerical integration on Friday
that I think are important.
And then I want to review briefly
implicit methods for ODE-IVPs because those
are going to be important for differential algebraic
equations.
So for numerical integration we talked
about quadrature of various integrals,
how to develop quadrature formulas.
But there's actually one type of integral
we didn't talk about there, which
is sort of improper integrals, where the limits of integration
are unbounded.
And these come up all the time in engineering problems.
We'd like to know what the, say, integrated
response of some function is over all possible times.
And we may have to evaluate that numerically.
So how do you do those sorts of integrals?
And it turns out the strategy that's most often used
is to divide this domain up into two pieces.
One piece, which is proper integral
over a finite domain, and another piece,
which is an improper integral over an infinite domain.
This first integral, you can handle it
with an ODE-IVP method or some higher order polynomial
interpolation.
But the second one we have to do differently.
And there are couple of ways to do this.
One is to transform variables.
So you map this domain onto a finite domain.
And the other option is to substitute
an asymptotic approximation.
So we say, when, say, tau is large, f of tau
has a characteristic asymptotic approximation.
Maybe it's exponentially decaying.
And we take the integral of that asymptotic approximation.
And the more accurate that approximation
is the better our approximation for this whole
integral will be.
There's another type of interval that we have to do,
where the same idea can be applied.
That's the integral over integrable singularities.
So here's an example.
So say we want to integrate the function cosine of tau
over square root of tau from 0 to some finite positive time
point.
As tau goes to 0, cosine goes to 1 and square root
of tau diverges.
But this integral actually has a finite value.
This 1 over square root of tau, that's
an integrable singularity.
So how do you do this?
Well, you want to split the domain again into two parts.
One part that contains the integrable singularity
and the other part, which excludes it.
And in the part that contains the singularity
we might do an asymptotic expansion of the function,
and integrate each of those asymptotically accurate terms
separately.
So the integral of 1 over square root of tau
is going to give us 2 tau 0 to the 1/2.
The integral of this ratio here is
going to give us minus 1/2 t0 to the 5/2.
And then we have an integral over a finite domain,
which doesn't include the singularity.
And the accuracy of this method can be dictated by two things
now.
One is how accurately can we do this integral.
And the other is how accurate is our asymptotic expansion here.
So if we want higher accuracy we need more terms.
We might need to be selective about how
we choose the end point for this domain
in order to minimize the error.
There are lots of methods to do this.
Sam.
AUDIENCE: Will you talk a little more
about what typical integral [INAUDIBLE] is?
WILLIAM GREEN, JR.: Yeah.
So you know these things.
If I try to integrate 1 over x, as x goes to 0
this thing will always diverge.
You know the antiderivative of 1 over x is going to be the log.
And the log diverges.
So power is weaker than 1 over x, like 1 over x to the 1/2.
Don't diverge like the log.
They actually go to a finite value.
The log is like the most weakly singular function there is.
And integrals over weaker powers of 1 over x actually
don't produce a singularity anymore.
They're what we call integrable singularities.
The function itself diverges but its integral does not.
So you usually have this asymptotic power law type
behavior.
OK.
So they come up all the time.
We see these things in places like, say,
squeezing a thin film of fluid between two plates.
How long does it take for these two plates to come together?
You might say, well, as the plates get closer and closer
together the pressure starts to diverge in the gap
and they'll never come together.
But maybe depending on the geometry of the plates
there can be a sort of finite time
at which they'll come together.
Depends on their geometry.
So that's one topic.
I think that's useful because these things come up
all the time.
And you can't actually handle those cases with quadrature
by itself.
The function diverges.
There's no polynomial interpolation
that's suitable to match it.
So the quadrature can't handle it.
And if you're trying to integrate over
an infinite domain, well, good luck.
You're never going to be able to fit enough points in
to know that you've accurately integrated out to infinity.
The other thing I want to recap here,
and this is a problem for you to try.
So here we have a simple first order ordinary differential
equation.
We want to use the implicit Euler method to solve it.
So can you give a closed form solution
for one step of the implicit Euler method,
or for n steps of the implicit Euler method?
What result does that produce?
Can you compute this?
Do you guys talk about backward Euler methods for for ODE-IVPs?
Backwards difference methods in general and the backwards Euler
method specifically?
No.
OK.
Well, clearly I should have been here.
So what's the strategy?
We have to approximate the derivative.
Yes?
So the approximation for the derivative we use
is a finite difference approximation.
So we write dx dt as x at a point k
plus 1 minus x at the point k divided by delta t.
Backwards difference.
We evaluate the right-hand side of the ODE-IVP at k
plus 1, xk plus 1.
And then we solve for xk plus 1.
We just substitute in here.
Backwards difference approximation
for the differential.
And xk plus 1 for the right-hand side.
And we find that xk plus 1 will be 1 over 1
minus delta t times this lambda times xk.
And if we iterate this from k equals 0 up to some finite k,
it's like multiplying by powers of 1 over 1
minus delta t times lambda.
Yes?
Does this makes sense?
So our backwards Euler solution to this fundamental problem
is going to look like this.
And a lot of times we ask about the stability
of these sorts of solutions.
So if for example, the solution to this equation
was supposed to decay to zero, we
would expect our numerical method
to also yield results which decay to zero.
We would call that stable.
So we all know under what conditions
does xk, for example, go to 0 as I change the product delta
t times lambda.
Lambda could be a real number or it can be a complex number too.
Doesn't really matter.
So what needs to be true for xk to the k to 0 as k gets bigger?
Well what needs to be true is this quantity in parentheses
needs to be smaller than 1 because I'm
taking lots of products of things
that are smaller than one.
That will continue to shrink xk from x 0 until it goes to zero.
If this is smaller than 1, then this quantity down here
needs to be bigger than 1 in absolute value.
And then I'll get solutions that slowly decay to zero.
So if I ask, for what values of delta t times lambda
is the absolute value of this thing bigger than 1,
and I allow lambda to be a complex number, for example?
I'll find out that's true when 1 minus delta t
times the real part of lambda squared
plus delta t times the imaginary part of lambda squared
is bigger than 1.
And if I plot the real part of delta t times
lambda versus the imaginary part of delta times lambda,
this inequality is represented by this pink zone
on the outside of the circle.
So any values of delta t times lambda
that live in this pink area will give stable integration.
And the solution xk, well, the k to zero as k goes to infinity.
So it makes sense?
This is a little funny though.
What does the solution to this equation
do when lambda is bigger than 1?
When it's bigger than 0, excuse me.
When lambda is bigger than 0, what
does the solution to this ODE-IVP do?
Does it decay?
AUDIENCE: [INAUDIBLE]
WILLIAM GREEN, JR.: It blows up.
It grows exponentially.
But the numerical method here, when
I have lambda is bigger than zero,
can be stable and produce decaying solutions.
So this is sort of peculiar.
So we're often concerned when we solve ODE-IVPs
and we try to use these sorts of backwards difference
methods, these implicit methods, to solve them.
We're often concerned with getting stability.
We want to get solutions that decay
when they're supposed to decay.
Sometimes we get solutions that decay even
when the solution is supposed to blow up.
So that's sort of funny.
So the exact solution should look like this.
But there are going to be circumstances where
lambda is bigger than zero.
In fact, lambda is bigger than--
lambda delta t is bigger than 2 along the real axis,
where the solution will actually blow up rather than decay away.
So stability and accuracy don't always
correlate with each other.
We really have to understand the underlying dynamics
of the problem we're trying to solve
before we choose a numerical method to apply to it.
And this is going to be true of differential algebraic
equations as well.
Are there any questions about this?
There's a discussion of these things in your textbook,
these sorts of stability regions associated
with different backwards difference methods.
The implicit Euler method is one backwards difference formula.
It's the simplest one.
It's the least accurate one that you can derive
but it's an easy one to use.
But now we want to move on to a different sort of problem.
And it's one that we come across in engineering all the time.
These problems are called differential
algebraic equations.
And they're problems of the sort a function of some state
variables, and the derivatives of those state variables
and time is equal to zero.
And there are some initial conditions.
The state variables at that time 0 is equal to some x0.
So this is a vector valued function
of a vector valued state and it's time derivatives.
And we want to understand the dynamics of this system.
How does x vary with time in a way
that's consistent with this governing equation?
Usually we call these well-posed problems
when the dimension of the state variable and the function
are the same.
They both live in the vector space rn.
So we have n different states.
I have n different functions that I have to satisfy.
And I want to understand how x varies with time.
Here's an example.
A stirred tank.
I call this stirred tank, example one.
So into the stirred tank comes some concentration of solute.
Out of the stirred tank comes some different concentration
of the same solute.
And we want to model the dynamics as a function of time.
So we know that a material balance
of the solute on the stirred tank, if it's well mixed,
will tell us that d c2 dt is related to the volumetric flow
rate into the tank divided by the volume of the tank
multiplied by the difference between c1 and c2,
the amount carried in and the amount carried out of the tank.
And let's put a constraint in here too.
Let's suppose we're trying to control this tank in some way.
So we say that c1 of time has got
to be equal to some control signal, gamma of t.
And then we want to understand the dynamics of this system.
At time 0, c2 is some c0, the initial concentration
of the tank.
And c1 is whatever the initial value of this control signal
is.
And the solution is c1 equals gamma of t.
And c2 is, well, the initial concentration's
got to decay out exponentially.
And then I've got to integrate over how
the input is changing in time.
Those also produce exponential decay.
So if the input is a little delta function spike
then I'll get an exponential decay leaving the tank.
If the input is changing dynamically in time,
I need this convolution of the input with this integral.
So you just to solve this system,
this ordinary differential equation
by substituting c1 in there.
But this system of equations here is not just
a differential equation.
It's a differential equation and an algebraic equation.
There's something peculiar about that that's
different than just solving systems
of ordinary differential equations.
Oftentimes in models, we write down all
the governing equations associated with the model.
Some of them may be differential, some of them
may be algebraic.
There are times when we can look at those equations
and see clever substitutions that we can make.
But if you have a hundred equations, or a thousand
equations, or a million equations,
there's no way to do that reliably.
We just kind of wind up with a system of ordinary differential
equations and algebraic equations mixed together.
And we have to solve these reliably.
Here's another example.
Oh, excuse me.
So let me write this in this form.
So we said there should be a vector valued
function of the states and the derivatives of the states.
The states are the concentration c1 and c2.
And the derivative here is dc 2 dt.
So here's my vector valued function.
The first element is this, and it's equal to 0.
The second element is this, and it's equal to 0.
Here's my initial condition over the states.
So I can just transform that list of equations
to the sort of canonical form for a differential
algebraic equation.
Here's a separate example.
Suppose instead, I'm trying to either control
or make a measurement of the outlet concentrations c2.
And we call that gamma of t.
And I want to know now what is c2 and what is c1?
So I've got to solve this system of equations for c2 and c1.
The solution for c2 is easy.
I know that.
c2 is gamma, by definition.
I got to substitute this into the first equation
and solve for c1.
And I see, c1 is gamma plus v over q gamma dot.
And I have to have initial conditions to go along
with this.
There's something funny here.
We can see the initial condition for c2 has to be gamma 0.
The initial condition for c1, what does that have to be?
Well, it needs to be consistent with the solution here.
There were no free variables when
I solved this equation for c1.
It isn't like solving a differential equation,
having a free coefficient to specify.
So it better be the case that this c0 here
is the same as gamma 0 plus v over q gamma dot 0.
Somehow I have to know this input signal to prescribe
the initial conditions for c1.
So that's peculiar and different from just ODE-IVPs.
You see this picture?
You see how this goes together?
The solution is funny too.
Suppose we are trying to use this system of equations
to do the following.
We are trying to measure c2 and use the measurement of c2
to predict c1.
So gamma is the measurement, and we're
trying to solve this system of differential algebraic
equations to predict what c1 is.
All measurements incur numerical error.
So even though our signal for gamma that we measure
may be continuous, it's bouncing around wildly.
It's not a constant signal.
It's going to move around a great deal.
And c1, our predicted value for c1,
depends on the derivatives of gamma,
which means c1 is going to be incredibly sensitive to how
that measurement is made.
So it's peculiar.
It means that there's going to be a lot of problems
with stably solving these equations because the solutions
can admit cases where they're not integrals of the input
but they're actually related to derivatives of the input.
Look at the solution back here again.
Oop, sorry.
The solution here was related to an integral of the input.
Suppose I was making a measurement of c1
in trying to predict c2 instead?
My prediction for c2 would be related to the integral
of that measured signal.
And we know integrals smooth signals out.
So it's not going to be so sensitive to the measurement.
So one way of doing this seems really sensible.
And the other way of doing it seems a little bizarre.
Well, we can construct these equations however we want.
So there's something about how we formulate models
that are going to make them either easily
solvable, as DAEs, or quite difficult to solve.
So these pop up all the time in engineering problems.
They pop up in dynamic simulations,
where we have some sort of underlying conservation
principle, or constraints, or equilibria.
So if we have to conserve something,
like total energy, total mass, total momentum, total particle
number, total atomic species, total charge,
there will always be an algebraic constraint that
tells us the total mass has to stay this,
or the total number of atoms of a certain type has to say this.
Yeah, Ian.
AUDIENCE: Just to be clear, were you saying, on the stirred tank
example, that those were the same problem
with different solution approaches?
WILLIAM GREEN, JR.: So good question.
Let me go back.
So they are fundamentally different problems.
So in one case I specified the value of c2 to be this gamma.
I tried to make a measurement of c2 and then predict c1.
In the other case, I specified c1.
I tried to measure c1 and predict c2.
So one case I measured the input and tried
to predict the output to the tank.
And in the other case I measured the output
and tried to use it to predict the input.
Yeah?
AUDIENCE: So physically the same system
with a choice by the operator.
WILLIAM GREEN, JR.: Yes.
That's right.
So I formulated the model differently.
The same underlying system but I formulated
the model differently.
And it led to completely different behavior
in the underlying solution to the differential algebraic
equations.
So you've already solved problems
that involve conservation.
You've done things like ODE-IVPs on batch reactors,
where the total amount of material in the reactor
had to remain constant.
Instead, you probably solved for the dynamics
of each of the components undergoing the reaction.
You might have checked to see if at each point in time
the total mass in the reactor stayed the same.
Because you solved all these differential equations,
they were interconnected with each other
but they all incurred some numerical error.
It's possible that numerical error
accrues in a way that actually has you losing mass
from your system.
So there may be benefits to actually trying
to solve these sorts of systems of equations
with, say, one less differential equation
for one of the components and instead an algebraic constraint
governing the total mass or total number of moles
in the system.
So these pop up all the time.
You see them in models of reaction networks,
where we try to use a pseudo steady state approximation.
We say some reactions go so fast that they equilibrate.
So they're not governed by differential equations
but by algebraic equations for equilibria.
They can pop up in models of control, where we neglect
the controller dynamics.
So you say I'm going to try to control this process
and I get to instantaneously turn on or off this valve
in response to a measurement.
Actually, controllers have inherent dynamics in them.
But if I don't put those dynamics in then
I may get an algebraic equation instead of a differential
equation for how the control process occurs.
There are different types of DAEs that we talk about.
So there's one type that's called semi-explicit DAEs.
And in these types of differential algebraic
equations we can write them in the form
M dx dt is equal to f of x and t, some initial condition.
You say that this looks like an ODE-IVP.
M is called the mass matrix.
And it may or may not be the case that M is invertible.
M may or may not be a full rank matrix.
So let's look at that stirred tank example one.
This was the underlying equation, f of x and dx dt
and t.
Can write it in semi-explicit form.
If c has two components, c1 and c2, then dc 2 dt-- oops,
this is a typo here.
I really apologize.
This should be a 0, and this element should be a 1.
I'll correct that online.
I apologize for that.
So this should be this matrix multiplied
with this differential should give me
dc 2 dt, this term here.
And a 0 for the second equation.
And then I've got q over v multiplying c1 and minus c2.
So that gives me the first line here.
I've got minus c1 coming on the second equation here.
And those balance with a 0 and a gamma.
So the lower line reads like this equation here
and the upper line--
fix this typo-- replace 0 with one,
reads like the upper equation here.
This is the semi-explicit form of the differential equation.
You can see this mass matrix is singular and has
a row that's all zeros.
There's no way to invert this matrix
and convert it into a system of ODE-IVPs, which you already
know how to solve.
If the mass matrix is full rank, which can happen,
you can formulate your model in such a way
that it takes on this semi-explicit form.
But if the mass matrix is full rank and invertible,
then you just invert the mass matrix.
And now you can apply all your ODE-IVP techniques
to solve this thing.
The mass matrix doesn't need to be constant.
Could depend on the concentration.
It can depend on the states as well in this form.
I just chose an example where that isn't the case.
But in general, it can depend on the states.
If it depends on the states and we invert it,
then we need to be sure that the mass matrix is
invertible for all values of the states, which may or may not
be easy to ensure.
So here's what I'd like you to try to do.
So here's that stirred tank example two.
The only difference between this and the previous one
is that c2 now is set equal to gamma instead of c1.
I want you to try to write it in semi-explicit form.
A sort of test of understanding.
Can you define the mass matrix?
Can you write out the right-hand side
of the equation in semi-explicit form?
Take a second, work with your neighbor.
See if you can figure out how to do that.
AUDIENCE: Is this clear?
OK to do this?
WILLIAM GREEN, JR.: You can look back on the previous slide.
Actually, the semi-exclusive form is going to be the same.
This is a 0, that's a 1.
The only difference is going to be c2 is now
balanced with gamma, the input.
So this minus 1 is got to shift over.
So I switched the 0 for the minus 1,
you got the second, the semi-explicit form.
There's another way that these equations can pop up.
So we have this mass matrix.
And it might be the case that M is
diagonal over many of the states and then
0 for all the other states.
So we might be able to write the equation in this form.
This comes up all the time.
So we might be able to write dx dt is
a function of x, y, and t.
And 0 is equal to another function of x, y, and t.
And we have some initial conditions
prescribed for the x's, of which we're
taking the differential in one of these equations.
And we've also got to satisfy those initial conditions
for these variables y associated with the second nonlinear
equation.
The x's are called the differential states
because they're involved in equations where they're
difference with respect to time, or the differential
with respect to time is taken.
y are called the algebraic states
because they only appear as themselves and not
as their time differential in any of these equations.
So we might want to solve for both the differential
and the algebraic states simultaneously.
We have to know these algebraic states to determine
the differential states.
We have to know the differential states to determine
the algebraic states.
AUDIENCE: Professor?
WILLIAM GREEN, JR.: Yes.
AUDIENCE: [INAUDIBLE]
WILLIAM GREEN, JR.: That's right.
So here the x that I have is not all of the states.
It's not all of the unknowns.
It's only the set of the unknowns of which I
have to take a differential.
That's right, that's right.
So we've sort of tooled down the kinds of equations
we might want to look at.
We've gone from the generic differential algebraic equation
to a semi-explicit form to the splitting, this differential
algebraic splitting.
Many problems can be formulated in this fashion.
This one's the easiest one to think
about so that's the one that we'll discuss
for most of the lecture.
So let's look at some examples.
I think examples are interesting to think about.
So here's one.
Think about a mass spring system.
There's an algebraic equation that
governs a conserved quantity, the total energy
of this mass spring system.
So the kinetic energy, 1/2 m times the velocity of the mass
squared.
Plus the energy stored in the spring, 1/2 kx squared
has got to be equal to the total energy.
And the total energy can change in time.
So we have a differential algebraic equation.
We have a differential state because we're taking
the time derivative of x.
And we'd like to be able to determine
x as a function of time.
So f of x, x dot and t.
That's this equation, 1/2 mv squared plus one half kx
squared minus E is equal to 0.
We'd like to solve this differential algebraic
equation.
It's well-posed.
We have one equation for one unknown, one state.
The equation has a solution, which
is a times cosine omega t.
It's not the only possible solution.
It's a solution to this equation.
We can substitute that solution into the equation
and determine, for example what this oscillation frequency
omega has to be.
It's square root of k over m.
Or what the amplitude has to be.
It's the square root of E over k.
So you've seen this from physics before.
We've already solved differential
algebraic equations.
The thing you saw before though, was actually
conservation of momentum.
You converted this entirely to a second order
differential equation in x and solved it.
But in principle, we could have tried to find various solutions
to this equation instead.
We would have said, well, we've observed masses and springs
and we know they oscillate in time.
So let's propose some oscillatory solutions
and see what values of the amplitude and omega we
need to satisfy this equation.
It would have been a perfectly acceptable way
of solving this problem.
The thing is, this mass spring problem,
the differential algebraic equation
contains nonlinearities with respect to the states.
So it's nonlinear in the velocity
and it's nonlinear in the position.
Nonlinear equations in general are hard to solve.
When you converted that to a momentum balance
you got linear equations.
You got the second derivative of position
was equal to stiffness over mass times position.
It was a linear differential equation to solve,
which is easy to solve.
So you've got higher order differentials,
but you linearize the equation so it made it easier to solve,
in a certain sense.
But we'd like to sometimes understand what
we can do with these equations.
There are certain limiting cases where
we can do interesting things with these equations.
I'm going to show you one now.
So in the case where the Jacobian of this
function f with respect to the differential of the state
variable.
So take the derivative of f with respect
to x dot and hold time in the state constant.
So this is a matrix, one that's full rank.
Then the DA be represented as an equivalent to ODE.
Here's how you can show that.
The total differential of f is this Jacobian with respect
to x dot times the change x dot plus the Jacobian with respect
to x times the change in x.
Plus the partial derivative of f with respect to time,
times the change in time.
And this total differential has to be equal to 0
because f is equal to 0 for all states in time.
So let's divide by dt.
We'd like to know how f changes with time.
So total differential of f with time.
We divide each of these little differentials by time.
And then we solve for dx dot dt.
That's going to involve moving this term
to the other side of the equation
and inverting this Jacobian.
So if I can invert the Jacobian then
I can convert my DAE system into a system of ODEs.
This equals 0 shouldn't be here.
I apologize.
That's another typo.
I move this term to the other side of the equation
to where the 0 is and I solve.
So I went from a first order equation, a first order
equation in x to a second order equation in x.
Two differentials of x with time instead of one differential
of x with time.
But I could convert it to an equivalent
system of ordinary differential equations.
Does that make sense?
Do you see how that works?
Any questions?
Let's go back and look at an example.
The mass spring system.
So here's our ODE that we're trying to solve.
I'm sorry, our DAE that we're trying to solve.
Can you convert this to an equivalent ordinary
differential equation using the approach I just described?
Think you guys can figure that out?
Try it.
You can work together.
[SIDE CONVERSATIONS]
WILLIAM GREEN, JR.: Guys are either tired
or friendships have been destroyed
over the mid-term season.
Maybe both.
So let's compute the things we need to know.
We needed to know the partials of this f with respect
to the states and with respect to time.
So partial f respect to x dot, holding all the other variables
constant, is mx dot.
Partial f with respect to x is kx.
Partial f with respect to t is t.
And I told you before that you could
write this says dx dot dt is equal to df
dx dot inverse multiplied by--
is equal to minus df dx dot inverse multiplied by df dx.
This was on the previous slide, which just gives you
conservation of momentum.
This gives you the acceleration is equal to the force minus kx
divided by the mass m.
So we can solve either, the DAE or the ODE.
Here's a more complicated example,
but it works, well, it's more complicated.
So a lot of times we do molecular dynamics simulations.
We try to model molecules moving around.
There we also want to conserve energy.
There's usually a conserved quantity, something
we're trying to hold constant.
So suppose we try to do that.
The total energy in this system of molecules,
maybe it's 1/2 mass times the velocity squared.
Now x is the position of all these molecules.
Ad x dot is the velocity of all these molecules
plus the potential energy, which is a function of f,
a function of their positions.
So the DAE representing this constraint is this.
It's the same as for the mass spring
system, where we've replaced the potential energy
with a generic one.
This actually isn't a well-posed DEA.
We have one equation for, let's see,
if we have n molecules in here we have three n states.
We can move around in three dimensions.
So this is a really ill-posed equation.
If we try to apply the same procedure, then we say,
well, this quantity still should be conserved.
Over time it needs to be equal to zero.
Then we can try to compute all these differentials.
But we'll find out that partial f, partial x dot,
the Jacobian of f with respect to x dot,
that's just the momentum.
That's not a matrix.
That's not invertible.
So we can't convert this equation
to an equivalent system of ordinary differential
equations.
So it's just something pathological
for a one-dimensional system.
Instead, what we find is zero is equal to the velocity
dotted with the mass times the acceleration plus the gradient
in the potential.
So minus the forces acting on the molecules.
We know that if we were to integrate
the equations of motion for the molecules exactly,
the momentum balances for those molecules exactly,
then this term here would always be equal to zero.
And we'd satisfy conservation of energy.
Of course, all solutions of ordinary differential equations
require approximations.
So as you move the molecules around,
their velocities and positions will
tend to drift away from the solution
where this is exactly in balance and equal to zero.
If you desire to do that in the molecular dynamic simulation
then you better have a method that
somehow respects this geometric property instead.
That the velocities are orthogonal to m
x dot plus grad v.
And there is a set of ODE-IVP methods that do that.
They're called symplectic integrators.
And they integrate the equations of motion
while exerting some control over the error in the total energy.
So if you were just the take a ODE45
and try to solve this system of equations,
over long times you would find that the energy drifts away
from the place where it started.
Even though your positions are being integrated
with reasonable accuracy, over time
the energy will drift away from where you started.
The system will heat up effectively or cool down.
That's undesirable if you're trying to do modeling.
So instead, one uses the symplectic integrators,
which are designed to control the error and the energy
in addition to integrate the equations of motion.
So you can look those up and read about them.
I think they're quite interesting.
If you're going to do molecular modeling
you'd like to understand better how they work.
But actually, conservation of energy
here doesn't give us a well-posed DAE.
We actually need the momentum equations
to formulate or to determine the trajectories of the molecules.
Let's talk about their numerical simulation.
So let's think about these semi-explicit DAEs,
the ones that can be split between differential
and algebraic states.
It's kind of useful to look at this problem in particular.
And you might say, well, let's do the simplest possible thing.
Let's do the forward Euler approximation.
So let's take our derivative here and represent it
as the state evaluated at x plus delta t minus the state at t
divided by delta t.
Let's make the right-hand side of our equation
be evaluated at time t.
And let's try to step in time, do time marching forward
to determine the trajectory of x and t.
How do the states change?
We see the first equation is an easy one
to solve because presumably I know x of t and y of t.
Well, do I know y of t? y of t has to satisfy
this second equation.
So first, if I know x of t I need to solve this nonlinear
equation determine y of t.
Then I know x of t and y of t and I
can step forward in time to determine x of t plus delta t.
So every iteration here with the simplest possible approximation
for the differential requires solution of a system
of nonlinear equations.
We already saw before that a forward Euler approximation
is not a great one to apply to solutions of IDE-IVPs
because it has very limited stability properties.
You need small time steps in order
to stably integrate in time.
So even though this seemed like it
was easy to get a solution for the differential equation,
we've always got to satisfy this nonlinear equation here.
So inherently, simulation of DAEs are implicit.
They require solution of nonlinear equations
in order to advance in time.
There's actually no point in using this weakly stable method
for any very small time steps.
It makes more sense to choose a naturally implicit method
for the differential equation, where I can take big time steps
and still have reasonable stability of the integration.
Does that make sense?
Consider now the fully implicit DAE.
So f of x, x dot, and t is equal to 0
and have some set of initial conditions. x of 0
is equal to x 0.
Here there is no way in general to avoid solving
systems of nonlinear equations.
And so one substitutes often backwards
difference approximations for x dot.
And then you solve for the next point in time.
So here's an example with the backward Euler,
implicit Euler approximations.
So I replaced the time derivative
with the difference between the state at point tk
and the state of point tk minus 1
divided by the difference between tk and tk minus 1.
The other state variables are evaluated at tk.
So x at tk.
And I evaluate time at tk.
And I try to satisfy this equation
and use it to determine x of tk.
So I solve this system of nonlinear equations
for x at tk.
And then I go to the next point in time, tk plus 1.
And I solve the same equation, but replacing tk with tk plus 1
and tk minus 1 with tk.
So at each iteration I solve a system of nonlinear equations.
No more painful than doing the forward Euler
approximation I showed you before.
So that's good.
But still a lot of work.
We might wonder how do we control
the error in these sorts of approximations.
And that's what we're going to talk about next.
So here's the equivalent time marching formula
applied to a semi-explicit DAE.
So I've got to satisfy the equation
0 is equal to the derivative of x with respect to time,
which I approximate with the backward Euler approximation,
minus f evaluated at tk.
And the equation g of x at tk, y at tk for this x of tk.
And actually y of tk as well.
I'm sorry, I left that off. y of tk as well.
Got to solve for the differential
and the algebraic states simultaneously.
And I use a backwards difference formula
to try to do this in a stable way.
One way to think about in particular these sorts of DAEs
is as exceptionally stiff ordinary differential
equations.
Did you guys discuss stiffness?
So imagine if I replace the 0 with, say, dy dt.
And I introduced a constant on the right-hand side,
a multiplicative constant, say, 1 over lambda here.
Let's say lambda, actually.
Say I introduce a multiplicative constant lambda
on the right-hand side here.
So as lambda becomes very large, the system
becomes increasingly stiff.
The dynamics of this equation take
place over very, very short time scales.
So short in fact, that eventually,
if I let lambda diverge to infinity,
I've just got to satisfy this algebraic equation
at each point in time.
So these are exceptionally stiff equations.
And the methods for solving these equations
share a lot in common with stiff ODE-IVP solvers.
So I've got a couple of minutes left.
I'm going to work through some examples.
So consider the stirred tank example one.
dc 2 dt.
It's related to q over v, flow rate over volume.
The difference in the concentrations
coming in and out of the tank.
And c1 is equal to this input signal gamma.
Maybe that's a measurement and I'm trying to predict c2.
So can you apply the backward Euler method to these equations
to determine c1 and c2?
What would one step of the backward Euler method look
like, say, going from time 0 to time delta t, or time t to t
plus delta t, or tk to tk plus 1, however you
want to write it?
So do you get something like this?
So this equation is easy.
c1 to tk is gamma of tk.
I got a substitute up here.
c2 of tk minus c2 of tk minus 1 over delta t.
T Just tk minus tk minus 1.
And then solve for c2 of dk.
And you get a formula like this.
This formula was based on an approximation
for the derivative, which had a leading order
error that was proportional to tk minus tk minus 1.
When I go through and I solve for c2 of tk, the leading order
error in c2, the local truncation error
is going to be order tk minus tk minus 1 squared.
So as I shrink the spacing between each of my two time
points, the local accuracy, the error
induced in one advance of this time stepping algorithm
is going to scale like the step size squared.
So I have the error--
I have the spacing.
The error will go down by a factor of four.
I reduce the spacing by an order of magnitude.
The error will go down by two orders of magnitude.
This is a really desirable sort of error control and a method.
And we can do better.
You've seen us do better with other methods.
Let's look at this equation now.
This is stirred tank example two.
I just replace c1 with c2.
You don't have to write this one down.
I'll just show it to you.
The notes are online so you're not missing anything.
All these things will pop up in the notes online.
So c2 at tk is equal to gamma of tk.
I got to come up to my first equation now
and I have to substitute an approximation
for the derivative.
It's this, the backwards difference formula.
So if find c1 of tk is c2 of tk plus v over q.
c2 of tk minus c2 of tk minus 1 divided by the time step.
Remember, this approximation for the derivative
has a leading order error that goes like tk minus tk minus 1.
So my approximation for c1, the local approximation for c1,
is not second order accurate.
It's only first order accurate.
I applied the same approximation method to both equations.
One equation I got something that was second order accurate.
The next equation I got something
that was first order accurate.
They look almost identical.
It might be hard to see from the start
why it should work out that way.
That's how it works out.
Let's do one more example.
Here's a system of DAEs.
c2 dot is equal to c1.
c3 dot is equal to c2.
And 0 is equal to c3 minus gamma.
So gamma, again, is some measurement.
We solve this equation to determine c3.
We substitute c3 up here to determine c2.
We substitute c2 up here to determine c1.
But we want to find the solution with the backward Euler method.
And the solution there is going to be c3 of tk is gamma of tk.
c2 is related to the derivative of c3.
So I use my backwards difference formula
for the derivative of c3.
c1 is related to the derivative of c2.
So I use my backward difference formula for c2.
But how accurate are each of these approximations?
So here I have to determine c2 from an approximation
for the derivative of c3.
We know this backwards difference approximation
has a leading order error proportional
to delta t, the time step.
Now I got to approximate c1 with the finite difference backwards
difference approximation in c2.
We know this should incur an error proportional to the time
step delta t, except I don't know c2 at tk exactly.
Actually, c2 tk, the exact value of c2 at tk
is this plus something proportional to the time step.
Which means the exact value of c1
is this plus something that's order one because I
carry over the error from my solution
to the previous equation.
And I divide it by the time step.
So I went from having something-- well, let's see.
The first equation had order delta t
squared local truncation error.
As I shrink the spacing, my solution
gets more and more locally accurate.
And it seems to do it in a fairly robust fashion.
The next system of equations we solved
had a local truncation error proportional to delta t.
Again, I can shrink delta t and I
can make my solution more and more accurate.
That seems like what you want in a numerical algorithm.
This one, there's no hope.
Change delta t, who knows what you're going to get?
And this looks fairly innocuous.
You just look at it and you say, well, OK.
Let's solve it and see what happens.
Well actually, you can't solve this problem accurately
with the backward Euler method.
There's something fundamentally different
about this problem from the previous one.
There's something fundamentally different about stirred
tank example two from stirred tank example one.
Right you might have already perceived that.
You can see what's going on here.
c3 is equal to gamma.
That's the exact solution.
c2 is equal to gamma dot.
And c1 is equal to two derivatives
of gamma, the exact solution to this equation.
So c1 is incredibly sensitive to changes in gamma.
And that's peculiar.
And it has important numerical consequences.
So I'm going to show you, our next lecture on Wednesday,
how to formalize an understanding
of the differences between these problems.
And we'll talk about some more subtleties associated
with differential algebraic equations.
You guys can pick up your quizzes.
We should do it outside rather than
inside because the next class has to start.