1
00:00:04,279 --> 00:00:21,599
cture we started talking about methods for
time integration of equation of motion and
2
00:00:21,599 --> 00:00:29,149
we introduced the basic terminologies and
the role of finite difference schemes in developing
3
00:00:29,149 --> 00:00:34,140
these methods, we started talking about the
forward difference method or the so
4
00:00:34,140 --> 00:00:41,260
called forward Euler method, so here we consider
the equilibrium equation, we are restricting
5
00:00:41,260 --> 00:00:47,760
our discussion to linear systems MU double
dot + CU dot + KU = F(t) with some prescribed
6
00:00:47,760 --> 00:00:54,999
initial conditions, and the time duration
of interest is between 0 to TF. Forward difference
7
00:00:54,999 --> 00:01:02,820
approximation to the derivatives lead to these
equations for displacements, velocities and
8
00:01:02,820 --> 00:01:12,820
I mean substituting that we will get this
equation for the evolution of the displacement.
9
00:01:12,820 --> 00:01:20,291
We were discussing the issue of how errors
grow as the time marching takes place, so
10
00:01:20,291 --> 00:01:28,640
to illustrate that we considered a single
degree freedom system which is driven by F(t),
11
00:01:28,640 --> 00:01:34,750
the idea is that even if you have a multi
degree freedom system moment you do normal
12
00:01:34,750 --> 00:01:40,159
mode analysis and transform the coordinates
to the natural coordinates, in the uncoupled
13
00:01:40,159 --> 00:01:44,850
state the equations of motion would be a set
of single degree freedom systems, so to understand
14
00:01:44,850 --> 00:01:52,550
growth of errors it is enough to focus on
how the errors would behave when we developed
15
00:01:52,550 --> 00:01:57,429
the scheme for a single degree freedom system,
so this is the equation based on forward difference
16
00:01:57,429 --> 00:02:06,060
scheme for this oscillator, and we can put
it in a matrix form capital Y N+1 is the vector
17
00:02:06,060 --> 00:02:13,320
consisting of a 2 cross 1 column consisting
of UN+1 and U dot N+1, and LN is this 0 delta
18
00:02:13,320 --> 00:02:19,000
T into FN, so A is this matrix.
Now what we do is at the nth step we assume
19
00:02:19,000 --> 00:02:25,209
that an error gamma N is introduced, because
of a truncation and round off and issues like
20
00:02:25,209 --> 00:02:32,930
that, now we substitute this equation into
the governing equation and we want to develop
21
00:02:32,930 --> 00:02:41,650
a scheme on understanding how the errors grow,
so from this equation we see that the governing
22
00:02:41,650 --> 00:02:50,490
equation for evolution of error is given by
gamma N+1 is A into gamma N. Now to be able
23
00:02:50,490 --> 00:02:54,470
to understand how the solutions to this type
of equations behave
24
00:02:54,470 --> 00:03:03,590
we considered, any equation of the form XN+1
is AXN, and we introduce the transformation
25
00:03:03,590 --> 00:03:09,659
on X which is phi Z, we introduce a new coordinate
system Z by introducing the transformation
26
00:03:09,659 --> 00:03:16,040
X = phi Z, and phi was selected to be the
matrix of eigenvectors of A matrix, and we
27
00:03:16,040 --> 00:03:22,780
showed in the last class that the condition
for the X the modulus of X, the J-th element
28
00:03:22,780 --> 00:03:28,420
of that to go to 0 as N tends to infinity
is that the maximum eigenvalue, the absolute
29
00:03:28,420 --> 00:03:35,099
value of maximum eigenvalue of matrix A must
be less than 1, so this quantity is known
30
00:03:35,099 --> 00:03:39,450
as the spectral radius of matrix A, so we
want that the spectral radius must be less
31
00:03:39,450 --> 00:03:49,040
than 1, an alternative way of looking at it
is if we plot the eigenvalues on the complex
32
00:03:49,040 --> 00:03:53,400
plane this is real part of the eigenvalue,
and imaginary part of the eigenvalue, and
33
00:03:53,400 --> 00:03:59,349
if this is a unit circle we want that for
stability the roots must lie within the unit
34
00:03:59,349 --> 00:04:05,489
circle, that is for asymptotic stability for
errors to go to 0, if errors not to grow it
35
00:04:05,489 --> 00:04:10,629
is sufficient if the roots lie on the unit
circle, and if they are outside that the solution
36
00:04:10,629 --> 00:04:15,840
would be unstable in the sense the errors
grow as we march in time.
37
00:04:15,840 --> 00:04:22,460
Now let us return to the analysis of forward
difference scheme and we had this equation
38
00:04:22,460 --> 00:04:31,270
YN+1 is AYN and from this I developed this
equation gamma N+1 = A gamma N. Now let us
39
00:04:31,270 --> 00:04:41,900
rewrite this in terms of gamma naught, so
we get that gamma N+ 1 will be A to the power
40
00:04:41,900 --> 00:04:50,139
of N into gamma naught, how do we see that?
If you write now gamma 1 this will be A into
41
00:04:50,139 --> 00:04:57,110
gamma naught, if I write gamma 2 it will be
A into gamma 1 which is A square into gamma
42
00:04:57,110 --> 00:05:06,419
naught, so continuing in this manner we get
this equation. Now this equation we can see
43
00:05:06,419 --> 00:05:13,710
that it describes how error introduced at
N = 0 grows in time, now what we do is we
44
00:05:13,710 --> 00:05:18,130
seek the solution of this equation in the
form gamma N is equal to alpha N gamma naught,
45
00:05:18,130 --> 00:05:23,520
where alpha is a scalar quantity which could
be complex valued, so upon substituting into
46
00:05:23,520 --> 00:05:32,020
this equation I get alpha N+1 gamma naught
is A alpha N gamma naught or in other words
47
00:05:32,020 --> 00:05:38,780
for non-trivial solutions we see that this
leads to an eigenvalue problem as stated here,
48
00:05:38,780 --> 00:05:44,259
and for non-trivial solution the determinant
of A - lambda I must be 0, so the roots of
49
00:05:44,259 --> 00:05:49,380
this equation if I write it as alpha naught
E raise to I theta, because we expect alpha
50
00:05:49,380 --> 00:05:54,740
to be typically complex valued, so I can write
the complex number in this form where alpha
51
00:05:54,740 --> 00:05:59,460
naught is the amplitude, and theta is the
phase, so gamma N will be alpha naught N exponential
52
00:05:59,460 --> 00:06:04,310
I theta into gamma naught.
So now if we are interested in the asymptotic
53
00:06:04,310 --> 00:06:14,400
behavior of gamma N, as N tends to infinity
we want that if this amplitude of these eigenvalues
54
00:06:14,400 --> 00:06:18,790
satisfy this requirement that the maximum
value of this absolute value of this eigenvalue
55
00:06:18,790 --> 00:06:27,110
is less than 1 then it guarantees that the
solution decays to 0, so the behavior of the
56
00:06:27,110 --> 00:06:31,479
growth of errors in forward difference scheme
is therefore governed by eigenvalues of this
57
00:06:31,479 --> 00:06:38,660
matrix A, what are the parameters that enter
this matrix A? Delta T omega and ETA, now
58
00:06:38,660 --> 00:06:46,889
ETA and omega are parameters dictated by the
problem on hand, on that we don't have any
59
00:06:46,889 --> 00:06:51,460
control when we are developing the numerical
scheme, so the only algorithmic parameter
60
00:06:51,460 --> 00:07:00,590
that we have is delta T. So now clearly for
this condition to be met it will place the
61
00:07:00,590 --> 00:07:07,129
restrictions on delta T, now if we were to
see that we will actually find the roots so
62
00:07:07,129 --> 00:07:10,960
if I
now use this equation determinant of A - alpha
63
00:07:10,960 --> 00:07:16,190
A is 0 I get this characteristic equation,
and this is a quadratic equation in this case
64
00:07:16,190 --> 00:07:21,810
we can quickly solve that and if we solve
this that we see that these two, the two roots
65
00:07:21,810 --> 00:07:29,880
are given by this. So upon simplification
we can show that the roots are given by 1
66
00:07:29,880 --> 00:07:35,400
- ETA omega T + - I Omega D delta T, where
omega D is the damped natural frequency, omega
67
00:07:35,400 --> 00:07:40,139
into 1 - ETA square, so bit of algebra will
lead to this.
68
00:07:40,139 --> 00:07:45,139
Now therefore what is the maximum value of
the amplitude of the eigenvalues, we can compute
69
00:07:45,139 --> 00:07:50,699
that and we get that to be this 1 + omega
square delta T square - 2 ETA omega delta
70
00:07:50,699 --> 00:07:55,600
T, for stability, asymptotic stability we
want that this maximum value should be less
71
00:07:55,600 --> 00:08:01,740
than 1, so if we impose this condition we
get the requirement that delta T should be
72
00:08:01,740 --> 00:08:09,270
less than 2 ETA /omega for the condition of
stability to be satisfied, asymptotic stability
73
00:08:09,270 --> 00:08:14,470
to be satisfied, so therefore we conclude
that the forward difference scheme is conditionally
74
00:08:14,470 --> 00:08:21,690
stable for damped systems, this critical step
size that is 2 ETA / Omega this one as we
75
00:08:21,690 --> 00:08:30,750
see here is a function of ETA and omega so
if system has damping then delta, the step
76
00:08:30,750 --> 00:08:37,950
size that you should be less than 2 ETA / omega
for to ensure that the errors don't grow in
77
00:08:37,950 --> 00:08:43,740
time, but if system is undamped we see that
there is no step size which guarantees stability,
78
00:08:43,740 --> 00:08:49,990
so this leads to the conclusion that the scheme
is unstable for all choices of step size for
79
00:08:49,990 --> 00:08:55,460
undamped systems. So in an undamped system
this would mean that the scheme introduces
80
00:08:55,460 --> 00:09:01,440
an artificial negative damping into the discretized
model, which alters the basic physics of the
81
00:09:01,440 --> 00:09:05,900
problem, so that alteration in the behavior
of the system is essentially due to the way
82
00:09:05,900 --> 00:09:12,040
we have discretized the governing differential
equation, so this would indicate that for
83
00:09:12,040 --> 00:09:19,290
all different scheme has many pitfalls and
we should be cautious if we want to use this.
84
00:09:19,290 --> 00:09:24,910
Now let's look at another first-order method
that is the backward Euler method, here instead
85
00:09:24,910 --> 00:09:28,960
of using the forward difference scheme for
approximating derivatives we use a backward
86
00:09:28,960 --> 00:09:37,000
difference scheme, as you would see this switch
you know results in a substantial change in
87
00:09:37,000 --> 00:09:42,220
the qualitative behavior of the integration
scheme, so again to see that let's again start
88
00:09:42,220 --> 00:09:50,070
with our multi degree damped forced oscillation
problem, MU double dot + CU dot + KU = F with
89
00:09:50,070 --> 00:09:56,250
prescribed initial conditions and time range
of interest from 0 to TF. Now let's consider
90
00:09:56,250 --> 00:10:02,930
the time instant T plus delta T, and we write
now the derivative at that point it will involve
91
00:10:02,930 --> 00:10:07,650
the displacement at T + delta T and at T,
so this is the approximation according to
92
00:10:07,650 --> 00:10:13,470
the backward difference scheme. Similarly
acceleration can be written in terms of velocities
93
00:10:13,470 --> 00:10:18,820
and suppose if you use now the expression
for U dot at T + delta T we see that if I
94
00:10:18,820 --> 00:10:29,670
solve for UN+1 it will be delta T into UN
dot N+1 +UN, similarly I get from the backward
95
00:10:29,670 --> 00:10:36,300
difference approximation to the acceleration
I get U dot N+1 is delta T, U double dot N+1
96
00:10:36,300 --> 00:10:42,250
+ UN dot. Now clearly you see that this is
an implicit scheme because the right hand
97
00:10:42,250 --> 00:10:51,820
side contains unknowns U dot N+1 and U double
dot N+1, so it's an implicit scheme.
98
00:10:51,820 --> 00:10:56,360
Now let's return to the governing equation
what we done, what we did just now was to
99
00:10:56,360 --> 00:11:00,000
obtain forward difference approximation to
the displacement and velocity, now let's written
100
00:11:00,000 --> 00:11:06,110
to the governing equation of motion and substitute
those approximations, so for UN dot I get
101
00:11:06,110 --> 00:11:12,740
this for ND1 I get this, that is I am considering
equilibrium equation at time TN+1, now this
102
00:11:12,740 --> 00:11:18,570
is the equilibrium equation as per the forward
difference scheme and by rearranging these
103
00:11:18,570 --> 00:11:27,170
terms and collecting coefficients of U double
dot N+1, I get the time evolution of U double
104
00:11:27,170 --> 00:11:33,750
dot N+1 is given by this equation, okay, so
in this let me again emphasize we have got,
105
00:11:33,750 --> 00:11:39,101
we have used backward difference scheme for
approximating velocity and displacement, velocity
106
00:11:39,101 --> 00:11:42,520
and acceleration terms. So how do we implement
this?
107
00:11:42,520 --> 00:11:49,010
In A version of this implementation is you
start with T = 0 and set N = 0, read U naught
108
00:11:49,010 --> 00:11:54,050
and U naught dot that is this information
we have at T = 0 and use this equation in
109
00:11:54,050 --> 00:11:58,770
this sequence, first we write the expression
for acceleration that is given by this, so
110
00:11:58,770 --> 00:12:05,550
when N = 0, I am writing U double dot 1, so
on the right hand side I have U naught dot
111
00:12:05,550 --> 00:12:13,050
which is known, U naught which is known, FN+1
is known, because F(t) is given to us. Then
112
00:12:13,050 --> 00:12:20,650
I go to the expression for U dot (N+1), I
would have computed this U double dot N+1
113
00:12:20,650 --> 00:12:25,550
in the previous step so I can use that, so
similarly when I come to UN+1, I would have
114
00:12:25,550 --> 00:12:30,640
computed U dot N+1 in the previous step I
will use that, then I increment N and I will
115
00:12:30,640 --> 00:12:36,010
stop this process if I cross TF, otherwise
I return to 2 and go through this.
116
00:12:36,010 --> 00:12:46,890
Now a more useful scheme up for implementation
is if we go back here we notice that at every
117
00:12:46,890 --> 00:12:52,410
time step we need not invert these matrices,
especially if delta T is constant for all
118
00:12:52,410 --> 00:12:58,740
time steps then this can be inverted only
once and that can be stored outside this loop,
119
00:12:58,740 --> 00:13:03,620
similarly the product of this matrix with
this C + delta TK, and K can be evaluated
120
00:13:03,620 --> 00:13:09,690
only once and kept outside the loop, this
is of course valid only if delta T is independent
121
00:13:09,690 --> 00:13:18,450
of steps N, for all N where I am using constant
delta T, so in this case I will first store
122
00:13:18,450 --> 00:13:22,380
all these matrices and initial conditions
and the forcing functions and then I will
123
00:13:22,380 --> 00:13:28,630
evaluate this matrix M + delta TC + delta
T square K inverse and product of that into
124
00:13:28,630 --> 00:13:35,690
C + delta TK and K, I will evaluate outside
the loop, then I will start my count N = 0,
125
00:13:35,690 --> 00:13:41,440
and then use these 3 relations, so here I
am using now BD and A, which have been computed
126
00:13:41,440 --> 00:13:48,060
outside the loop. Now increment N, and I will
stop if I reach the final time instant otherwise
127
00:13:48,060 --> 00:13:53,670
I return to this, so that means as I am running
through these steps 4, 5, 6 I am not evaluating
128
00:13:53,670 --> 00:14:03,890
A, B, and D, repeatedly so this will obviously
enhance the speed of computation.
129
00:14:03,890 --> 00:14:12,270
Now let's see how, what are the stability
characteristics of the Euler backward Euler
130
00:14:12,270 --> 00:14:20,140
method, now in the previous study on stability
of forward difference scheme we notice that
131
00:14:20,140 --> 00:14:26,460
the presence of forcing function doesn't affect
the matter solidity to stability, so we need
132
00:14:26,460 --> 00:14:31,080
not consider a forced oscillator but we can
consider now in single degree freedom system
133
00:14:31,080 --> 00:14:36,950
which is unforced, so the finite difference
schemes that we are using are listed here
134
00:14:36,950 --> 00:14:46,220
and by reorganizing these terms I can write
the evolution equation as some matrix which
135
00:14:46,220 --> 00:14:53,790
is A in our earlier case, and UN+1, U dot
N+1 is a previous step, so the next step if
136
00:14:53,790 --> 00:15:00,800
you want you have to invert this matrix and
get this, okay, now this matrix we have to
137
00:15:00,800 --> 00:15:07,380
know find the eigenvalues of this matrix,
that means this 1 - delta T this inverse,
138
00:15:07,380 --> 00:15:13,380
so this amplification matrix here is inverse
of this, this is A.
139
00:15:13,380 --> 00:15:23,740
Now this matrix
we already encountered in the application
140
00:15:23,740 --> 00:15:29,330
of forward difference scheme, so we have found
eigenvalues of A there, now we want the eigenvalues
141
00:15:29,330 --> 00:15:35,830
of inverse of that so we can see that eigenvalues
of A can be obtained as reciprocals of eigenvalues
142
00:15:35,830 --> 00:15:45,000
of A inverse so that you can verify why that
is true. Then the characteristic equation
143
00:15:45,000 --> 00:15:50,930
we have already solved and these were the
roots that we obtained for A matrix in the
144
00:15:50,930 --> 00:15:56,470
earlier case, so this is this, and the reciprocals
of that would be this now.
145
00:15:56,470 --> 00:16:03,910
So now if you separate real and imaginary
parts I get this expression for their 2 eigenvalues,
146
00:16:03,910 --> 00:16:10,500
so the spectral radius which is maximum value
of lambda bar is given by this expression.
147
00:16:10,500 --> 00:16:16,020
Now we want that this spectral radius should
be less than 1 for stability, now you observe
148
00:16:16,020 --> 00:16:22,360
this there is 1 here and there are other quantities
in the denominator which are all positive,
149
00:16:22,360 --> 00:16:28,610
so that would mean the denominator here is
always greater than 1 for any choice of delta
150
00:16:28,610 --> 00:16:35,660
T, so that also means that the spectral radius
is less than 1 for any choice of delta T,
151
00:16:35,660 --> 00:16:41,080
so this means the integration scheme is unconditionally
stable so you can select delta T in whichever
152
00:16:41,080 --> 00:16:49,030
manner you want the errors won't grow.
Now let me show some numerical results, so
153
00:16:49,030 --> 00:16:54,770
let us consider an undamped single degree
freedom system with frequency, a natural frequency
154
00:16:54,770 --> 00:16:59,470
2 phi, the period will be 1 second, and we
are starting from 0 displacement and unit
155
00:16:59,470 --> 00:17:04,400
velocity that means we are basically computing
the impulse response function of the system,
156
00:17:04,400 --> 00:17:08,350
so let us start with a forward difference
scheme and I will take about a 100 points
157
00:17:08,350 --> 00:17:13,680
in one cycle of oscillation, so notwithstanding
that you see that the blue line that I am
158
00:17:13,680 --> 00:17:19,600
showing here is the result obtained from numerical
simulations and we already seen that forward
159
00:17:19,600 --> 00:17:25,430
difference scheme is always unstable for undamped
system and that you see here the error is
160
00:17:25,430 --> 00:17:33,740
growing, okay, so the system, the discretized
system has a negative damping.
161
00:17:33,740 --> 00:17:38,660
Now let us introduce damping, suppose about
3% damping is introduced and we know the delta
162
00:17:38,660 --> 00:17:45,430
T critical we can find out 2 ETA / omega is
a critical step size, and if I now take for
163
00:17:45,430 --> 00:17:51,670
illustration delta T to be 2 times the delta
T critical I expect solutions to be, solutions
164
00:17:51,670 --> 00:17:56,430
to grow as time passes and that's what happens
here, because forward difference scheme is
165
00:17:56,430 --> 00:18:02,190
unstable if delta TX is delta T critical.
On the other hand if it is exactly equal to
166
00:18:02,190 --> 00:18:08,630
delta T critical you see that the blue line
amplitude remains constant, but still it's
167
00:18:08,630 --> 00:18:17,630
not a good approximation for the red line
which is the exact solution, okay. Now the
168
00:18:17,630 --> 00:18:24,850
physical system here has a natural damping
but my numerical solution appears as if the
169
00:18:24,850 --> 00:18:32,730
system is undamped, so this again an influence
of discretization which is undesirable. Now
170
00:18:32,730 --> 00:18:36,430
how do you remedy this if you are using the
same method? You reduce step
171
00:18:36,430 --> 00:18:43,120
size, so I have halved the step size delta
T critical by 0.5, so again red is exact and
172
00:18:43,120 --> 00:18:49,180
blue is the approximate solution. Now the
blue lines start showing the DK that we expect
173
00:18:49,180 --> 00:18:57,650
but still the results are not accurate enough.
Now I take delta T critical by 20, now I see
174
00:18:57,650 --> 00:19:05,900
that the two solutions start comparing quite
well so the forward difference scheme in this
175
00:19:05,900 --> 00:19:14,970
case with the right choice of step size produces
stable and acceptable solutions.
176
00:19:14,970 --> 00:19:22,500
Now let me take the backward difference scheme,
so it is a damped system, this system the
177
00:19:22,500 --> 00:19:28,740
backward difference scheme is stable for any
delta T, so now if I take 25 points in a cycle
178
00:19:28,740 --> 00:19:34,060
it looks adequate from an engineering point
of view, but the numerical solutions show
179
00:19:34,060 --> 00:19:39,820
that that's not quite adequate because the
red line is the exact solution and blue is
180
00:19:39,820 --> 00:19:44,650
this, so it's not showing, it's showing where
the comparison is pretty bad.
181
00:19:44,650 --> 00:19:50,190
Now instead of 25 points, suppose if I take
50 points in a second there is some improvement
182
00:19:50,190 --> 00:19:59,050
but still not impressive 100 points still
not good enough, 200 points there is improvement
183
00:19:59,050 --> 00:20:05,160
but still
not acceptable, 400 points we're nearing acceptance
184
00:20:05,160 --> 00:20:12,800
but still not quite there, 800 points seems
okay but still there are differences, 1600
185
00:20:12,800 --> 00:20:17,630
points the comparison is pretty good so that
means if you want to use backward difference
186
00:20:17,630 --> 00:20:24,250
scheme you have to take about 1600 points
in a cycle for this case to get a good approximation
187
00:20:24,250 --> 00:20:28,770
to the impulse response of the system, this
also points towards the fact that these methods
188
00:20:28,770 --> 00:20:38,720
are a first-order method so you need very
small step size to get acceptable solutions.
189
00:20:38,720 --> 00:20:43,980
Now I have plotted the norm of the error for
the backward difference scheme and the norm
190
00:20:43,980 --> 00:20:49,220
that I have used is defined here, since we
know the exact solution it is X exact - X
191
00:20:49,220 --> 00:20:56,550
numerical whole square + X dot exact - X dot
numerical, and there is a to make units consistent
192
00:20:56,550 --> 00:21:01,240
I have divided by omega, so this is the norm
I have shown and this is with 200 points in
193
00:21:01,240 --> 00:21:08,930
a cycle the value of this error is about 0.035
and with 1600 points the error comes to about
194
00:21:08,930 --> 00:21:20,390
you know 5 into 10 to the power of - 3.
So this is how the backward and forward difference
195
00:21:20,390 --> 00:21:24,140
schemes perform with respect to a single degree
freedom system, so but our objective is really
196
00:21:24,140 --> 00:21:28,890
to apply these methods to multi degree forced
oscillation problems. Now we can make few
197
00:21:28,890 --> 00:21:35,840
remarks now, what are the modes of unsatisfactory
performance of a given scheme? Lack of stability
198
00:21:35,840 --> 00:21:41,300
that means errors grow, the numerical dissipation
and amplitude distortion even if errors don't
199
00:21:41,300 --> 00:21:48,530
grow there may be undesirable damping induced
by the integration scheme which may distort
200
00:21:48,530 --> 00:21:54,250
the amplitude of the response and it also
could lead to distortion of the frequency,
201
00:21:54,250 --> 00:22:01,830
then low order of accuracy resulting in requiring
smaller step size to get acceptable solutions,
202
00:22:01,830 --> 00:22:06,850
so these are some unsatisfactory features
of the integration scheme that we have seen
203
00:22:06,850 --> 00:22:13,130
so far. And one more thing that we have to
notice is stability of the scheme and does
204
00:22:13,130 --> 00:22:17,851
not guarantee that solutions obtain would
possess acceptable accuracy, okay, we have
205
00:22:17,851 --> 00:22:23,110
seen now we have used step sizes which are
quite satisfactory from the point of view
206
00:22:23,110 --> 00:22:29,270
of stability requirement for both forward
and backward cases for difference schemes,
207
00:22:29,270 --> 00:22:38,230
but accuracy was not always acceptable.
Now when we apply the integration schemes
208
00:22:38,230 --> 00:22:43,930
to multi degree freedom system we need to
make you know take into account a few facts
209
00:22:43,930 --> 00:22:51,090
so I will try to start discussing that, suppose
we have equation MU double dot + CU dot +
210
00:22:51,090 --> 00:22:58,110
K = F(t) with specified initial conditions,
and if we are integrating this equation numerically
211
00:22:58,110 --> 00:23:03,550
that is directly that is without introducing
transformation on the dependent variable that
212
00:23:03,550 --> 00:23:09,940
is known as direct integration, so here we
would be integrating a set of N coupled ordinary
213
00:23:09,940 --> 00:23:16,030
differential equation, that is direct integration
an alternative would be we will consider this
214
00:23:16,030 --> 00:23:23,230
equation but we will make this transformation
U = phi Z, where phi is the matrix of eigenvectors
215
00:23:23,230 --> 00:23:29,640
of the underlying stiffness and mass matrices,
and phi transpose M phi is identity, phi transpose
216
00:23:29,640 --> 00:23:34,490
C phi is diagonal, we assume damping to be
classical and phi transpose K phi is diagonal
217
00:23:34,490 --> 00:23:38,800
with elements being the squares of the natural
frequencies, with this normalization scheme
218
00:23:38,800 --> 00:23:45,690
the governing equations for Z will result
in a set of uncoupled second order ordinary
219
00:23:45,690 --> 00:23:52,320
differential equation this we have seen this.
So now P(t) is the generalized force obtained
220
00:23:52,320 --> 00:23:58,360
by using the relation phi transpose P(t) = Phi
transpose F, and similarly we can derive the
221
00:23:58,360 --> 00:24:06,620
initial conditions on Z using the orthogonality
relations and the mass matrix of the system.
222
00:24:06,620 --> 00:24:11,790
Now what we could do is we can integrate each
one of these equations numerically find out
223
00:24:11,790 --> 00:24:17,500
ZN(t) as time history and our interest would
be to find elements of U(t) so we had to go
224
00:24:17,500 --> 00:24:23,700
back to this equation and find elements of
U(t), so in a mode superposition method that
225
00:24:23,700 --> 00:24:29,610
the numerical solution of a differential equation
is carried out in the Z domain, and to get
226
00:24:29,610 --> 00:24:38,010
back to U domain you have to again go back
to U = phi Z, the relation U = phi Z.
227
00:24:38,010 --> 00:24:45,100
Now suppose we are doing modal superposition
and we consider this set of equations N = 1
228
00:24:45,100 --> 00:24:51,670
to N, in a N degree, capital N degree of freedom
system there will be N such oscillators and
229
00:24:51,670 --> 00:24:55,811
suppose if we integrate all these equations
with a common step size of delta T, delta
230
00:24:55,811 --> 00:25:03,050
T is common not only for each mode for at
all times, but also the same delta T is used
231
00:25:03,050 --> 00:25:11,200
to integrate Z1, Z2, Z3 and Z capital N, this
would be equivalent to in evaluating U(t)
232
00:25:11,200 --> 00:25:15,480
by integrating directly the governing equation
without doing mode superposition with the
233
00:25:15,480 --> 00:25:23,100
same step size delta T, so this is something
that is very important to take cognizance
234
00:25:23,100 --> 00:25:29,240
of. Suppose we are using forward difference
scheme which is conditionally stable to integrate
235
00:25:29,240 --> 00:25:36,190
this equation say MU double dot + CU dot KU
= F(t), in order that we get stable solutions
236
00:25:36,190 --> 00:25:42,650
we need to use a step size of delta T which
is minimum of 2 ETA and omega N, where these
237
00:25:42,650 --> 00:25:49,460
N runs from the first mode to the last mode.
Now suppose if ETA N is constant for all modes
238
00:25:49,460 --> 00:25:55,600
then the critical step size is given by 2
ETA omega capital N, that means the highest
239
00:25:55,600 --> 00:26:03,270
natural frequency plays a crucial role in
determining the critical step size if you
240
00:26:03,270 --> 00:26:07,580
are using the forward difference scheme, so
this places severe demand on the step size,
241
00:26:07,580 --> 00:26:14,460
because as you go higher up in the frequency
the time periods become small and our requirements
242
00:26:14,460 --> 00:26:20,460
on step size become more and more stringent.
On the other hand if you use backward difference
243
00:26:20,460 --> 00:26:24,010
scheme which is unconditionally stable to
integrate this equation
244
00:26:24,010 --> 00:26:29,210
the scheme would provide stable solution for
any step size, so there is no requirement
245
00:26:29,210 --> 00:26:36,260
on step size from point of view of stability,
so the choice of step size now, step size
246
00:26:36,260 --> 00:26:42,430
to be used in this case is governed by the
considerations on accuracy, so here the analyst
247
00:26:42,430 --> 00:26:50,260
has to exercise his or her judgment in selecting
delta T, so this calls for engineering judgment
248
00:26:50,260 --> 00:26:58,100
that mean we should have an idea on how to
select step size because any step size that
249
00:26:58,100 --> 00:27:05,390
we use will result in stable solutions. Now
from the numerical schemes and the illustrations
250
00:27:05,390 --> 00:27:13,390
that we have shown so far, we see that about
800 to 1600 steps within a cycle of oscillations
251
00:27:13,390 --> 00:27:20,550
at the highest frequency expected in the response
has to be used to get acceptable accuracy.
252
00:27:20,550 --> 00:27:26,920
Now when an conditionally stable explicit
scheme is used for analyzing multi degree
253
00:27:26,920 --> 00:27:31,990
freedom systems the requirement on step size
to ensure stable solution typically results
254
00:27:31,990 --> 00:27:42,110
in stringer requirements than requirements
on step size imposed by consideration of accuracy,
255
00:27:42,110 --> 00:27:47,720
that means if you ensure that the solution
scheme is stable in a multi degree freedom
256
00:27:47,720 --> 00:27:51,390
system that means you have to go to the highest
natural frequency and find out the critical
257
00:27:51,390 --> 00:27:58,150
step size and use the step size which is less
than the delta T critical then by enlarge
258
00:27:58,150 --> 00:28:06,000
the requirement on accuracy would be met,
because in the response of multi degree freedom
259
00:28:06,000 --> 00:28:14,380
systems if we use modal superposition often
we see that the majority of the contribution
260
00:28:14,380 --> 00:28:19,591
to the response would be made by first few
modes, so if you use now the highest mode
261
00:28:19,591 --> 00:28:25,380
in determining step size then you will have
adequate number of points within a cycle of
262
00:28:25,380 --> 00:28:32,090
lower order modes, so that is likely to lead
to acceptable solutions, so in conditionally
263
00:28:32,090 --> 00:28:40,460
stable schemes therefore accuracy is likely
to be guaranteed if you use step size which
264
00:28:40,460 --> 00:28:45,880
is less than critical step size, but it is
not a universally valid statement but one
265
00:28:45,880 --> 00:28:53,670
could typically expect this to happen.
Now in the methods that we discussed so far
266
00:28:53,670 --> 00:29:02,820
we have used first order approximations and
that led to you know for getting acceptable
267
00:29:02,820 --> 00:29:08,750
solutions from point of view of accuracy we
saw that we need very fine step size to be
268
00:29:08,750 --> 00:29:16,010
used in the scheme, now that requirement can
be mitigated if we use now, second order finite
269
00:29:16,010 --> 00:29:24,179
difference schemes to approximate the derivatives,
so we start discussing about or what is known
270
00:29:24,179 --> 00:29:28,170
as central difference method, so first let
us review quickly what is the central difference
271
00:29:28,170 --> 00:29:37,580
approximation to a derivative. Suppose I consider
U dot(t), I can write this in terms of value
272
00:29:37,580 --> 00:29:44,550
of displacement at T + delta T/2 and T - delta
T/2, so delta T is the time difference between
273
00:29:44,550 --> 00:29:50,160
these two time instants, and this will have
accuracy order delta T square. Similarly U
274
00:29:50,160 --> 00:29:59,340
double dot (t) I can write in terms of velocities
as this. Now for each of these U dot and U
275
00:29:59,340 --> 00:30:05,950
dot T + delta T/2 and T - delta T/2 I can
further use the finite difference approximation,
276
00:30:05,950 --> 00:30:16,910
so that means what we are doing here is we
have time instant T here, and we have T +
277
00:30:16,910 --> 00:30:27,700
delta T/2, and we have T - delta T/2, so at
this point we are writing U dot, so that will
278
00:30:27,700 --> 00:30:39,160
be the displacement here at T + delta T/2
- T minus delta T/2 by the difference delta
279
00:30:39,160 --> 00:30:45,210
T, similarly acceleration can be written,
but now when it comes to approximating velocities
280
00:30:45,210 --> 00:30:48,920
I need
velocities at this point as well as at this
281
00:30:48,920 --> 00:30:59,620
point, what I do? I consider other time instant
T + delta-T and one more time instant T - delta
282
00:30:59,620 --> 00:31:08,290
T, so they will start appearing in our approximation
because we need to get velocity here again
283
00:31:08,290 --> 00:31:13,620
central difference method requires the displacement
here minus displacement here, so that is what
284
00:31:13,620 --> 00:31:18,910
we are doing here, so this is the approximation
to acceleration, this is approximation to
285
00:31:18,910 --> 00:31:23,411
velocity.
Now if I now return to the governing equation
286
00:31:23,411 --> 00:31:30,020
of motion this is MU double dot + CU dot +
KU = F(t), now we will consider equilibrium
287
00:31:30,020 --> 00:31:34,750
at time T and we will use for velocity this
approximation, and for acceleration this approximation,
288
00:31:34,750 --> 00:31:41,750
this is what we have derived just now, substitute
for U dot and U double dot into this equation
289
00:31:41,750 --> 00:31:51,150
we get this equation. Now if I rearrange these
terms and collect terms which multiply U(t)
290
00:31:51,150 --> 00:31:58,000
+ delta T and so on and so forth I get this
equation. Now the equation for U(t) + delta
291
00:31:58,000 --> 00:32:08,470
T can be deduced from this and we get this
equation, therefore UN + 1 in this approximation
292
00:32:08,470 --> 00:32:18,180
will be this matrix inverse into this column
vector, so to implement this scheme what we
293
00:32:18,180 --> 00:32:24,140
should do we will see that. Now if we now,
there
294
00:32:24,140 --> 00:32:33,870
is one point that we should notice so this
is the expression for UN+1, now if N = 0,
295
00:32:33,870 --> 00:32:46,210
U1 will have if you see carefully U naught
and U(-1) okay, this will be U(-1). Now what
296
00:32:46,210 --> 00:32:53,290
is U(-1)? How do we get U(-1)? U(-1) is a
hypothetical quantity, so this also means
297
00:32:53,290 --> 00:33:00,580
that this scheme is not a self-starting scheme,
what we do is we manipulate these expressions
298
00:33:00,580 --> 00:33:07,350
and get in approximation to U(-1) as follows,
so we have U naught dot is U(1) - U(-1)/2
299
00:33:07,350 --> 00:33:14,170
delta T, from this I get U1 as this, similarly
U naught double dot will be this, from this
300
00:33:14,170 --> 00:33:21,820
I get this expression. So now by solving for
U(-1) I get this, so this is what I will be
301
00:33:21,820 --> 00:33:27,380
doing and we should notice that U double dot
of 0 is reducible from the governing equation,
302
00:33:27,380 --> 00:33:33,210
so this will be simply M inverse F(0) - CU
naught dot - KU naught, so U double dot of
303
00:33:33,210 --> 00:33:37,610
0 is known from the definition of the problem
therefore what the terms that lie on the right
304
00:33:37,610 --> 00:33:46,910
hand side of this equation are also known.
Now how do we implement this method? So we
305
00:33:46,910 --> 00:33:56,030
can write now the steps so we set N = 0 and
obtain U double dot of 0 from the governing
306
00:33:56,030 --> 00:34:01,210
equation because we know F and F(t) for all
T, then from this I reduce U(-1) which is
307
00:34:01,210 --> 00:34:08,809
this, then I evaluate these matrices I am
assuming that delta T is constant for all
308
00:34:08,809 --> 00:34:14,419
times that would mean the inversion of this
matrix need not be done within the time loop,
309
00:34:14,419 --> 00:34:19,250
within the loops that through which the time
advances so that can be done outside, so I
310
00:34:19,250 --> 00:34:25,350
evaluate all these matrices which I need subsequently
outside the loop. Now in this step I evaluate
311
00:34:25,350 --> 00:34:30,980
E1 +1, UN dot and UN double dot, for which
we already reduce the equation. Now we advance
312
00:34:30,980 --> 00:34:36,159
N and we will stop if it crosses the final
time instant, otherwise we'll go back and
313
00:34:36,159 --> 00:34:42,110
continue advancing the time, so this is the
scheme so again you should notice that the
314
00:34:42,110 --> 00:34:48,869
matrices A, B and D are evaluated only once
outside this looping, of course this is possible
315
00:34:48,869 --> 00:34:54,760
only when if we are dealing with linear systems,
if the system is nonlinear then at every time
316
00:34:54,760 --> 00:35:01,270
T we need to solve an algebraic non-linear
equation and we will deal with that sometime
317
00:35:01,270 --> 00:35:08,970
later.
Now before we consider the questions about
318
00:35:08,970 --> 00:35:15,130
stability of the central difference scheme,
will digress a bit and we discussed what is
319
00:35:15,130 --> 00:35:21,320
known as Jury's criterion. In our stability
analysis we have already saying if A is the
320
00:35:21,320 --> 00:35:25,480
amplification matrix the question that we
need to answer is whether the spectral radius
321
00:35:25,480 --> 00:35:32,070
of A is less than 1 or not, we don't really
need to know the numerical values of eigenvalues
322
00:35:32,070 --> 00:35:39,070
of A, we need simply to check an inequality,
so that means the question we are asking is
323
00:35:39,070 --> 00:35:43,490
very simple, far simpler than asking what
are the eigenvalues we are simply asking whether
324
00:35:43,490 --> 00:35:48,130
spectral radius is magnitude of spectral radius,
and the magnitude of the highest eigenvalue
325
00:35:48,130 --> 00:35:55,300
is less than 1 or not, to answer that question
there is a simpler strategy which avoids the
326
00:35:55,300 --> 00:36:02,000
actual determination of eigenvalues and that
strategy uses what is known as Jury's criterion.
327
00:36:02,000 --> 00:36:07,109
So to illustrate that we will consider the
polynomial P(z) as A naught Z to the power
328
00:36:07,109 --> 00:36:16,510
of N A1 ZN-1 and this, this is the polynomial
where A naught is greater than 0.
329
00:36:16,510 --> 00:36:23,220
Now Jury's criterion provides a necessary
and sufficient conditions, necessary and sufficient
330
00:36:23,220 --> 00:36:27,240
condition for the roots of the polynomial
to lie within the unit circle in the complex
331
00:36:27,240 --> 00:36:34,390
plane, so the necessary conditions are P(1)
must be greater than 0 and - 1 to the power
332
00:36:34,390 --> 00:36:40,980
of N, P(-1) should be greater than 0. Now
to construct sufficient conditions we construct
333
00:36:40,980 --> 00:36:45,150
what is known as Jury's
table, the table has this form so from the
334
00:36:45,150 --> 00:36:50,520
given polynomial we start constructing these
rows and the number of rows to be constructed
335
00:36:50,520 --> 00:36:55,470
is 2N - 3, suppose this is the second order
polynomial there will be only one row, so
336
00:36:55,470 --> 00:37:00,350
the rows are shown here, the first row clearly
provides the coefficients of Z to the power
337
00:37:00,350 --> 00:37:09,060
of 0, Z to the power of 1, etcetera. The second
row is also you know derived from that and
338
00:37:09,060 --> 00:37:17,400
subsequent rows have B's and C's etcetera
where the BCK up to Q this continues, so this
339
00:37:17,400 --> 00:37:25,960
is a sequence continues, is given by this
set of relations so from the given polynomial
340
00:37:25,960 --> 00:37:32,340
the first step is to construct this table,
then the necessary conditions we already seen
341
00:37:32,340 --> 00:37:38,829
P(1) greater than 0, P(-1) to be greater than
0, for N even otherwise less than 0 for N
342
00:37:38,829 --> 00:37:44,750
= R, the set of sufficient conditions is given
by terms contained in this table, and they
343
00:37:44,750 --> 00:37:53,850
are listed here. So the number of sufficient
conditions is actually N - 1. Now the criteria
344
00:37:53,850 --> 00:37:57,850
do not provide information of the value of
the roots that means we are not finding the
345
00:37:57,850 --> 00:38:04,720
eigenvalues, but we are only checking whether
the roots lie within unit circle.
346
00:38:04,720 --> 00:38:11,300
Now so let's try to now return to the question
of stability analysis of the central difference
347
00:38:11,300 --> 00:38:18,410
method, now let's consider this equation X
double dot + 2 ETA omega X dot + omega square
348
00:38:18,410 --> 00:38:25,400
= R(t) and if we apply now the central difference
method I get for XN double dot this relation
349
00:38:25,400 --> 00:38:33,450
and for XN dot this, so substituting that
this we have done we get this equation. For
350
00:38:33,450 --> 00:38:40,930
reasons that will become clear we'll combine
this equation with the identity XN = XN and
351
00:38:40,930 --> 00:38:50,630
we get this equation, so this is a vector
of system states at N and N+1 and this is
352
00:38:50,630 --> 00:38:56,720
the amplification matrix given by this, so
the stability of the scheme depends on eigenvalues
353
00:38:56,720 --> 00:39:02,450
of A, and eigenvalues of A intern depends
on the natural frequency, damping and delta
354
00:39:02,450 --> 00:39:08,750
T, natural frequency and damping are parts
of this problem specification so it is not,
355
00:39:08,750 --> 00:39:13,130
they are not the parameters that we can choose,
so the only parameter that we can choose is
356
00:39:13,130 --> 00:39:18,690
delta T and by examining eigenvalues of this
we will be able to tell what is the critical
357
00:39:18,690 --> 00:39:24,850
delta T, so we get by formulating A - lambda
I
358
00:39:24,850 --> 00:39:31,240
determinant of that equal to 0 this polynomial
equation, this is a characteristic equation.
359
00:39:31,240 --> 00:39:36,500
Now of course I can solve for this lambda
and find out what the eigenvalues are, but
360
00:39:36,500 --> 00:39:42,290
let us now use Jury's criteria and see what
we get. Now P(lambda) therefore is lambda
361
00:39:42,290 --> 00:39:47,330
square minus lambda this, this is the polynomial
associated with a characteristic equation,
362
00:39:47,330 --> 00:39:53,170
the necessary condition is P(1) must be greater
than 0, so that leads to the requirement omega
363
00:39:53,170 --> 00:39:58,810
square delta T squared / 1 + ETA omega delta
T must be greater than 0, numerator is positive,
364
00:39:58,810 --> 00:40:04,860
ETA is positive, omega is positive, therefore
this requirement is met. Similarly - 1 whole
365
00:40:04,860 --> 00:40:13,290
square P (-1) if we manipulate we get this
expression and this has to be greater than
366
00:40:13,290 --> 00:40:19,310
0, so for this to happen we get a condition
on delta T and that turns out to be delta
367
00:40:19,310 --> 00:40:25,859
T less than 2 / Omega which is this 2 by omega
is nothing but TN / phi, where TN is the period,
368
00:40:25,859 --> 00:40:32,280
2 phi / omega N.
Now the sufficient conditions how many are
369
00:40:32,280 --> 00:40:38,490
there? 2 - 1 which is 1, so that is modulus
of A naught must be less than modulus of A2,
370
00:40:38,490 --> 00:40:43,920
that means 1 should be less than 1 - ETA omega,
delta T / 1 + ETA omega delta, so this is
371
00:40:43,920 --> 00:40:49,900
always satisfied we can see this, examine
this and we can verify that it is always satisfied.
372
00:40:49,900 --> 00:41:03,610
So
now delta T critical is obtained as 2/omega,
373
00:41:03,610 --> 00:41:10,010
so the step size that we select must satisfy
the requirement that delta T should be less
374
00:41:10,010 --> 00:41:15,530
than 2/omega and we can immediately notice
that this critical step size which is 2/omega
375
00:41:15,530 --> 00:41:19,730
is independent of damping in the system, so
in the forward difference scheme if you recall
376
00:41:19,730 --> 00:41:24,900
we got 2 ETA/omega is the critical step size,
backward difference scheme of course was unconditionally
377
00:41:24,900 --> 00:41:33,770
stable but here we get the requirement that
delta T should be less than 2/omega, so therefore
378
00:41:33,770 --> 00:41:38,990
we can observe few properties of this method,
the method is an explicit method with second-order
379
00:41:38,990 --> 00:41:45,840
accuracy, the method is conditionally stable
with critical step size of 2/omega, the critical
380
00:41:45,840 --> 00:41:52,010
size step size is independent of damping.
Now the method requires a special starting
381
00:41:52,010 --> 00:41:56,841
procedure, it is not a self-starting scheme
and we have discussed, we have assumed in
382
00:41:56,841 --> 00:42:03,240
our discussion that the step size is constant.
Now if M is diagonal and C is proportional
383
00:42:03,240 --> 00:42:14,350
to M, then this matrix A, if M is diagonal
and C is proportional to M this matrix will
384
00:42:14,350 --> 00:42:20,870
be a diagonal matrix therefore the implementation
of the method does not require inversion of
385
00:42:20,870 --> 00:42:26,370
any matrix, so that is when you are dealing
with large system that is something that we
386
00:42:26,370 --> 00:42:34,970
should take note of, and it is also important
to note that this requirement may, that M
387
00:42:34,970 --> 00:42:39,280
is diagonal and C is proportional the observation
that the method does not require inversion
388
00:42:39,280 --> 00:42:45,070
of the matrix is not particularly significant
if you are dealing with linear systems, if
389
00:42:45,070 --> 00:42:49,790
we are dealing with either multi-step methods
or non-linear systems this becomes an important
390
00:42:49,790 --> 00:42:56,930
observation. If this condition is not satisfied
then we need to invert this at least once
391
00:42:56,930 --> 00:43:02,370
if you are using same delta T, otherwise you
have to invert this at every delta T that
392
00:43:02,370 --> 00:43:06,290
is if delta T is not constant then the matrix
needs to be inverted at every step, so this
393
00:43:06,290 --> 00:43:12,700
can place severe demands on computation if
you are dealing with large-scale problems.
394
00:43:12,700 --> 00:43:17,190
The method is easily implementable even if
the system is nonlinear, because it's an explicit
395
00:43:17,190 --> 00:43:24,040
method there is no solution of nonlinear algebraic
equations as time progresses. Now delta T
396
00:43:24,040 --> 00:43:31,230
is the only algorithmic parameter to be selected
and that has to be less than 2/omega.
397
00:43:31,230 --> 00:43:37,380
So some numerical illustrations, so let's
start with an undamped free vibration of a
398
00:43:37,380 --> 00:43:43,860
system with natural frequency 2 phi, ETA is
0 and will select a step size which is about
399
00:43:43,860 --> 00:43:48,850
5% more than the critical step size, so we
see that the red line is the finite difference
400
00:43:48,850 --> 00:43:53,690
approximation so the error is growing and
it swamps the true solution, so on this scale
401
00:43:53,690 --> 00:43:59,210
the true solution appears like a flat line.
Now on the other hand if delta T is about
402
00:43:59,210 --> 00:44:05,250
95% of the delta T critical, so things look
the at least the errors don't grow but the
403
00:44:05,250 --> 00:44:11,230
solution obtained is quite unsatisfactory,
okay, the blue is the exact solution, red
404
00:44:11,230 --> 00:44:21,420
is the solution from central difference scheme.
Now if I now take step size to be delta T
405
00:44:21,420 --> 00:44:26,030
critical/5, the solution starts looking much
better, but still as you see, as you progress
406
00:44:26,030 --> 00:44:32,400
in time there is still a distortion in phase
and there is some observable difference between
407
00:44:32,400 --> 00:44:39,859
blue and the red lines, if you take now delta
T, delta T to be delta T critical / 10 then
408
00:44:39,859 --> 00:44:47,170
the solution looks quite good, the red and
blue look almost identical, of course one
409
00:44:47,170 --> 00:44:52,520
can plot the error, norms and examine what
exactly are these differences, but through
410
00:44:52,520 --> 00:44:59,590
these graphs we can see at least qualitatively
that the solutions are agreeing quite well.
411
00:44:59,590 --> 00:45:06,420
Now let's consider damped system, this critical
step size is independent of damping and that
412
00:45:06,420 --> 00:45:11,119
remains constant, suppose if I take now a
delta T to be one point naught 5 times delta
413
00:45:11,119 --> 00:45:15,250
T critical which is not a right choice to
make the solutions will become unstable and
414
00:45:15,250 --> 00:45:20,590
that is what we see here.
Now on the other hand if I take 95% of delta
415
00:45:20,590 --> 00:45:27,530
T critical, the errors there has no growth
of errors but still the solution is, numerical
416
00:45:27,530 --> 00:45:34,580
solution is very poor, so with factor 5 that
means delta T is delta T critical by 5, we
417
00:45:34,580 --> 00:45:42,140
start getting good solutions, and with delta
T critical by 10 I get quite pretty good solution.
418
00:45:42,140 --> 00:45:55,060
Now
we want to now study the stability of other
419
00:45:55,060 --> 00:46:00,190
integration schemes that I will be introducing
shortly, central difference method is an explicit
420
00:46:00,190 --> 00:46:06,570
second-order method it is conditionally stable,
so the next question that we need to ask is
421
00:46:06,570 --> 00:46:14,740
are there any schemes which are implicit and
unconditionally stable and have second-order
422
00:46:14,740 --> 00:46:20,300
accuracy, so that discussion will take up
in the next class, as a prelude to that we
423
00:46:20,300 --> 00:46:26,220
need to have some idea on what are known as
Z transforms, that we'll be using in analyzing
424
00:46:26,220 --> 00:46:32,790
those other schemes and we can also as well
analyze what we have already done using this
425
00:46:32,790 --> 00:46:40,080
approach, so let's quickly recall some preliminaries
about that, so let F(t) be a function of time
426
00:46:40,080 --> 00:46:46,869
such that F(t) is 0 for T less than 0, and
let S be a complex number, we define the Laplace
427
00:46:46,869 --> 00:46:53,650
transform of F(t) as F(s) 0 to infinity, F(t)
E raise to - ST DT.
428
00:46:53,650 --> 00:47:00,200
Now if F(t) is a sampled, that means it is
multiplied by a comb like function, a series
429
00:47:00,200 --> 00:47:10,010
of direct delta functions, so how does this
function look like? It is a series of direct
430
00:47:10,010 --> 00:47:22,119
delta functions place that a common interval
a capital T, now if F(t) is written like this
431
00:47:22,119 --> 00:47:26,980
that means, F(t) is sampled at all these time
instants, the sampled waveform can be written
432
00:47:26,980 --> 00:47:31,290
in this form. Now let us consider Laplace
transform of that so I substitute this into
433
00:47:31,290 --> 00:47:36,730
this so upon carrying out the integration
I get Laplace transform F(s) as N = 0, F(nT)
434
00:47:36,730 --> 00:47:46,170
E raise to - SN capital T. Now this can be
written as F(nT) exponential ST to the power
435
00:47:46,170 --> 00:47:53,670
of - N, suppose if I call this quantity E
raise to ST as Z, I can write this as F(nT),
436
00:47:53,670 --> 00:48:00,950
Z to the power of -N, so this quantity is
known as Z transform of F(t) and it is written
437
00:48:00,950 --> 00:48:10,830
as F(z) as this, so this is a discrete version
of Laplace transform, so this is the definition.
438
00:48:10,830 --> 00:48:18,480
Now we need to recall some result from geometric
progression a series like 1 + R + R square
439
00:48:18,480 --> 00:48:26,820
up to RN is called a geometric progression,
and if I multiply now this by R, I get R +
440
00:48:26,820 --> 00:48:33,420
R square, so on and so forth +R to the power
of N+1, suppose if I subtract this I get 1
441
00:48:33,420 --> 00:48:38,490
- R into SN which is
equal to this, from which I get sum of first
442
00:48:38,490 --> 00:48:44,849
N terms in a geometric progression to be given
by this, so if modulus of R is less than 1
443
00:48:44,849 --> 00:48:51,640
then as N tends to infinity SN goes to S which
is 1/1-R, now we can quickly see few examples
444
00:48:51,640 --> 00:48:58,359
if F(t) is a step function, F(nT) will be
1 for all N, otherwise it is 0 now using that
445
00:48:58,359 --> 00:49:03,570
and using this result we can show that the
Laplace, the Z transform of a unit step function
446
00:49:03,570 --> 00:49:11,430
is given by Z /Z - 1.
Similarly if I have F(t) is KT that means
447
00:49:11,430 --> 00:49:16,770
a straight line, then F(nT) will be n capital
T and the Laplace transform of that we can
448
00:49:16,770 --> 00:49:25,970
show through some manipulation that it is
given by TZ/Z-1 whole square, this you can
449
00:49:25,970 --> 00:49:31,320
verify. There are few more examples, for example
we having
450
00:49:31,320 --> 00:49:36,510
exponential, if you are having a F(t) as E
raise to - AT, then F(nT) will be exponential
451
00:49:36,510 --> 00:49:45,710
- anT. Now the Laplace transform of this if
I do I will get this as Z/Z - exponential
452
00:49:45,710 --> 00:49:54,770
- AT, similarly we can verify that the Z transform
of sin omega T and cos omega T are obtained,
453
00:49:54,770 --> 00:50:01,369
through these expressions, so I'll leave this
as an exercise for you to verify.
454
00:50:01,369 --> 00:50:09,700
Now a few properties of Z transform can be
quickly recalled, the Z transform of AF(nT)
455
00:50:09,700 --> 00:50:20,200
is AF(z), similarly the Z transform of A to
the power of N, F(nT) is FA - 1 Z, Z transform
456
00:50:20,200 --> 00:50:32,310
F(t) -NT is Z - N F(z), Z transform of F(t)
+ NT is given by this, okay, now if we now,
457
00:50:32,310 --> 00:50:40,270
I mean this can be used for K+1, K+2, K+N,
so on and so forth, suppose if I now consider
458
00:50:40,270 --> 00:50:50,650
Y(k) to be summation N = 0 to KX(n), if we
have this Y(k) I can write Y(K) as Y(K-1)
459
00:50:50,650 --> 00:50:59,510
+ X(K). Now if I take Z transform, Z transform
of Y(k) - Y (k-1) is Z transform of X(k) which
460
00:50:59,510 --> 00:51:08,380
is X(z), so from this I get Y(z) to be X(z)/
1 - Z, so if there is a dynamical system in
461
00:51:08,380 --> 00:51:18,020
which evolves as YK = YK-1 + XK, then the
X(z)/1-Z would be its transfer function.
462
00:51:18,020 --> 00:51:24,771
Now to make that clear we'll consider a dynamical
system with input X(t) and output Y(t), let
463
00:51:24,771 --> 00:51:32,600
the input-output relation for the sample data
be given by this equation, let be Y(n) + A1
464
00:51:32,600 --> 00:51:39,220
Y(n-1) + so few terms up to AK YN-K this is
equal to B naught X(n) + B1 XN - 1 and so
465
00:51:39,220 --> 00:51:45,110
on and so forth, up to N - M. Now if I take
now Z transforms on both sides and use identities
466
00:51:45,110 --> 00:51:51,320
that I briefly mentioned we can show that
the frequency response function in the Z domain
467
00:51:51,320 --> 00:51:57,130
is given by this ratio and we get this as
the frequency response function. Now the idea
468
00:51:57,130 --> 00:52:03,790
of introducing this Z transform is basically
to visualize the finite difference schemes
469
00:52:03,790 --> 00:52:10,630
that we are going to develop as a kind of
input-output relation of this form, so consequently
470
00:52:10,630 --> 00:52:15,520
associated with the finite difference scheme
that we use we can define a frequency response
471
00:52:15,520 --> 00:52:19,250
function.
Now this frequency response function will
472
00:52:19,250 --> 00:52:29,330
be in the form of a ratio like this and the
zeroes of the denominator zeroes, and the
473
00:52:29,330 --> 00:52:34,440
zeroes of the numerator are known as poles.
So the zeroes of the frequency response function
474
00:52:34,440 --> 00:52:40,170
must lie within unit circle is a requirement
that we need to verify, so here if you are
475
00:52:40,170 --> 00:52:44,330
following this approach we need not formulate
the A matrix, the amplification matrix instead
476
00:52:44,330 --> 00:52:50,380
we can formulate the frequency response function
in the Z domain and look at the zeros and
477
00:52:50,380 --> 00:52:57,390
poles of the frequency response function.
Now to quickly see that let's say we revisit
478
00:52:57,390 --> 00:53:03,960
the central difference scheme and we got this
as evolution equation, which is this, now
479
00:53:03,960 --> 00:53:09,599
if I write this for N I get this equation,
on the right hand side there will be N - 1,
480
00:53:09,599 --> 00:53:20,840
now if I take Z transform I get this equation,
so for non-trivial solutions you know this
481
00:53:20,840 --> 00:53:26,240
bracket, terms inside this bracket must be
equal to 0 and this leads to the polynomial
482
00:53:26,240 --> 00:53:32,810
equation that we are looking for, and this
is identical to the characteristic equation
483
00:53:32,810 --> 00:53:37,770
we obtain by considering eigenvalues of the
amplification matrix.
484
00:53:37,770 --> 00:53:42,480
So now what we can do is instead of applying
Jury's criteria on this, criterion on this
485
00:53:42,480 --> 00:53:47,619
we can apply the Jury's criterion on this
equation, so in this case will get the same
486
00:53:47,619 --> 00:53:51,590
solution as we got before because we apply
Jury's criterion on this equation to derive
487
00:53:51,590 --> 00:53:57,230
the critical step size, and if we had done
this using this logic would have got the same
488
00:53:57,230 --> 00:54:02,460
answer.
So now let's conclude this lecture at this
489
00:54:02,460 --> 00:54:08,760
stage, in the next class I will ask the question,
are there any second order schemes which are
490
00:54:08,760 --> 00:54:14,850
implicit and which are unconditionally stable,
so that takes us to discussion of what are
491
00:54:14,850 --> 00:54:20,069
known as Newmark family of methods and we'll
take up that discussion in the next lecture.