WEBVTT
00:00:01.580 --> 00:00:03.920
The following content is
provided under a Creative
00:00:03.920 --> 00:00:05.340
Commons license.
00:00:05.340 --> 00:00:07.550
Your support will help
MIT OpenCourseWare
00:00:07.550 --> 00:00:11.640
continue to offer high quality
educational resources for free.
00:00:11.640 --> 00:00:14.180
To make a donation or to
view additional materials
00:00:14.180 --> 00:00:18.110
from hundreds of MIT courses,
visit MIT OpenCourseWare
00:00:18.110 --> 00:00:19.340
at ocw.mit.edu.
00:00:23.562 --> 00:00:25.270
WILLIAM GREEN, JR.:
Shall we get started?
00:00:25.270 --> 00:00:27.340
Yes.
00:00:27.340 --> 00:00:29.710
Apologies for the notes
being posted so late.
00:00:29.710 --> 00:00:33.020
I had trouble with the new two
factor authentication system,
00:00:33.020 --> 00:00:35.980
which I need to get into
Stellar and put things up there.
00:00:35.980 --> 00:00:39.610
My phone wasn't working so I
had to find my iPad someplace
00:00:39.610 --> 00:00:42.542
and use that to get
credentials to log in.
00:00:42.542 --> 00:00:44.500
So if you don't have them
don't worry about it.
00:00:44.500 --> 00:00:47.540
The full notes are posted
online, no missing pieces.
00:00:47.540 --> 00:00:49.480
So you can utilize
those afterwards
00:00:49.480 --> 00:00:51.060
if something is unclear.
00:00:51.060 --> 00:00:51.560
All right?
00:00:56.250 --> 00:00:58.070
So your quizzes are graded.
00:00:58.070 --> 00:01:00.390
The TAs have them.
00:01:00.390 --> 00:01:02.610
I think overall, we
were really pleased
00:01:02.610 --> 00:01:04.349
with the outcome of this quiz.
00:01:04.349 --> 00:01:06.360
It is very conceptual
in nature and we
00:01:06.360 --> 00:01:08.790
thought people really
understood a lot of the material
00:01:08.790 --> 00:01:09.540
quite well.
00:01:09.540 --> 00:01:12.440
And you can see that from
the distribution of scores.
00:01:12.440 --> 00:01:15.650
So the mean on the
quiz was about 77,
00:01:15.650 --> 00:01:17.250
a standard deviation of 12.
00:01:17.250 --> 00:01:19.930
I think everyone did
really quite well.
00:01:19.930 --> 00:01:23.580
And on a number of the problems,
perfect solutions were posted.
00:01:23.580 --> 00:01:26.670
So I'll take a second here.
00:01:26.670 --> 00:01:28.800
If there are any questions
about the material that
00:01:28.800 --> 00:01:31.464
was on the quiz,
or the format of it
00:01:31.464 --> 00:01:32.880
that you'd like
to have addressed,
00:01:32.880 --> 00:01:35.351
I'm happy to answer them.
00:01:35.351 --> 00:01:35.850
No?
00:01:38.580 --> 00:01:41.440
Does it seem like it
took too much time
00:01:41.440 --> 00:01:44.080
to complete all three questions
in the two hours allotted?
00:01:44.080 --> 00:01:49.510
Or was it about the right length
of material for the time slot?
00:01:49.510 --> 00:01:51.030
So we vetted it
with the two TAs.
00:01:51.030 --> 00:01:53.720
And it took each of them less
than an hour to complete it.
00:01:53.720 --> 00:01:55.969
So that was a good
sign that we thought
00:01:55.969 --> 00:01:57.760
you guys should be able
to get through most
00:01:57.760 --> 00:01:59.140
of the material in two hours.
00:02:05.560 --> 00:02:08.012
So let's recap.
00:02:08.012 --> 00:02:10.470
We're going to move on today
to a topic called differential
00:02:10.470 --> 00:02:11.650
algebraic equations.
00:02:11.650 --> 00:02:14.310
But before doing that, there
are a couple of concepts
00:02:14.310 --> 00:02:16.260
I want to review from
numerical integration,
00:02:16.260 --> 00:02:19.669
or actually concepts that
weren't covered when we talked
00:02:19.669 --> 00:02:21.210
about numerical
integration on Friday
00:02:21.210 --> 00:02:22.480
that I think are important.
00:02:22.480 --> 00:02:23.855
And then I want
to review briefly
00:02:23.855 --> 00:02:26.190
implicit methods for
ODE-IVPs because those
00:02:26.190 --> 00:02:28.710
are going to be important
for differential algebraic
00:02:28.710 --> 00:02:29.250
equations.
00:02:29.250 --> 00:02:32.550
So for numerical
integration we talked
00:02:32.550 --> 00:02:35.250
about quadrature of
various integrals,
00:02:35.250 --> 00:02:36.960
how to develop
quadrature formulas.
00:02:36.960 --> 00:02:38.130
But there's actually
one type of integral
00:02:38.130 --> 00:02:39.150
we didn't talk
about there, which
00:02:39.150 --> 00:02:41.910
is sort of improper integrals,
where the limits of integration
00:02:41.910 --> 00:02:43.440
are unbounded.
00:02:43.440 --> 00:02:46.670
And these come up all the
time in engineering problems.
00:02:46.670 --> 00:02:49.560
We'd like to know what
the, say, integrated
00:02:49.560 --> 00:02:53.576
response of some function
is over all possible times.
00:02:53.576 --> 00:02:55.450
And we may have to
evaluate that numerically.
00:02:55.450 --> 00:02:57.860
So how do you do those
sorts of integrals?
00:02:57.860 --> 00:03:00.090
And it turns out the strategy
that's most often used
00:03:00.090 --> 00:03:03.420
is to divide this domain
up into two pieces.
00:03:03.420 --> 00:03:05.370
One piece, which
is proper integral
00:03:05.370 --> 00:03:07.680
over a finite domain,
and another piece,
00:03:07.680 --> 00:03:09.930
which is an improper integral
over an infinite domain.
00:03:09.930 --> 00:03:11.513
This first integral,
you can handle it
00:03:11.513 --> 00:03:15.120
with an ODE-IVP method or
some higher order polynomial
00:03:15.120 --> 00:03:16.560
interpolation.
00:03:16.560 --> 00:03:19.494
But the second one we
have to do differently.
00:03:19.494 --> 00:03:21.160
And there are couple
of ways to do this.
00:03:21.160 --> 00:03:23.370
One is to transform variables.
00:03:23.370 --> 00:03:27.070
So you map this domain
onto a finite domain.
00:03:27.070 --> 00:03:28.920
And the other option
is to substitute
00:03:28.920 --> 00:03:30.420
an asymptotic approximation.
00:03:30.420 --> 00:03:34.740
So we say, when, say,
tau is large, f of tau
00:03:34.740 --> 00:03:37.410
has a characteristic
asymptotic approximation.
00:03:37.410 --> 00:03:40.200
Maybe it's
exponentially decaying.
00:03:40.200 --> 00:03:43.080
And we take the integral of
that asymptotic approximation.
00:03:43.080 --> 00:03:44.910
And the more accurate
that approximation
00:03:44.910 --> 00:03:48.150
is the better our
approximation for this whole
00:03:48.150 --> 00:03:49.800
integral will be.
00:03:49.800 --> 00:03:52.370
There's another type of
interval that we have to do,
00:03:52.370 --> 00:03:54.960
where the same idea
can be applied.
00:03:54.960 --> 00:03:58.485
That's the integral over
integrable singularities.
00:03:58.485 --> 00:03:59.360
So here's an example.
00:04:02.810 --> 00:04:06.300
So say we want to integrate
the function cosine of tau
00:04:06.300 --> 00:04:10.185
over square root of tau from
0 to some finite positive time
00:04:10.185 --> 00:04:12.800
point.
00:04:12.800 --> 00:04:15.260
As tau goes to 0, cosine
goes to 1 and square root
00:04:15.260 --> 00:04:17.430
of tau diverges.
00:04:17.430 --> 00:04:20.390
But this integral actually
has a finite value.
00:04:20.390 --> 00:04:21.980
This 1 over square
root of tau, that's
00:04:21.980 --> 00:04:25.650
an integrable singularity.
00:04:25.650 --> 00:04:26.670
So how do you do this?
00:04:26.670 --> 00:04:29.730
Well, you want to split the
domain again into two parts.
00:04:29.730 --> 00:04:32.840
One part that contains
the integrable singularity
00:04:32.840 --> 00:04:35.270
and the other part,
which excludes it.
00:04:35.270 --> 00:04:38.210
And in the part that
contains the singularity
00:04:38.210 --> 00:04:41.540
we might do an asymptotic
expansion of the function,
00:04:41.540 --> 00:04:44.510
and integrate each of those
asymptotically accurate terms
00:04:44.510 --> 00:04:45.740
separately.
00:04:45.740 --> 00:04:48.530
So the integral of 1
over square root of tau
00:04:48.530 --> 00:04:51.680
is going to give us
2 tau 0 to the 1/2.
00:04:51.680 --> 00:04:54.230
The integral of
this ratio here is
00:04:54.230 --> 00:04:57.840
going to give us minus
1/2 t0 to the 5/2.
00:04:57.840 --> 00:05:00.600
And then we have an integral
over a finite domain,
00:05:00.600 --> 00:05:02.580
which doesn't include
the singularity.
00:05:02.580 --> 00:05:05.390
And the accuracy of this method
can be dictated by two things
00:05:05.390 --> 00:05:06.140
now.
00:05:06.140 --> 00:05:08.370
One is how accurately
can we do this integral.
00:05:08.370 --> 00:05:13.190
And the other is how accurate is
our asymptotic expansion here.
00:05:13.190 --> 00:05:16.130
So if we want higher
accuracy we need more terms.
00:05:16.130 --> 00:05:17.810
We might need to be
selective about how
00:05:17.810 --> 00:05:20.000
we choose the end
point for this domain
00:05:20.000 --> 00:05:21.489
in order to minimize the error.
00:05:21.489 --> 00:05:23.030
There are lots of
methods to do this.
00:05:23.030 --> 00:05:23.490
Sam.
00:05:23.490 --> 00:05:25.031
AUDIENCE: Will you
talk a little more
00:05:25.031 --> 00:05:28.360
about what typical
integral [INAUDIBLE] is?
00:05:28.360 --> 00:05:31.210
WILLIAM GREEN, JR.: Yeah.
00:05:31.210 --> 00:05:33.850
So you know these things.
00:05:33.850 --> 00:05:40.260
If I try to integrate 1
over x, as x goes to 0
00:05:40.260 --> 00:05:42.310
this thing will always diverge.
00:05:42.310 --> 00:05:46.180
You know the antiderivative of
1 over x is going to be the log.
00:05:46.180 --> 00:05:47.890
And the log diverges.
00:05:47.890 --> 00:05:53.950
So power is weaker than 1 over
x, like 1 over x to the 1/2.
00:05:53.950 --> 00:05:55.290
Don't diverge like the log.
00:05:55.290 --> 00:05:56.830
They actually go
to a finite value.
00:05:56.830 --> 00:05:59.770
The log is like the most weakly
singular function there is.
00:05:59.770 --> 00:06:02.800
And integrals over weaker
powers of 1 over x actually
00:06:02.800 --> 00:06:04.994
don't produce a
singularity anymore.
00:06:04.994 --> 00:06:06.910
They're what we call
integrable singularities.
00:06:06.910 --> 00:06:11.200
The function itself diverges
but its integral does not.
00:06:11.200 --> 00:06:15.100
So you usually have this
asymptotic power law type
00:06:15.100 --> 00:06:16.270
behavior.
00:06:16.270 --> 00:06:18.140
OK.
00:06:18.140 --> 00:06:20.100
So they come up all the time.
00:06:20.100 --> 00:06:23.910
We see these things
in places like, say,
00:06:23.910 --> 00:06:27.460
squeezing a thin film of
fluid between two plates.
00:06:27.460 --> 00:06:31.567
How long does it take for these
two plates to come together?
00:06:31.567 --> 00:06:33.900
You might say, well, as the
plates get closer and closer
00:06:33.900 --> 00:06:36.737
together the pressure
starts to diverge in the gap
00:06:36.737 --> 00:06:38.070
and they'll never come together.
00:06:38.070 --> 00:06:40.950
But maybe depending on
the geometry of the plates
00:06:40.950 --> 00:06:42.750
there can be a
sort of finite time
00:06:42.750 --> 00:06:44.190
at which they'll come together.
00:06:44.190 --> 00:06:45.273
Depends on their geometry.
00:06:48.557 --> 00:06:49.390
So that's one topic.
00:06:49.390 --> 00:06:51.473
I think that's useful
because these things come up
00:06:51.473 --> 00:06:52.030
all the time.
00:06:52.030 --> 00:06:55.330
And you can't actually handle
those cases with quadrature
00:06:55.330 --> 00:06:56.800
by itself.
00:06:56.800 --> 00:06:57.910
The function diverges.
00:06:57.910 --> 00:06:59.680
There's no polynomial
interpolation
00:06:59.680 --> 00:07:01.150
that's suitable to match it.
00:07:01.150 --> 00:07:02.777
So the quadrature
can't handle it.
00:07:02.777 --> 00:07:04.360
And if you're trying
to integrate over
00:07:04.360 --> 00:07:07.210
an infinite domain,
well, good luck.
00:07:07.210 --> 00:07:09.580
You're never going to be
able to fit enough points in
00:07:09.580 --> 00:07:14.440
to know that you've accurately
integrated out to infinity.
00:07:14.440 --> 00:07:17.260
The other thing I
want to recap here,
00:07:17.260 --> 00:07:19.750
and this is a problem
for you to try.
00:07:19.750 --> 00:07:23.890
So here we have a simple first
order ordinary differential
00:07:23.890 --> 00:07:24.390
equation.
00:07:24.390 --> 00:07:27.080
We want to use the implicit
Euler method to solve it.
00:07:27.080 --> 00:07:29.840
So can you give a
closed form solution
00:07:29.840 --> 00:07:33.260
for one step of the
implicit Euler method,
00:07:33.260 --> 00:07:36.620
or for n steps of the
implicit Euler method?
00:07:36.620 --> 00:07:38.300
What result does that produce?
00:07:38.300 --> 00:07:40.520
Can you compute this?
00:07:40.520 --> 00:07:48.556
Do you guys talk about backward
Euler methods for for ODE-IVPs?
00:07:48.556 --> 00:07:51.180
Backwards difference methods in
general and the backwards Euler
00:07:51.180 --> 00:07:54.600
method specifically?
00:07:54.600 --> 00:07:55.940
No.
00:07:55.940 --> 00:07:57.150
OK.
00:07:57.150 --> 00:08:01.770
Well, clearly I
should have been here.
00:08:01.770 --> 00:08:03.150
So what's the strategy?
00:08:03.150 --> 00:08:05.100
We have to approximate
the derivative.
00:08:05.100 --> 00:08:06.960
Yes?
00:08:06.960 --> 00:08:08.880
So the approximation for
the derivative we use
00:08:08.880 --> 00:08:11.160
is a finite difference
approximation.
00:08:11.160 --> 00:08:14.940
So we write dx dt
as x at a point k
00:08:14.940 --> 00:08:18.870
plus 1 minus x at the
point k divided by delta t.
00:08:18.870 --> 00:08:20.700
Backwards difference.
00:08:20.700 --> 00:08:24.930
We evaluate the right-hand
side of the ODE-IVP at k
00:08:24.930 --> 00:08:27.940
plus 1, xk plus 1.
00:08:27.940 --> 00:08:30.360
And then we solve for xk plus 1.
00:08:30.360 --> 00:08:32.159
We just substitute in here.
00:08:32.159 --> 00:08:33.690
Backwards difference
approximation
00:08:33.690 --> 00:08:35.320
for the differential.
00:08:35.320 --> 00:08:37.979
And xk plus 1 for
the right-hand side.
00:08:37.979 --> 00:08:41.010
And we find that xk
plus 1 will be 1 over 1
00:08:41.010 --> 00:08:46.950
minus delta t times
this lambda times xk.
00:08:46.950 --> 00:08:51.680
And if we iterate this from k
equals 0 up to some finite k,
00:08:51.680 --> 00:08:54.240
it's like multiplying
by powers of 1 over 1
00:08:54.240 --> 00:08:55.980
minus delta t times lambda.
00:08:55.980 --> 00:08:57.044
Yes?
00:08:57.044 --> 00:08:57.960
Does this makes sense?
00:09:01.240 --> 00:09:05.910
So our backwards Euler solution
to this fundamental problem
00:09:05.910 --> 00:09:08.386
is going to look like this.
00:09:08.386 --> 00:09:10.260
And a lot of times we
ask about the stability
00:09:10.260 --> 00:09:12.240
of these sorts of solutions.
00:09:12.240 --> 00:09:16.850
So if for example, the
solution to this equation
00:09:16.850 --> 00:09:20.000
was supposed to
decay to zero, we
00:09:20.000 --> 00:09:21.770
would expect our
numerical method
00:09:21.770 --> 00:09:23.890
to also yield results
which decay to zero.
00:09:23.890 --> 00:09:27.110
We would call that stable.
00:09:27.110 --> 00:09:29.660
So we all know under
what conditions
00:09:29.660 --> 00:09:34.662
does xk, for example, go to 0
as I change the product delta
00:09:34.662 --> 00:09:36.350
t times lambda.
00:09:36.350 --> 00:09:40.850
Lambda could be a real number or
it can be a complex number too.
00:09:40.850 --> 00:09:42.360
Doesn't really matter.
00:09:42.360 --> 00:09:46.830
So what needs to be true for xk
to the k to 0 as k gets bigger?
00:09:46.830 --> 00:09:49.340
Well what needs to be true is
this quantity in parentheses
00:09:49.340 --> 00:09:50.960
needs to be smaller
than 1 because I'm
00:09:50.960 --> 00:09:52.400
taking lots of
products of things
00:09:52.400 --> 00:09:53.483
that are smaller than one.
00:09:53.483 --> 00:09:57.770
That will continue to shrink xk
from x 0 until it goes to zero.
00:09:57.770 --> 00:10:01.790
If this is smaller than 1,
then this quantity down here
00:10:01.790 --> 00:10:04.610
needs to be bigger than
1 in absolute value.
00:10:04.610 --> 00:10:07.010
And then I'll get solutions
that slowly decay to zero.
00:10:07.010 --> 00:10:10.490
So if I ask, for what values
of delta t times lambda
00:10:10.490 --> 00:10:13.770
is the absolute value of
this thing bigger than 1,
00:10:13.770 --> 00:10:17.245
and I allow lambda to be a
complex number, for example?
00:10:17.245 --> 00:10:20.110
I'll find out that's
true when 1 minus delta t
00:10:20.110 --> 00:10:22.490
times the real part
of lambda squared
00:10:22.490 --> 00:10:26.630
plus delta t times the
imaginary part of lambda squared
00:10:26.630 --> 00:10:27.530
is bigger than 1.
00:10:30.100 --> 00:10:33.440
And if I plot the real
part of delta t times
00:10:33.440 --> 00:10:37.630
lambda versus the imaginary
part of delta times lambda,
00:10:37.630 --> 00:10:41.160
this inequality is
represented by this pink zone
00:10:41.160 --> 00:10:42.470
on the outside of the circle.
00:10:42.470 --> 00:10:44.510
So any values of
delta t times lambda
00:10:44.510 --> 00:10:49.460
that live in this pink area
will give stable integration.
00:10:49.460 --> 00:10:54.500
And the solution xk, well, the
k to zero as k goes to infinity.
00:10:57.020 --> 00:10:59.530
So it makes sense?
00:10:59.530 --> 00:11:02.010
This is a little funny though.
00:11:02.010 --> 00:11:03.780
What does the solution
to this equation
00:11:03.780 --> 00:11:08.012
do when lambda is bigger than 1?
00:11:08.012 --> 00:11:09.470
When it's bigger
than 0, excuse me.
00:11:09.470 --> 00:11:10.886
When lambda is
bigger than 0, what
00:11:10.886 --> 00:11:14.860
does the solution
to this ODE-IVP do?
00:11:14.860 --> 00:11:15.560
Does it decay?
00:11:15.560 --> 00:11:16.640
AUDIENCE: [INAUDIBLE]
00:11:16.640 --> 00:11:17.973
WILLIAM GREEN, JR.: It blows up.
00:11:17.973 --> 00:11:19.100
It grows exponentially.
00:11:19.100 --> 00:11:21.920
But the numerical
method here, when
00:11:21.920 --> 00:11:25.070
I have lambda is
bigger than zero,
00:11:25.070 --> 00:11:28.659
can be stable and produce
decaying solutions.
00:11:28.659 --> 00:11:29.825
So this is sort of peculiar.
00:11:32.890 --> 00:11:35.290
So we're often concerned
when we solve ODE-IVPs
00:11:35.290 --> 00:11:37.640
and we try to use these
sorts of backwards difference
00:11:37.640 --> 00:11:40.880
methods, these implicit
methods, to solve them.
00:11:40.880 --> 00:11:44.390
We're often concerned
with getting stability.
00:11:44.390 --> 00:11:46.299
We want to get
solutions that decay
00:11:46.299 --> 00:11:47.590
when they're supposed to decay.
00:11:47.590 --> 00:11:49.340
Sometimes we get
solutions that decay even
00:11:49.340 --> 00:11:52.690
when the solution is
supposed to blow up.
00:11:52.690 --> 00:11:55.240
So that's sort of funny.
00:11:55.240 --> 00:11:57.835
So the exact solution
should look like this.
00:11:57.835 --> 00:11:59.710
But there are going to
be circumstances where
00:11:59.710 --> 00:12:00.940
lambda is bigger than zero.
00:12:00.940 --> 00:12:02.530
In fact, lambda is bigger than--
00:12:02.530 --> 00:12:06.340
lambda delta t is bigger
than 2 along the real axis,
00:12:06.340 --> 00:12:11.060
where the solution will actually
blow up rather than decay away.
00:12:11.060 --> 00:12:12.970
So stability and
accuracy don't always
00:12:12.970 --> 00:12:15.760
correlate with each other.
00:12:15.760 --> 00:12:18.760
We really have to understand
the underlying dynamics
00:12:18.760 --> 00:12:21.040
of the problem we're
trying to solve
00:12:21.040 --> 00:12:24.110
before we choose a numerical
method to apply to it.
00:12:24.110 --> 00:12:27.280
And this is going to be true
of differential algebraic
00:12:27.280 --> 00:12:29.740
equations as well.
00:12:29.740 --> 00:12:33.310
Are there any
questions about this?
00:12:33.310 --> 00:12:37.030
There's a discussion of these
things in your textbook,
00:12:37.030 --> 00:12:39.160
these sorts of stability
regions associated
00:12:39.160 --> 00:12:41.290
with different backwards
difference methods.
00:12:41.290 --> 00:12:46.450
The implicit Euler method is one
backwards difference formula.
00:12:46.450 --> 00:12:48.760
It's the simplest one.
00:12:48.760 --> 00:12:52.160
It's the least accurate
one that you can derive
00:12:52.160 --> 00:12:53.420
but it's an easy one to use.
00:12:56.074 --> 00:12:58.490
But now we want to move on to
a different sort of problem.
00:12:58.490 --> 00:13:01.737
And it's one that we come across
in engineering all the time.
00:13:01.737 --> 00:13:03.320
These problems are
called differential
00:13:03.320 --> 00:13:04.670
algebraic equations.
00:13:04.670 --> 00:13:08.270
And they're problems of the
sort a function of some state
00:13:08.270 --> 00:13:12.500
variables, and the derivatives
of those state variables
00:13:12.500 --> 00:13:15.147
and time is equal to zero.
00:13:15.147 --> 00:13:16.730
And there are some
initial conditions.
00:13:16.730 --> 00:13:19.825
The state variables at that
time 0 is equal to some x0.
00:13:24.840 --> 00:13:26.490
So this is a vector
valued function
00:13:26.490 --> 00:13:29.320
of a vector valued state
and it's time derivatives.
00:13:29.320 --> 00:13:31.650
And we want to understand
the dynamics of this system.
00:13:31.650 --> 00:13:33.660
How does x vary
with time in a way
00:13:33.660 --> 00:13:37.500
that's consistent with
this governing equation?
00:13:37.500 --> 00:13:39.840
Usually we call these
well-posed problems
00:13:39.840 --> 00:13:42.690
when the dimension of the
state variable and the function
00:13:42.690 --> 00:13:43.590
are the same.
00:13:43.590 --> 00:13:46.730
They both live in
the vector space rn.
00:13:46.730 --> 00:13:48.375
So we have n different states.
00:13:48.375 --> 00:13:50.550
I have n different functions
that I have to satisfy.
00:13:50.550 --> 00:13:55.090
And I want to understand
how x varies with time.
00:13:55.090 --> 00:13:57.320
Here's an example.
00:13:57.320 --> 00:13:58.200
A stirred tank.
00:13:58.200 --> 00:14:00.460
I call this stirred
tank, example one.
00:14:00.460 --> 00:14:06.684
So into the stirred tank comes
some concentration of solute.
00:14:06.684 --> 00:14:09.100
Out of the stirred tank comes
some different concentration
00:14:09.100 --> 00:14:09.891
of the same solute.
00:14:09.891 --> 00:14:14.120
And we want to model the
dynamics as a function of time.
00:14:14.120 --> 00:14:16.900
So we know that a
material balance
00:14:16.900 --> 00:14:19.390
of the solute on the stirred
tank, if it's well mixed,
00:14:19.390 --> 00:14:26.320
will tell us that d c2 dt is
related to the volumetric flow
00:14:26.320 --> 00:14:28.960
rate into the tank divided
by the volume of the tank
00:14:28.960 --> 00:14:32.140
multiplied by the difference
between c1 and c2,
00:14:32.140 --> 00:14:37.112
the amount carried in and the
amount carried out of the tank.
00:14:37.112 --> 00:14:38.890
And let's put a
constraint in here too.
00:14:38.890 --> 00:14:42.430
Let's suppose we're trying to
control this tank in some way.
00:14:42.430 --> 00:14:45.820
So we say that c1
of time has got
00:14:45.820 --> 00:14:49.190
to be equal to some
control signal, gamma of t.
00:14:52.390 --> 00:14:55.810
And then we want to understand
the dynamics of this system.
00:14:55.810 --> 00:14:59.710
At time 0, c2 is some c0,
the initial concentration
00:14:59.710 --> 00:15:00.570
of the tank.
00:15:00.570 --> 00:15:03.490
And c1 is whatever the initial
value of this control signal
00:15:03.490 --> 00:15:04.960
is.
00:15:04.960 --> 00:15:08.350
And the solution is
c1 equals gamma of t.
00:15:08.350 --> 00:15:11.470
And c2 is, well, the
initial concentration's
00:15:11.470 --> 00:15:13.090
got to decay out exponentially.
00:15:13.090 --> 00:15:15.220
And then I've got to
integrate over how
00:15:15.220 --> 00:15:17.770
the input is changing in time.
00:15:17.770 --> 00:15:20.060
Those also produce
exponential decay.
00:15:20.060 --> 00:15:22.870
So if the input is a
little delta function spike
00:15:22.870 --> 00:15:25.330
then I'll get an exponential
decay leaving the tank.
00:15:25.330 --> 00:15:27.310
If the input is changing
dynamically in time,
00:15:27.310 --> 00:15:30.550
I need this convolution of
the input with this integral.
00:15:30.550 --> 00:15:33.130
So you just to
solve this system,
00:15:33.130 --> 00:15:34.900
this ordinary
differential equation
00:15:34.900 --> 00:15:37.310
by substituting c1 in there.
00:15:37.310 --> 00:15:40.280
But this system of
equations here is not just
00:15:40.280 --> 00:15:42.040
a differential equation.
00:15:42.040 --> 00:15:45.847
It's a differential equation
and an algebraic equation.
00:15:45.847 --> 00:15:47.680
There's something
peculiar about that that's
00:15:47.680 --> 00:15:50.140
different than just
solving systems
00:15:50.140 --> 00:15:52.440
of ordinary
differential equations.
00:15:52.440 --> 00:15:57.400
Oftentimes in models,
we write down all
00:15:57.400 --> 00:15:59.720
the governing equations
associated with the model.
00:15:59.720 --> 00:16:01.636
Some of them may be
differential, some of them
00:16:01.636 --> 00:16:02.800
may be algebraic.
00:16:02.800 --> 00:16:04.960
There are times when we
can look at those equations
00:16:04.960 --> 00:16:07.480
and see clever substitutions
that we can make.
00:16:07.480 --> 00:16:10.840
But if you have a hundred
equations, or a thousand
00:16:10.840 --> 00:16:12.900
equations, or a
million equations,
00:16:12.900 --> 00:16:14.820
there's no way to
do that reliably.
00:16:14.820 --> 00:16:18.310
We just kind of wind up with a
system of ordinary differential
00:16:18.310 --> 00:16:21.730
equations and algebraic
equations mixed together.
00:16:21.730 --> 00:16:24.777
And we have to solve
these reliably.
00:16:24.777 --> 00:16:25.735
Here's another example.
00:16:28.650 --> 00:16:29.720
Oh, excuse me.
00:16:29.720 --> 00:16:31.800
So let me write
this in this form.
00:16:31.800 --> 00:16:34.110
So we said there should
be a vector valued
00:16:34.110 --> 00:16:38.460
function of the states and
the derivatives of the states.
00:16:38.460 --> 00:16:41.490
The states are the
concentration c1 and c2.
00:16:41.490 --> 00:16:45.090
And the derivative
here is dc 2 dt.
00:16:45.090 --> 00:16:46.720
So here's my vector
valued function.
00:16:46.720 --> 00:16:48.900
The first element is
this, and it's equal to 0.
00:16:48.900 --> 00:16:51.660
The second element is
this, and it's equal to 0.
00:16:51.660 --> 00:16:54.900
Here's my initial
condition over the states.
00:16:54.900 --> 00:16:57.930
So I can just transform
that list of equations
00:16:57.930 --> 00:17:00.450
to the sort of canonical
form for a differential
00:17:00.450 --> 00:17:01.420
algebraic equation.
00:17:04.690 --> 00:17:07.619
Here's a separate example.
00:17:07.619 --> 00:17:12.790
Suppose instead, I'm
trying to either control
00:17:12.790 --> 00:17:16.450
or make a measurement of the
outlet concentrations c2.
00:17:16.450 --> 00:17:18.220
And we call that gamma of t.
00:17:22.060 --> 00:17:26.329
And I want to know now
what is c2 and what is c1?
00:17:26.329 --> 00:17:29.340
So I've got to solve this system
of equations for c2 and c1.
00:17:29.340 --> 00:17:31.740
The solution for c2 is easy.
00:17:31.740 --> 00:17:32.360
I know that.
00:17:32.360 --> 00:17:35.060
c2 is gamma, by definition.
00:17:35.060 --> 00:17:37.430
I got to substitute this
into the first equation
00:17:37.430 --> 00:17:38.420
and solve for c1.
00:17:38.420 --> 00:17:43.199
And I see, c1 is gamma
plus v over q gamma dot.
00:17:43.199 --> 00:17:45.240
And I have to have initial
conditions to go along
00:17:45.240 --> 00:17:47.992
with this.
00:17:47.992 --> 00:17:49.200
There's something funny here.
00:17:49.200 --> 00:17:52.800
We can see the initial condition
for c2 has to be gamma 0.
00:17:52.800 --> 00:17:57.710
The initial condition for c1,
what does that have to be?
00:17:57.710 --> 00:18:01.430
Well, it needs to be consistent
with the solution here.
00:18:01.430 --> 00:18:03.800
There were no free
variables when
00:18:03.800 --> 00:18:05.484
I solved this equation for c1.
00:18:05.484 --> 00:18:07.400
It isn't like solving a
differential equation,
00:18:07.400 --> 00:18:10.220
having a free
coefficient to specify.
00:18:10.220 --> 00:18:13.130
So it better be the
case that this c0 here
00:18:13.130 --> 00:18:18.680
is the same as gamma 0
plus v over q gamma dot 0.
00:18:18.680 --> 00:18:21.410
Somehow I have to know this
input signal to prescribe
00:18:21.410 --> 00:18:23.030
the initial conditions for c1.
00:18:23.030 --> 00:18:28.805
So that's peculiar and
different from just ODE-IVPs.
00:18:28.805 --> 00:18:29.680
You see this picture?
00:18:29.680 --> 00:18:32.940
You see how this goes together?
00:18:32.940 --> 00:18:34.630
The solution is funny too.
00:18:34.630 --> 00:18:36.880
Suppose we are trying to
use this system of equations
00:18:36.880 --> 00:18:39.040
to do the following.
00:18:39.040 --> 00:18:43.720
We are trying to measure c2
and use the measurement of c2
00:18:43.720 --> 00:18:46.030
to predict c1.
00:18:46.030 --> 00:18:48.580
So gamma is the
measurement, and we're
00:18:48.580 --> 00:18:52.060
trying to solve this system
of differential algebraic
00:18:52.060 --> 00:18:55.120
equations to predict what c1 is.
00:18:55.120 --> 00:18:58.300
All measurements
incur numerical error.
00:18:58.300 --> 00:19:01.560
So even though our signal
for gamma that we measure
00:19:01.560 --> 00:19:06.670
may be continuous, it's
bouncing around wildly.
00:19:06.670 --> 00:19:08.290
It's not a constant signal.
00:19:08.290 --> 00:19:10.900
It's going to move
around a great deal.
00:19:10.900 --> 00:19:13.960
And c1, our predicted
value for c1,
00:19:13.960 --> 00:19:15.670
depends on the
derivatives of gamma,
00:19:15.670 --> 00:19:18.490
which means c1 is going to be
incredibly sensitive to how
00:19:18.490 --> 00:19:21.440
that measurement is made.
00:19:21.440 --> 00:19:22.395
So it's peculiar.
00:19:22.395 --> 00:19:24.520
It means that there's going
to be a lot of problems
00:19:24.520 --> 00:19:28.840
with stably solving these
equations because the solutions
00:19:28.840 --> 00:19:34.390
can admit cases where they're
not integrals of the input
00:19:34.390 --> 00:19:37.990
but they're actually related
to derivatives of the input.
00:19:37.990 --> 00:19:39.960
Look at the solution
back here again.
00:19:39.960 --> 00:19:43.140
Oop, sorry.
00:19:43.140 --> 00:19:45.960
The solution here was related
to an integral of the input.
00:19:45.960 --> 00:19:48.030
Suppose I was making
a measurement of c1
00:19:48.030 --> 00:19:50.940
in trying to predict c2 instead?
00:19:50.940 --> 00:19:53.480
My prediction for c2 would
be related to the integral
00:19:53.480 --> 00:19:54.480
of that measured signal.
00:19:54.480 --> 00:19:56.850
And we know integrals
smooth signals out.
00:19:56.850 --> 00:19:59.340
So it's not going to be so
sensitive to the measurement.
00:19:59.340 --> 00:20:01.740
So one way of doing this
seems really sensible.
00:20:01.740 --> 00:20:06.540
And the other way of doing
it seems a little bizarre.
00:20:06.540 --> 00:20:09.317
Well, we can construct these
equations however we want.
00:20:09.317 --> 00:20:11.400
So there's something about
how we formulate models
00:20:11.400 --> 00:20:13.140
that are going to make
them either easily
00:20:13.140 --> 00:20:18.850
solvable, as DAEs, or
quite difficult to solve.
00:20:18.850 --> 00:20:23.190
So these pop up all the time
in engineering problems.
00:20:23.190 --> 00:20:25.530
They pop up in
dynamic simulations,
00:20:25.530 --> 00:20:27.900
where we have some sort
of underlying conservation
00:20:27.900 --> 00:20:31.950
principle, or constraints,
or equilibria.
00:20:31.950 --> 00:20:33.450
So if we have to
conserve something,
00:20:33.450 --> 00:20:38.550
like total energy, total mass,
total momentum, total particle
00:20:38.550 --> 00:20:41.760
number, total atomic
species, total charge,
00:20:41.760 --> 00:20:44.040
there will always be an
algebraic constraint that
00:20:44.040 --> 00:20:48.210
tells us the total
mass has to stay this,
00:20:48.210 --> 00:20:51.500
or the total number of atoms of
a certain type has to say this.
00:20:51.500 --> 00:20:52.230
Yeah, Ian.
00:20:52.230 --> 00:20:56.950
AUDIENCE: Just to be clear, were
you saying, on the stirred tank
00:20:56.950 --> 00:20:59.722
example, that those
were the same problem
00:20:59.722 --> 00:21:02.119
with different
solution approaches?
00:21:02.119 --> 00:21:03.660
WILLIAM GREEN, JR.:
So good question.
00:21:03.660 --> 00:21:04.285
Let me go back.
00:21:07.260 --> 00:21:10.960
So they are fundamentally
different problems.
00:21:10.960 --> 00:21:17.450
So in one case I specified the
value of c2 to be this gamma.
00:21:17.450 --> 00:21:21.450
I tried to make a measurement
of c2 and then predict c1.
00:21:21.450 --> 00:21:24.180
In the other case,
I specified c1.
00:21:24.180 --> 00:21:28.100
I tried to measure
c1 and predict c2.
00:21:28.100 --> 00:21:29.850
So one case I measured
the input and tried
00:21:29.850 --> 00:21:31.869
to predict the
output to the tank.
00:21:31.869 --> 00:21:33.660
And in the other case
I measured the output
00:21:33.660 --> 00:21:36.600
and tried to use it
to predict the input.
00:21:36.600 --> 00:21:37.273
Yeah?
00:21:37.273 --> 00:21:39.125
AUDIENCE: So physically
the same system
00:21:39.125 --> 00:21:40.790
with a choice by the operator.
00:21:40.790 --> 00:21:41.790
WILLIAM GREEN, JR.: Yes.
00:21:41.790 --> 00:21:42.569
That's right.
00:21:42.569 --> 00:21:43.110
That's right.
00:21:43.110 --> 00:21:45.270
So I formulated the
model differently.
00:21:45.270 --> 00:21:49.080
The same underlying
system but I formulated
00:21:49.080 --> 00:21:50.700
the model differently.
00:21:50.700 --> 00:21:52.950
And it led to completely
different behavior
00:21:52.950 --> 00:21:55.770
in the underlying solution
to the differential algebraic
00:21:55.770 --> 00:21:57.520
equations.
00:21:57.520 --> 00:21:59.160
So you've already
solved problems
00:21:59.160 --> 00:22:00.360
that involve conservation.
00:22:00.360 --> 00:22:04.770
You've done things like
ODE-IVPs on batch reactors,
00:22:04.770 --> 00:22:07.320
where the total amount of
material in the reactor
00:22:07.320 --> 00:22:09.015
had to remain constant.
00:22:09.015 --> 00:22:10.890
Instead, you probably
solved for the dynamics
00:22:10.890 --> 00:22:14.130
of each of the components
undergoing the reaction.
00:22:14.130 --> 00:22:18.540
You might have checked to
see if at each point in time
00:22:18.540 --> 00:22:21.200
the total mass in the
reactor stayed the same.
00:22:21.200 --> 00:22:23.754
Because you solved all these
differential equations,
00:22:23.754 --> 00:22:25.420
they were interconnected
with each other
00:22:25.420 --> 00:22:27.450
but they all incurred
some numerical error.
00:22:27.450 --> 00:22:29.370
It's possible that
numerical error
00:22:29.370 --> 00:22:33.150
accrues in a way that
actually has you losing mass
00:22:33.150 --> 00:22:35.430
from your system.
00:22:35.430 --> 00:22:37.230
So there may be benefits
to actually trying
00:22:37.230 --> 00:22:39.360
to solve these sorts
of systems of equations
00:22:39.360 --> 00:22:42.150
with, say, one less
differential equation
00:22:42.150 --> 00:22:45.390
for one of the components and
instead an algebraic constraint
00:22:45.390 --> 00:22:48.330
governing the total mass
or total number of moles
00:22:48.330 --> 00:22:49.050
in the system.
00:22:53.010 --> 00:22:55.320
So these pop up all the time.
00:22:55.320 --> 00:22:58.080
You see them in models
of reaction networks,
00:22:58.080 --> 00:23:00.540
where we try to use a pseudo
steady state approximation.
00:23:00.540 --> 00:23:03.960
We say some reactions go so
fast that they equilibrate.
00:23:03.960 --> 00:23:06.120
So they're not governed
by differential equations
00:23:06.120 --> 00:23:11.790
but by algebraic
equations for equilibria.
00:23:11.790 --> 00:23:14.520
They can pop up in models
of control, where we neglect
00:23:14.520 --> 00:23:16.680
the controller dynamics.
00:23:16.680 --> 00:23:19.440
So you say I'm going to
try to control this process
00:23:19.440 --> 00:23:23.010
and I get to instantaneously
turn on or off this valve
00:23:23.010 --> 00:23:24.930
in response to a measurement.
00:23:24.930 --> 00:23:28.472
Actually, controllers have
inherent dynamics in them.
00:23:28.472 --> 00:23:30.180
But if I don't put
those dynamics in then
00:23:30.180 --> 00:23:32.760
I may get an algebraic equation
instead of a differential
00:23:32.760 --> 00:23:34.755
equation for how the
control process occurs.
00:23:39.690 --> 00:23:44.250
There are different types
of DAEs that we talk about.
00:23:44.250 --> 00:23:48.640
So there's one type that's
called semi-explicit DAEs.
00:23:48.640 --> 00:23:50.980
And in these types of
differential algebraic
00:23:50.980 --> 00:23:53.710
equations we can
write them in the form
00:23:53.710 --> 00:24:00.310
M dx dt is equal to f of x
and t, some initial condition.
00:24:00.310 --> 00:24:03.410
You say that this
looks like an ODE-IVP.
00:24:03.410 --> 00:24:06.430
M is called the mass matrix.
00:24:06.430 --> 00:24:10.990
And it may or may not be the
case that M is invertible.
00:24:10.990 --> 00:24:16.150
M may or may not be
a full rank matrix.
00:24:16.150 --> 00:24:20.600
So let's look at that
stirred tank example one.
00:24:20.600 --> 00:24:24.890
This was the underlying
equation, f of x and dx dt
00:24:24.890 --> 00:24:25.790
and t.
00:24:25.790 --> 00:24:27.860
Can write it in
semi-explicit form.
00:24:27.860 --> 00:24:35.590
If c has two components, c1
and c2, then dc 2 dt-- oops,
00:24:35.590 --> 00:24:36.640
this is a typo here.
00:24:36.640 --> 00:24:37.810
I really apologize.
00:24:37.810 --> 00:24:42.180
This should be a 0, and
this element should be a 1.
00:24:42.180 --> 00:24:43.230
I'll correct that online.
00:24:43.230 --> 00:24:46.200
I apologize for that.
00:24:46.200 --> 00:24:48.709
So this should be
this matrix multiplied
00:24:48.709 --> 00:24:50.250
with this differential
should give me
00:24:50.250 --> 00:24:52.750
dc 2 dt, this term here.
00:24:52.750 --> 00:24:56.410
And a 0 for the second equation.
00:24:56.410 --> 00:25:02.450
And then I've got q over v
multiplying c1 and minus c2.
00:25:02.450 --> 00:25:05.010
So that gives me
the first line here.
00:25:05.010 --> 00:25:09.850
I've got minus c1 coming on
the second equation here.
00:25:09.850 --> 00:25:12.820
And those balance
with a 0 and a gamma.
00:25:12.820 --> 00:25:16.000
So the lower line reads
like this equation here
00:25:16.000 --> 00:25:17.650
and the upper line--
00:25:17.650 --> 00:25:19.660
fix this typo--
replace 0 with one,
00:25:19.660 --> 00:25:21.790
reads like the
upper equation here.
00:25:21.790 --> 00:25:24.950
This is the semi-explicit form
of the differential equation.
00:25:24.950 --> 00:25:27.970
You can see this mass
matrix is singular and has
00:25:27.970 --> 00:25:29.200
a row that's all zeros.
00:25:29.200 --> 00:25:31.000
There's no way to
invert this matrix
00:25:31.000 --> 00:25:34.540
and convert it into a system
of ODE-IVPs, which you already
00:25:34.540 --> 00:25:37.360
know how to solve.
00:25:37.360 --> 00:25:41.080
If the mass matrix is full
rank, which can happen,
00:25:41.080 --> 00:25:43.570
you can formulate your
model in such a way
00:25:43.570 --> 00:25:45.980
that it takes on this
semi-explicit form.
00:25:45.980 --> 00:25:48.530
But if the mass matrix is
full rank and invertible,
00:25:48.530 --> 00:25:50.730
then you just invert
the mass matrix.
00:25:50.730 --> 00:25:53.260
And now you can apply all
your ODE-IVP techniques
00:25:53.260 --> 00:25:56.710
to solve this thing.
00:25:56.710 --> 00:25:58.960
The mass matrix doesn't
need to be constant.
00:25:58.960 --> 00:26:01.720
Could depend on
the concentration.
00:26:01.720 --> 00:26:06.397
It can depend on the states
as well in this form.
00:26:06.397 --> 00:26:08.480
I just chose an example
where that isn't the case.
00:26:08.480 --> 00:26:10.345
But in general, it can
depend on the states.
00:26:10.345 --> 00:26:12.662
If it depends on the
states and we invert it,
00:26:12.662 --> 00:26:14.620
then we need to be sure
that the mass matrix is
00:26:14.620 --> 00:26:18.100
invertible for all values of
the states, which may or may not
00:26:18.100 --> 00:26:19.395
be easy to ensure.
00:26:24.700 --> 00:26:28.010
So here's what I'd
like you to try to do.
00:26:28.010 --> 00:26:30.850
So here's that stirred
tank example two.
00:26:30.850 --> 00:26:34.240
The only difference between
this and the previous one
00:26:34.240 --> 00:26:38.130
is that c2 now is set equal
to gamma instead of c1.
00:26:38.130 --> 00:26:40.430
I want you to try to write
it in semi-explicit form.
00:26:40.430 --> 00:26:42.790
A sort of test of understanding.
00:26:42.790 --> 00:26:46.700
Can you define the mass matrix?
00:26:46.700 --> 00:26:48.380
Can you write out
the right-hand side
00:26:48.380 --> 00:26:50.175
of the equation in
semi-explicit form?
00:26:50.175 --> 00:26:51.800
Take a second, work
with your neighbor.
00:26:51.800 --> 00:26:53.508
See if you can figure
out how to do that.
00:28:04.658 --> 00:28:06.020
AUDIENCE: Is this clear?
00:28:06.020 --> 00:28:07.296
OK to do this?
00:28:12.920 --> 00:28:15.420
WILLIAM GREEN, JR.: You can
look back on the previous slide.
00:28:15.420 --> 00:28:18.690
Actually, the semi-exclusive
form is going to be the same.
00:28:18.690 --> 00:28:20.910
This is a 0, that's a 1.
00:28:20.910 --> 00:28:23.850
The only difference is
going to be c2 is now
00:28:23.850 --> 00:28:25.970
balanced with gamma, the input.
00:28:25.970 --> 00:28:27.880
So this minus 1 is
got to shift over.
00:28:27.880 --> 00:28:29.460
So I switched the
0 for the minus 1,
00:28:29.460 --> 00:28:32.040
you got the second,
the semi-explicit form.
00:28:38.270 --> 00:28:43.140
There's another way that
these equations can pop up.
00:28:43.140 --> 00:28:44.880
So we have this mass matrix.
00:28:44.880 --> 00:28:47.790
And it might be
the case that M is
00:28:47.790 --> 00:28:51.450
diagonal over many of
the states and then
00:28:51.450 --> 00:28:52.789
0 for all the other states.
00:28:52.789 --> 00:28:55.080
So we might be able to write
the equation in this form.
00:28:55.080 --> 00:28:56.204
This comes up all the time.
00:28:56.204 --> 00:28:58.260
So we might be able
to write dx dt is
00:28:58.260 --> 00:29:01.110
a function of x, y, and t.
00:29:01.110 --> 00:29:06.950
And 0 is equal to another
function of x, y, and t.
00:29:06.950 --> 00:29:09.300
And we have some
initial conditions
00:29:09.300 --> 00:29:13.410
prescribed for the
x's, of which we're
00:29:13.410 --> 00:29:16.560
taking the differential
in one of these equations.
00:29:16.560 --> 00:29:19.890
And we've also got to satisfy
those initial conditions
00:29:19.890 --> 00:29:23.400
for these variables y associated
with the second nonlinear
00:29:23.400 --> 00:29:24.390
equation.
00:29:24.390 --> 00:29:27.750
The x's are called the
differential states
00:29:27.750 --> 00:29:30.144
because they're involved
in equations where they're
00:29:30.144 --> 00:29:32.310
difference with respect to
time, or the differential
00:29:32.310 --> 00:29:34.260
with respect to time is taken.
00:29:34.260 --> 00:29:37.740
y are called the
algebraic states
00:29:37.740 --> 00:29:41.280
because they only appear
as themselves and not
00:29:41.280 --> 00:29:44.470
as their time differential
in any of these equations.
00:29:44.470 --> 00:29:47.590
So we might want to solve
for both the differential
00:29:47.590 --> 00:29:50.040
and the algebraic
states simultaneously.
00:29:50.040 --> 00:29:52.710
We have to know these
algebraic states to determine
00:29:52.710 --> 00:29:54.024
the differential states.
00:29:54.024 --> 00:29:56.190
We have to know the
differential states to determine
00:29:56.190 --> 00:29:57.120
the algebraic states.
00:30:00.771 --> 00:30:01.740
AUDIENCE: Professor?
00:30:01.740 --> 00:30:02.739
WILLIAM GREEN, JR.: Yes.
00:30:02.739 --> 00:30:07.695
AUDIENCE: [INAUDIBLE]
00:30:07.695 --> 00:30:09.070
WILLIAM GREEN,
JR.: That's right.
00:30:09.070 --> 00:30:13.370
So here the x that I have
is not all of the states.
00:30:13.370 --> 00:30:14.690
It's not all of the unknowns.
00:30:14.690 --> 00:30:17.390
It's only the set of
the unknowns of which I
00:30:17.390 --> 00:30:18.576
have to take a differential.
00:30:18.576 --> 00:30:19.700
That's right, that's right.
00:30:25.040 --> 00:30:29.030
So we've sort of tooled
down the kinds of equations
00:30:29.030 --> 00:30:30.810
we might want to look at.
00:30:30.810 --> 00:30:33.860
We've gone from the generic
differential algebraic equation
00:30:33.860 --> 00:30:38.390
to a semi-explicit form to the
splitting, this differential
00:30:38.390 --> 00:30:39.320
algebraic splitting.
00:30:39.320 --> 00:30:43.400
Many problems can be
formulated in this fashion.
00:30:43.400 --> 00:30:45.251
This one's the
easiest one to think
00:30:45.251 --> 00:30:47.000
about so that's the
one that we'll discuss
00:30:47.000 --> 00:30:48.530
for most of the lecture.
00:30:51.740 --> 00:30:53.720
So let's look at some examples.
00:30:53.720 --> 00:30:57.340
I think examples are
interesting to think about.
00:30:57.340 --> 00:30:58.420
So here's one.
00:30:58.420 --> 00:31:02.140
Think about a mass
spring system.
00:31:02.140 --> 00:31:04.360
There's an algebraic
equation that
00:31:04.360 --> 00:31:07.270
governs a conserved
quantity, the total energy
00:31:07.270 --> 00:31:09.370
of this mass spring system.
00:31:09.370 --> 00:31:13.060
So the kinetic energy, 1/2 m
times the velocity of the mass
00:31:13.060 --> 00:31:14.200
squared.
00:31:14.200 --> 00:31:18.444
Plus the energy stored in
the spring, 1/2 kx squared
00:31:18.444 --> 00:31:20.110
has got to be equal
to the total energy.
00:31:20.110 --> 00:31:22.670
And the total energy
can change in time.
00:31:22.670 --> 00:31:27.010
So we have a differential
algebraic equation.
00:31:27.010 --> 00:31:29.860
We have a differential
state because we're taking
00:31:29.860 --> 00:31:31.420
the time derivative of x.
00:31:31.420 --> 00:31:33.610
And we'd like to be
able to determine
00:31:33.610 --> 00:31:36.040
x as a function of time.
00:31:36.040 --> 00:31:38.920
So f of x, x dot and t.
00:31:38.920 --> 00:31:42.250
That's this equation, 1/2
mv squared plus one half kx
00:31:42.250 --> 00:31:45.254
squared minus E is equal to 0.
00:31:45.254 --> 00:31:47.170
We'd like to solve this
differential algebraic
00:31:47.170 --> 00:31:48.070
equation.
00:31:48.070 --> 00:31:49.300
It's well-posed.
00:31:49.300 --> 00:31:54.580
We have one equation for
one unknown, one state.
00:31:54.580 --> 00:31:58.260
The equation has
a solution, which
00:31:58.260 --> 00:32:00.520
is a times cosine omega t.
00:32:00.520 --> 00:32:02.020
It's not the only
possible solution.
00:32:02.020 --> 00:32:04.290
It's a solution
to this equation.
00:32:04.290 --> 00:32:06.660
We can substitute that
solution into the equation
00:32:06.660 --> 00:32:10.560
and determine, for example
what this oscillation frequency
00:32:10.560 --> 00:32:11.490
omega has to be.
00:32:11.490 --> 00:32:12.810
It's square root of k over m.
00:32:12.810 --> 00:32:14.250
Or what the amplitude has to be.
00:32:14.250 --> 00:32:17.130
It's the square
root of E over k.
00:32:17.130 --> 00:32:20.496
So you've seen this
from physics before.
00:32:20.496 --> 00:32:21.870
We've already
solved differential
00:32:21.870 --> 00:32:24.030
algebraic equations.
00:32:24.030 --> 00:32:26.520
The thing you saw before
though, was actually
00:32:26.520 --> 00:32:28.500
conservation of momentum.
00:32:28.500 --> 00:32:32.190
You converted this
entirely to a second order
00:32:32.190 --> 00:32:34.380
differential equation
in x and solved it.
00:32:34.380 --> 00:32:37.980
But in principle, we could have
tried to find various solutions
00:32:37.980 --> 00:32:39.150
to this equation instead.
00:32:39.150 --> 00:32:41.802
We would have said, well, we've
observed masses and springs
00:32:41.802 --> 00:32:43.260
and we know they
oscillate in time.
00:32:43.260 --> 00:32:45.450
So let's propose some
oscillatory solutions
00:32:45.450 --> 00:32:48.300
and see what values of
the amplitude and omega we
00:32:48.300 --> 00:32:50.356
need to satisfy this equation.
00:32:50.356 --> 00:32:52.230
It would have been a
perfectly acceptable way
00:32:52.230 --> 00:32:53.280
of solving this problem.
00:33:01.010 --> 00:33:03.940
The thing is, this
mass spring problem,
00:33:03.940 --> 00:33:05.980
the differential
algebraic equation
00:33:05.980 --> 00:33:09.880
contains nonlinearities
with respect to the states.
00:33:09.880 --> 00:33:11.470
So it's nonlinear
in the velocity
00:33:11.470 --> 00:33:12.940
and it's nonlinear
in the position.
00:33:12.940 --> 00:33:15.865
Nonlinear equations in
general are hard to solve.
00:33:15.865 --> 00:33:17.740
When you converted that
to a momentum balance
00:33:17.740 --> 00:33:20.110
you got linear equations.
00:33:20.110 --> 00:33:23.590
You got the second
derivative of position
00:33:23.590 --> 00:33:27.160
was equal to stiffness
over mass times position.
00:33:27.160 --> 00:33:29.290
It was a linear differential
equation to solve,
00:33:29.290 --> 00:33:30.300
which is easy to solve.
00:33:30.300 --> 00:33:32.950
So you've got higher
order differentials,
00:33:32.950 --> 00:33:36.310
but you linearize the equation
so it made it easier to solve,
00:33:36.310 --> 00:33:38.740
in a certain sense.
00:33:38.740 --> 00:33:42.790
But we'd like to
sometimes understand what
00:33:42.790 --> 00:33:44.245
we can do with these equations.
00:33:44.245 --> 00:33:45.970
There are certain
limiting cases where
00:33:45.970 --> 00:33:48.053
we can do interesting
things with these equations.
00:33:48.053 --> 00:33:49.600
I'm going to show you one now.
00:33:49.600 --> 00:33:53.650
So in the case where
the Jacobian of this
00:33:53.650 --> 00:33:58.021
function f with respect to
the differential of the state
00:33:58.021 --> 00:33:58.520
variable.
00:33:58.520 --> 00:34:00.340
So take the derivative
of f with respect
00:34:00.340 --> 00:34:04.080
to x dot and hold time
in the state constant.
00:34:04.080 --> 00:34:06.830
So this is a matrix,
one that's full rank.
00:34:06.830 --> 00:34:10.330
Then the DA be represented
as an equivalent to ODE.
00:34:10.330 --> 00:34:12.190
Here's how you can show that.
00:34:12.190 --> 00:34:16.389
The total differential of f
is this Jacobian with respect
00:34:16.389 --> 00:34:21.969
to x dot times the change x dot
plus the Jacobian with respect
00:34:21.969 --> 00:34:24.760
to x times the change in x.
00:34:24.760 --> 00:34:28.400
Plus the partial derivative
of f with respect to time,
00:34:28.400 --> 00:34:30.230
times the change in time.
00:34:30.230 --> 00:34:32.230
And this total differential
has to be equal to 0
00:34:32.230 --> 00:34:37.080
because f is equal to 0
for all states in time.
00:34:37.080 --> 00:34:39.949
So let's divide by dt.
00:34:39.949 --> 00:34:43.050
We'd like to know how
f changes with time.
00:34:43.050 --> 00:34:45.139
So total differential
of f with time.
00:34:45.139 --> 00:34:48.300
We divide each of these
little differentials by time.
00:34:48.300 --> 00:34:53.060
And then we solve for dx dot dt.
00:34:53.060 --> 00:34:55.060
That's going to involve
moving this term
00:34:55.060 --> 00:34:56.469
to the other side
of the equation
00:34:56.469 --> 00:34:58.340
and inverting this Jacobian.
00:34:58.340 --> 00:35:00.910
So if I can invert
the Jacobian then
00:35:00.910 --> 00:35:05.125
I can convert my DAE system
into a system of ODEs.
00:35:08.197 --> 00:35:09.530
This equals 0 shouldn't be here.
00:35:09.530 --> 00:35:10.050
I apologize.
00:35:10.050 --> 00:35:13.120
That's another typo.
00:35:13.120 --> 00:35:15.920
I move this term to the
other side of the equation
00:35:15.920 --> 00:35:19.580
to where the 0 is and I solve.
00:35:19.580 --> 00:35:24.180
So I went from a first order
equation, a first order
00:35:24.180 --> 00:35:27.330
equation in x to a second
order equation in x.
00:35:27.330 --> 00:35:30.120
Two differentials of x with
time instead of one differential
00:35:30.120 --> 00:35:30.900
of x with time.
00:35:30.900 --> 00:35:32.525
But I could convert
it to an equivalent
00:35:32.525 --> 00:35:34.990
system of ordinary
differential equations.
00:35:34.990 --> 00:35:36.836
Does that make sense?
00:35:36.836 --> 00:35:39.220
Do you see how that works?
00:35:39.220 --> 00:35:40.776
Any questions?
00:35:40.776 --> 00:35:44.860
Let's go back and
look at an example.
00:35:44.860 --> 00:35:45.900
The mass spring system.
00:35:49.430 --> 00:35:53.810
So here's our ODE that
we're trying to solve.
00:35:53.810 --> 00:35:56.360
I'm sorry, our DAE that
we're trying to solve.
00:35:56.360 --> 00:36:00.440
Can you convert this to
an equivalent ordinary
00:36:00.440 --> 00:36:04.030
differential equation using
the approach I just described?
00:36:04.030 --> 00:36:06.900
Think you guys can
figure that out?
00:36:06.900 --> 00:36:07.400
Try it.
00:36:07.400 --> 00:36:10.555
You can work together.
00:36:10.555 --> 00:37:07.000
[SIDE CONVERSATIONS]
00:37:07.000 --> 00:37:10.070
WILLIAM GREEN, JR.:
Guys are either tired
00:37:10.070 --> 00:37:11.880
or friendships
have been destroyed
00:37:11.880 --> 00:37:15.150
over the mid-term season.
00:37:15.150 --> 00:37:15.770
Maybe both.
00:37:38.757 --> 00:37:40.590
So let's compute the
things we need to know.
00:37:40.590 --> 00:37:44.400
We needed to know the partials
of this f with respect
00:37:44.400 --> 00:37:46.180
to the states and
with respect to time.
00:37:46.180 --> 00:37:50.190
So partial f respect to x dot,
holding all the other variables
00:37:50.190 --> 00:37:51.840
constant, is mx dot.
00:37:51.840 --> 00:37:55.120
Partial f with
respect to x is kx.
00:37:55.120 --> 00:37:58.591
Partial f with
respect to t is t.
00:37:58.591 --> 00:38:00.090
And I told you
before that you could
00:38:00.090 --> 00:38:05.645
write this says dx
dot dt is equal to df
00:38:05.645 --> 00:38:10.670
dx dot inverse multiplied by--
00:38:10.670 --> 00:38:14.264
is equal to minus df dx dot
inverse multiplied by df dx.
00:38:14.264 --> 00:38:16.430
This was on the previous
slide, which just gives you
00:38:16.430 --> 00:38:19.240
conservation of momentum.
00:38:19.240 --> 00:38:25.390
This gives you the acceleration
is equal to the force minus kx
00:38:25.390 --> 00:38:26.870
divided by the mass m.
00:38:31.300 --> 00:38:33.875
So we can solve either,
the DAE or the ODE.
00:38:37.890 --> 00:38:40.220
Here's a more
complicated example,
00:38:40.220 --> 00:38:43.320
but it works, well,
it's more complicated.
00:38:43.320 --> 00:38:46.010
So a lot of times we do
molecular dynamics simulations.
00:38:46.010 --> 00:38:48.130
We try to model
molecules moving around.
00:38:48.130 --> 00:38:50.232
There we also want
to conserve energy.
00:38:50.232 --> 00:38:52.190
There's usually a conserved
quantity, something
00:38:52.190 --> 00:38:53.870
we're trying to hold constant.
00:38:53.870 --> 00:38:55.160
So suppose we try to do that.
00:38:55.160 --> 00:38:57.770
The total energy in this
system of molecules,
00:38:57.770 --> 00:39:01.610
maybe it's 1/2 mass times
the velocity squared.
00:39:01.610 --> 00:39:04.010
Now x is the position
of all these molecules.
00:39:04.010 --> 00:39:05.990
Ad x dot is the velocity
of all these molecules
00:39:05.990 --> 00:39:09.650
plus the potential energy,
which is a function of f,
00:39:09.650 --> 00:39:12.510
a function of their positions.
00:39:12.510 --> 00:39:17.391
So the DAE representing
this constraint is this.
00:39:17.391 --> 00:39:18.890
It's the same as
for the mass spring
00:39:18.890 --> 00:39:20.931
system, where we've replaced
the potential energy
00:39:20.931 --> 00:39:21.980
with a generic one.
00:39:21.980 --> 00:39:24.520
This actually isn't
a well-posed DEA.
00:39:24.520 --> 00:39:28.430
We have one equation
for, let's see,
00:39:28.430 --> 00:39:31.146
if we have n molecules in
here we have three n states.
00:39:31.146 --> 00:39:32.770
We can move around
in three dimensions.
00:39:32.770 --> 00:39:35.605
So this is a really
ill-posed equation.
00:39:35.605 --> 00:39:37.730
If we try to apply the same
procedure, then we say,
00:39:37.730 --> 00:39:40.040
well, this quantity still
should be conserved.
00:39:40.040 --> 00:39:43.180
Over time it needs
to be equal to zero.
00:39:43.180 --> 00:39:47.630
Then we can try to compute
all these differentials.
00:39:47.630 --> 00:39:50.630
But we'll find out that
partial f, partial x dot,
00:39:50.630 --> 00:39:53.630
the Jacobian of f
with respect to x dot,
00:39:53.630 --> 00:39:56.320
that's just the momentum.
00:39:56.320 --> 00:39:57.605
That's not a matrix.
00:39:57.605 --> 00:40:00.120
That's not invertible.
00:40:00.120 --> 00:40:02.390
So we can't convert
this equation
00:40:02.390 --> 00:40:04.600
to an equivalent system
of ordinary differential
00:40:04.600 --> 00:40:05.100
equations.
00:40:05.100 --> 00:40:06.558
So it's just
something pathological
00:40:06.558 --> 00:40:08.930
for a one-dimensional system.
00:40:08.930 --> 00:40:13.010
Instead, what we find is
zero is equal to the velocity
00:40:13.010 --> 00:40:17.740
dotted with the mass times the
acceleration plus the gradient
00:40:17.740 --> 00:40:19.040
in the potential.
00:40:19.040 --> 00:40:22.880
So minus the forces
acting on the molecules.
00:40:22.880 --> 00:40:24.590
We know that if we
were to integrate
00:40:24.590 --> 00:40:27.590
the equations of motion
for the molecules exactly,
00:40:27.590 --> 00:40:30.710
the momentum balances for
those molecules exactly,
00:40:30.710 --> 00:40:32.900
then this term here would
always be equal to zero.
00:40:32.900 --> 00:40:37.090
And we'd satisfy
conservation of energy.
00:40:37.090 --> 00:40:40.130
Of course, all solutions of
ordinary differential equations
00:40:40.130 --> 00:40:42.170
require approximations.
00:40:42.170 --> 00:40:44.312
So as you move the
molecules around,
00:40:44.312 --> 00:40:45.770
their velocities
and positions will
00:40:45.770 --> 00:40:48.200
tend to drift away
from the solution
00:40:48.200 --> 00:40:51.950
where this is exactly in
balance and equal to zero.
00:40:51.950 --> 00:40:54.950
If you desire to do that in the
molecular dynamic simulation
00:40:54.950 --> 00:40:56.630
then you better
have a method that
00:40:56.630 --> 00:40:59.960
somehow respects this
geometric property instead.
00:40:59.960 --> 00:41:02.950
That the velocities
are orthogonal to m
00:41:02.950 --> 00:41:05.000
x dot plus grad v.
00:41:05.000 --> 00:41:08.180
And there is a set of
ODE-IVP methods that do that.
00:41:08.180 --> 00:41:10.841
They're called
symplectic integrators.
00:41:10.841 --> 00:41:12.590
And they integrate the
equations of motion
00:41:12.590 --> 00:41:16.220
while exerting some control over
the error in the total energy.
00:41:16.220 --> 00:41:19.040
So if you were just
the take a ODE45
00:41:19.040 --> 00:41:21.080
and try to solve this
system of equations,
00:41:21.080 --> 00:41:24.530
over long times you would find
that the energy drifts away
00:41:24.530 --> 00:41:26.340
from the place where it started.
00:41:26.340 --> 00:41:28.670
Even though your positions
are being integrated
00:41:28.670 --> 00:41:31.297
with reasonable
accuracy, over time
00:41:31.297 --> 00:41:33.380
the energy will drift away
from where you started.
00:41:33.380 --> 00:41:36.590
The system will heat up
effectively or cool down.
00:41:36.590 --> 00:41:39.200
That's undesirable if you're
trying to do modeling.
00:41:39.200 --> 00:41:41.630
So instead, one uses the
symplectic integrators,
00:41:41.630 --> 00:41:45.810
which are designed to control
the error and the energy
00:41:45.810 --> 00:41:48.336
in addition to integrate
the equations of motion.
00:41:48.336 --> 00:41:50.210
So you can look those
up and read about them.
00:41:50.210 --> 00:41:51.694
I think they're
quite interesting.
00:41:51.694 --> 00:41:53.360
If you're going to
do molecular modeling
00:41:53.360 --> 00:41:55.430
you'd like to understand
better how they work.
00:41:55.430 --> 00:41:56.929
But actually,
conservation of energy
00:41:56.929 --> 00:41:59.870
here doesn't give
us a well-posed DAE.
00:41:59.870 --> 00:42:02.060
We actually need the
momentum equations
00:42:02.060 --> 00:42:07.680
to formulate or to determine the
trajectories of the molecules.
00:42:07.680 --> 00:42:12.120
Let's talk about their
numerical simulation.
00:42:12.120 --> 00:42:14.540
So let's think about
these semi-explicit DAEs,
00:42:14.540 --> 00:42:16.640
the ones that can be
split between differential
00:42:16.640 --> 00:42:17.790
and algebraic states.
00:42:17.790 --> 00:42:21.740
It's kind of useful to look
at this problem in particular.
00:42:21.740 --> 00:42:25.700
And you might say, well, let's
do the simplest possible thing.
00:42:25.700 --> 00:42:29.000
Let's do the forward
Euler approximation.
00:42:29.000 --> 00:42:31.880
So let's take our derivative
here and represent it
00:42:31.880 --> 00:42:36.500
as the state evaluated at x plus
delta t minus the state at t
00:42:36.500 --> 00:42:38.080
divided by delta t.
00:42:38.080 --> 00:42:40.070
Let's make the right-hand
side of our equation
00:42:40.070 --> 00:42:43.340
be evaluated at time t.
00:42:43.340 --> 00:42:46.550
And let's try to step in
time, do time marching forward
00:42:46.550 --> 00:42:49.550
to determine the
trajectory of x and t.
00:42:49.550 --> 00:42:51.440
How do the states change?
00:42:51.440 --> 00:42:53.210
We see the first
equation is an easy one
00:42:53.210 --> 00:42:58.870
to solve because presumably
I know x of t and y of t.
00:42:58.870 --> 00:43:01.550
Well, do I know y of t?
y of t has to satisfy
00:43:01.550 --> 00:43:02.720
this second equation.
00:43:02.720 --> 00:43:06.200
So first, if I know x of t I
need to solve this nonlinear
00:43:06.200 --> 00:43:08.190
equation determine y of t.
00:43:08.190 --> 00:43:09.830
Then I know x of
t and y of t and I
00:43:09.830 --> 00:43:14.390
can step forward in time to
determine x of t plus delta t.
00:43:14.390 --> 00:43:17.720
So every iteration here with the
simplest possible approximation
00:43:17.720 --> 00:43:20.030
for the differential
requires solution of a system
00:43:20.030 --> 00:43:21.330
of nonlinear equations.
00:43:24.230 --> 00:43:27.540
We already saw before that a
forward Euler approximation
00:43:27.540 --> 00:43:32.640
is not a great one to apply
to solutions of IDE-IVPs
00:43:32.640 --> 00:43:35.660
because it has very limited
stability properties.
00:43:35.660 --> 00:43:37.630
You need small
time steps in order
00:43:37.630 --> 00:43:39.560
to stably integrate in time.
00:43:42.334 --> 00:43:43.750
So even though
this seemed like it
00:43:43.750 --> 00:43:46.480
was easy to get a solution
for the differential equation,
00:43:46.480 --> 00:43:51.620
we've always got to satisfy
this nonlinear equation here.
00:43:51.620 --> 00:43:54.970
So inherently, simulation
of DAEs are implicit.
00:43:54.970 --> 00:43:57.700
They require solution
of nonlinear equations
00:43:57.700 --> 00:43:59.980
in order to advance in time.
00:43:59.980 --> 00:44:03.590
There's actually no point in
using this weakly stable method
00:44:03.590 --> 00:44:04.990
for any very small time steps.
00:44:04.990 --> 00:44:08.200
It makes more sense to choose
a naturally implicit method
00:44:08.200 --> 00:44:11.830
for the differential equation,
where I can take big time steps
00:44:11.830 --> 00:44:16.376
and still have reasonable
stability of the integration.
00:44:16.376 --> 00:44:17.250
Does that make sense?
00:44:21.040 --> 00:44:25.180
Consider now the
fully implicit DAE.
00:44:25.180 --> 00:44:28.480
So f of x, x dot,
and t is equal to 0
00:44:28.480 --> 00:44:30.565
and have some set of
initial conditions. x of 0
00:44:30.565 --> 00:44:32.590
is equal to x 0.
00:44:32.590 --> 00:44:34.960
Here there is no way in
general to avoid solving
00:44:34.960 --> 00:44:37.330
systems of nonlinear equations.
00:44:37.330 --> 00:44:39.850
And so one substitutes
often backwards
00:44:39.850 --> 00:44:42.670
difference
approximations for x dot.
00:44:42.670 --> 00:44:45.355
And then you solve for
the next point in time.
00:44:45.355 --> 00:44:47.230
So here's an example
with the backward Euler,
00:44:47.230 --> 00:44:48.730
implicit Euler approximations.
00:44:48.730 --> 00:44:51.070
So I replaced the
time derivative
00:44:51.070 --> 00:44:55.000
with the difference between
the state at point tk
00:44:55.000 --> 00:44:57.160
and the state of
point tk minus 1
00:44:57.160 --> 00:45:00.880
divided by the difference
between tk and tk minus 1.
00:45:00.880 --> 00:45:03.670
The other state variables
are evaluated at tk.
00:45:03.670 --> 00:45:05.130
So x at tk.
00:45:05.130 --> 00:45:06.430
And I evaluate time at tk.
00:45:06.430 --> 00:45:08.650
And I try to satisfy
this equation
00:45:08.650 --> 00:45:10.460
and use it to determine x of tk.
00:45:10.460 --> 00:45:13.690
So I solve this system
of nonlinear equations
00:45:13.690 --> 00:45:15.160
for x at tk.
00:45:15.160 --> 00:45:18.580
And then I go to the next
point in time, tk plus 1.
00:45:18.580 --> 00:45:22.440
And I solve the same equation,
but replacing tk with tk plus 1
00:45:22.440 --> 00:45:25.060
and tk minus 1 with tk.
00:45:25.060 --> 00:45:28.060
So at each iteration I solve a
system of nonlinear equations.
00:45:28.060 --> 00:45:30.400
No more painful than
doing the forward Euler
00:45:30.400 --> 00:45:32.860
approximation I
showed you before.
00:45:32.860 --> 00:45:35.590
So that's good.
00:45:35.590 --> 00:45:38.865
But still a lot of work.
00:45:38.865 --> 00:45:40.240
We might wonder
how do we control
00:45:40.240 --> 00:45:42.872
the error in these
sorts of approximations.
00:45:42.872 --> 00:45:44.830
And that's what we're
going to talk about next.
00:45:48.394 --> 00:45:50.310
So here's the equivalent
time marching formula
00:45:50.310 --> 00:45:52.770
applied to a semi-explicit DAE.
00:45:55.910 --> 00:46:02.790
So I've got to
satisfy the equation
00:46:02.790 --> 00:46:05.850
0 is equal to the derivative
of x with respect to time,
00:46:05.850 --> 00:46:09.060
which I approximate with the
backward Euler approximation,
00:46:09.060 --> 00:46:12.300
minus f evaluated at tk.
00:46:12.300 --> 00:46:19.410
And the equation g of x at
tk, y at tk for this x of tk.
00:46:19.410 --> 00:46:20.670
And actually y of tk as well.
00:46:20.670 --> 00:46:22.800
I'm sorry, I left that
off. y of tk as well.
00:46:22.800 --> 00:46:24.890
Got to solve for
the differential
00:46:24.890 --> 00:46:28.340
and the algebraic
states simultaneously.
00:46:28.340 --> 00:46:30.030
And I use a backwards
difference formula
00:46:30.030 --> 00:46:31.877
to try to do this
in a stable way.
00:46:36.850 --> 00:46:43.180
One way to think about in
particular these sorts of DAEs
00:46:43.180 --> 00:46:46.990
is as exceptionally stiff
ordinary differential
00:46:46.990 --> 00:46:47.490
equations.
00:46:47.490 --> 00:46:50.320
Did you guys discuss stiffness?
00:46:50.320 --> 00:46:55.750
So imagine if I replace
the 0 with, say, dy dt.
00:46:55.750 --> 00:46:58.090
And I introduced a constant
on the right-hand side,
00:46:58.090 --> 00:47:05.485
a multiplicative constant,
say, 1 over lambda here.
00:47:05.485 --> 00:47:06.610
Let's say lambda, actually.
00:47:06.610 --> 00:47:08.610
Say I introduce a
multiplicative constant lambda
00:47:08.610 --> 00:47:10.940
on the right-hand side here.
00:47:10.940 --> 00:47:15.010
So as lambda becomes
very large, the system
00:47:15.010 --> 00:47:17.620
becomes increasingly stiff.
00:47:17.620 --> 00:47:19.600
The dynamics of
this equation take
00:47:19.600 --> 00:47:22.790
place over very, very
short time scales.
00:47:22.790 --> 00:47:24.820
So short in fact,
that eventually,
00:47:24.820 --> 00:47:26.500
if I let lambda
diverge to infinity,
00:47:26.500 --> 00:47:29.500
I've just got to satisfy
this algebraic equation
00:47:29.500 --> 00:47:31.400
at each point in time.
00:47:31.400 --> 00:47:34.240
So these are exceptionally
stiff equations.
00:47:34.240 --> 00:47:36.270
And the methods for
solving these equations
00:47:36.270 --> 00:47:39.880
share a lot in common with
stiff ODE-IVP solvers.
00:47:45.069 --> 00:47:46.610
So I've got a couple
of minutes left.
00:47:46.610 --> 00:47:48.880
I'm going to work
through some examples.
00:47:48.880 --> 00:47:51.540
So consider the stirred
tank example one.
00:47:51.540 --> 00:47:53.250
dc 2 dt.
00:47:53.250 --> 00:47:57.210
It's related to q over
v, flow rate over volume.
00:47:57.210 --> 00:47:58.770
The difference in
the concentrations
00:47:58.770 --> 00:48:00.570
coming in and out of the tank.
00:48:00.570 --> 00:48:03.060
And c1 is equal to this
input signal gamma.
00:48:03.060 --> 00:48:06.030
Maybe that's a measurement
and I'm trying to predict c2.
00:48:06.030 --> 00:48:08.850
So can you apply the backward
Euler method to these equations
00:48:08.850 --> 00:48:12.810
to determine c1 and c2?
00:48:12.810 --> 00:48:15.670
What would one step of the
backward Euler method look
00:48:15.670 --> 00:48:24.500
like, say, going from time 0
to time delta t, or time t to t
00:48:24.500 --> 00:48:27.260
plus delta t, or tk to
tk plus 1, however you
00:48:27.260 --> 00:48:28.748
want to write it?
00:49:22.314 --> 00:49:23.730
So do you get
something like this?
00:49:26.480 --> 00:49:28.530
So this equation is easy.
00:49:28.530 --> 00:49:32.360
c1 to tk is gamma of tk.
00:49:32.360 --> 00:49:34.390
I got a substitute up here.
00:49:34.390 --> 00:49:38.570
c2 of tk minus c2 of tk
minus 1 over delta t.
00:49:38.570 --> 00:49:41.210
T Just tk minus tk minus 1.
00:49:41.210 --> 00:49:42.800
And then solve for c2 of dk.
00:49:42.800 --> 00:49:44.390
And you get a formula like this.
00:49:49.010 --> 00:49:52.010
This formula was based
on an approximation
00:49:52.010 --> 00:49:56.510
for the derivative,
which had a leading order
00:49:56.510 --> 00:50:03.140
error that was proportional
to tk minus tk minus 1.
00:50:03.140 --> 00:50:07.070
When I go through and I solve
for c2 of tk, the leading order
00:50:07.070 --> 00:50:10.600
error in c2, the
local truncation error
00:50:10.600 --> 00:50:16.250
is going to be order tk
minus tk minus 1 squared.
00:50:16.250 --> 00:50:20.840
So as I shrink the spacing
between each of my two time
00:50:20.840 --> 00:50:24.050
points, the local
accuracy, the error
00:50:24.050 --> 00:50:27.830
induced in one advance of
this time stepping algorithm
00:50:27.830 --> 00:50:29.790
is going to scale like
the step size squared.
00:50:29.790 --> 00:50:31.550
So I have the error--
00:50:31.550 --> 00:50:32.369
I have the spacing.
00:50:32.369 --> 00:50:34.160
The error will go down
by a factor of four.
00:50:34.160 --> 00:50:37.667
I reduce the spacing by
an order of magnitude.
00:50:37.667 --> 00:50:39.750
The error will go down by
two orders of magnitude.
00:50:39.750 --> 00:50:43.605
This is a really desirable sort
of error control and a method.
00:50:43.605 --> 00:50:44.480
And we can do better.
00:50:44.480 --> 00:50:48.460
You've seen us do better
with other methods.
00:50:48.460 --> 00:50:51.180
Let's look at this equation now.
00:50:51.180 --> 00:50:52.590
This is stirred
tank example two.
00:50:52.590 --> 00:50:53.797
I just replace c1 with c2.
00:50:53.797 --> 00:50:55.380
You don't have to
write this one down.
00:50:55.380 --> 00:50:56.484
I'll just show it to you.
00:50:56.484 --> 00:50:58.650
The notes are online so
you're not missing anything.
00:50:58.650 --> 00:51:00.900
All these things will pop
up in the notes online.
00:51:00.900 --> 00:51:05.220
So c2 at tk is equal
to gamma of tk.
00:51:05.220 --> 00:51:06.930
I got to come up to
my first equation now
00:51:06.930 --> 00:51:08.700
and I have to substitute
an approximation
00:51:08.700 --> 00:51:09.810
for the derivative.
00:51:09.810 --> 00:51:13.390
It's this, the backwards
difference formula.
00:51:13.390 --> 00:51:17.610
So if find c1 of tk is
c2 of tk plus v over q.
00:51:17.610 --> 00:51:21.660
c2 of tk minus c2 of tk minus
1 divided by the time step.
00:51:21.660 --> 00:51:23.640
Remember, this approximation
for the derivative
00:51:23.640 --> 00:51:27.292
has a leading order error that
goes like tk minus tk minus 1.
00:51:27.292 --> 00:51:32.235
So my approximation for c1,
the local approximation for c1,
00:51:32.235 --> 00:51:34.040
is not second order accurate.
00:51:34.040 --> 00:51:38.010
It's only first order accurate.
00:51:38.010 --> 00:51:42.450
I applied the same approximation
method to both equations.
00:51:42.450 --> 00:51:45.050
One equation I got something
that was second order accurate.
00:51:45.050 --> 00:51:46.425
The next equation
I got something
00:51:46.425 --> 00:51:47.840
that was first order accurate.
00:51:47.840 --> 00:51:49.520
They look almost identical.
00:51:49.520 --> 00:51:51.530
It might be hard to
see from the start
00:51:51.530 --> 00:51:53.420
why it should work out that way.
00:51:53.420 --> 00:51:56.190
That's how it works out.
00:51:56.190 --> 00:51:57.410
Let's do one more example.
00:52:00.640 --> 00:52:03.680
Here's a system of DAEs.
00:52:03.680 --> 00:52:05.760
c2 dot is equal to c1.
00:52:05.760 --> 00:52:07.900
c3 dot is equal to c2.
00:52:07.900 --> 00:52:10.240
And 0 is equal to
c3 minus gamma.
00:52:10.240 --> 00:52:12.250
So gamma, again, is
some measurement.
00:52:12.250 --> 00:52:14.660
We solve this equation
to determine c3.
00:52:14.660 --> 00:52:17.640
We substitute c3 up
here to determine c2.
00:52:17.640 --> 00:52:20.600
We substitute c2 up
here to determine c1.
00:52:20.600 --> 00:52:23.965
But we want to find the solution
with the backward Euler method.
00:52:27.100 --> 00:52:32.680
And the solution there is going
to be c3 of tk is gamma of tk.
00:52:32.680 --> 00:52:35.410
c2 is related to the
derivative of c3.
00:52:35.410 --> 00:52:37.180
So I use my backwards
difference formula
00:52:37.180 --> 00:52:39.060
for the derivative of c3.
00:52:39.060 --> 00:52:41.620
c1 is related to the
derivative of c2.
00:52:41.620 --> 00:52:44.710
So I use my backward
difference formula for c2.
00:52:44.710 --> 00:52:46.930
But how accurate are each
of these approximations?
00:52:49.840 --> 00:52:52.780
So here I have to determine
c2 from an approximation
00:52:52.780 --> 00:52:54.880
for the derivative of c3.
00:52:54.880 --> 00:52:56.920
We know this backwards
difference approximation
00:52:56.920 --> 00:52:58.510
has a leading order
error proportional
00:52:58.510 --> 00:53:02.350
to delta t, the time step.
00:53:02.350 --> 00:53:07.750
Now I got to approximate c1 with
the finite difference backwards
00:53:07.750 --> 00:53:09.922
difference approximation in c2.
00:53:09.922 --> 00:53:12.380
We know this should incur an
error proportional to the time
00:53:12.380 --> 00:53:17.350
step delta t, except I
don't know c2 at tk exactly.
00:53:17.350 --> 00:53:23.530
Actually, c2 tk, the
exact value of c2 at tk
00:53:23.530 --> 00:53:27.830
is this plus something
proportional to the time step.
00:53:27.830 --> 00:53:30.560
Which means the
exact value of c1
00:53:30.560 --> 00:53:34.430
is this plus something
that's order one because I
00:53:34.430 --> 00:53:36.299
carry over the error
from my solution
00:53:36.299 --> 00:53:37.340
to the previous equation.
00:53:37.340 --> 00:53:40.477
And I divide it
by the time step.
00:53:40.477 --> 00:53:42.560
So I went from having
something-- well, let's see.
00:53:42.560 --> 00:53:44.340
The first equation
had order delta t
00:53:44.340 --> 00:53:47.040
squared local truncation error.
00:53:47.040 --> 00:53:50.010
As I shrink the
spacing, my solution
00:53:50.010 --> 00:53:52.440
gets more and more
locally accurate.
00:53:52.440 --> 00:53:57.300
And it seems to do it in
a fairly robust fashion.
00:53:57.300 --> 00:53:59.910
The next system of
equations we solved
00:53:59.910 --> 00:54:02.170
had a local truncation error
proportional to delta t.
00:54:02.170 --> 00:54:03.545
Again, I can shrink
delta t and I
00:54:03.545 --> 00:54:05.800
can make my solution
more and more accurate.
00:54:05.800 --> 00:54:08.220
That seems like what you want
in a numerical algorithm.
00:54:08.220 --> 00:54:11.580
This one, there's no hope.
00:54:11.580 --> 00:54:15.980
Change delta t, who knows
what you're going to get?
00:54:15.980 --> 00:54:18.855
And this looks fairly innocuous.
00:54:18.855 --> 00:54:20.990
You just look at it
and you say, well, OK.
00:54:20.990 --> 00:54:22.980
Let's solve it and
see what happens.
00:54:22.980 --> 00:54:26.040
Well actually, you can't
solve this problem accurately
00:54:26.040 --> 00:54:28.322
with the backward Euler method.
00:54:28.322 --> 00:54:30.030
There's something
fundamentally different
00:54:30.030 --> 00:54:32.359
about this problem
from the previous one.
00:54:32.359 --> 00:54:34.650
There's something fundamentally
different about stirred
00:54:34.650 --> 00:54:37.500
tank example two from
stirred tank example one.
00:54:37.500 --> 00:54:40.170
Right you might have
already perceived that.
00:54:40.170 --> 00:54:42.160
You can see what's
going on here.
00:54:42.160 --> 00:54:44.520
c3 is equal to gamma.
00:54:44.520 --> 00:54:46.025
That's the exact solution.
00:54:46.025 --> 00:54:48.390
c2 is equal to gamma dot.
00:54:48.390 --> 00:54:50.700
And c1 is equal
to two derivatives
00:54:50.700 --> 00:54:53.860
of gamma, the exact
solution to this equation.
00:54:53.860 --> 00:54:59.310
So c1 is incredibly sensitive
to changes in gamma.
00:54:59.310 --> 00:55:00.620
And that's peculiar.
00:55:00.620 --> 00:55:03.180
And it has important
numerical consequences.
00:55:03.180 --> 00:55:06.380
So I'm going to show you, our
next lecture on Wednesday,
00:55:06.380 --> 00:55:09.620
how to formalize
an understanding
00:55:09.620 --> 00:55:11.450
of the differences
between these problems.
00:55:11.450 --> 00:55:14.990
And we'll talk about some
more subtleties associated
00:55:14.990 --> 00:55:17.150
with differential
algebraic equations.
00:55:17.150 --> 00:55:18.650
You guys can pick
up your quizzes.
00:55:18.650 --> 00:55:23.000
We should do it
outside rather than
00:55:23.000 --> 00:55:25.596
inside because the next
class has to start.