1
00:00:14,440 --> 00:00:21,580
We'll begin discussion on a new topic today,
this is on time integration of equation of
2
00:00:21,580 --> 00:00:28,090
motion,
so at the end of our previous discussions
3
00:00:28,090 --> 00:00:35,399
we have seen that the governing equation of
motion typically is formulated in this form
4
00:00:35,399 --> 00:00:41,830
MU double dot + U dot + KU, for sake of generality
I am now introducing a nonlinear term we have
5
00:00:41,830 --> 00:00:47,329
not discussed how this originates, but we
can believe that this type of terms would
6
00:00:47,329 --> 00:00:53,460
arise if we include nonlinear behavior either
in strain displacement relations or in stress
7
00:00:53,460 --> 00:00:59,859
strain relations, so this is the equation
of motion, and these are the specified initial
8
00:00:59,859 --> 00:01:05,500
condition, so this equation constitutes a
set of coupled second-order nonlinear ordinary
9
00:01:05,500 --> 00:01:12,909
differential equations, we say that these
equations are semi discretized equations because
10
00:01:12,909 --> 00:01:18,610
in arriving at this equation we have discretized
space but the time is still a continuous parameter,
11
00:01:18,610 --> 00:01:27,540
so they also constitute a set of initial value
problem, so these are a set of ordinary differential
12
00:01:27,540 --> 00:01:33,290
equations, semi discretized ordinary differential
equations which constitute an initial value
13
00:01:33,290 --> 00:01:41,310
problem and these equations are coupled.
Now we consider the problem of discretizing
14
00:01:41,310 --> 00:01:48,670
in time so what we aim to do is to we consider
solution of this equation at a set of discrete
15
00:01:48,670 --> 00:01:57,000
time instants ordered as T naught, T1, T2,
TN with delta TN being the step size, TN +
16
00:01:57,000 --> 00:02:02,280
1 - TN. The basic idea here is to replace
the derivatives appearing in the equation
17
00:02:02,280 --> 00:02:08,350
of motion by finite difference approximations
and then solve the resulting algebraic equation
18
00:02:08,350 --> 00:02:12,540
that is a main idea. Now to be able to do
that what we will do is we
19
00:02:12,540 --> 00:02:17,340
will reorganize the equation of motion in
a slightly different form, so this is the
20
00:02:17,340 --> 00:02:25,000
equation of motion that we have MU double
dot + CU dot + KU + nonlinear terms is F(t)
21
00:02:25,000 --> 00:02:31,530
and specified initial conditions. Now I pre
multiply by M inverse so this equation now
22
00:02:31,530 --> 00:02:40,310
becomes, takes this form. I will introduce
now a new set of variables X1 and X2, X1 is
23
00:02:40,310 --> 00:02:47,510
U(t), X2 is U dot(t), I can define another
2N cross 1 vector in which I assemble X1 and
24
00:02:47,510 --> 00:02:53,090
X2 as shown here, we call this as system state
vector, so this is a state vector consisting
25
00:02:53,090 --> 00:02:56,870
of displacements and velocities at the end
degrees of freedom.
26
00:02:56,870 --> 00:03:05,880
Now what is X1 dot? X1 dot is U dot, U dot
is X2, so this is X1 dot is X2, X2 dot is
27
00:03:05,880 --> 00:03:13,319
U double dot so that I obtained from the governing
equation and we write in this form. Now I
28
00:03:13,319 --> 00:03:17,480
will assemble this equation in the matrix
form, I will have X1 dot, X2 dot is equal
29
00:03:17,480 --> 00:03:24,860
to this matrix into X1, X2 plus the non-linear
term, plus the excitation term. Now I can
30
00:03:24,860 --> 00:03:33,110
rewrite this set of equations in a general
form as X dot is some function A, X(t,t) where
31
00:03:33,110 --> 00:03:44,000
A, this vector A is 2N cross 1, we say that
this equation is in state space form, by state
32
00:03:44,000 --> 00:03:49,450
I mean a displacement and velocity vector
together they had constitute the state, and
33
00:03:49,450 --> 00:03:54,349
this equation is said to be written in the
configuration space where we have only displacements
34
00:03:54,349 --> 00:04:01,239
and we get second order differential equations,
in state space we have a set of 2N first order
35
00:04:01,239 --> 00:04:05,950
coupled ordinary differential equations which
are initial value problems and we say that
36
00:04:05,950 --> 00:04:12,200
the equation is in the state space form.
So what we have done is the set of N coupled
37
00:04:12,200 --> 00:04:16,850
second order ordinary differential equations
given by this in the configuration space have
38
00:04:16,850 --> 00:04:22,470
been thus recast as the following set of 2N
first order coupled ordinary differential
39
00:04:22,470 --> 00:04:31,120
equation. Now what I do is I introduce a notation
X(t,x naught, t naught) I denote solution
40
00:04:31,120 --> 00:04:37,160
of this differential equation with X(t naught)
= X naught, so there is T naught, X naught,
41
00:04:37,160 --> 00:04:43,420
and T of course is the independent variable.
In discrete time approximation, time will
42
00:04:43,420 --> 00:04:50,840
be discretized as increasing sequence T naught,
T1, T2, TN with steps as delta TN, as I mentioned
43
00:04:50,840 --> 00:05:01,560
before and I denote the value of the system
state at T = TN by YN, so YN is X(TN) where
44
00:05:01,560 --> 00:05:08,350
initial conditions at T = T naught is X naught.
Now the main question is how to approximate
45
00:05:08,350 --> 00:05:16,690
X dot in terms of YN, okay, that is the problem.
Now we could use finite difference approximation
46
00:05:16,690 --> 00:05:22,650
for example X dot(t) can be written using
forward difference scheme, X dot(t) is approximately
47
00:05:22,650 --> 00:05:29,320
equal to X(t) + delta T - X(t)/delta T, so
consequently I can write this as YN dot = YN
48
00:05:29,320 --> 00:05:36,120
+ 1 - YN / delta TN, so from this I can derive
now the equation for YN + 1, which will be
49
00:05:36,120 --> 00:05:47,030
YN + delta TN into YN dot, YN is nothing but
A(YN,TN), where N runs from 0, 1, 2, 3, with
50
00:05:47,030 --> 00:05:57,910
Y naught = X naught at T = T naught. Now delta
TN are the algorithmic parameter delta T1,
51
00:05:57,910 --> 00:06:02,380
delta T2, delta T3 are the algorithmic parameters,
if all of them are equal then delta T is the
52
00:06:02,380 --> 00:06:07,460
only algorithmic parameter. Now the accuracy
of this solution depends upon this choice
53
00:06:07,460 --> 00:06:13,080
of these algorithmic parameters delta TN,
if delta TN are all equal the accuracy depends
54
00:06:13,080 --> 00:06:21,840
basically on the choice of delta T.
Now we talk about errors, we call something
55
00:06:21,840 --> 00:06:32,800
on a local discretization error, this is X(TN+1,
XN, TN) - YN+1 that means this is, I have
56
00:06:32,800 --> 00:06:41,620
moved from T = TN to T = TN+1 to get X(TN+1)
and that is approximated as YN + 1, so this
57
00:06:41,620 --> 00:06:45,680
error is called local discretization scheme,
that means we are basically moving from TN
58
00:06:45,680 --> 00:06:50,800
to TN+1. On the other hand if you move from
T naught to TN+1 that means all the preceding
59
00:06:50,800 --> 00:06:56,980
steps up to N+1 step if you cover then we
call this error as the global discretization
60
00:06:56,980 --> 00:07:03,790
error, there is other error called round off
error which is due to finite precision calculations,
61
00:07:03,790 --> 00:07:07,440
so you do either double precision or single
precision calculation, typically we do double
62
00:07:07,440 --> 00:07:12,810
precision calculations, so on a computer the
digits get terminated so it is inevitable
63
00:07:12,810 --> 00:07:19,220
so there will be errors, that we call as RN+1.
Of course there are other errors which we
64
00:07:19,220 --> 00:07:24,300
call as dynamical system modeling errors,
for example does algorithm correctly capture
65
00:07:24,300 --> 00:07:29,800
dynamical properties of system like natural
frequencies, frequency response functions,
66
00:07:29,800 --> 00:07:34,680
free vibration amplitudes, is the energy conserved
suppose I am dealing with an undamped free
67
00:07:34,680 --> 00:07:40,580
vibration problem, is the energy conserved?
So if there is any compromise on any of these
68
00:07:40,580 --> 00:07:49,110
issues, these can be grouped as dynamical
system modeling errors.
69
00:07:49,110 --> 00:07:53,850
Now the basic question that we are looking
for is how to formulate the integration algorithms
70
00:07:53,850 --> 00:07:58,870
to achieve acceptable accuracy, there's so
many errors of different kinds so associated
71
00:07:58,870 --> 00:08:04,460
with that there is a definition of an accuracy
so how to understand the nature of different
72
00:08:04,460 --> 00:08:11,321
errors, okay, and does this error as N tends
to infinity norm of EN, does it go to 0, or
73
00:08:11,321 --> 00:08:17,430
does it remain finite, or does it become unbounded,
so at every time step certain error is made
74
00:08:17,430 --> 00:08:23,360
either because of truncating of, we are approximating
the derivative by a finite difference approximation
75
00:08:23,360 --> 00:08:29,210
and there is a round off error, so consequently
at every time step there will be an error,
76
00:08:29,210 --> 00:08:34,200
so how does this error committed at one step,
how does it propagate to the future steps,
77
00:08:34,200 --> 00:08:39,389
so does it go to 0 or it remains bounded or
does it become unbounded. So now how to select
78
00:08:39,389 --> 00:08:44,039
algorithmic parameters to achieve acceptable
performance of the integrators, so what are
79
00:08:44,039 --> 00:08:47,630
the expectations from a good integrator, so
these are the type of questions we want to
80
00:08:47,630 --> 00:08:53,420
now address.
Now let's spend some time on dynamical system
81
00:08:53,420 --> 00:08:59,090
modeling errors so that we understand what
exactly is meant, so now let us consider the
82
00:08:59,090 --> 00:09:04,260
scalar equation of motion U double dot + omega
square U = 0, so it is a single degree freedom
83
00:09:04,260 --> 00:09:07,940
system, natural frequency is omega and it
starts with initial condition U naught and
84
00:09:07,940 --> 00:09:14,200
U naught dot, we can solve this problem exactly
and we know that this is the solution, we
85
00:09:14,200 --> 00:09:18,580
have an amplitude and FS that they are given
by this, they are functions of initial conditions
86
00:09:18,580 --> 00:09:23,890
and system natural frequency. Now the question
is suppose if I solve this equation using
87
00:09:23,890 --> 00:09:30,270
the discrete approximation, do we get amplitude
of response to be constant, for example do
88
00:09:30,270 --> 00:09:35,220
we get the solution in this form where R is
a constant, otherwise we get numerical dissipation
89
00:09:35,220 --> 00:09:41,180
or growth, if amplitude is not constant either
it will decay or it will grow, so then that
90
00:09:41,180 --> 00:09:46,470
would mean that the algorithm has introduced
some artificial dissipation mechanism into
91
00:09:46,470 --> 00:09:51,260
the system which is not there in the original
equation of motion, then do we get the impulse
92
00:09:51,260 --> 00:09:57,920
response function correctly, is the distribution
distortions acceptable, okay, what is acceptable
93
00:09:57,920 --> 00:10:04,180
distortion? Now is the vibration energy conserved,
some of these questions are important, so
94
00:10:04,180 --> 00:10:10,890
answering these questions would enable us
to discuss the dynamical system modeling errors.
95
00:10:10,890 --> 00:10:17,020
We could also consider for example a damped
harmonically driven oscillator and we know
96
00:10:17,020 --> 00:10:21,740
that in steady state this type of systems
have certain well-known features, there is
97
00:10:21,740 --> 00:10:26,101
a resonance and where the resonance occurs,
what is resonance amplitude, now what is the
98
00:10:26,101 --> 00:10:30,380
shape of the frequency response function so
on and so forth, there are well known qualitative
99
00:10:30,380 --> 00:10:34,580
features associated with the response of the
system. Now the question we can ask is does
100
00:10:34,580 --> 00:10:41,510
the discretized version of the solution process
these well-known features or not, okay, so
101
00:10:41,510 --> 00:10:45,430
that is one. Similarly if there is a nonlinear
term we are not discussed about nonlinear
102
00:10:45,430 --> 00:10:52,830
term, but we can think of adding say for example
alpha U cube here, now such systems display
103
00:10:52,830 --> 00:10:59,240
a pattern of bifurcations we will come to
that later but at this stage it is suffice
104
00:10:59,240 --> 00:11:05,640
to observe that nonlinear systems have certain
qualitative features in the associated with
105
00:11:05,640 --> 00:11:11,930
their behavior, and we can ask the question
does the discretized version of the equation
106
00:11:11,930 --> 00:11:17,880
display the same type of qualitative behavior
as the continuous version, okay, this type
107
00:11:17,880 --> 00:11:24,519
of issues can be grouped under the heading
of dynamical system modeling errors.
108
00:11:24,519 --> 00:11:29,500
So what we are going to do is, there will
be a few mathematical preliminaries that are
109
00:11:29,500 --> 00:11:34,830
needed to understand the questions and the
answers that we will be discussing related
110
00:11:34,830 --> 00:11:39,680
to the properties of these integration schemes,
and there are some set of terminologies we
111
00:11:39,680 --> 00:11:46,120
talked about implicit schemes, explicit scheme,
single step method, multi-step methods, self-starting,
112
00:11:46,120 --> 00:11:50,910
and not self-starting, we talk about consistency
of this, stability of these methods and we
113
00:11:50,910 --> 00:11:55,990
discussed inner accuracy, energy conservation
properties, and so on and so forth, so what
114
00:11:55,990 --> 00:12:02,060
we will do is we will try to discuss these
issues as we go along, so we'll start with
115
00:12:02,060 --> 00:12:05,450
some simple mathematical
preliminaries, we discussed this earlier in
116
00:12:05,450 --> 00:12:12,120
one of the lectures, we will be using this
so-called gauge notations the O, and capital
117
00:12:12,120 --> 00:12:24,620
O and lowercase o notations, we say that F(x)
is of order G(x) as X tends to A, if as X
118
00:12:24,620 --> 00:12:34,490
tends to A this ratio F(x)/G(x) it remains
finite, okay, and if this goes to 0 then we
119
00:12:34,490 --> 00:12:39,600
say that it is lowercase order G(x) okay,
the mathematical description is given here
120
00:12:39,600 --> 00:12:44,290
I will not get into these details I have provided
it for sake of completion, but I will list
121
00:12:44,290 --> 00:12:51,880
it with few examples. Now let's consider this
function A into X to the power of 7, BX cube
122
00:12:51,880 --> 00:12:57,910
+ CX + D, we make a statement that this function
is order X to the power of 7, as X tends to
123
00:12:57,910 --> 00:13:04,140
infinity, so how do you verify, you would
divide this by X to the power of 7, okay and
124
00:13:04,140 --> 00:13:14,490
you can see that this goes to as X tends to
infinity, this goes to A, right the first
125
00:13:14,490 --> 00:13:19,910
term will be A, second term will be B/X to
the power of 4, C/X to the power of 6, D/X
126
00:13:19,910 --> 00:13:25,050
to the power of 7, so in the denominator as
X tends to infinity all the terms go to 0
127
00:13:25,050 --> 00:13:31,279
except the first term which is A, so we say
that this statement is verified by checking
128
00:13:31,279 --> 00:13:37,269
this, similarly this function is order X to
the power of 0 as X goes to 0, how do I check?
129
00:13:37,269 --> 00:13:47,990
You take extra 0, as X goes to 0 only D will
remain, so this is finite, right, so this
130
00:13:47,990 --> 00:13:52,300
polynomial is order X to the power of 0.
Now A to the power of, X to the power of,
131
00:13:52,300 --> 00:13:57,149
A into X to the power of 7 is order X to the
power of 7 as X goes to 0, how do you verify,
132
00:13:57,149 --> 00:14:03,130
you divide by X to the power of 7 and go take
X to power of 0, I get A, which is finite,
133
00:14:03,130 --> 00:14:10,470
so there are a few more examples you can verify,
you know, get a feel for what these terminologies
134
00:14:10,470 --> 00:14:14,709
mean. So I have a few more examples, let's
quickly run
135
00:14:14,709 --> 00:14:19,950
through this sin X is order X, as X goes to
0, why? You divide sin X/X and take X to 0
136
00:14:19,950 --> 00:14:26,070
we know that this limit is 1, which is finite.
Similarly sin X square is order X square by
137
00:14:26,070 --> 00:14:32,760
the same logic we come to this conclusion,
cos X is order 1 as X goes to 0, because cos
138
00:14:32,760 --> 00:14:38,210
X/1, as X goes to 0 is 1, and so on and so
forth, so there are few more examples I'll
139
00:14:38,210 --> 00:14:45,339
leave it for you to verify.
There is one another mathematical preliminary
140
00:14:45,339 --> 00:14:51,611
that we would be needing, we should quickly
recall what is the mean value theorem, now
141
00:14:51,611 --> 00:14:58,430
consider a function F(x) to be continuous
for X lying between A and B, and differentiable
142
00:14:58,430 --> 00:15:03,980
for, this is the closed interval, this is
the open interval, it is differentiable in
143
00:15:03,980 --> 00:15:09,600
the open interval, then according to the mean
value theorem there exists at least one C
144
00:15:09,600 --> 00:15:16,710
satisfying the condition that C lies between
A and B such that DF/DX at C is exactly given
145
00:15:16,710 --> 00:15:25,139
by this F(B) - F(A) / B - A, that you can
see here, this is my F(x) and this is A and
146
00:15:25,139 --> 00:15:34,720
this is B, and you can see that at this point,
and at this point, if you draw the curve DF/DX
147
00:15:34,720 --> 00:15:43,060
which is given by this, this will be parallel
to the 2, you know this is the approximation
148
00:15:43,060 --> 00:15:49,260
so this will be parallel to this, which is
a tangent to F(x) at those points, so at C1
149
00:15:49,260 --> 00:15:54,190
and C2 this is an exact representation, but
of course for a given F(x) we would not know
150
00:15:54,190 --> 00:15:58,990
where is that C, how many of them are there,
whether there is you know where it is located
151
00:15:58,990 --> 00:16:05,209
we would not know.
Now Taylor's series is something that we are
152
00:16:05,209 --> 00:16:10,810
all familiar with, so let's quickly see the
statement of some of the results associated
153
00:16:10,810 --> 00:16:15,400
with Taylor's series expansions, so again
let F be a function of X with derivatives
154
00:16:15,400 --> 00:16:20,769
of all orders throughout an interval containing
point A. The Taylor's series generated by
155
00:16:20,769 --> 00:16:28,470
F at X = A is given by F(A) + F prime A into
X - A, F double prime A by 2 factorial X - A
156
00:16:28,470 --> 00:16:32,540
whole square and so on and so forth, this
is the well-known Taylor's series. Now
157
00:16:32,540 --> 00:16:37,720
if you truncate the Taylor's series at say
the nth term we get what is known as the Taylor's
158
00:16:37,720 --> 00:16:49,089
polynomial, so that is other conditions remaining
the same I write PN(x) as this, so this is
159
00:16:49,089 --> 00:16:58,040
known as Taylor's polynomial, okay.
Now the corollary to the Taylor's theorem,
160
00:16:58,040 --> 00:17:03,540
if F has derivatives of all orders in an open
interval I containing A, then for each positive
161
00:17:03,540 --> 00:17:10,262
integer N and for each X in the interval I,
F(x) can be written as F(A) + F prime A, X
162
00:17:10,262 --> 00:17:16,970
- A so on and so forth after nth term there
is a reminder term, and this reminder is given
163
00:17:16,970 --> 00:17:24,989
by this for some C between A and X, this again
we are using mean value theorem in writing
164
00:17:24,989 --> 00:17:31,479
this, now if this limit of this remainder
term goes to 0, right for all X in I, we say
165
00:17:31,479 --> 00:17:37,889
that Taylor series generated by F at X = A
converges to F on I, and we write F(x) is
166
00:17:37,889 --> 00:17:44,809
this, okay.
Now we talked about forward difference, backward
167
00:17:44,809 --> 00:17:47,620
difference, and central difference, where
I talked about finite difference, so some
168
00:17:47,620 --> 00:17:51,109
of the examples are forward difference, backward
difference, central difference, so let's see
169
00:17:51,109 --> 00:17:57,500
quickly what it means, so let's consider a
function F(x) and I write the Taylor's expansion
170
00:17:57,500 --> 00:18:04,309
F(x) + H is F(x) + HF prime X + H square/2
factorial F double prime FX and so on and
171
00:18:04,309 --> 00:18:11,350
so forth. Now from this I can write, solve
for F prime(x), F prime(x) will be given by
172
00:18:11,350 --> 00:18:18,809
F(x) + H - F (x) / H, so you are taking these,
other term to the other side you are dividing
173
00:18:18,809 --> 00:18:24,690
by H therefore this becomes H, H/2 factorial
+ H square by 3 factorial and so on and so
174
00:18:24,690 --> 00:18:32,370
forth. So this term is clearly of the order
H, because if you divide by H and take H20
175
00:18:32,370 --> 00:18:38,539
this will be a finite quantity, right, so
this we say that this is a forward difference
176
00:18:38,539 --> 00:18:47,549
approximation to F(x) at X, the backward difference
approximation to derive that we consider instead
177
00:18:47,549 --> 00:18:54,419
of F(x) + H, I consider F(x) - H, so this
will be F(x) - HF prime(x) and so on and so
178
00:18:54,419 --> 00:18:58,539
forth.
Again if I solve for F prime(x) I get this
179
00:18:58,539 --> 00:19:03,940
expression, and divided by, we have divided
by H and we can see that these terms are of
180
00:19:03,940 --> 00:19:10,450
order H, as H goes to 0, so this approximation
is known as backward difference approximation
181
00:19:10,450 --> 00:19:17,399
to F prime(x). To derive the central difference
approximation what we do is we consider F(x)
182
00:19:17,399 --> 00:19:24,859
+ H which is this expansion, and also F(x)
- H which is this expansion, now I subtract
183
00:19:24,859 --> 00:19:31,600
these two, so if moment I subtract these 2,
these 2 get cancelled, these 2 add up, I get
184
00:19:31,600 --> 00:19:37,320
2H F prime(x), and similarly this H square
term will get cancelled, and I will get the
185
00:19:37,320 --> 00:19:42,999
next term will be of the 2H cube / 3 factorial
so on and so forth. Now if I solve for F prime(x)
186
00:19:42,999 --> 00:19:53,190
from this, I get F prime(x) of X is F(x) +
H - F(x) - H / 2H + this term will be divided
187
00:19:53,190 --> 00:20:02,750
by H, so this will be order H square, so as
H goes to 0, okay, so this is the central
188
00:20:02,750 --> 00:20:11,190
difference approximation for derivative of
F at X. Intuitively we can see that if error
189
00:20:11,190 --> 00:20:18,879
is of the order H, and if we have H, that
is the step size H the error gets halved,
190
00:20:18,879 --> 00:20:24,450
similarly if the order, error is of the order
H is square and if we half the step size the
191
00:20:24,450 --> 00:20:34,749
error will get quartered, okay so that is
the advantage of order H square, okay.
192
00:20:34,749 --> 00:20:40,749
Now let's return to the equilibrium equation
in the standard form, we already seen that
193
00:20:40,749 --> 00:20:45,220
this equation can be written in this standard
form in the state space, now what I will do
194
00:20:45,220 --> 00:20:51,039
is, I will replace this X dot by finite difference
approximation and see what we get, so we will
195
00:20:51,039 --> 00:20:57,059
consider a sequence of increasing time instant
0, T1, T2, T capital N, capital TN is TF which
196
00:20:57,059 --> 00:21:02,429
is a final time instant up to which we want
to integrate this equation, and we denote
197
00:21:02,429 --> 00:21:08,160
YNS value of the system state at T = TN by
starting at T = 0 with an initial condition
198
00:21:08,160 --> 00:21:15,239
X naught, now for X dot (t) if I use forward
difference approximation, X dot(t) will be
199
00:21:15,239 --> 00:21:22,790
written as X(t) + D delta T - x of T / delta
T, and this must be equal to A(x(t), t), so
200
00:21:22,790 --> 00:21:30,919
now if I solve for X(t + delta t) which is
YN+1 will be equal to YN, that is X(t) + delta
201
00:21:30,919 --> 00:21:41,070
T is on the other side, delta T YN dot, YN
dot is A(YN), so I'll get YN+1 as YN + delta
202
00:21:41,070 --> 00:21:47,999
TA YN and this order of approximation is order
delta T, okay.
203
00:21:47,999 --> 00:21:55,249
Now if we want now backward difference approximation
we consider time instant T + delta T, and
204
00:21:55,249 --> 00:22:03,009
go back in time and I get X dot(t) + delta
T is X(t) + delta T - X(t) / delta T, so this
205
00:22:03,009 --> 00:22:07,559
looks similar to this but here you must notice
that I am writing this at time instant T,
206
00:22:07,559 --> 00:22:14,669
okay, T + delta T is, I heard of the time
at which I am writing this whereas here the
207
00:22:14,669 --> 00:22:23,000
T is lagging, okay T + delta T is the current
time and whereas T is the previous time instant.
208
00:22:23,000 --> 00:22:28,559
Now here again I can write YN + 1 is YN +
delta T this is important, this is Y dot(N+1)
209
00:22:28,559 --> 00:22:35,570
whereas here it is Y dot(n), so this will
be YN + delta T, A(YN+1), whereas this is
210
00:22:35,570 --> 00:22:44,479
A(YN) and the order of accuracy is still order
of delta T, the central difference scheme
211
00:22:44,479 --> 00:22:52,340
I can write this so again at time T, I write
X dot(t) is
212
00:22:52,340 --> 00:22:58,779
X(T+delta T) - X(t) - delta T/2 delta T, this
must be equal to A(x(t),t), so I want to write
213
00:22:58,779 --> 00:23:06,169
for X(T+delta T) I will get YN + 1, this is
YN - 1, 2 delta T into YN dot, this is YN
214
00:23:06,169 --> 00:23:13,149
dot, so that is YN - 1 + 2 delta TA(YN) and
this approximation is order delta T square,
215
00:23:13,149 --> 00:23:18,240
okay.
Now we can also develop another scheme which
216
00:23:18,240 --> 00:23:23,830
is slightly different from, logic is slightly
different suppose X(t) I write it as 0 to
217
00:23:23,830 --> 00:23:31,169
T, X dot (s) DS, that's a definition of a
derivative, so if you consider T = TN+1, this
218
00:23:31,169 --> 00:23:40,320
will be 0 to TN+1, X dot(s) DS. This 0 to
TN+1 I can write it as 0 to TN + TN to TN+1,
219
00:23:40,320 --> 00:23:48,720
so this first term is nothing but X(TN) and
this is the increment. For this second term
220
00:23:48,720 --> 00:23:55,200
I will use the trapezoidal rule of integration,
so I will therefore I will get this as YN+1
221
00:23:55,200 --> 00:24:04,999
as YN + 1/2 of delta T, YN dot N+1 +YN dot
, okay, so from this I get this as the approximation.
222
00:24:04,999 --> 00:24:13,039
Now for YN + 1, Y dot N+1, I will write A(YN+1),
and YN dot I'll write it as YN + delta TA
223
00:24:13,039 --> 00:24:19,950
YN, so this approximation is known as the
approximation based on trapezoidal rule.
224
00:24:19,950 --> 00:24:25,710
Now we can summarize all this, so in the forward
difference scheme I had this representation
225
00:24:25,710 --> 00:24:32,299
in backward different scheme I had this representation,
in central difference I had this, in trapezoidal
226
00:24:32,299 --> 00:24:38,429
rule I had this. Now a generic form based
on this we can see that a generic form can
227
00:24:38,429 --> 00:24:46,570
be written as YN+1 is alpha 1YN + alpha 2YN
- 1, so on and so forth, alpha M YN+1 - M
228
00:24:46,570 --> 00:24:56,119
+ delta T this is on Y(t) and these are derivatives,
some beta naught Y dot N+1, beta 1 Y dot N,
229
00:24:56,119 --> 00:25:01,669
so on and so forth up to the K-th term. This
is a general form into the, all these schemes
230
00:25:01,669 --> 00:25:09,759
will fit into these by selecting alpha 1,
alpha 2, etcetera,, beta naught, beta 1 etcetera
231
00:25:09,759 --> 00:25:17,321
in a suitable manner I can recover these schemes.
Let's consider the generic form again, you
232
00:25:17,321 --> 00:25:26,299
look at the term
beta naught, so I am writing YN + Y at T = TN+1
233
00:25:26,299 --> 00:25:32,989
and this involves the derivative of the state
at N + 1, if B naught is not 0, okay if B
234
00:25:32,989 --> 00:25:38,179
naught is not 0 therefore YN plus depends
upon derivatives of the state at TN+1 that
235
00:25:38,179 --> 00:25:46,580
is upon Y dot N+1, so such schemes are known
as implicit schemes, okay, otherwise if beta
236
00:25:46,580 --> 00:25:50,779
naught is 0 then the scheme is said to be
explicit.
237
00:25:50,779 --> 00:25:59,100
Now if we now consider alpha 2, alpha 3, and
alpha M as 0, and beta 2 up to beta M are
238
00:25:59,100 --> 00:26:08,570
0, then we say that the scheme is a single
step scheme, otherwise it is called multi-step
239
00:26:08,570 --> 00:26:15,169
scheme. Now the scheme is set to be self-starting
if YN for N less than 0 does not enter the
240
00:26:15,169 --> 00:26:24,659
calculations, okay, right, so if I write this
for N this term is a previous step, okay,
241
00:26:24,659 --> 00:26:33,169
suppose if N = 0, this will be YN and this
is alpha 1 Y naught, there is no problem,
242
00:26:33,169 --> 00:26:39,359
but in certain schemes you will soon see that
I want to start with N+1 = 0 and this I will
243
00:26:39,359 --> 00:26:45,130
be needing Y of - M, and such schemes are
said to be not self-starting okay we need
244
00:26:45,130 --> 00:26:49,509
to see that, so let's quickly see in this
scheme of things how does
245
00:26:49,509 --> 00:26:54,440
this different classification schemes can
be viewed as, suppose if you now consider
246
00:26:54,440 --> 00:26:59,720
forward difference approximation which is
given by this, now you can see that this is
247
00:26:59,720 --> 00:27:10,710
a single step explicit self-starting scheme
okay, explicit because YN plus, we have YN+1
248
00:27:10,710 --> 00:27:15,979
here and there is no derivative of YN+1 on
the other side, backward difference scheme
249
00:27:15,979 --> 00:27:21,820
is again single step implicit, it is implicit
self-starting okay, because you see here YN+1
250
00:27:21,820 --> 00:27:30,929
is here A(YN+1) is nothing but YN+1, derivative
of YN+1 so this is implicit, backward scheme
251
00:27:30,929 --> 00:27:37,159
is implicit. Similarly central difference
approximation is, it is a multi-step method
252
00:27:37,159 --> 00:27:44,719
it is explicit, and it is not a self-starting
scheme because here YN - 1 you will see that
253
00:27:44,719 --> 00:27:50,470
in the due course, for example for N = 0,
this will be YN and I will be needing Y of
254
00:27:50,470 --> 00:27:56,499
- 1, and Y of - 1 is a hypothetical quantity
which is not there as a part of definition
255
00:27:56,499 --> 00:28:00,940
of the original problems of mechanics, problem
of mechanics, so this is not a self-starting
256
00:28:00,940 --> 00:28:06,019
scheme and we need to devise some special
methods to start the solution.
257
00:28:06,019 --> 00:28:10,759
Similarly trapezoidal rule is single step
implicit and it is self-starting, okay, why
258
00:28:10,759 --> 00:28:18,499
implicit, because I have A(YN+1) here which
is YN+1 dot, okay so if system is non-linear
259
00:28:18,499 --> 00:28:23,980
you can see that in implicit schemes you need
to solve a times every time step a algebraic
260
00:28:23,980 --> 00:28:31,049
non-linear equation okay, so the implicit
schemes that way demand computational efforts,
261
00:28:31,049 --> 00:28:35,330
especially if system is non-linear, if system
is linear we can take these terms to the left
262
00:28:35,330 --> 00:28:39,580
side and rearrange the terms which is simply
straight forward, so the difference between
263
00:28:39,580 --> 00:28:44,729
implicit scheme and explicit scheme will be
strongly felt in computational effort for
264
00:28:44,729 --> 00:28:51,889
non-linear systems.
Okay, let's make some observations and remarks
265
00:28:51,889 --> 00:28:59,950
so we have discretized time T naught to TF
in terms of 0, T1, T2, TN which is a non-decreasing
266
00:28:59,950 --> 00:29:05,500
sequence of, actually increasing sequence
of time instance, if we now take TN to be
267
00:29:05,500 --> 00:29:12,219
N delta T, we say that we have a constant
step size, okay, now system states are U(t),
268
00:29:12,219 --> 00:29:19,460
U dot(t) and U double dot(t) in our calculations,
now if we assume that system states at T = TN
269
00:29:19,460 --> 00:29:24,049
are known and that is what we assume then
the problem on hand is to find system states
270
00:29:24,049 --> 00:29:31,609
at T = TN+1, the integration algorithms essentially
achieve this, so these are known at time marching
271
00:29:31,609 --> 00:29:39,349
techniques, so we move from TN to TN+1 and
we want to take state from TN to TN+1, so
272
00:29:39,349 --> 00:29:43,919
these are called time marching techniques.
Now the computation effort is proportional
273
00:29:43,919 --> 00:29:52,119
to the number of time steps, okay to advance
the solution from T to 0 to T to TF, therefore
274
00:29:52,119 --> 00:29:58,320
the step size need not be smaller than what
is essential they say what is the meaning
275
00:29:58,320 --> 00:30:03,899
of something being essential, it is something
with accuracy basically, so what we need to
276
00:30:03,899 --> 00:30:09,450
be concerned about, we need to be concerned
about growth of errors, does error committed
277
00:30:09,450 --> 00:30:17,049
at a time steps AT and grow as time advances,
if it grows we say that the solutions are
278
00:30:17,049 --> 00:30:23,389
unstable, the scheme is unstable otherwise
it is, if it remains constant we say it is
279
00:30:23,389 --> 00:30:29,219
stable, if it goes to 0 we say that it is
as impractically stable, the accuracy of the
280
00:30:29,219 --> 00:30:35,019
solution step, if the errors do not grow,
it does not mean that we get accurate answers,
281
00:30:35,019 --> 00:30:42,070
okay, how to judge accuracy? That is a different
issue, now for example for undamped free vibration
282
00:30:42,070 --> 00:30:46,769
of a single degree freedom system we can check
for amplitude, amplitude should remain constant
283
00:30:46,769 --> 00:30:50,630
and frequency distortions must not be there,
we know that frequency let's say for example
284
00:30:50,630 --> 00:30:59,999
square root K/M do we get that square root
K/M in numerical simulations or not.
285
00:30:59,999 --> 00:31:05,190
The time integration method is the most generally
applicable method to solve equilibrium equations,
286
00:31:05,190 --> 00:31:08,820
we have talked about frequency response function,
based methods, a dynamic stiffness, and so
287
00:31:08,820 --> 00:31:14,159
on and so forth, they are valid only for linear
systems and of course there are additional
288
00:31:14,159 --> 00:31:22,749
requirements that system should reach harmonic
steady-state, but here the system can be nonlinear,
289
00:31:22,749 --> 00:31:30,749
the excitations could be transient, this method
remains applicable. Now for linear systems
290
00:31:30,749 --> 00:31:35,390
the solution actually does not require transformation
to natural coordinates, see for linear system
291
00:31:35,390 --> 00:31:39,379
we can do, we can find natural frequencies,
mode shapes, and uncouple the equation we
292
00:31:39,379 --> 00:31:44,869
have that option, but if you are trying to
integrate the equation numerically it is not
293
00:31:44,869 --> 00:31:50,139
needed to actually perform the uncoupling
of equations of motion, if you can do it,
294
00:31:50,139 --> 00:31:57,059
it helps but it is not essential, so this
would mean that we have greater flexibility
295
00:31:57,059 --> 00:32:03,130
in modeling damping, see we have seen that
if damping is classical there are many advantages
296
00:32:03,130 --> 00:32:10,230
in doing a more analysis based on mode superposition,
but if you are doing direct integration there
297
00:32:10,230 --> 00:32:15,139
is no special requirement on what should be
the damping matrix, because you can directly
298
00:32:15,139 --> 00:32:20,929
integrate.
Now but in practice we prefer to specify damping
299
00:32:20,929 --> 00:32:30,139
in terms of modal ratios, in that case damping
will be specified in terms of as formulation
300
00:32:30,139 --> 00:32:35,940
which employs normal modes, natural coordinates
and normal modes, so if you intend to use
301
00:32:35,940 --> 00:32:42,129
direct integration in such situations you
have to formulate the C matrix from the known
302
00:32:42,129 --> 00:32:47,690
information about modal damping ratios, so
we have discussed how to do that in the previous
303
00:32:47,690 --> 00:32:54,450
class. As I already said the method remains
applicable even when equations of motions
304
00:32:54,450 --> 00:33:02,379
are nonlinear.
Now let's look at modeling of linear multi
305
00:33:02,379 --> 00:33:07,979
degree freedom systems, so how many degrees
of freedom should be included in our model,
306
00:33:07,979 --> 00:33:13,169
we will discuss this in again at a later stage
but we can start asking these questions at
307
00:33:13,169 --> 00:33:19,409
this stage, now the size of the model that
is number of degrees of freedom actually depend
308
00:33:19,409 --> 00:33:28,299
on details of mesh used in the spatial discretization,
now these details need to be chosen such that
309
00:33:28,299 --> 00:33:34,389
the system behavior is represented with acceptable
accuracy over a given frequency range, what
310
00:33:34,389 --> 00:33:39,221
is the frequency range that we should select?
We should look at excitation, suppose you
311
00:33:39,221 --> 00:33:43,169
are dealing with say earthquake like excitation
the frequency range can be up to say about
312
00:33:43,169 --> 00:33:49,970
20 to 30 Hertz from low frequency to up to
that frequency, wind it could be 0 to 2 Hertz
313
00:33:49,970 --> 00:33:55,619
or 3 Hertz, wave also in the same region,
so if you are performing say seismic response
314
00:33:55,619 --> 00:34:00,049
analysis of engineering structure you should
ensure that say within the frequency range
315
00:34:00,049 --> 00:34:05,649
of say 0 to 20 Hertz whatever are the natural
modes that the structure possess all of them
316
00:34:05,649 --> 00:34:14,730
should be modeled accurately, so we have seen
that if you want to capture first say 5 modes
317
00:34:14,730 --> 00:34:19,679
in a model, you should have about 50 degrees
of freedom in your model, that means the spatial
318
00:34:19,679 --> 00:34:25,960
discretization should be fine enough so that
you get a model with 50 degrees of freedom
319
00:34:25,960 --> 00:34:32,639
so that at least one tenth of the natural
frequencies will be computationally trustworthy,
320
00:34:32,639 --> 00:34:39,200
so that would mean if omega max is the highest
frequency we will say that the discretization
321
00:34:39,200 --> 00:34:45,119
scheme should be such that all the modes that
lie in frequency range up to 1.25 Omega max
322
00:34:45,119 --> 00:34:52,080
need to be captured well, so the mesh details
should be such that the model has at least
323
00:34:52,080 --> 00:34:56,950
10R degrees of freedom, where R is the number
of degrees of freedom, the R is the number
324
00:34:56,950 --> 00:35:02,779
of modes that you expect the structure to
have in frequency range up to 1.25 omega max,
325
00:35:02,779 --> 00:35:08,080
a consequence of this is that when you are
using direct integration you will always
326
00:35:08,080 --> 00:35:14,369
be dealing with models which have, which has
spurious higher-order modes, suppose you are
327
00:35:14,369 --> 00:35:18,609
making a model with 100 degrees of freedom
we know that the first 10 modes are likely
328
00:35:18,609 --> 00:35:26,049
to contribute to the response, and modes from
10 to 100 are the higher order spurious modes.
329
00:35:26,049 --> 00:35:30,619
Now it is desirable that the discretization
scheme that is the time discretization scheme
330
00:35:30,619 --> 00:35:38,359
should have some inherent capability to numerically
dissipate the higher order spurious modes,
331
00:35:38,359 --> 00:35:42,849
but not affecting the genuine lower order
modes which contributes significantly to the
332
00:35:42,849 --> 00:35:48,569
response, so this would mean that we look
for certain numerical dissipation characteristics
333
00:35:48,569 --> 00:35:56,480
in our integration schemes when we formulate
the time integration schemes, so what is desirable
334
00:35:56,480 --> 00:36:02,589
is that the lower modes should not be numerically
dissipated, but the higher spurious mode should
335
00:36:02,589 --> 00:36:07,230
be numerically dissipated, okay, we will see
more about that, but this is one of the concern
336
00:36:07,230 --> 00:36:12,900
that we need to appreciate at the outside.
So what we will do now is will take up discussions
337
00:36:12,900 --> 00:36:19,150
on few methods I propose to discuss the following
method, the forward Euler method, the backward
338
00:36:19,150 --> 00:36:24,150
Euler method, the central difference method,
and there are methods known as Newmark's family
339
00:36:24,150 --> 00:36:29,660
of methods, and HHT alpha method, and so on
and so forth, so we will see the development
340
00:36:29,660 --> 00:36:34,670
of these methods and as we go along we will
see what motivates, what are the motivations
341
00:36:34,670 --> 00:36:42,670
to develop different schemes after getting,
for example backward Euler scheme what is
342
00:36:42,670 --> 00:36:46,609
the need to go for central difference method,
and what is the need to go to Newmark's family
343
00:36:46,609 --> 00:36:50,490
of methods when you already have done central
difference method and so on and so forth.
344
00:36:50,490 --> 00:36:55,410
So the questions will be on stability, accuracy,
dissipation of higher or spurious modes and
345
00:36:55,410 --> 00:37:02,539
so on and so forth, methods being implicit
or explicit, these issues also will come up.
346
00:37:02,539 --> 00:37:06,329
So what we, the strategy we will take is,
we will develop the basic formulary for each
347
00:37:06,329 --> 00:37:12,010
of these schemes, and I will provide a Pseudocode
for implementation, then we will analyze the
348
00:37:12,010 --> 00:37:19,200
method for its, how the errors behave in each,
the specific method and then we will illustrate
349
00:37:19,200 --> 00:37:23,710
each of these methods with a few examples
and we'll conclude by discussing the relative
350
00:37:23,710 --> 00:37:29,000
merits of different methods, so this is a
scheme of things that we will follow.
351
00:37:29,000 --> 00:37:34,910
So let me start with the discussion on forward
Euler method, so this is a time axis T naught
352
00:37:34,910 --> 00:37:41,069
is here, TF is a final time instant which
is here, so let's the equation of equilibrium
353
00:37:41,069 --> 00:37:48,500
is MU double dot + U dot + KU = F(t) at time
T, these are the initial conditions. Now at
354
00:37:48,500 --> 00:37:53,569
time T, I will approximate the velocity U
dot(t) using forward difference scheme, which
355
00:37:53,569 --> 00:38:00,570
is U(t) + delta T - U(t) delta T, so from
this I can get UN+ 1 which is U(t) + delta
356
00:38:00,570 --> 00:38:08,819
T as delta T into U dot, that is U dot N +
UN, a similar approximation I can make for
357
00:38:08,819 --> 00:38:17,210
acceleration also, so that will be U dot(t+delta
t) - U dot(t) / delta T, where U dot (N+1)
358
00:38:17,210 --> 00:38:23,200
is obtained as delta T, UN double dot + UN
dot, so you can see here this is an explicit
359
00:38:23,200 --> 00:38:31,039
scheme, because when I am writing the state
at N+1 I do not need acceleration at N+1,
360
00:38:31,039 --> 00:38:39,290
okay, now so let us consider the condition
for equilibrium at TN+1, so this is written
361
00:38:39,290 --> 00:38:49,780
as MU double dot N+1 plus for UN dot UN+1
dot I will use this, okay this is this, then
362
00:38:49,780 --> 00:38:57,290
for UN+1 I will use this, so this is this,
this is equal to FN where F1 is F(tn). Now
363
00:38:57,290 --> 00:39:03,510
I will rearrange these terms what is not known
is UN+1 double dot and that is obtained in
364
00:39:03,510 --> 00:39:10,200
terms of the remaining terms like this.
So how do we implement the method you can
365
00:39:10,200 --> 00:39:14,890
see here the scheme is very clear, I have
to bank on this equation, this equation, and
366
00:39:14,890 --> 00:39:20,381
the final equation here, so I will start with
T = 0, N = 0, we will input the initial condition
367
00:39:20,381 --> 00:39:25,440
U naught and U naught dot, I will also need
the initial acceleration, I will obtain it
368
00:39:25,440 --> 00:39:32,039
from the equation of motion, then I start
with this equation, U double dot N +1 which
369
00:39:32,039 --> 00:39:38,519
is given by this, so here UN double dot, UN
dot etcetera are known, for example when N
370
00:39:38,519 --> 00:39:44,770
= 0 this will be U double dot 1, this will
be U naught double dot, U naught dot and U
371
00:39:44,770 --> 00:39:49,410
naught, and they are already specified here,
U naught and U naught dot are given initial
372
00:39:49,410 --> 00:39:53,319
conditions and U naught dot is obtained from
the governing equilibrium equation, so this
373
00:39:53,319 --> 00:40:00,019
is right. Then I'll find the velocity and
the displacement, then I'll increment N and
374
00:40:00,019 --> 00:40:08,289
if N delta T is more than TF I'll stop, otherwise
I'll go to 2, so this is the simple-minded
375
00:40:08,289 --> 00:40:13,000
outline of how to implement the forward difference
scheme, but if you carefully look at this
376
00:40:13,000 --> 00:40:18,481
it is not necessary to invert M matrix at
every time T, you need not have to multiply
377
00:40:18,481 --> 00:40:24,019
M inverse and C at every time T and so on
and so forth, so we can refine this a bit
378
00:40:24,019 --> 00:40:29,460
so what I will do is I will calculate these
matrices inverse and this product outside
379
00:40:29,460 --> 00:40:33,619
the loop, they need to be done only once,
we need not do every time T, so if you do
380
00:40:33,619 --> 00:40:39,089
that I get a more usable implementation scheme,
so I will input
381
00:40:39,089 --> 00:40:45,650
KMC matrices and delta T there is some algorithmic
parameter and we will also store all the excitation
382
00:40:45,650 --> 00:40:50,720
and then the initial conditions, so I will
define A as M inverse, B as A into C, C is
383
00:40:50,720 --> 00:40:55,859
damping, D is A into K.
Now we will start with N = 0, and we will
384
00:40:55,859 --> 00:41:00,119
accept this U naught and U naught dot form
U naught double dot, so I am not inverting
385
00:41:00,119 --> 00:41:06,309
any of these matrices nor I am multiplying
any matrix now, so UN+1 double dot is given
386
00:41:06,309 --> 00:41:14,410
by this, I get similarly velocity and displacement,
and I will increment time and if I cross the
387
00:41:14,410 --> 00:41:18,559
final time instant I will stop, otherwise
I will restart with this time, so the most
388
00:41:18,559 --> 00:41:25,010
of the calculation gets done in the steps
4 to 6, and I am not repeating any calculation,
389
00:41:25,010 --> 00:41:31,450
if I can avoid I'll avoided all that, okay,
so this is the forward Euler forward difference
390
00:41:31,450 --> 00:41:40,420
scheme, okay, it is logically simple, so what
we can do now is we can think of asking
391
00:41:40,420 --> 00:41:47,960
how does errors behave, when I apply this
scheme to analysis of simple systems.
392
00:41:47,960 --> 00:41:51,569
Now for purpose of illustration I will consider
a single degree freedom system under some
393
00:41:51,569 --> 00:41:57,180
excitation F(t), now this is the scheme I
will get using forward difference scheme for
394
00:41:57,180 --> 00:42:04,599
this problem I will get this is the scheme,
okay, now this can be written as capital YN+
395
00:42:04,599 --> 00:42:14,651
1 as AYN+LN okay, where A is this matrix,
okay, now suppose at the N-th step the state
396
00:42:14,651 --> 00:42:22,289
that is UN+1 and UN+1 dot is contaminated
by noise gamma N, now the question is how
397
00:42:22,289 --> 00:42:28,869
does gamma N behave? So I will substitute
this into this equation so YN+ 1 will be YN+1
398
00:42:28,869 --> 00:42:37,839
+ gamma N+1 = AYN + A gamma N + LN, so now
YN+1 is solution to this equation, therefore
399
00:42:37,839 --> 00:42:42,880
this, this and this cancel out, they are equal,
the sum of these two is equal to this, so
400
00:42:42,880 --> 00:42:50,170
that gets out of reckoning, so the error is
gamma N+1 is A gamma N, so this is the equation
401
00:42:50,170 --> 00:42:57,349
for evolution of the error gamma, so we will
digress a bit now, will return to that shortly,
402
00:42:57,349 --> 00:43:01,190
but let us examine the
nature of this kind of finite difference equations
403
00:43:01,190 --> 00:43:07,040
suppose you consider a scalar equation XN+1
is AXN, and we start with some nonzero initial
404
00:43:07,040 --> 00:43:12,910
conditions, so X1 will be AX naught, X2 will
be AX 1 which is A square X naught, and XN
405
00:43:12,910 --> 00:43:17,510
will be A to the power of N X naught, now
if I am interested in knowing what happens
406
00:43:17,510 --> 00:43:24,450
to XN as N tends to infinity, we can easily
see here that this function will go to 0,
407
00:43:24,450 --> 00:43:29,140
if modulus of A is less than 1, suppose A
is 0.1, A square will be 0.01, A to the power
408
00:43:29,140 --> 00:43:36,089
of 3 will be 0.001 and so on and so forth,
A to the power of 10 will be 10 to, you know
409
00:43:36,089 --> 00:43:42,030
1 into 10 to the power of - 10, so on and
so forth, so this is the condition.
410
00:43:42,030 --> 00:43:52,809
Now if A is 1, + 1 or - 1, XN will stay at
X naught, as N tends to infinity, okay, and
411
00:43:52,809 --> 00:43:58,320
similarly if A is greater than 1, suppose
A is 2, initially it will be 2X naught, then
412
00:43:58,320 --> 00:44:02,440
4X naught, then 16 X naught and so on and
so forth, it goes to infinity as N tends to
413
00:44:02,440 --> 00:44:07,710
infinity, so the behavior of this XN as N
tends to infinity is crucially governed by
414
00:44:07,710 --> 00:44:13,670
the value of A, if the absolute value of A
is less than 1, XN goes to 0 as N tends to
415
00:44:13,670 --> 00:44:17,339
infinity, this is scalar equation. Suppose
now you
416
00:44:17,339 --> 00:44:24,109
consider S cross 1 vector equation, suppose
this is XN + 1 = AXN, where A is S cross S
417
00:44:24,109 --> 00:44:31,940
matrix, okay. Now X1 is AX naught, X2 is AX1
which is A square X naught, similarly XN will
418
00:44:31,940 --> 00:44:39,721
be A to the power of N, X naught. Now let's
do the following, let us introduce a transformation
419
00:44:39,721 --> 00:44:51,710
XN = phi ZN, so this equation for XN+1 will
be phi ZN +1 will be equal to A phi ZN.
420
00:44:51,710 --> 00:44:56,349
Now let phi be such that phi transpose phi
is identity, and phi transpose A phi is a
421
00:44:56,349 --> 00:45:01,779
diagonal matrix, say if A is such that this
is possible, how do we select, we can select
422
00:45:01,779 --> 00:45:06,240
phi by solving the eigenvalue problem associated
with A, that we have seen in few vacations
423
00:45:06,240 --> 00:45:12,329
how to do that, now suppose A, B such that
this is possible, then I can write phi ZN
424
00:45:12,329 --> 00:45:17,460
+ 1 = A phi ZN and pre
multiply by phi transpose I get this equation,
425
00:45:17,460 --> 00:45:24,799
now since phi transpose phi is I, I get ZN+1
and phi transpose A phi is a diagonal matrix
426
00:45:24,799 --> 00:45:32,460
of lambda I, I get this, okay, so Z is a S
cross 1 vector therefore I can write the K-th
427
00:45:32,460 --> 00:45:41,190
element of that, ZN +1K as lambda K ZNK where
K runs from 1 to S, so this is the equation
428
00:45:41,190 --> 00:45:51,750
I get, ZN+1 K, now the behavior of this individual
Z as K tends to infinity, depends on absolute
429
00:45:51,750 --> 00:45:56,819
value of lambda K, the K-th eigenvalue, right,
now if you are interested in J-th component
430
00:45:56,819 --> 00:46:06,630
of your vector X, that is given by this summation,
okay, so the behavior of this as N tends to
431
00:46:06,630 --> 00:46:13,010
infinity crucially depends on the highest
eigenvalue that if you rank order this lambda
432
00:46:13,010 --> 00:46:20,290
K depending on the absolute value of its,
that is the absolute value of lambda K, we
433
00:46:20,290 --> 00:46:26,059
are interested in, if you are interested in
knowing what happens to this modulus of XNJ
434
00:46:26,059 --> 00:46:32,490
as N tends to infinity, for this to go to
0 we require that the maximum value of this
435
00:46:32,490 --> 00:46:38,640
lambda K must less than 1, okay, because each
1 is a scalar equation, for each scalar equation
436
00:46:38,640 --> 00:46:43,700
to go to 0, the each of the eigenvalue, absolute
value should be less than 1, so if the highest
437
00:46:43,700 --> 00:46:50,420
eigenvalue is a modulus is less than 1, then
XNJ will go to 0, as N tends to infinity,
438
00:46:50,420 --> 00:46:57,750
so therefore if you are interested in behavior
of XNJ as N tends to infinity, that is controlled
439
00:46:57,750 --> 00:47:01,849
by the highest modulus of the eigenvalues
of A, okay.
440
00:47:01,849 --> 00:47:09,050
Now alternatively we can also consider XN
as AXN + 1 into A into XN-1, I can write it
441
00:47:09,050 --> 00:47:16,000
as A to the power of N, X naught, now if we
seek the solution now of the form XN is alpha
442
00:47:16,000 --> 00:47:21,650
N into X naught, suppose if you seek this
solution, for some value of alpha the such
443
00:47:21,650 --> 00:47:28,150
a solution may be possible, for which value
of alpha this is possible we have to see,
444
00:47:28,150 --> 00:47:32,820
this alpha could be scalar, the scalar which
can, in general be complex value, so now you
445
00:47:32,820 --> 00:47:38,010
substitute this into this equation, I get
XN is alpha N X naught, therefore XN+1 is
446
00:47:38,010 --> 00:47:44,029
alpha N+1 X naught. Now if you consider this
equation XN + 1 is AXN, A to the power of
447
00:47:44,029 --> 00:47:50,509
X naught and make these substitutions I get
alpha N + 1 X naught is A alpha N X naught,
448
00:47:50,509 --> 00:47:55,869
so from this I get AX naught is alpha X naught,
so this would mean for non-trivial solutions
449
00:47:55,869 --> 00:48:02,359
the determinant of A - alpha I must be equal
to 0, and if we represent the roots of this
450
00:48:02,359 --> 00:48:08,289
equation which are the eigenvalues of A, suppose
I-th route is written as alpha 0I exponential
451
00:48:08,289 --> 00:48:13,180
I theta I, this is a complex number so I can
always represent like this.
452
00:48:13,180 --> 00:48:21,210
The condition for XNJ modulus of that as N
tends to infinity to go to 0 is given by the
453
00:48:21,210 --> 00:48:29,279
condition that the maximum value of, absolute
value of alpha I must be less than 1, so that
454
00:48:29,279 --> 00:48:35,319
would mean actually this quantity which is
a maximum value of modulus of the eigenvalue
455
00:48:35,319 --> 00:48:42,079
for I running from 1 to S is known as spectral
radius of A, okay so what we are asking is
456
00:48:42,079 --> 00:48:46,670
the spectral radius
of A must be less than 1 so that would mean
457
00:48:46,670 --> 00:48:53,670
if you look at the eigenvalue in the complex
plane we have real part here, imaginary part
458
00:48:53,670 --> 00:48:59,509
here, and this is the so-called unit circle,
okay it has radius equal to 1, what we want
459
00:48:59,509 --> 00:49:04,670
is the root should be inside this unit circle
for stability asymptotic stability, if it
460
00:49:04,670 --> 00:49:11,720
is right on the unit circle the errors don't
grow, but it is stable, but if it is outside
461
00:49:11,720 --> 00:49:20,849
the unit circle the errors would grow, okay,
so this is what we need to verify if you are
462
00:49:20,849 --> 00:49:26,079
interested in studying the growth of errors.
One more thing that we should notice here
463
00:49:26,079 --> 00:49:33,220
is the only question we need to answer is
whether roots lie within unit circle or not,
464
00:49:33,220 --> 00:49:37,560
we are not so much interested in knowing the
absolute value of the route, we are simply
465
00:49:37,560 --> 00:49:44,099
interested in knowing a qualitative feature
associated with the roots. So we do not need
466
00:49:44,099 --> 00:49:48,559
the value of eigenvalues, we need to only
ascertain if all the roots of the characteristic
467
00:49:48,559 --> 00:49:55,500
equation lie within the unit circle in the
complex plane, given that we are asking a
468
00:49:55,500 --> 00:50:03,730
simple question then actually not asking,
we are not asking, we are not interested in
469
00:50:03,730 --> 00:50:08,079
determining the value of all the eigenvalues,
we are simply interested in knowing whether
470
00:50:08,079 --> 00:50:13,010
all eigenvalues lie within the unit circle
or not, to answer this simpler question we
471
00:50:13,010 --> 00:50:19,119
can develop an alternative method which is
computationally simpler than finding all the
472
00:50:19,119 --> 00:50:24,410
roots and finding out whether they lie within
unit circle or not, so that takes us to discussion
473
00:50:24,410 --> 00:50:31,089
on what is known as Jury's criterion where
given a polynomial by using Jury's criterion
474
00:50:31,089 --> 00:50:37,619
we can verify whether roots of a polynomial
lie within a unit circle or not, so that has
475
00:50:37,619 --> 00:50:45,430
bearing on discussion on stability of the
finite difference time integration schemes,
476
00:50:45,430 --> 00:50:50,079
so what we will do is, we will take up that
issue in the next lecture, and we will conclude
477
00:50:50,079 --> 00:50:51,709
this lecture at this stage.