1
00:00:01,580 --> 00:00:03,920
The following content is
provided under a Creative
2
00:00:03,920 --> 00:00:05,340
Commons license.
3
00:00:05,340 --> 00:00:07,550
Your support will help
MIT OpenCourseWare
4
00:00:07,550 --> 00:00:11,640
continue to offer high quality
educational resources for free.
5
00:00:11,640 --> 00:00:14,180
To make a donation or to
view additional materials
6
00:00:14,180 --> 00:00:18,110
from hundreds of MIT courses,
visit MIT OpenCourseWare
7
00:00:18,110 --> 00:00:19,340
at ocw.mit.edu.
8
00:00:23,562 --> 00:00:25,270
WILLIAM GREEN, JR.:
Shall we get started?
9
00:00:25,270 --> 00:00:27,340
Yes.
10
00:00:27,340 --> 00:00:29,710
Apologies for the notes
being posted so late.
11
00:00:29,710 --> 00:00:33,020
I had trouble with the new two
factor authentication system,
12
00:00:33,020 --> 00:00:35,980
which I need to get into
Stellar and put things up there.
13
00:00:35,980 --> 00:00:39,610
My phone wasn't working so I
had to find my iPad someplace
14
00:00:39,610 --> 00:00:42,542
and use that to get
credentials to log in.
15
00:00:42,542 --> 00:00:44,500
So if you don't have them
don't worry about it.
16
00:00:44,500 --> 00:00:47,540
The full notes are posted
online, no missing pieces.
17
00:00:47,540 --> 00:00:49,480
So you can utilize
those afterwards
18
00:00:49,480 --> 00:00:51,060
if something is unclear.
19
00:00:51,060 --> 00:00:51,560
All right?
20
00:00:56,250 --> 00:00:58,070
So your quizzes are graded.
21
00:00:58,070 --> 00:01:00,390
The TAs have them.
22
00:01:00,390 --> 00:01:02,610
I think overall, we
were really pleased
23
00:01:02,610 --> 00:01:04,349
with the outcome of this quiz.
24
00:01:04,349 --> 00:01:06,360
It is very conceptual
in nature and we
25
00:01:06,360 --> 00:01:08,790
thought people really
understood a lot of the material
26
00:01:08,790 --> 00:01:09,540
quite well.
27
00:01:09,540 --> 00:01:12,440
And you can see that from
the distribution of scores.
28
00:01:12,440 --> 00:01:15,650
So the mean on the
quiz was about 77,
29
00:01:15,650 --> 00:01:17,250
a standard deviation of 12.
30
00:01:17,250 --> 00:01:19,930
I think everyone did
really quite well.
31
00:01:19,930 --> 00:01:23,580
And on a number of the problems,
perfect solutions were posted.
32
00:01:23,580 --> 00:01:26,670
So I'll take a second here.
33
00:01:26,670 --> 00:01:28,800
If there are any questions
about the material that
34
00:01:28,800 --> 00:01:31,464
was on the quiz,
or the format of it
35
00:01:31,464 --> 00:01:32,880
that you'd like
to have addressed,
36
00:01:32,880 --> 00:01:35,351
I'm happy to answer them.
37
00:01:35,351 --> 00:01:35,850
No?
38
00:01:38,580 --> 00:01:41,440
Does it seem like it
took too much time
39
00:01:41,440 --> 00:01:44,080
to complete all three questions
in the two hours allotted?
40
00:01:44,080 --> 00:01:49,510
Or was it about the right length
of material for the time slot?
41
00:01:49,510 --> 00:01:51,030
So we vetted it
with the two TAs.
42
00:01:51,030 --> 00:01:53,720
And it took each of them less
than an hour to complete it.
43
00:01:53,720 --> 00:01:55,969
So that was a good
sign that we thought
44
00:01:55,969 --> 00:01:57,760
you guys should be able
to get through most
45
00:01:57,760 --> 00:01:59,140
of the material in two hours.
46
00:02:05,560 --> 00:02:08,012
So let's recap.
47
00:02:08,012 --> 00:02:10,470
We're going to move on today
to a topic called differential
48
00:02:10,470 --> 00:02:11,650
algebraic equations.
49
00:02:11,650 --> 00:02:14,310
But before doing that, there
are a couple of concepts
50
00:02:14,310 --> 00:02:16,260
I want to review from
numerical integration,
51
00:02:16,260 --> 00:02:19,669
or actually concepts that
weren't covered when we talked
52
00:02:19,669 --> 00:02:21,210
about numerical
integration on Friday
53
00:02:21,210 --> 00:02:22,480
that I think are important.
54
00:02:22,480 --> 00:02:23,855
And then I want
to review briefly
55
00:02:23,855 --> 00:02:26,190
implicit methods for
ODE-IVPs because those
56
00:02:26,190 --> 00:02:28,710
are going to be important
for differential algebraic
57
00:02:28,710 --> 00:02:29,250
equations.
58
00:02:29,250 --> 00:02:32,550
So for numerical
integration we talked
59
00:02:32,550 --> 00:02:35,250
about quadrature of
various integrals,
60
00:02:35,250 --> 00:02:36,960
how to develop
quadrature formulas.
61
00:02:36,960 --> 00:02:38,130
But there's actually
one type of integral
62
00:02:38,130 --> 00:02:39,150
we didn't talk
about there, which
63
00:02:39,150 --> 00:02:41,910
is sort of improper integrals,
where the limits of integration
64
00:02:41,910 --> 00:02:43,440
are unbounded.
65
00:02:43,440 --> 00:02:46,670
And these come up all the
time in engineering problems.
66
00:02:46,670 --> 00:02:49,560
We'd like to know what
the, say, integrated
67
00:02:49,560 --> 00:02:53,576
response of some function
is over all possible times.
68
00:02:53,576 --> 00:02:55,450
And we may have to
evaluate that numerically.
69
00:02:55,450 --> 00:02:57,860
So how do you do those
sorts of integrals?
70
00:02:57,860 --> 00:03:00,090
And it turns out the strategy
that's most often used
71
00:03:00,090 --> 00:03:03,420
is to divide this domain
up into two pieces.
72
00:03:03,420 --> 00:03:05,370
One piece, which
is proper integral
73
00:03:05,370 --> 00:03:07,680
over a finite domain,
and another piece,
74
00:03:07,680 --> 00:03:09,930
which is an improper integral
over an infinite domain.
75
00:03:09,930 --> 00:03:11,513
This first integral,
you can handle it
76
00:03:11,513 --> 00:03:15,120
with an ODE-IVP method or
some higher order polynomial
77
00:03:15,120 --> 00:03:16,560
interpolation.
78
00:03:16,560 --> 00:03:19,494
But the second one we
have to do differently.
79
00:03:19,494 --> 00:03:21,160
And there are couple
of ways to do this.
80
00:03:21,160 --> 00:03:23,370
One is to transform variables.
81
00:03:23,370 --> 00:03:27,070
So you map this domain
onto a finite domain.
82
00:03:27,070 --> 00:03:28,920
And the other option
is to substitute
83
00:03:28,920 --> 00:03:30,420
an asymptotic approximation.
84
00:03:30,420 --> 00:03:34,740
So we say, when, say,
tau is large, f of tau
85
00:03:34,740 --> 00:03:37,410
has a characteristic
asymptotic approximation.
86
00:03:37,410 --> 00:03:40,200
Maybe it's
exponentially decaying.
87
00:03:40,200 --> 00:03:43,080
And we take the integral of
that asymptotic approximation.
88
00:03:43,080 --> 00:03:44,910
And the more accurate
that approximation
89
00:03:44,910 --> 00:03:48,150
is the better our
approximation for this whole
90
00:03:48,150 --> 00:03:49,800
integral will be.
91
00:03:49,800 --> 00:03:52,370
There's another type of
interval that we have to do,
92
00:03:52,370 --> 00:03:54,960
where the same idea
can be applied.
93
00:03:54,960 --> 00:03:58,485
That's the integral over
integrable singularities.
94
00:03:58,485 --> 00:03:59,360
So here's an example.
95
00:04:02,810 --> 00:04:06,300
So say we want to integrate
the function cosine of tau
96
00:04:06,300 --> 00:04:10,185
over square root of tau from
0 to some finite positive time
97
00:04:10,185 --> 00:04:12,800
point.
98
00:04:12,800 --> 00:04:15,260
As tau goes to 0, cosine
goes to 1 and square root
99
00:04:15,260 --> 00:04:17,430
of tau diverges.
100
00:04:17,430 --> 00:04:20,390
But this integral actually
has a finite value.
101
00:04:20,390 --> 00:04:21,980
This 1 over square
root of tau, that's
102
00:04:21,980 --> 00:04:25,650
an integrable singularity.
103
00:04:25,650 --> 00:04:26,670
So how do you do this?
104
00:04:26,670 --> 00:04:29,730
Well, you want to split the
domain again into two parts.
105
00:04:29,730 --> 00:04:32,840
One part that contains
the integrable singularity
106
00:04:32,840 --> 00:04:35,270
and the other part,
which excludes it.
107
00:04:35,270 --> 00:04:38,210
And in the part that
contains the singularity
108
00:04:38,210 --> 00:04:41,540
we might do an asymptotic
expansion of the function,
109
00:04:41,540 --> 00:04:44,510
and integrate each of those
asymptotically accurate terms
110
00:04:44,510 --> 00:04:45,740
separately.
111
00:04:45,740 --> 00:04:48,530
So the integral of 1
over square root of tau
112
00:04:48,530 --> 00:04:51,680
is going to give us
2 tau 0 to the 1/2.
113
00:04:51,680 --> 00:04:54,230
The integral of
this ratio here is
114
00:04:54,230 --> 00:04:57,840
going to give us minus
1/2 t0 to the 5/2.
115
00:04:57,840 --> 00:05:00,600
And then we have an integral
over a finite domain,
116
00:05:00,600 --> 00:05:02,580
which doesn't include
the singularity.
117
00:05:02,580 --> 00:05:05,390
And the accuracy of this method
can be dictated by two things
118
00:05:05,390 --> 00:05:06,140
now.
119
00:05:06,140 --> 00:05:08,370
One is how accurately
can we do this integral.
120
00:05:08,370 --> 00:05:13,190
And the other is how accurate is
our asymptotic expansion here.
121
00:05:13,190 --> 00:05:16,130
So if we want higher
accuracy we need more terms.
122
00:05:16,130 --> 00:05:17,810
We might need to be
selective about how
123
00:05:17,810 --> 00:05:20,000
we choose the end
point for this domain
124
00:05:20,000 --> 00:05:21,489
in order to minimize the error.
125
00:05:21,489 --> 00:05:23,030
There are lots of
methods to do this.
126
00:05:23,030 --> 00:05:23,490
Sam.
127
00:05:23,490 --> 00:05:25,031
AUDIENCE: Will you
talk a little more
128
00:05:25,031 --> 00:05:28,360
about what typical
integral [INAUDIBLE] is?
129
00:05:28,360 --> 00:05:31,210
WILLIAM GREEN, JR.: Yeah.
130
00:05:31,210 --> 00:05:33,850
So you know these things.
131
00:05:33,850 --> 00:05:40,260
If I try to integrate 1
over x, as x goes to 0
132
00:05:40,260 --> 00:05:42,310
this thing will always diverge.
133
00:05:42,310 --> 00:05:46,180
You know the antiderivative of
1 over x is going to be the log.
134
00:05:46,180 --> 00:05:47,890
And the log diverges.
135
00:05:47,890 --> 00:05:53,950
So power is weaker than 1 over
x, like 1 over x to the 1/2.
136
00:05:53,950 --> 00:05:55,290
Don't diverge like the log.
137
00:05:55,290 --> 00:05:56,830
They actually go
to a finite value.
138
00:05:56,830 --> 00:05:59,770
The log is like the most weakly
singular function there is.
139
00:05:59,770 --> 00:06:02,800
And integrals over weaker
powers of 1 over x actually
140
00:06:02,800 --> 00:06:04,994
don't produce a
singularity anymore.
141
00:06:04,994 --> 00:06:06,910
They're what we call
integrable singularities.
142
00:06:06,910 --> 00:06:11,200
The function itself diverges
but its integral does not.
143
00:06:11,200 --> 00:06:15,100
So you usually have this
asymptotic power law type
144
00:06:15,100 --> 00:06:16,270
behavior.
145
00:06:16,270 --> 00:06:18,140
OK.
146
00:06:18,140 --> 00:06:20,100
So they come up all the time.
147
00:06:20,100 --> 00:06:23,910
We see these things
in places like, say,
148
00:06:23,910 --> 00:06:27,460
squeezing a thin film of
fluid between two plates.
149
00:06:27,460 --> 00:06:31,567
How long does it take for these
two plates to come together?
150
00:06:31,567 --> 00:06:33,900
You might say, well, as the
plates get closer and closer
151
00:06:33,900 --> 00:06:36,737
together the pressure
starts to diverge in the gap
152
00:06:36,737 --> 00:06:38,070
and they'll never come together.
153
00:06:38,070 --> 00:06:40,950
But maybe depending on
the geometry of the plates
154
00:06:40,950 --> 00:06:42,750
there can be a
sort of finite time
155
00:06:42,750 --> 00:06:44,190
at which they'll come together.
156
00:06:44,190 --> 00:06:45,273
Depends on their geometry.
157
00:06:48,557 --> 00:06:49,390
So that's one topic.
158
00:06:49,390 --> 00:06:51,473
I think that's useful
because these things come up
159
00:06:51,473 --> 00:06:52,030
all the time.
160
00:06:52,030 --> 00:06:55,330
And you can't actually handle
those cases with quadrature
161
00:06:55,330 --> 00:06:56,800
by itself.
162
00:06:56,800 --> 00:06:57,910
The function diverges.
163
00:06:57,910 --> 00:06:59,680
There's no polynomial
interpolation
164
00:06:59,680 --> 00:07:01,150
that's suitable to match it.
165
00:07:01,150 --> 00:07:02,777
So the quadrature
can't handle it.
166
00:07:02,777 --> 00:07:04,360
And if you're trying
to integrate over
167
00:07:04,360 --> 00:07:07,210
an infinite domain,
well, good luck.
168
00:07:07,210 --> 00:07:09,580
You're never going to be
able to fit enough points in
169
00:07:09,580 --> 00:07:14,440
to know that you've accurately
integrated out to infinity.
170
00:07:14,440 --> 00:07:17,260
The other thing I
want to recap here,
171
00:07:17,260 --> 00:07:19,750
and this is a problem
for you to try.
172
00:07:19,750 --> 00:07:23,890
So here we have a simple first
order ordinary differential
173
00:07:23,890 --> 00:07:24,390
equation.
174
00:07:24,390 --> 00:07:27,080
We want to use the implicit
Euler method to solve it.
175
00:07:27,080 --> 00:07:29,840
So can you give a
closed form solution
176
00:07:29,840 --> 00:07:33,260
for one step of the
implicit Euler method,
177
00:07:33,260 --> 00:07:36,620
or for n steps of the
implicit Euler method?
178
00:07:36,620 --> 00:07:38,300
What result does that produce?
179
00:07:38,300 --> 00:07:40,520
Can you compute this?
180
00:07:40,520 --> 00:07:48,556
Do you guys talk about backward
Euler methods for for ODE-IVPs?
181
00:07:48,556 --> 00:07:51,180
Backwards difference methods in
general and the backwards Euler
182
00:07:51,180 --> 00:07:54,600
method specifically?
183
00:07:54,600 --> 00:07:55,940
No.
184
00:07:55,940 --> 00:07:57,150
OK.
185
00:07:57,150 --> 00:08:01,770
Well, clearly I
should have been here.
186
00:08:01,770 --> 00:08:03,150
So what's the strategy?
187
00:08:03,150 --> 00:08:05,100
We have to approximate
the derivative.
188
00:08:05,100 --> 00:08:06,960
Yes?
189
00:08:06,960 --> 00:08:08,880
So the approximation for
the derivative we use
190
00:08:08,880 --> 00:08:11,160
is a finite difference
approximation.
191
00:08:11,160 --> 00:08:14,940
So we write dx dt
as x at a point k
192
00:08:14,940 --> 00:08:18,870
plus 1 minus x at the
point k divided by delta t.
193
00:08:18,870 --> 00:08:20,700
Backwards difference.
194
00:08:20,700 --> 00:08:24,930
We evaluate the right-hand
side of the ODE-IVP at k
195
00:08:24,930 --> 00:08:27,940
plus 1, xk plus 1.
196
00:08:27,940 --> 00:08:30,360
And then we solve for xk plus 1.
197
00:08:30,360 --> 00:08:32,159
We just substitute in here.
198
00:08:32,159 --> 00:08:33,690
Backwards difference
approximation
199
00:08:33,690 --> 00:08:35,320
for the differential.
200
00:08:35,320 --> 00:08:37,979
And xk plus 1 for
the right-hand side.
201
00:08:37,979 --> 00:08:41,010
And we find that xk
plus 1 will be 1 over 1
202
00:08:41,010 --> 00:08:46,950
minus delta t times
this lambda times xk.
203
00:08:46,950 --> 00:08:51,680
And if we iterate this from k
equals 0 up to some finite k,
204
00:08:51,680 --> 00:08:54,240
it's like multiplying
by powers of 1 over 1
205
00:08:54,240 --> 00:08:55,980
minus delta t times lambda.
206
00:08:55,980 --> 00:08:57,044
Yes?
207
00:08:57,044 --> 00:08:57,960
Does this makes sense?
208
00:09:01,240 --> 00:09:05,910
So our backwards Euler solution
to this fundamental problem
209
00:09:05,910 --> 00:09:08,386
is going to look like this.
210
00:09:08,386 --> 00:09:10,260
And a lot of times we
ask about the stability
211
00:09:10,260 --> 00:09:12,240
of these sorts of solutions.
212
00:09:12,240 --> 00:09:16,850
So if for example, the
solution to this equation
213
00:09:16,850 --> 00:09:20,000
was supposed to
decay to zero, we
214
00:09:20,000 --> 00:09:21,770
would expect our
numerical method
215
00:09:21,770 --> 00:09:23,890
to also yield results
which decay to zero.
216
00:09:23,890 --> 00:09:27,110
We would call that stable.
217
00:09:27,110 --> 00:09:29,660
So we all know under
what conditions
218
00:09:29,660 --> 00:09:34,662
does xk, for example, go to 0
as I change the product delta
219
00:09:34,662 --> 00:09:36,350
t times lambda.
220
00:09:36,350 --> 00:09:40,850
Lambda could be a real number or
it can be a complex number too.
221
00:09:40,850 --> 00:09:42,360
Doesn't really matter.
222
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?
223
00:09:46,830 --> 00:09:49,340
Well what needs to be true is
this quantity in parentheses
224
00:09:49,340 --> 00:09:50,960
needs to be smaller
than 1 because I'm
225
00:09:50,960 --> 00:09:52,400
taking lots of
products of things
226
00:09:52,400 --> 00:09:53,483
that are smaller than one.
227
00:09:53,483 --> 00:09:57,770
That will continue to shrink xk
from x 0 until it goes to zero.
228
00:09:57,770 --> 00:10:01,790
If this is smaller than 1,
then this quantity down here
229
00:10:01,790 --> 00:10:04,610
needs to be bigger than
1 in absolute value.
230
00:10:04,610 --> 00:10:07,010
And then I'll get solutions
that slowly decay to zero.
231
00:10:07,010 --> 00:10:10,490
So if I ask, for what values
of delta t times lambda
232
00:10:10,490 --> 00:10:13,770
is the absolute value of
this thing bigger than 1,
233
00:10:13,770 --> 00:10:17,245
and I allow lambda to be a
complex number, for example?
234
00:10:17,245 --> 00:10:20,110
I'll find out that's
true when 1 minus delta t
235
00:10:20,110 --> 00:10:22,490
times the real part
of lambda squared
236
00:10:22,490 --> 00:10:26,630
plus delta t times the
imaginary part of lambda squared
237
00:10:26,630 --> 00:10:27,530
is bigger than 1.
238
00:10:30,100 --> 00:10:33,440
And if I plot the real
part of delta t times
239
00:10:33,440 --> 00:10:37,630
lambda versus the imaginary
part of delta times lambda,
240
00:10:37,630 --> 00:10:41,160
this inequality is
represented by this pink zone
241
00:10:41,160 --> 00:10:42,470
on the outside of the circle.
242
00:10:42,470 --> 00:10:44,510
So any values of
delta t times lambda
243
00:10:44,510 --> 00:10:49,460
that live in this pink area
will give stable integration.
244
00:10:49,460 --> 00:10:54,500
And the solution xk, well, the
k to zero as k goes to infinity.
245
00:10:57,020 --> 00:10:59,530
So it makes sense?
246
00:10:59,530 --> 00:11:02,010
This is a little funny though.
247
00:11:02,010 --> 00:11:03,780
What does the solution
to this equation
248
00:11:03,780 --> 00:11:08,012
do when lambda is bigger than 1?
249
00:11:08,012 --> 00:11:09,470
When it's bigger
than 0, excuse me.
250
00:11:09,470 --> 00:11:10,886
When lambda is
bigger than 0, what
251
00:11:10,886 --> 00:11:14,860
does the solution
to this ODE-IVP do?
252
00:11:14,860 --> 00:11:15,560
Does it decay?
253
00:11:15,560 --> 00:11:16,640
AUDIENCE: [INAUDIBLE]
254
00:11:16,640 --> 00:11:17,973
WILLIAM GREEN, JR.: It blows up.
255
00:11:17,973 --> 00:11:19,100
It grows exponentially.
256
00:11:19,100 --> 00:11:21,920
But the numerical
method here, when
257
00:11:21,920 --> 00:11:25,070
I have lambda is
bigger than zero,
258
00:11:25,070 --> 00:11:28,659
can be stable and produce
decaying solutions.
259
00:11:28,659 --> 00:11:29,825
So this is sort of peculiar.
260
00:11:32,890 --> 00:11:35,290
So we're often concerned
when we solve ODE-IVPs
261
00:11:35,290 --> 00:11:37,640
and we try to use these
sorts of backwards difference
262
00:11:37,640 --> 00:11:40,880
methods, these implicit
methods, to solve them.
263
00:11:40,880 --> 00:11:44,390
We're often concerned
with getting stability.
264
00:11:44,390 --> 00:11:46,299
We want to get
solutions that decay
265
00:11:46,299 --> 00:11:47,590
when they're supposed to decay.
266
00:11:47,590 --> 00:11:49,340
Sometimes we get
solutions that decay even
267
00:11:49,340 --> 00:11:52,690
when the solution is
supposed to blow up.
268
00:11:52,690 --> 00:11:55,240
So that's sort of funny.
269
00:11:55,240 --> 00:11:57,835
So the exact solution
should look like this.
270
00:11:57,835 --> 00:11:59,710
But there are going to
be circumstances where
271
00:11:59,710 --> 00:12:00,940
lambda is bigger than zero.
272
00:12:00,940 --> 00:12:02,530
In fact, lambda is bigger than--
273
00:12:02,530 --> 00:12:06,340
lambda delta t is bigger
than 2 along the real axis,
274
00:12:06,340 --> 00:12:11,060
where the solution will actually
blow up rather than decay away.
275
00:12:11,060 --> 00:12:12,970
So stability and
accuracy don't always
276
00:12:12,970 --> 00:12:15,760
correlate with each other.
277
00:12:15,760 --> 00:12:18,760
We really have to understand
the underlying dynamics
278
00:12:18,760 --> 00:12:21,040
of the problem we're
trying to solve
279
00:12:21,040 --> 00:12:24,110
before we choose a numerical
method to apply to it.
280
00:12:24,110 --> 00:12:27,280
And this is going to be true
of differential algebraic
281
00:12:27,280 --> 00:12:29,740
equations as well.
282
00:12:29,740 --> 00:12:33,310
Are there any
questions about this?
283
00:12:33,310 --> 00:12:37,030
There's a discussion of these
things in your textbook,
284
00:12:37,030 --> 00:12:39,160
these sorts of stability
regions associated
285
00:12:39,160 --> 00:12:41,290
with different backwards
difference methods.
286
00:12:41,290 --> 00:12:46,450
The implicit Euler method is one
backwards difference formula.
287
00:12:46,450 --> 00:12:48,760
It's the simplest one.
288
00:12:48,760 --> 00:12:52,160
It's the least accurate
one that you can derive
289
00:12:52,160 --> 00:12:53,420
but it's an easy one to use.
290
00:12:56,074 --> 00:12:58,490
But now we want to move on to
a different sort of problem.
291
00:12:58,490 --> 00:13:01,737
And it's one that we come across
in engineering all the time.
292
00:13:01,737 --> 00:13:03,320
These problems are
called differential
293
00:13:03,320 --> 00:13:04,670
algebraic equations.
294
00:13:04,670 --> 00:13:08,270
And they're problems of the
sort a function of some state
295
00:13:08,270 --> 00:13:12,500
variables, and the derivatives
of those state variables
296
00:13:12,500 --> 00:13:15,147
and time is equal to zero.
297
00:13:15,147 --> 00:13:16,730
And there are some
initial conditions.
298
00:13:16,730 --> 00:13:19,825
The state variables at that
time 0 is equal to some x0.
299
00:13:24,840 --> 00:13:26,490
So this is a vector
valued function
300
00:13:26,490 --> 00:13:29,320
of a vector valued state
and it's time derivatives.
301
00:13:29,320 --> 00:13:31,650
And we want to understand
the dynamics of this system.
302
00:13:31,650 --> 00:13:33,660
How does x vary
with time in a way
303
00:13:33,660 --> 00:13:37,500
that's consistent with
this governing equation?
304
00:13:37,500 --> 00:13:39,840
Usually we call these
well-posed problems
305
00:13:39,840 --> 00:13:42,690
when the dimension of the
state variable and the function
306
00:13:42,690 --> 00:13:43,590
are the same.
307
00:13:43,590 --> 00:13:46,730
They both live in
the vector space rn.
308
00:13:46,730 --> 00:13:48,375
So we have n different states.
309
00:13:48,375 --> 00:13:50,550
I have n different functions
that I have to satisfy.
310
00:13:50,550 --> 00:13:55,090
And I want to understand
how x varies with time.
311
00:13:55,090 --> 00:13:57,320
Here's an example.
312
00:13:57,320 --> 00:13:58,200
A stirred tank.
313
00:13:58,200 --> 00:14:00,460
I call this stirred
tank, example one.
314
00:14:00,460 --> 00:14:06,684
So into the stirred tank comes
some concentration of solute.
315
00:14:06,684 --> 00:14:09,100
Out of the stirred tank comes
some different concentration
316
00:14:09,100 --> 00:14:09,891
of the same solute.
317
00:14:09,891 --> 00:14:14,120
And we want to model the
dynamics as a function of time.
318
00:14:14,120 --> 00:14:16,900
So we know that a
material balance
319
00:14:16,900 --> 00:14:19,390
of the solute on the stirred
tank, if it's well mixed,
320
00:14:19,390 --> 00:14:26,320
will tell us that d c2 dt is
related to the volumetric flow
321
00:14:26,320 --> 00:14:28,960
rate into the tank divided
by the volume of the tank
322
00:14:28,960 --> 00:14:32,140
multiplied by the difference
between c1 and c2,
323
00:14:32,140 --> 00:14:37,112
the amount carried in and the
amount carried out of the tank.
324
00:14:37,112 --> 00:14:38,890
And let's put a
constraint in here too.
325
00:14:38,890 --> 00:14:42,430
Let's suppose we're trying to
control this tank in some way.
326
00:14:42,430 --> 00:14:45,820
So we say that c1
of time has got
327
00:14:45,820 --> 00:14:49,190
to be equal to some
control signal, gamma of t.
328
00:14:52,390 --> 00:14:55,810
And then we want to understand
the dynamics of this system.
329
00:14:55,810 --> 00:14:59,710
At time 0, c2 is some c0,
the initial concentration
330
00:14:59,710 --> 00:15:00,570
of the tank.
331
00:15:00,570 --> 00:15:03,490
And c1 is whatever the initial
value of this control signal
332
00:15:03,490 --> 00:15:04,960
is.
333
00:15:04,960 --> 00:15:08,350
And the solution is
c1 equals gamma of t.
334
00:15:08,350 --> 00:15:11,470
And c2 is, well, the
initial concentration's
335
00:15:11,470 --> 00:15:13,090
got to decay out exponentially.
336
00:15:13,090 --> 00:15:15,220
And then I've got to
integrate over how
337
00:15:15,220 --> 00:15:17,770
the input is changing in time.
338
00:15:17,770 --> 00:15:20,060
Those also produce
exponential decay.
339
00:15:20,060 --> 00:15:22,870
So if the input is a
little delta function spike
340
00:15:22,870 --> 00:15:25,330
then I'll get an exponential
decay leaving the tank.
341
00:15:25,330 --> 00:15:27,310
If the input is changing
dynamically in time,
342
00:15:27,310 --> 00:15:30,550
I need this convolution of
the input with this integral.
343
00:15:30,550 --> 00:15:33,130
So you just to
solve this system,
344
00:15:33,130 --> 00:15:34,900
this ordinary
differential equation
345
00:15:34,900 --> 00:15:37,310
by substituting c1 in there.
346
00:15:37,310 --> 00:15:40,280
But this system of
equations here is not just
347
00:15:40,280 --> 00:15:42,040
a differential equation.
348
00:15:42,040 --> 00:15:45,847
It's a differential equation
and an algebraic equation.
349
00:15:45,847 --> 00:15:47,680
There's something
peculiar about that that's
350
00:15:47,680 --> 00:15:50,140
different than just
solving systems
351
00:15:50,140 --> 00:15:52,440
of ordinary
differential equations.
352
00:15:52,440 --> 00:15:57,400
Oftentimes in models,
we write down all
353
00:15:57,400 --> 00:15:59,720
the governing equations
associated with the model.
354
00:15:59,720 --> 00:16:01,636
Some of them may be
differential, some of them
355
00:16:01,636 --> 00:16:02,800
may be algebraic.
356
00:16:02,800 --> 00:16:04,960
There are times when we
can look at those equations
357
00:16:04,960 --> 00:16:07,480
and see clever substitutions
that we can make.
358
00:16:07,480 --> 00:16:10,840
But if you have a hundred
equations, or a thousand
359
00:16:10,840 --> 00:16:12,900
equations, or a
million equations,
360
00:16:12,900 --> 00:16:14,820
there's no way to
do that reliably.
361
00:16:14,820 --> 00:16:18,310
We just kind of wind up with a
system of ordinary differential
362
00:16:18,310 --> 00:16:21,730
equations and algebraic
equations mixed together.
363
00:16:21,730 --> 00:16:24,777
And we have to solve
these reliably.
364
00:16:24,777 --> 00:16:25,735
Here's another example.
365
00:16:28,650 --> 00:16:29,720
Oh, excuse me.
366
00:16:29,720 --> 00:16:31,800
So let me write
this in this form.
367
00:16:31,800 --> 00:16:34,110
So we said there should
be a vector valued
368
00:16:34,110 --> 00:16:38,460
function of the states and
the derivatives of the states.
369
00:16:38,460 --> 00:16:41,490
The states are the
concentration c1 and c2.
370
00:16:41,490 --> 00:16:45,090
And the derivative
here is dc 2 dt.
371
00:16:45,090 --> 00:16:46,720
So here's my vector
valued function.
372
00:16:46,720 --> 00:16:48,900
The first element is
this, and it's equal to 0.
373
00:16:48,900 --> 00:16:51,660
The second element is
this, and it's equal to 0.
374
00:16:51,660 --> 00:16:54,900
Here's my initial
condition over the states.
375
00:16:54,900 --> 00:16:57,930
So I can just transform
that list of equations
376
00:16:57,930 --> 00:17:00,450
to the sort of canonical
form for a differential
377
00:17:00,450 --> 00:17:01,420
algebraic equation.
378
00:17:04,690 --> 00:17:07,619
Here's a separate example.
379
00:17:07,619 --> 00:17:12,790
Suppose instead, I'm
trying to either control
380
00:17:12,790 --> 00:17:16,450
or make a measurement of the
outlet concentrations c2.
381
00:17:16,450 --> 00:17:18,220
And we call that gamma of t.
382
00:17:22,060 --> 00:17:26,329
And I want to know now
what is c2 and what is c1?
383
00:17:26,329 --> 00:17:29,340
So I've got to solve this system
of equations for c2 and c1.
384
00:17:29,340 --> 00:17:31,740
The solution for c2 is easy.
385
00:17:31,740 --> 00:17:32,360
I know that.
386
00:17:32,360 --> 00:17:35,060
c2 is gamma, by definition.
387
00:17:35,060 --> 00:17:37,430
I got to substitute this
into the first equation
388
00:17:37,430 --> 00:17:38,420
and solve for c1.
389
00:17:38,420 --> 00:17:43,199
And I see, c1 is gamma
plus v over q gamma dot.
390
00:17:43,199 --> 00:17:45,240
And I have to have initial
conditions to go along
391
00:17:45,240 --> 00:17:47,992
with this.
392
00:17:47,992 --> 00:17:49,200
There's something funny here.
393
00:17:49,200 --> 00:17:52,800
We can see the initial condition
for c2 has to be gamma 0.
394
00:17:52,800 --> 00:17:57,710
The initial condition for c1,
what does that have to be?
395
00:17:57,710 --> 00:18:01,430
Well, it needs to be consistent
with the solution here.
396
00:18:01,430 --> 00:18:03,800
There were no free
variables when
397
00:18:03,800 --> 00:18:05,484
I solved this equation for c1.
398
00:18:05,484 --> 00:18:07,400
It isn't like solving a
differential equation,
399
00:18:07,400 --> 00:18:10,220
having a free
coefficient to specify.
400
00:18:10,220 --> 00:18:13,130
So it better be the
case that this c0 here
401
00:18:13,130 --> 00:18:18,680
is the same as gamma 0
plus v over q gamma dot 0.
402
00:18:18,680 --> 00:18:21,410
Somehow I have to know this
input signal to prescribe
403
00:18:21,410 --> 00:18:23,030
the initial conditions for c1.
404
00:18:23,030 --> 00:18:28,805
So that's peculiar and
different from just ODE-IVPs.
405
00:18:28,805 --> 00:18:29,680
You see this picture?
406
00:18:29,680 --> 00:18:32,940
You see how this goes together?
407
00:18:32,940 --> 00:18:34,630
The solution is funny too.
408
00:18:34,630 --> 00:18:36,880
Suppose we are trying to
use this system of equations
409
00:18:36,880 --> 00:18:39,040
to do the following.
410
00:18:39,040 --> 00:18:43,720
We are trying to measure c2
and use the measurement of c2
411
00:18:43,720 --> 00:18:46,030
to predict c1.
412
00:18:46,030 --> 00:18:48,580
So gamma is the
measurement, and we're
413
00:18:48,580 --> 00:18:52,060
trying to solve this system
of differential algebraic
414
00:18:52,060 --> 00:18:55,120
equations to predict what c1 is.
415
00:18:55,120 --> 00:18:58,300
All measurements
incur numerical error.
416
00:18:58,300 --> 00:19:01,560
So even though our signal
for gamma that we measure
417
00:19:01,560 --> 00:19:06,670
may be continuous, it's
bouncing around wildly.
418
00:19:06,670 --> 00:19:08,290
It's not a constant signal.
419
00:19:08,290 --> 00:19:10,900
It's going to move
around a great deal.
420
00:19:10,900 --> 00:19:13,960
And c1, our predicted
value for c1,
421
00:19:13,960 --> 00:19:15,670
depends on the
derivatives of gamma,
422
00:19:15,670 --> 00:19:18,490
which means c1 is going to be
incredibly sensitive to how
423
00:19:18,490 --> 00:19:21,440
that measurement is made.
424
00:19:21,440 --> 00:19:22,395
So it's peculiar.
425
00:19:22,395 --> 00:19:24,520
It means that there's going
to be a lot of problems
426
00:19:24,520 --> 00:19:28,840
with stably solving these
equations because the solutions
427
00:19:28,840 --> 00:19:34,390
can admit cases where they're
not integrals of the input
428
00:19:34,390 --> 00:19:37,990
but they're actually related
to derivatives of the input.
429
00:19:37,990 --> 00:19:39,960
Look at the solution
back here again.
430
00:19:39,960 --> 00:19:43,140
Oop, sorry.
431
00:19:43,140 --> 00:19:45,960
The solution here was related
to an integral of the input.
432
00:19:45,960 --> 00:19:48,030
Suppose I was making
a measurement of c1
433
00:19:48,030 --> 00:19:50,940
in trying to predict c2 instead?
434
00:19:50,940 --> 00:19:53,480
My prediction for c2 would
be related to the integral
435
00:19:53,480 --> 00:19:54,480
of that measured signal.
436
00:19:54,480 --> 00:19:56,850
And we know integrals
smooth signals out.
437
00:19:56,850 --> 00:19:59,340
So it's not going to be so
sensitive to the measurement.
438
00:19:59,340 --> 00:20:01,740
So one way of doing this
seems really sensible.
439
00:20:01,740 --> 00:20:06,540
And the other way of doing
it seems a little bizarre.
440
00:20:06,540 --> 00:20:09,317
Well, we can construct these
equations however we want.
441
00:20:09,317 --> 00:20:11,400
So there's something about
how we formulate models
442
00:20:11,400 --> 00:20:13,140
that are going to make
them either easily
443
00:20:13,140 --> 00:20:18,850
solvable, as DAEs, or
quite difficult to solve.
444
00:20:18,850 --> 00:20:23,190
So these pop up all the time
in engineering problems.
445
00:20:23,190 --> 00:20:25,530
They pop up in
dynamic simulations,
446
00:20:25,530 --> 00:20:27,900
where we have some sort
of underlying conservation
447
00:20:27,900 --> 00:20:31,950
principle, or constraints,
or equilibria.
448
00:20:31,950 --> 00:20:33,450
So if we have to
conserve something,
449
00:20:33,450 --> 00:20:38,550
like total energy, total mass,
total momentum, total particle
450
00:20:38,550 --> 00:20:41,760
number, total atomic
species, total charge,
451
00:20:41,760 --> 00:20:44,040
there will always be an
algebraic constraint that
452
00:20:44,040 --> 00:20:48,210
tells us the total
mass has to stay this,
453
00:20:48,210 --> 00:20:51,500
or the total number of atoms of
a certain type has to say this.
454
00:20:51,500 --> 00:20:52,230
Yeah, Ian.
455
00:20:52,230 --> 00:20:56,950
AUDIENCE: Just to be clear, were
you saying, on the stirred tank
456
00:20:56,950 --> 00:20:59,722
example, that those
were the same problem
457
00:20:59,722 --> 00:21:02,119
with different
solution approaches?
458
00:21:02,119 --> 00:21:03,660
WILLIAM GREEN, JR.:
So good question.
459
00:21:03,660 --> 00:21:04,285
Let me go back.
460
00:21:07,260 --> 00:21:10,960
So they are fundamentally
different problems.
461
00:21:10,960 --> 00:21:17,450
So in one case I specified the
value of c2 to be this gamma.
462
00:21:17,450 --> 00:21:21,450
I tried to make a measurement
of c2 and then predict c1.
463
00:21:21,450 --> 00:21:24,180
In the other case,
I specified c1.
464
00:21:24,180 --> 00:21:28,100
I tried to measure
c1 and predict c2.
465
00:21:28,100 --> 00:21:29,850
So one case I measured
the input and tried
466
00:21:29,850 --> 00:21:31,869
to predict the
output to the tank.
467
00:21:31,869 --> 00:21:33,660
And in the other case
I measured the output
468
00:21:33,660 --> 00:21:36,600
and tried to use it
to predict the input.
469
00:21:36,600 --> 00:21:37,273
Yeah?
470
00:21:37,273 --> 00:21:39,125
AUDIENCE: So physically
the same system
471
00:21:39,125 --> 00:21:40,790
with a choice by the operator.
472
00:21:40,790 --> 00:21:41,790
WILLIAM GREEN, JR.: Yes.
473
00:21:41,790 --> 00:21:42,569
That's right.
474
00:21:42,569 --> 00:21:43,110
That's right.
475
00:21:43,110 --> 00:21:45,270
So I formulated the
model differently.
476
00:21:45,270 --> 00:21:49,080
The same underlying
system but I formulated
477
00:21:49,080 --> 00:21:50,700
the model differently.
478
00:21:50,700 --> 00:21:52,950
And it led to completely
different behavior
479
00:21:52,950 --> 00:21:55,770
in the underlying solution
to the differential algebraic
480
00:21:55,770 --> 00:21:57,520
equations.
481
00:21:57,520 --> 00:21:59,160
So you've already
solved problems
482
00:21:59,160 --> 00:22:00,360
that involve conservation.
483
00:22:00,360 --> 00:22:04,770
You've done things like
ODE-IVPs on batch reactors,
484
00:22:04,770 --> 00:22:07,320
where the total amount of
material in the reactor
485
00:22:07,320 --> 00:22:09,015
had to remain constant.
486
00:22:09,015 --> 00:22:10,890
Instead, you probably
solved for the dynamics
487
00:22:10,890 --> 00:22:14,130
of each of the components
undergoing the reaction.
488
00:22:14,130 --> 00:22:18,540
You might have checked to
see if at each point in time
489
00:22:18,540 --> 00:22:21,200
the total mass in the
reactor stayed the same.
490
00:22:21,200 --> 00:22:23,754
Because you solved all these
differential equations,
491
00:22:23,754 --> 00:22:25,420
they were interconnected
with each other
492
00:22:25,420 --> 00:22:27,450
but they all incurred
some numerical error.
493
00:22:27,450 --> 00:22:29,370
It's possible that
numerical error
494
00:22:29,370 --> 00:22:33,150
accrues in a way that
actually has you losing mass
495
00:22:33,150 --> 00:22:35,430
from your system.
496
00:22:35,430 --> 00:22:37,230
So there may be benefits
to actually trying
497
00:22:37,230 --> 00:22:39,360
to solve these sorts
of systems of equations
498
00:22:39,360 --> 00:22:42,150
with, say, one less
differential equation
499
00:22:42,150 --> 00:22:45,390
for one of the components and
instead an algebraic constraint
500
00:22:45,390 --> 00:22:48,330
governing the total mass
or total number of moles
501
00:22:48,330 --> 00:22:49,050
in the system.
502
00:22:53,010 --> 00:22:55,320
So these pop up all the time.
503
00:22:55,320 --> 00:22:58,080
You see them in models
of reaction networks,
504
00:22:58,080 --> 00:23:00,540
where we try to use a pseudo
steady state approximation.
505
00:23:00,540 --> 00:23:03,960
We say some reactions go so
fast that they equilibrate.
506
00:23:03,960 --> 00:23:06,120
So they're not governed
by differential equations
507
00:23:06,120 --> 00:23:11,790
but by algebraic
equations for equilibria.
508
00:23:11,790 --> 00:23:14,520
They can pop up in models
of control, where we neglect
509
00:23:14,520 --> 00:23:16,680
the controller dynamics.
510
00:23:16,680 --> 00:23:19,440
So you say I'm going to
try to control this process
511
00:23:19,440 --> 00:23:23,010
and I get to instantaneously
turn on or off this valve
512
00:23:23,010 --> 00:23:24,930
in response to a measurement.
513
00:23:24,930 --> 00:23:28,472
Actually, controllers have
inherent dynamics in them.
514
00:23:28,472 --> 00:23:30,180
But if I don't put
those dynamics in then
515
00:23:30,180 --> 00:23:32,760
I may get an algebraic equation
instead of a differential
516
00:23:32,760 --> 00:23:34,755
equation for how the
control process occurs.
517
00:23:39,690 --> 00:23:44,250
There are different types
of DAEs that we talk about.
518
00:23:44,250 --> 00:23:48,640
So there's one type that's
called semi-explicit DAEs.
519
00:23:48,640 --> 00:23:50,980
And in these types of
differential algebraic
520
00:23:50,980 --> 00:23:53,710
equations we can
write them in the form
521
00:23:53,710 --> 00:24:00,310
M dx dt is equal to f of x
and t, some initial condition.
522
00:24:00,310 --> 00:24:03,410
You say that this
looks like an ODE-IVP.
523
00:24:03,410 --> 00:24:06,430
M is called the mass matrix.
524
00:24:06,430 --> 00:24:10,990
And it may or may not be the
case that M is invertible.
525
00:24:10,990 --> 00:24:16,150
M may or may not be
a full rank matrix.
526
00:24:16,150 --> 00:24:20,600
So let's look at that
stirred tank example one.
527
00:24:20,600 --> 00:24:24,890
This was the underlying
equation, f of x and dx dt
528
00:24:24,890 --> 00:24:25,790
and t.
529
00:24:25,790 --> 00:24:27,860
Can write it in
semi-explicit form.
530
00:24:27,860 --> 00:24:35,590
If c has two components, c1
and c2, then dc 2 dt-- oops,
531
00:24:35,590 --> 00:24:36,640
this is a typo here.
532
00:24:36,640 --> 00:24:37,810
I really apologize.
533
00:24:37,810 --> 00:24:42,180
This should be a 0, and
this element should be a 1.
534
00:24:42,180 --> 00:24:43,230
I'll correct that online.
535
00:24:43,230 --> 00:24:46,200
I apologize for that.
536
00:24:46,200 --> 00:24:48,709
So this should be
this matrix multiplied
537
00:24:48,709 --> 00:24:50,250
with this differential
should give me
538
00:24:50,250 --> 00:24:52,750
dc 2 dt, this term here.
539
00:24:52,750 --> 00:24:56,410
And a 0 for the second equation.
540
00:24:56,410 --> 00:25:02,450
And then I've got q over v
multiplying c1 and minus c2.
541
00:25:02,450 --> 00:25:05,010
So that gives me
the first line here.
542
00:25:05,010 --> 00:25:09,850
I've got minus c1 coming on
the second equation here.
543
00:25:09,850 --> 00:25:12,820
And those balance
with a 0 and a gamma.
544
00:25:12,820 --> 00:25:16,000
So the lower line reads
like this equation here
545
00:25:16,000 --> 00:25:17,650
and the upper line--
546
00:25:17,650 --> 00:25:19,660
fix this typo--
replace 0 with one,
547
00:25:19,660 --> 00:25:21,790
reads like the
upper equation here.
548
00:25:21,790 --> 00:25:24,950
This is the semi-explicit form
of the differential equation.
549
00:25:24,950 --> 00:25:27,970
You can see this mass
matrix is singular and has
550
00:25:27,970 --> 00:25:29,200
a row that's all zeros.
551
00:25:29,200 --> 00:25:31,000
There's no way to
invert this matrix
552
00:25:31,000 --> 00:25:34,540
and convert it into a system
of ODE-IVPs, which you already
553
00:25:34,540 --> 00:25:37,360
know how to solve.
554
00:25:37,360 --> 00:25:41,080
If the mass matrix is full
rank, which can happen,
555
00:25:41,080 --> 00:25:43,570
you can formulate your
model in such a way
556
00:25:43,570 --> 00:25:45,980
that it takes on this
semi-explicit form.
557
00:25:45,980 --> 00:25:48,530
But if the mass matrix is
full rank and invertible,
558
00:25:48,530 --> 00:25:50,730
then you just invert
the mass matrix.
559
00:25:50,730 --> 00:25:53,260
And now you can apply all
your ODE-IVP techniques
560
00:25:53,260 --> 00:25:56,710
to solve this thing.
561
00:25:56,710 --> 00:25:58,960
The mass matrix doesn't
need to be constant.
562
00:25:58,960 --> 00:26:01,720
Could depend on
the concentration.
563
00:26:01,720 --> 00:26:06,397
It can depend on the states
as well in this form.
564
00:26:06,397 --> 00:26:08,480
I just chose an example
where that isn't the case.
565
00:26:08,480 --> 00:26:10,345
But in general, it can
depend on the states.
566
00:26:10,345 --> 00:26:12,662
If it depends on the
states and we invert it,
567
00:26:12,662 --> 00:26:14,620
then we need to be sure
that the mass matrix is
568
00:26:14,620 --> 00:26:18,100
invertible for all values of
the states, which may or may not
569
00:26:18,100 --> 00:26:19,395
be easy to ensure.
570
00:26:24,700 --> 00:26:28,010
So here's what I'd
like you to try to do.
571
00:26:28,010 --> 00:26:30,850
So here's that stirred
tank example two.
572
00:26:30,850 --> 00:26:34,240
The only difference between
this and the previous one
573
00:26:34,240 --> 00:26:38,130
is that c2 now is set equal
to gamma instead of c1.
574
00:26:38,130 --> 00:26:40,430
I want you to try to write
it in semi-explicit form.
575
00:26:40,430 --> 00:26:42,790
A sort of test of understanding.
576
00:26:42,790 --> 00:26:46,700
Can you define the mass matrix?
577
00:26:46,700 --> 00:26:48,380
Can you write out
the right-hand side
578
00:26:48,380 --> 00:26:50,175
of the equation in
semi-explicit form?
579
00:26:50,175 --> 00:26:51,800
Take a second, work
with your neighbor.
580
00:26:51,800 --> 00:26:53,508
See if you can figure
out how to do that.
581
00:28:04,658 --> 00:28:06,020
AUDIENCE: Is this clear?
582
00:28:06,020 --> 00:28:07,296
OK to do this?
583
00:28:12,920 --> 00:28:15,420
WILLIAM GREEN, JR.: You can
look back on the previous slide.
584
00:28:15,420 --> 00:28:18,690
Actually, the semi-exclusive
form is going to be the same.
585
00:28:18,690 --> 00:28:20,910
This is a 0, that's a 1.
586
00:28:20,910 --> 00:28:23,850
The only difference is
going to be c2 is now
587
00:28:23,850 --> 00:28:25,970
balanced with gamma, the input.
588
00:28:25,970 --> 00:28:27,880
So this minus 1 is
got to shift over.
589
00:28:27,880 --> 00:28:29,460
So I switched the
0 for the minus 1,
590
00:28:29,460 --> 00:28:32,040
you got the second,
the semi-explicit form.
591
00:28:38,270 --> 00:28:43,140
There's another way that
these equations can pop up.
592
00:28:43,140 --> 00:28:44,880
So we have this mass matrix.
593
00:28:44,880 --> 00:28:47,790
And it might be
the case that M is
594
00:28:47,790 --> 00:28:51,450
diagonal over many of
the states and then
595
00:28:51,450 --> 00:28:52,789
0 for all the other states.
596
00:28:52,789 --> 00:28:55,080
So we might be able to write
the equation in this form.
597
00:28:55,080 --> 00:28:56,204
This comes up all the time.
598
00:28:56,204 --> 00:28:58,260
So we might be able
to write dx dt is
599
00:28:58,260 --> 00:29:01,110
a function of x, y, and t.
600
00:29:01,110 --> 00:29:06,950
And 0 is equal to another
function of x, y, and t.
601
00:29:06,950 --> 00:29:09,300
And we have some
initial conditions
602
00:29:09,300 --> 00:29:13,410
prescribed for the
x's, of which we're
603
00:29:13,410 --> 00:29:16,560
taking the differential
in one of these equations.
604
00:29:16,560 --> 00:29:19,890
And we've also got to satisfy
those initial conditions
605
00:29:19,890 --> 00:29:23,400
for these variables y associated
with the second nonlinear
606
00:29:23,400 --> 00:29:24,390
equation.
607
00:29:24,390 --> 00:29:27,750
The x's are called the
differential states
608
00:29:27,750 --> 00:29:30,144
because they're involved
in equations where they're
609
00:29:30,144 --> 00:29:32,310
difference with respect to
time, or the differential
610
00:29:32,310 --> 00:29:34,260
with respect to time is taken.
611
00:29:34,260 --> 00:29:37,740
y are called the
algebraic states
612
00:29:37,740 --> 00:29:41,280
because they only appear
as themselves and not
613
00:29:41,280 --> 00:29:44,470
as their time differential
in any of these equations.
614
00:29:44,470 --> 00:29:47,590
So we might want to solve
for both the differential
615
00:29:47,590 --> 00:29:50,040
and the algebraic
states simultaneously.
616
00:29:50,040 --> 00:29:52,710
We have to know these
algebraic states to determine
617
00:29:52,710 --> 00:29:54,024
the differential states.
618
00:29:54,024 --> 00:29:56,190
We have to know the
differential states to determine
619
00:29:56,190 --> 00:29:57,120
the algebraic states.
620
00:30:00,771 --> 00:30:01,740
AUDIENCE: Professor?
621
00:30:01,740 --> 00:30:02,739
WILLIAM GREEN, JR.: Yes.
622
00:30:02,739 --> 00:30:07,695
AUDIENCE: [INAUDIBLE]
623
00:30:07,695 --> 00:30:09,070
WILLIAM GREEN,
JR.: That's right.
624
00:30:09,070 --> 00:30:13,370
So here the x that I have
is not all of the states.
625
00:30:13,370 --> 00:30:14,690
It's not all of the unknowns.
626
00:30:14,690 --> 00:30:17,390
It's only the set of
the unknowns of which I
627
00:30:17,390 --> 00:30:18,576
have to take a differential.
628
00:30:18,576 --> 00:30:19,700
That's right, that's right.
629
00:30:25,040 --> 00:30:29,030
So we've sort of tooled
down the kinds of equations
630
00:30:29,030 --> 00:30:30,810
we might want to look at.
631
00:30:30,810 --> 00:30:33,860
We've gone from the generic
differential algebraic equation
632
00:30:33,860 --> 00:30:38,390
to a semi-explicit form to the
splitting, this differential
633
00:30:38,390 --> 00:30:39,320
algebraic splitting.
634
00:30:39,320 --> 00:30:43,400
Many problems can be
formulated in this fashion.
635
00:30:43,400 --> 00:30:45,251
This one's the
easiest one to think
636
00:30:45,251 --> 00:30:47,000
about so that's the
one that we'll discuss
637
00:30:47,000 --> 00:30:48,530
for most of the lecture.
638
00:30:51,740 --> 00:30:53,720
So let's look at some examples.
639
00:30:53,720 --> 00:30:57,340
I think examples are
interesting to think about.
640
00:30:57,340 --> 00:30:58,420
So here's one.
641
00:30:58,420 --> 00:31:02,140
Think about a mass
spring system.
642
00:31:02,140 --> 00:31:04,360
There's an algebraic
equation that
643
00:31:04,360 --> 00:31:07,270
governs a conserved
quantity, the total energy
644
00:31:07,270 --> 00:31:09,370
of this mass spring system.
645
00:31:09,370 --> 00:31:13,060
So the kinetic energy, 1/2 m
times the velocity of the mass
646
00:31:13,060 --> 00:31:14,200
squared.
647
00:31:14,200 --> 00:31:18,444
Plus the energy stored in
the spring, 1/2 kx squared
648
00:31:18,444 --> 00:31:20,110
has got to be equal
to the total energy.
649
00:31:20,110 --> 00:31:22,670
And the total energy
can change in time.
650
00:31:22,670 --> 00:31:27,010
So we have a differential
algebraic equation.
651
00:31:27,010 --> 00:31:29,860
We have a differential
state because we're taking
652
00:31:29,860 --> 00:31:31,420
the time derivative of x.
653
00:31:31,420 --> 00:31:33,610
And we'd like to be
able to determine
654
00:31:33,610 --> 00:31:36,040
x as a function of time.
655
00:31:36,040 --> 00:31:38,920
So f of x, x dot and t.
656
00:31:38,920 --> 00:31:42,250
That's this equation, 1/2
mv squared plus one half kx
657
00:31:42,250 --> 00:31:45,254
squared minus E is equal to 0.
658
00:31:45,254 --> 00:31:47,170
We'd like to solve this
differential algebraic
659
00:31:47,170 --> 00:31:48,070
equation.
660
00:31:48,070 --> 00:31:49,300
It's well-posed.
661
00:31:49,300 --> 00:31:54,580
We have one equation for
one unknown, one state.
662
00:31:54,580 --> 00:31:58,260
The equation has
a solution, which
663
00:31:58,260 --> 00:32:00,520
is a times cosine omega t.
664
00:32:00,520 --> 00:32:02,020
It's not the only
possible solution.
665
00:32:02,020 --> 00:32:04,290
It's a solution
to this equation.
666
00:32:04,290 --> 00:32:06,660
We can substitute that
solution into the equation
667
00:32:06,660 --> 00:32:10,560
and determine, for example
what this oscillation frequency
668
00:32:10,560 --> 00:32:11,490
omega has to be.
669
00:32:11,490 --> 00:32:12,810
It's square root of k over m.
670
00:32:12,810 --> 00:32:14,250
Or what the amplitude has to be.
671
00:32:14,250 --> 00:32:17,130
It's the square
root of E over k.
672
00:32:17,130 --> 00:32:20,496
So you've seen this
from physics before.
673
00:32:20,496 --> 00:32:21,870
We've already
solved differential
674
00:32:21,870 --> 00:32:24,030
algebraic equations.
675
00:32:24,030 --> 00:32:26,520
The thing you saw before
though, was actually
676
00:32:26,520 --> 00:32:28,500
conservation of momentum.
677
00:32:28,500 --> 00:32:32,190
You converted this
entirely to a second order
678
00:32:32,190 --> 00:32:34,380
differential equation
in x and solved it.
679
00:32:34,380 --> 00:32:37,980
But in principle, we could have
tried to find various solutions
680
00:32:37,980 --> 00:32:39,150
to this equation instead.
681
00:32:39,150 --> 00:32:41,802
We would have said, well, we've
observed masses and springs
682
00:32:41,802 --> 00:32:43,260
and we know they
oscillate in time.
683
00:32:43,260 --> 00:32:45,450
So let's propose some
oscillatory solutions
684
00:32:45,450 --> 00:32:48,300
and see what values of
the amplitude and omega we
685
00:32:48,300 --> 00:32:50,356
need to satisfy this equation.
686
00:32:50,356 --> 00:32:52,230
It would have been a
perfectly acceptable way
687
00:32:52,230 --> 00:32:53,280
of solving this problem.
688
00:33:01,010 --> 00:33:03,940
The thing is, this
mass spring problem,
689
00:33:03,940 --> 00:33:05,980
the differential
algebraic equation
690
00:33:05,980 --> 00:33:09,880
contains nonlinearities
with respect to the states.
691
00:33:09,880 --> 00:33:11,470
So it's nonlinear
in the velocity
692
00:33:11,470 --> 00:33:12,940
and it's nonlinear
in the position.
693
00:33:12,940 --> 00:33:15,865
Nonlinear equations in
general are hard to solve.
694
00:33:15,865 --> 00:33:17,740
When you converted that
to a momentum balance
695
00:33:17,740 --> 00:33:20,110
you got linear equations.
696
00:33:20,110 --> 00:33:23,590
You got the second
derivative of position
697
00:33:23,590 --> 00:33:27,160
was equal to stiffness
over mass times position.
698
00:33:27,160 --> 00:33:29,290
It was a linear differential
equation to solve,
699
00:33:29,290 --> 00:33:30,300
which is easy to solve.
700
00:33:30,300 --> 00:33:32,950
So you've got higher
order differentials,
701
00:33:32,950 --> 00:33:36,310
but you linearize the equation
so it made it easier to solve,
702
00:33:36,310 --> 00:33:38,740
in a certain sense.
703
00:33:38,740 --> 00:33:42,790
But we'd like to
sometimes understand what
704
00:33:42,790 --> 00:33:44,245
we can do with these equations.
705
00:33:44,245 --> 00:33:45,970
There are certain
limiting cases where
706
00:33:45,970 --> 00:33:48,053
we can do interesting
things with these equations.
707
00:33:48,053 --> 00:33:49,600
I'm going to show you one now.
708
00:33:49,600 --> 00:33:53,650
So in the case where
the Jacobian of this
709
00:33:53,650 --> 00:33:58,021
function f with respect to
the differential of the state
710
00:33:58,021 --> 00:33:58,520
variable.
711
00:33:58,520 --> 00:34:00,340
So take the derivative
of f with respect
712
00:34:00,340 --> 00:34:04,080
to x dot and hold time
in the state constant.
713
00:34:04,080 --> 00:34:06,830
So this is a matrix,
one that's full rank.
714
00:34:06,830 --> 00:34:10,330
Then the DA be represented
as an equivalent to ODE.
715
00:34:10,330 --> 00:34:12,190
Here's how you can show that.
716
00:34:12,190 --> 00:34:16,389
The total differential of f
is this Jacobian with respect
717
00:34:16,389 --> 00:34:21,969
to x dot times the change x dot
plus the Jacobian with respect
718
00:34:21,969 --> 00:34:24,760
to x times the change in x.
719
00:34:24,760 --> 00:34:28,400
Plus the partial derivative
of f with respect to time,
720
00:34:28,400 --> 00:34:30,230
times the change in time.
721
00:34:30,230 --> 00:34:32,230
And this total differential
has to be equal to 0
722
00:34:32,230 --> 00:34:37,080
because f is equal to 0
for all states in time.
723
00:34:37,080 --> 00:34:39,949
So let's divide by dt.
724
00:34:39,949 --> 00:34:43,050
We'd like to know how
f changes with time.
725
00:34:43,050 --> 00:34:45,139
So total differential
of f with time.
726
00:34:45,139 --> 00:34:48,300
We divide each of these
little differentials by time.
727
00:34:48,300 --> 00:34:53,060
And then we solve for dx dot dt.
728
00:34:53,060 --> 00:34:55,060
That's going to involve
moving this term
729
00:34:55,060 --> 00:34:56,469
to the other side
of the equation
730
00:34:56,469 --> 00:34:58,340
and inverting this Jacobian.
731
00:34:58,340 --> 00:35:00,910
So if I can invert
the Jacobian then
732
00:35:00,910 --> 00:35:05,125
I can convert my DAE system
into a system of ODEs.
733
00:35:08,197 --> 00:35:09,530
This equals 0 shouldn't be here.
734
00:35:09,530 --> 00:35:10,050
I apologize.
735
00:35:10,050 --> 00:35:13,120
That's another typo.
736
00:35:13,120 --> 00:35:15,920
I move this term to the
other side of the equation
737
00:35:15,920 --> 00:35:19,580
to where the 0 is and I solve.
738
00:35:19,580 --> 00:35:24,180
So I went from a first order
equation, a first order
739
00:35:24,180 --> 00:35:27,330
equation in x to a second
order equation in x.
740
00:35:27,330 --> 00:35:30,120
Two differentials of x with
time instead of one differential
741
00:35:30,120 --> 00:35:30,900
of x with time.
742
00:35:30,900 --> 00:35:32,525
But I could convert
it to an equivalent
743
00:35:32,525 --> 00:35:34,990
system of ordinary
differential equations.
744
00:35:34,990 --> 00:35:36,836
Does that make sense?
745
00:35:36,836 --> 00:35:39,220
Do you see how that works?
746
00:35:39,220 --> 00:35:40,776
Any questions?
747
00:35:40,776 --> 00:35:44,860
Let's go back and
look at an example.
748
00:35:44,860 --> 00:35:45,900
The mass spring system.
749
00:35:49,430 --> 00:35:53,810
So here's our ODE that
we're trying to solve.
750
00:35:53,810 --> 00:35:56,360
I'm sorry, our DAE that
we're trying to solve.
751
00:35:56,360 --> 00:36:00,440
Can you convert this to
an equivalent ordinary
752
00:36:00,440 --> 00:36:04,030
differential equation using
the approach I just described?
753
00:36:04,030 --> 00:36:06,900
Think you guys can
figure that out?
754
00:36:06,900 --> 00:36:07,400
Try it.
755
00:36:07,400 --> 00:36:10,555
You can work together.
756
00:36:10,555 --> 00:37:07,000
[SIDE CONVERSATIONS]
757
00:37:07,000 --> 00:37:10,070
WILLIAM GREEN, JR.:
Guys are either tired
758
00:37:10,070 --> 00:37:11,880
or friendships
have been destroyed
759
00:37:11,880 --> 00:37:15,150
over the mid-term season.
760
00:37:15,150 --> 00:37:15,770
Maybe both.
761
00:37:38,757 --> 00:37:40,590
So let's compute the
things we need to know.
762
00:37:40,590 --> 00:37:44,400
We needed to know the partials
of this f with respect
763
00:37:44,400 --> 00:37:46,180
to the states and
with respect to time.
764
00:37:46,180 --> 00:37:50,190
So partial f respect to x dot,
holding all the other variables
765
00:37:50,190 --> 00:37:51,840
constant, is mx dot.
766
00:37:51,840 --> 00:37:55,120
Partial f with
respect to x is kx.
767
00:37:55,120 --> 00:37:58,591
Partial f with
respect to t is t.
768
00:37:58,591 --> 00:38:00,090
And I told you
before that you could
769
00:38:00,090 --> 00:38:05,645
write this says dx
dot dt is equal to df
770
00:38:05,645 --> 00:38:10,670
dx dot inverse multiplied by--
771
00:38:10,670 --> 00:38:14,264
is equal to minus df dx dot
inverse multiplied by df dx.
772
00:38:14,264 --> 00:38:16,430
This was on the previous
slide, which just gives you
773
00:38:16,430 --> 00:38:19,240
conservation of momentum.
774
00:38:19,240 --> 00:38:25,390
This gives you the acceleration
is equal to the force minus kx
775
00:38:25,390 --> 00:38:26,870
divided by the mass m.
776
00:38:31,300 --> 00:38:33,875
So we can solve either,
the DAE or the ODE.
777
00:38:37,890 --> 00:38:40,220
Here's a more
complicated example,
778
00:38:40,220 --> 00:38:43,320
but it works, well,
it's more complicated.
779
00:38:43,320 --> 00:38:46,010
So a lot of times we do
molecular dynamics simulations.
780
00:38:46,010 --> 00:38:48,130
We try to model
molecules moving around.
781
00:38:48,130 --> 00:38:50,232
There we also want
to conserve energy.
782
00:38:50,232 --> 00:38:52,190
There's usually a conserved
quantity, something
783
00:38:52,190 --> 00:38:53,870
we're trying to hold constant.
784
00:38:53,870 --> 00:38:55,160
So suppose we try to do that.
785
00:38:55,160 --> 00:38:57,770
The total energy in this
system of molecules,
786
00:38:57,770 --> 00:39:01,610
maybe it's 1/2 mass times
the velocity squared.
787
00:39:01,610 --> 00:39:04,010
Now x is the position
of all these molecules.
788
00:39:04,010 --> 00:39:05,990
Ad x dot is the velocity
of all these molecules
789
00:39:05,990 --> 00:39:09,650
plus the potential energy,
which is a function of f,
790
00:39:09,650 --> 00:39:12,510
a function of their positions.
791
00:39:12,510 --> 00:39:17,391
So the DAE representing
this constraint is this.
792
00:39:17,391 --> 00:39:18,890
It's the same as
for the mass spring
793
00:39:18,890 --> 00:39:20,931
system, where we've replaced
the potential energy
794
00:39:20,931 --> 00:39:21,980
with a generic one.
795
00:39:21,980 --> 00:39:24,520
This actually isn't
a well-posed DEA.
796
00:39:24,520 --> 00:39:28,430
We have one equation
for, let's see,
797
00:39:28,430 --> 00:39:31,146
if we have n molecules in
here we have three n states.
798
00:39:31,146 --> 00:39:32,770
We can move around
in three dimensions.
799
00:39:32,770 --> 00:39:35,605
So this is a really
ill-posed equation.
800
00:39:35,605 --> 00:39:37,730
If we try to apply the same
procedure, then we say,
801
00:39:37,730 --> 00:39:40,040
well, this quantity still
should be conserved.
802
00:39:40,040 --> 00:39:43,180
Over time it needs
to be equal to zero.
803
00:39:43,180 --> 00:39:47,630
Then we can try to compute
all these differentials.
804
00:39:47,630 --> 00:39:50,630
But we'll find out that
partial f, partial x dot,
805
00:39:50,630 --> 00:39:53,630
the Jacobian of f
with respect to x dot,
806
00:39:53,630 --> 00:39:56,320
that's just the momentum.
807
00:39:56,320 --> 00:39:57,605
That's not a matrix.
808
00:39:57,605 --> 00:40:00,120
That's not invertible.
809
00:40:00,120 --> 00:40:02,390
So we can't convert
this equation
810
00:40:02,390 --> 00:40:04,600
to an equivalent system
of ordinary differential
811
00:40:04,600 --> 00:40:05,100
equations.
812
00:40:05,100 --> 00:40:06,558
So it's just
something pathological
813
00:40:06,558 --> 00:40:08,930
for a one-dimensional system.
814
00:40:08,930 --> 00:40:13,010
Instead, what we find is
zero is equal to the velocity
815
00:40:13,010 --> 00:40:17,740
dotted with the mass times the
acceleration plus the gradient
816
00:40:17,740 --> 00:40:19,040
in the potential.
817
00:40:19,040 --> 00:40:22,880
So minus the forces
acting on the molecules.
818
00:40:22,880 --> 00:40:24,590
We know that if we
were to integrate
819
00:40:24,590 --> 00:40:27,590
the equations of motion
for the molecules exactly,
820
00:40:27,590 --> 00:40:30,710
the momentum balances for
those molecules exactly,
821
00:40:30,710 --> 00:40:32,900
then this term here would
always be equal to zero.
822
00:40:32,900 --> 00:40:37,090
And we'd satisfy
conservation of energy.
823
00:40:37,090 --> 00:40:40,130
Of course, all solutions of
ordinary differential equations
824
00:40:40,130 --> 00:40:42,170
require approximations.
825
00:40:42,170 --> 00:40:44,312
So as you move the
molecules around,
826
00:40:44,312 --> 00:40:45,770
their velocities
and positions will
827
00:40:45,770 --> 00:40:48,200
tend to drift away
from the solution
828
00:40:48,200 --> 00:40:51,950
where this is exactly in
balance and equal to zero.
829
00:40:51,950 --> 00:40:54,950
If you desire to do that in the
molecular dynamic simulation
830
00:40:54,950 --> 00:40:56,630
then you better
have a method that
831
00:40:56,630 --> 00:40:59,960
somehow respects this
geometric property instead.
832
00:40:59,960 --> 00:41:02,950
That the velocities
are orthogonal to m
833
00:41:02,950 --> 00:41:05,000
x dot plus grad v.
834
00:41:05,000 --> 00:41:08,180
And there is a set of
ODE-IVP methods that do that.
835
00:41:08,180 --> 00:41:10,841
They're called
symplectic integrators.
836
00:41:10,841 --> 00:41:12,590
And they integrate the
equations of motion
837
00:41:12,590 --> 00:41:16,220
while exerting some control over
the error in the total energy.
838
00:41:16,220 --> 00:41:19,040
So if you were just
the take a ODE45
839
00:41:19,040 --> 00:41:21,080
and try to solve this
system of equations,
840
00:41:21,080 --> 00:41:24,530
over long times you would find
that the energy drifts away
841
00:41:24,530 --> 00:41:26,340
from the place where it started.
842
00:41:26,340 --> 00:41:28,670
Even though your positions
are being integrated
843
00:41:28,670 --> 00:41:31,297
with reasonable
accuracy, over time
844
00:41:31,297 --> 00:41:33,380
the energy will drift away
from where you started.
845
00:41:33,380 --> 00:41:36,590
The system will heat up
effectively or cool down.
846
00:41:36,590 --> 00:41:39,200
That's undesirable if you're
trying to do modeling.
847
00:41:39,200 --> 00:41:41,630
So instead, one uses the
symplectic integrators,
848
00:41:41,630 --> 00:41:45,810
which are designed to control
the error and the energy
849
00:41:45,810 --> 00:41:48,336
in addition to integrate
the equations of motion.
850
00:41:48,336 --> 00:41:50,210
So you can look those
up and read about them.
851
00:41:50,210 --> 00:41:51,694
I think they're
quite interesting.
852
00:41:51,694 --> 00:41:53,360
If you're going to
do molecular modeling
853
00:41:53,360 --> 00:41:55,430
you'd like to understand
better how they work.
854
00:41:55,430 --> 00:41:56,929
But actually,
conservation of energy
855
00:41:56,929 --> 00:41:59,870
here doesn't give
us a well-posed DAE.
856
00:41:59,870 --> 00:42:02,060
We actually need the
momentum equations
857
00:42:02,060 --> 00:42:07,680
to formulate or to determine the
trajectories of the molecules.
858
00:42:07,680 --> 00:42:12,120
Let's talk about their
numerical simulation.
859
00:42:12,120 --> 00:42:14,540
So let's think about
these semi-explicit DAEs,
860
00:42:14,540 --> 00:42:16,640
the ones that can be
split between differential
861
00:42:16,640 --> 00:42:17,790
and algebraic states.
862
00:42:17,790 --> 00:42:21,740
It's kind of useful to look
at this problem in particular.
863
00:42:21,740 --> 00:42:25,700
And you might say, well, let's
do the simplest possible thing.
864
00:42:25,700 --> 00:42:29,000
Let's do the forward
Euler approximation.
865
00:42:29,000 --> 00:42:31,880
So let's take our derivative
here and represent it
866
00:42:31,880 --> 00:42:36,500
as the state evaluated at x plus
delta t minus the state at t
867
00:42:36,500 --> 00:42:38,080
divided by delta t.
868
00:42:38,080 --> 00:42:40,070
Let's make the right-hand
side of our equation
869
00:42:40,070 --> 00:42:43,340
be evaluated at time t.
870
00:42:43,340 --> 00:42:46,550
And let's try to step in
time, do time marching forward
871
00:42:46,550 --> 00:42:49,550
to determine the
trajectory of x and t.
872
00:42:49,550 --> 00:42:51,440
How do the states change?
873
00:42:51,440 --> 00:42:53,210
We see the first
equation is an easy one
874
00:42:53,210 --> 00:42:58,870
to solve because presumably
I know x of t and y of t.
875
00:42:58,870 --> 00:43:01,550
Well, do I know y of t?
y of t has to satisfy
876
00:43:01,550 --> 00:43:02,720
this second equation.
877
00:43:02,720 --> 00:43:06,200
So first, if I know x of t I
need to solve this nonlinear
878
00:43:06,200 --> 00:43:08,190
equation determine y of t.
879
00:43:08,190 --> 00:43:09,830
Then I know x of
t and y of t and I
880
00:43:09,830 --> 00:43:14,390
can step forward in time to
determine x of t plus delta t.
881
00:43:14,390 --> 00:43:17,720
So every iteration here with the
simplest possible approximation
882
00:43:17,720 --> 00:43:20,030
for the differential
requires solution of a system
883
00:43:20,030 --> 00:43:21,330
of nonlinear equations.
884
00:43:24,230 --> 00:43:27,540
We already saw before that a
forward Euler approximation
885
00:43:27,540 --> 00:43:32,640
is not a great one to apply
to solutions of IDE-IVPs
886
00:43:32,640 --> 00:43:35,660
because it has very limited
stability properties.
887
00:43:35,660 --> 00:43:37,630
You need small
time steps in order
888
00:43:37,630 --> 00:43:39,560
to stably integrate in time.
889
00:43:42,334 --> 00:43:43,750
So even though
this seemed like it
890
00:43:43,750 --> 00:43:46,480
was easy to get a solution
for the differential equation,
891
00:43:46,480 --> 00:43:51,620
we've always got to satisfy
this nonlinear equation here.
892
00:43:51,620 --> 00:43:54,970
So inherently, simulation
of DAEs are implicit.
893
00:43:54,970 --> 00:43:57,700
They require solution
of nonlinear equations
894
00:43:57,700 --> 00:43:59,980
in order to advance in time.
895
00:43:59,980 --> 00:44:03,590
There's actually no point in
using this weakly stable method
896
00:44:03,590 --> 00:44:04,990
for any very small time steps.
897
00:44:04,990 --> 00:44:08,200
It makes more sense to choose
a naturally implicit method
898
00:44:08,200 --> 00:44:11,830
for the differential equation,
where I can take big time steps
899
00:44:11,830 --> 00:44:16,376
and still have reasonable
stability of the integration.
900
00:44:16,376 --> 00:44:17,250
Does that make sense?
901
00:44:21,040 --> 00:44:25,180
Consider now the
fully implicit DAE.
902
00:44:25,180 --> 00:44:28,480
So f of x, x dot,
and t is equal to 0
903
00:44:28,480 --> 00:44:30,565
and have some set of
initial conditions. x of 0
904
00:44:30,565 --> 00:44:32,590
is equal to x 0.
905
00:44:32,590 --> 00:44:34,960
Here there is no way in
general to avoid solving
906
00:44:34,960 --> 00:44:37,330
systems of nonlinear equations.
907
00:44:37,330 --> 00:44:39,850
And so one substitutes
often backwards
908
00:44:39,850 --> 00:44:42,670
difference
approximations for x dot.
909
00:44:42,670 --> 00:44:45,355
And then you solve for
the next point in time.
910
00:44:45,355 --> 00:44:47,230
So here's an example
with the backward Euler,
911
00:44:47,230 --> 00:44:48,730
implicit Euler approximations.
912
00:44:48,730 --> 00:44:51,070
So I replaced the
time derivative
913
00:44:51,070 --> 00:44:55,000
with the difference between
the state at point tk
914
00:44:55,000 --> 00:44:57,160
and the state of
point tk minus 1
915
00:44:57,160 --> 00:45:00,880
divided by the difference
between tk and tk minus 1.
916
00:45:00,880 --> 00:45:03,670
The other state variables
are evaluated at tk.
917
00:45:03,670 --> 00:45:05,130
So x at tk.
918
00:45:05,130 --> 00:45:06,430
And I evaluate time at tk.
919
00:45:06,430 --> 00:45:08,650
And I try to satisfy
this equation
920
00:45:08,650 --> 00:45:10,460
and use it to determine x of tk.
921
00:45:10,460 --> 00:45:13,690
So I solve this system
of nonlinear equations
922
00:45:13,690 --> 00:45:15,160
for x at tk.
923
00:45:15,160 --> 00:45:18,580
And then I go to the next
point in time, tk plus 1.
924
00:45:18,580 --> 00:45:22,440
And I solve the same equation,
but replacing tk with tk plus 1
925
00:45:22,440 --> 00:45:25,060
and tk minus 1 with tk.
926
00:45:25,060 --> 00:45:28,060
So at each iteration I solve a
system of nonlinear equations.
927
00:45:28,060 --> 00:45:30,400
No more painful than
doing the forward Euler
928
00:45:30,400 --> 00:45:32,860
approximation I
showed you before.
929
00:45:32,860 --> 00:45:35,590
So that's good.
930
00:45:35,590 --> 00:45:38,865
But still a lot of work.
931
00:45:38,865 --> 00:45:40,240
We might wonder
how do we control
932
00:45:40,240 --> 00:45:42,872
the error in these
sorts of approximations.
933
00:45:42,872 --> 00:45:44,830
And that's what we're
going to talk about next.
934
00:45:48,394 --> 00:45:50,310
So here's the equivalent
time marching formula
935
00:45:50,310 --> 00:45:52,770
applied to a semi-explicit DAE.
936
00:45:55,910 --> 00:46:02,790
So I've got to
satisfy the equation
937
00:46:02,790 --> 00:46:05,850
0 is equal to the derivative
of x with respect to time,
938
00:46:05,850 --> 00:46:09,060
which I approximate with the
backward Euler approximation,
939
00:46:09,060 --> 00:46:12,300
minus f evaluated at tk.
940
00:46:12,300 --> 00:46:19,410
And the equation g of x at
tk, y at tk for this x of tk.
941
00:46:19,410 --> 00:46:20,670
And actually y of tk as well.
942
00:46:20,670 --> 00:46:22,800
I'm sorry, I left that
off. y of tk as well.
943
00:46:22,800 --> 00:46:24,890
Got to solve for
the differential
944
00:46:24,890 --> 00:46:28,340
and the algebraic
states simultaneously.
945
00:46:28,340 --> 00:46:30,030
And I use a backwards
difference formula
946
00:46:30,030 --> 00:46:31,877
to try to do this
in a stable way.
947
00:46:36,850 --> 00:46:43,180
One way to think about in
particular these sorts of DAEs
948
00:46:43,180 --> 00:46:46,990
is as exceptionally stiff
ordinary differential
949
00:46:46,990 --> 00:46:47,490
equations.
950
00:46:47,490 --> 00:46:50,320
Did you guys discuss stiffness?
951
00:46:50,320 --> 00:46:55,750
So imagine if I replace
the 0 with, say, dy dt.
952
00:46:55,750 --> 00:46:58,090
And I introduced a constant
on the right-hand side,
953
00:46:58,090 --> 00:47:05,485
a multiplicative constant,
say, 1 over lambda here.
954
00:47:05,485 --> 00:47:06,610
Let's say lambda, actually.
955
00:47:06,610 --> 00:47:08,610
Say I introduce a
multiplicative constant lambda
956
00:47:08,610 --> 00:47:10,940
on the right-hand side here.
957
00:47:10,940 --> 00:47:15,010
So as lambda becomes
very large, the system
958
00:47:15,010 --> 00:47:17,620
becomes increasingly stiff.
959
00:47:17,620 --> 00:47:19,600
The dynamics of
this equation take
960
00:47:19,600 --> 00:47:22,790
place over very, very
short time scales.
961
00:47:22,790 --> 00:47:24,820
So short in fact,
that eventually,
962
00:47:24,820 --> 00:47:26,500
if I let lambda
diverge to infinity,
963
00:47:26,500 --> 00:47:29,500
I've just got to satisfy
this algebraic equation
964
00:47:29,500 --> 00:47:31,400
at each point in time.
965
00:47:31,400 --> 00:47:34,240
So these are exceptionally
stiff equations.
966
00:47:34,240 --> 00:47:36,270
And the methods for
solving these equations
967
00:47:36,270 --> 00:47:39,880
share a lot in common with
stiff ODE-IVP solvers.
968
00:47:45,069 --> 00:47:46,610
So I've got a couple
of minutes left.
969
00:47:46,610 --> 00:47:48,880
I'm going to work
through some examples.
970
00:47:48,880 --> 00:47:51,540
So consider the stirred
tank example one.
971
00:47:51,540 --> 00:47:53,250
dc 2 dt.
972
00:47:53,250 --> 00:47:57,210
It's related to q over
v, flow rate over volume.
973
00:47:57,210 --> 00:47:58,770
The difference in
the concentrations
974
00:47:58,770 --> 00:48:00,570
coming in and out of the tank.
975
00:48:00,570 --> 00:48:03,060
And c1 is equal to this
input signal gamma.
976
00:48:03,060 --> 00:48:06,030
Maybe that's a measurement
and I'm trying to predict c2.
977
00:48:06,030 --> 00:48:08,850
So can you apply the backward
Euler method to these equations
978
00:48:08,850 --> 00:48:12,810
to determine c1 and c2?
979
00:48:12,810 --> 00:48:15,670
What would one step of the
backward Euler method look
980
00:48:15,670 --> 00:48:24,500
like, say, going from time 0
to time delta t, or time t to t
981
00:48:24,500 --> 00:48:27,260
plus delta t, or tk to
tk plus 1, however you
982
00:48:27,260 --> 00:48:28,748
want to write it?
983
00:49:22,314 --> 00:49:23,730
So do you get
something like this?
984
00:49:26,480 --> 00:49:28,530
So this equation is easy.
985
00:49:28,530 --> 00:49:32,360
c1 to tk is gamma of tk.
986
00:49:32,360 --> 00:49:34,390
I got a substitute up here.
987
00:49:34,390 --> 00:49:38,570
c2 of tk minus c2 of tk
minus 1 over delta t.
988
00:49:38,570 --> 00:49:41,210
T Just tk minus tk minus 1.
989
00:49:41,210 --> 00:49:42,800
And then solve for c2 of dk.
990
00:49:42,800 --> 00:49:44,390
And you get a formula like this.
991
00:49:49,010 --> 00:49:52,010
This formula was based
on an approximation
992
00:49:52,010 --> 00:49:56,510
for the derivative,
which had a leading order
993
00:49:56,510 --> 00:50:03,140
error that was proportional
to tk minus tk minus 1.
994
00:50:03,140 --> 00:50:07,070
When I go through and I solve
for c2 of tk, the leading order
995
00:50:07,070 --> 00:50:10,600
error in c2, the
local truncation error
996
00:50:10,600 --> 00:50:16,250
is going to be order tk
minus tk minus 1 squared.
997
00:50:16,250 --> 00:50:20,840
So as I shrink the spacing
between each of my two time
998
00:50:20,840 --> 00:50:24,050
points, the local
accuracy, the error
999
00:50:24,050 --> 00:50:27,830
induced in one advance of
this time stepping algorithm
1000
00:50:27,830 --> 00:50:29,790
is going to scale like
the step size squared.
1001
00:50:29,790 --> 00:50:31,550
So I have the error--
1002
00:50:31,550 --> 00:50:32,369
I have the spacing.
1003
00:50:32,369 --> 00:50:34,160
The error will go down
by a factor of four.
1004
00:50:34,160 --> 00:50:37,667
I reduce the spacing by
an order of magnitude.
1005
00:50:37,667 --> 00:50:39,750
The error will go down by
two orders of magnitude.
1006
00:50:39,750 --> 00:50:43,605
This is a really desirable sort
of error control and a method.
1007
00:50:43,605 --> 00:50:44,480
And we can do better.
1008
00:50:44,480 --> 00:50:48,460
You've seen us do better
with other methods.
1009
00:50:48,460 --> 00:50:51,180
Let's look at this equation now.
1010
00:50:51,180 --> 00:50:52,590
This is stirred
tank example two.
1011
00:50:52,590 --> 00:50:53,797
I just replace c1 with c2.
1012
00:50:53,797 --> 00:50:55,380
You don't have to
write this one down.
1013
00:50:55,380 --> 00:50:56,484
I'll just show it to you.
1014
00:50:56,484 --> 00:50:58,650
The notes are online so
you're not missing anything.
1015
00:50:58,650 --> 00:51:00,900
All these things will pop
up in the notes online.
1016
00:51:00,900 --> 00:51:05,220
So c2 at tk is equal
to gamma of tk.
1017
00:51:05,220 --> 00:51:06,930
I got to come up to
my first equation now
1018
00:51:06,930 --> 00:51:08,700
and I have to substitute
an approximation
1019
00:51:08,700 --> 00:51:09,810
for the derivative.
1020
00:51:09,810 --> 00:51:13,390
It's this, the backwards
difference formula.
1021
00:51:13,390 --> 00:51:17,610
So if find c1 of tk is
c2 of tk plus v over q.
1022
00:51:17,610 --> 00:51:21,660
c2 of tk minus c2 of tk minus
1 divided by the time step.
1023
00:51:21,660 --> 00:51:23,640
Remember, this approximation
for the derivative
1024
00:51:23,640 --> 00:51:27,292
has a leading order error that
goes like tk minus tk minus 1.
1025
00:51:27,292 --> 00:51:32,235
So my approximation for c1,
the local approximation for c1,
1026
00:51:32,235 --> 00:51:34,040
is not second order accurate.
1027
00:51:34,040 --> 00:51:38,010
It's only first order accurate.
1028
00:51:38,010 --> 00:51:42,450
I applied the same approximation
method to both equations.
1029
00:51:42,450 --> 00:51:45,050
One equation I got something
that was second order accurate.
1030
00:51:45,050 --> 00:51:46,425
The next equation
I got something
1031
00:51:46,425 --> 00:51:47,840
that was first order accurate.
1032
00:51:47,840 --> 00:51:49,520
They look almost identical.
1033
00:51:49,520 --> 00:51:51,530
It might be hard to
see from the start
1034
00:51:51,530 --> 00:51:53,420
why it should work out that way.
1035
00:51:53,420 --> 00:51:56,190
That's how it works out.
1036
00:51:56,190 --> 00:51:57,410
Let's do one more example.
1037
00:52:00,640 --> 00:52:03,680
Here's a system of DAEs.
1038
00:52:03,680 --> 00:52:05,760
c2 dot is equal to c1.
1039
00:52:05,760 --> 00:52:07,900
c3 dot is equal to c2.
1040
00:52:07,900 --> 00:52:10,240
And 0 is equal to
c3 minus gamma.
1041
00:52:10,240 --> 00:52:12,250
So gamma, again, is
some measurement.
1042
00:52:12,250 --> 00:52:14,660
We solve this equation
to determine c3.
1043
00:52:14,660 --> 00:52:17,640
We substitute c3 up
here to determine c2.
1044
00:52:17,640 --> 00:52:20,600
We substitute c2 up
here to determine c1.
1045
00:52:20,600 --> 00:52:23,965
But we want to find the solution
with the backward Euler method.
1046
00:52:27,100 --> 00:52:32,680
And the solution there is going
to be c3 of tk is gamma of tk.
1047
00:52:32,680 --> 00:52:35,410
c2 is related to the
derivative of c3.
1048
00:52:35,410 --> 00:52:37,180
So I use my backwards
difference formula
1049
00:52:37,180 --> 00:52:39,060
for the derivative of c3.
1050
00:52:39,060 --> 00:52:41,620
c1 is related to the
derivative of c2.
1051
00:52:41,620 --> 00:52:44,710
So I use my backward
difference formula for c2.
1052
00:52:44,710 --> 00:52:46,930
But how accurate are each
of these approximations?
1053
00:52:49,840 --> 00:52:52,780
So here I have to determine
c2 from an approximation
1054
00:52:52,780 --> 00:52:54,880
for the derivative of c3.
1055
00:52:54,880 --> 00:52:56,920
We know this backwards
difference approximation
1056
00:52:56,920 --> 00:52:58,510
has a leading order
error proportional
1057
00:52:58,510 --> 00:53:02,350
to delta t, the time step.
1058
00:53:02,350 --> 00:53:07,750
Now I got to approximate c1 with
the finite difference backwards
1059
00:53:07,750 --> 00:53:09,922
difference approximation in c2.
1060
00:53:09,922 --> 00:53:12,380
We know this should incur an
error proportional to the time
1061
00:53:12,380 --> 00:53:17,350
step delta t, except I
don't know c2 at tk exactly.
1062
00:53:17,350 --> 00:53:23,530
Actually, c2 tk, the
exact value of c2 at tk
1063
00:53:23,530 --> 00:53:27,830
is this plus something
proportional to the time step.
1064
00:53:27,830 --> 00:53:30,560
Which means the
exact value of c1
1065
00:53:30,560 --> 00:53:34,430
is this plus something
that's order one because I
1066
00:53:34,430 --> 00:53:36,299
carry over the error
from my solution
1067
00:53:36,299 --> 00:53:37,340
to the previous equation.
1068
00:53:37,340 --> 00:53:40,477
And I divide it
by the time step.
1069
00:53:40,477 --> 00:53:42,560
So I went from having
something-- well, let's see.
1070
00:53:42,560 --> 00:53:44,340
The first equation
had order delta t
1071
00:53:44,340 --> 00:53:47,040
squared local truncation error.
1072
00:53:47,040 --> 00:53:50,010
As I shrink the
spacing, my solution
1073
00:53:50,010 --> 00:53:52,440
gets more and more
locally accurate.
1074
00:53:52,440 --> 00:53:57,300
And it seems to do it in
a fairly robust fashion.
1075
00:53:57,300 --> 00:53:59,910
The next system of
equations we solved
1076
00:53:59,910 --> 00:54:02,170
had a local truncation error
proportional to delta t.
1077
00:54:02,170 --> 00:54:03,545
Again, I can shrink
delta t and I
1078
00:54:03,545 --> 00:54:05,800
can make my solution
more and more accurate.
1079
00:54:05,800 --> 00:54:08,220
That seems like what you want
in a numerical algorithm.
1080
00:54:08,220 --> 00:54:11,580
This one, there's no hope.
1081
00:54:11,580 --> 00:54:15,980
Change delta t, who knows
what you're going to get?
1082
00:54:15,980 --> 00:54:18,855
And this looks fairly innocuous.
1083
00:54:18,855 --> 00:54:20,990
You just look at it
and you say, well, OK.
1084
00:54:20,990 --> 00:54:22,980
Let's solve it and
see what happens.
1085
00:54:22,980 --> 00:54:26,040
Well actually, you can't
solve this problem accurately
1086
00:54:26,040 --> 00:54:28,322
with the backward Euler method.
1087
00:54:28,322 --> 00:54:30,030
There's something
fundamentally different
1088
00:54:30,030 --> 00:54:32,359
about this problem
from the previous one.
1089
00:54:32,359 --> 00:54:34,650
There's something fundamentally
different about stirred
1090
00:54:34,650 --> 00:54:37,500
tank example two from
stirred tank example one.
1091
00:54:37,500 --> 00:54:40,170
Right you might have
already perceived that.
1092
00:54:40,170 --> 00:54:42,160
You can see what's
going on here.
1093
00:54:42,160 --> 00:54:44,520
c3 is equal to gamma.
1094
00:54:44,520 --> 00:54:46,025
That's the exact solution.
1095
00:54:46,025 --> 00:54:48,390
c2 is equal to gamma dot.
1096
00:54:48,390 --> 00:54:50,700
And c1 is equal
to two derivatives
1097
00:54:50,700 --> 00:54:53,860
of gamma, the exact
solution to this equation.
1098
00:54:53,860 --> 00:54:59,310
So c1 is incredibly sensitive
to changes in gamma.
1099
00:54:59,310 --> 00:55:00,620
And that's peculiar.
1100
00:55:00,620 --> 00:55:03,180
And it has important
numerical consequences.
1101
00:55:03,180 --> 00:55:06,380
So I'm going to show you, our
next lecture on Wednesday,
1102
00:55:06,380 --> 00:55:09,620
how to formalize
an understanding
1103
00:55:09,620 --> 00:55:11,450
of the differences
between these problems.
1104
00:55:11,450 --> 00:55:14,990
And we'll talk about some
more subtleties associated
1105
00:55:14,990 --> 00:55:17,150
with differential
algebraic equations.
1106
00:55:17,150 --> 00:55:18,650
You guys can pick
up your quizzes.
1107
00:55:18,650 --> 00:55:23,000
We should do it
outside rather than
1108
00:55:23,000 --> 00:55:25,596
inside because the next
class has to start.