1
00:00:18,050 --> 00:00:25,289
In the previous class, I have introduced you
to numerical methods of solving differential
2
00:00:25,289 --> 00:00:32,720
equations. In this particular lecture, we
will continue that strength and we will do
3
00:00:32,720 --> 00:00:38,820
a few examples on psi lab and to illustrate
the basic properties of most of the commonly
4
00:00:38,820 --> 00:00:45,239
used numerical methods. The examples which
I will chose today are examples which we have
5
00:00:45,239 --> 00:00:52,329
done before we have tried to analyze them
using some intuition as well as linear analysis
6
00:00:52,329 --> 00:00:59,000
in some cases.
Today, let us just first revise what we have
7
00:00:59,000 --> 00:01:04,539
done in the previous class. In the previous
class we learnt a few numerical methods like
8
00:01:04,539 --> 00:01:05,539
Euler method.
9
00:01:05,539 --> 00:01:12,770
So, let us just quickly go through this and
then begin with our examples. Euler method
10
00:01:12,770 --> 00:01:18,610
or is also known as forward Euler method is
the first is a first method which is taught
11
00:01:18,610 --> 00:01:25,640
to any student studying numerical methods
in for integration of differential equations.
12
00:01:25,640 --> 00:01:29,860
It is a very simple method. In fact, if you
have got a differential equation x dot is
13
00:01:29,860 --> 00:01:37,691
equal to x comma t, this comma looks a bit
big, so I will just write it down again x
14
00:01:37,691 --> 00:01:44,410
of t.
In that case, we evaluate the function using
15
00:01:44,410 --> 00:01:54,500
this algorithm at discrete time steps spaced
by h. So, we have got x k 1 is equal to x
16
00:01:54,500 --> 00:02:07,020
k this is x x k plus h into f into x k and
that is the function f evaluated at the point
17
00:02:07,020 --> 00:02:17,590
x k and t k. Forward Euler method or simply
Euler method is having order 1 that is, if
18
00:02:17,590 --> 00:02:24,379
your numerical solution can be expressed simply
as a first order polynomial in t. Then, Euler
19
00:02:24,379 --> 00:02:30,950
method will you exact answers. Typically of
course, you know that our responses are not
20
00:02:30,950 --> 00:02:36,280
of first order are not first order polynomials
in t, we have in fact come across sinusoids
21
00:02:36,280 --> 00:02:40,970
and exponential functions which are not first
order polynomial. So, Euler method is likely
22
00:02:40,970 --> 00:02:48,420
to have some error whenever you integrate
most common differential equations.
23
00:02:48,420 --> 00:02:54,190
It is also an explicit method in the sense
to evaluate x k plus 1; you really do not
24
00:02:54,190 --> 00:02:57,730
have to solve any differential equation. You
have to just write the, if you know the previous
25
00:02:57,730 --> 00:03:05,850
value of x k. In one shot you can get x k
plus 1. So, you get essentially in explicit
26
00:03:05,850 --> 00:03:13,870
algebraic equation where in x k plus 1 is
written explicitly in terms of x k. Forward
27
00:03:13,870 --> 00:03:21,980
Euler is also a single step method. So, you
just require one previous, the value at the
28
00:03:21,980 --> 00:03:28,700
previous time step in order to evaluate the
value at the next time step.
29
00:03:28,700 --> 00:03:36,540
If you look at a special class of systems
that are linear systems then, one thing we
30
00:03:36,540 --> 00:03:41,510
can attempt to do since we know the solution
of linear equations, we can compare what the
31
00:03:41,510 --> 00:03:47,099
numerical method gives us and what is the
actual solution of a linear system. Now, we
32
00:03:47,099 --> 00:03:52,650
cannot do this in general for non-linear system
because, you cannot write down the response
33
00:03:52,650 --> 00:03:58,680
of a non-linear system in terms of simple
functions. In contrast, the linear systems
34
00:03:58,680 --> 00:04:03,950
for example, x dot is equal to x we know,
it is stable when a is less than 0. If a is
35
00:04:03,950 --> 00:04:10,350
a real number then, if a is less than 0 then
it is a stable.
36
00:04:10,350 --> 00:04:15,990
Euler method, if you apply Euler method, if
you apply this algorithm where f of x is nothing
37
00:04:15,990 --> 00:04:21,989
but, a of A x. In that case, Euler method
tells you that the system is stable if mode
38
00:04:21,989 --> 00:04:29,940
of the absolute value of 1 plus A h is less
than 1. So obviously, the Euler method may
39
00:04:29,940 --> 00:04:35,639
give these two conditions in fact are not
equivalent. You know, if you look at, you
40
00:04:35,639 --> 00:04:47,539
know, for example, if h is greater than 2
by a, in that case even though a is less than
41
00:04:47,539 --> 00:04:51,809
0 you may still have Euler method showing
the system to be unstable.
42
00:04:51,809 --> 00:04:57,139
So, Euler method when you integrate the system
may be give you wrong answer. Even, not only
43
00:04:57,139 --> 00:05:02,059
wrong answers, it will give you a qualitatively
wrong result. It may show a stable to system
44
00:05:02,059 --> 00:05:09,690
to be a unstable one. So, this is basically
1 one problem with Euler method that, you
45
00:05:09,690 --> 00:05:13,389
know, you will not get a accurate answer.
Then, you may even get qualitatively wrong
46
00:05:13,389 --> 00:05:14,560
answers.
47
00:05:14,560 --> 00:05:20,610
Look at backward Euler method, backward Euler
method; the only difference is that the function
48
00:05:20,610 --> 00:05:29,990
evaluation is done at the new value of x.
It is in fact done at the point which you
49
00:05:29,990 --> 00:05:34,240
want to evaluate.
So although, backward Euler is also order
50
00:05:34,240 --> 00:05:40,300
1 this is something you can prove. I am not
proved this in this class. It is implicit,
51
00:05:40,300 --> 00:05:46,300
it is an implicit method in the sense to evaluate
x k plus 1. You have to solve this equation
52
00:05:46,300 --> 00:05:52,270
x k plus 1 appears in within this function.
Therefore, you need a numerical method at
53
00:05:52,270 --> 00:06:01,629
every time step in order to get x a plus 1.
So, this is; if for, if x is f is a non-linear
54
00:06:01,629 --> 00:06:05,919
function, you will have to use some numerical
method to solve algebraic equations like Newton
55
00:06:05,919 --> 00:06:10,889
Raphson or Godse did.
Backward Euler also is a single step method.
56
00:06:10,889 --> 00:06:18,659
For a linear system, Backward Euler is stable
if, this condition is satisfied, that is,
57
00:06:18,659 --> 00:06:27,080
1 upon the absolute value of 1 minus A h less
than 1. This condition also is not equivalent
58
00:06:27,080 --> 00:06:31,280
to this condition. In fact, an interesting
point which I will show you sometime later
59
00:06:31,280 --> 00:06:38,069
1 one of those examples is that, you can chose
values of h so that, an unstable system is
60
00:06:38,069 --> 00:06:42,409
shown as a stable system if you use backward
Euler method. So, backward Euler method also
61
00:06:42,409 --> 00:06:49,949
can in sometimes, if depending on the choice
of the value h, give you qualitatively wrong
62
00:06:49,949 --> 00:06:58,669
answers as far as stability is concerned of
a linear system of this kind.
63
00:06:58,669 --> 00:07:05,330
Trapezoidal method on the other hand is an
order 2 two method. The basic algorithm is
64
00:07:05,330 --> 00:07:10,860
this; it takes the average of the function
evaluation here, at the previous point and
65
00:07:10,860 --> 00:07:16,389
the new point. It is an implicit method again
and it is a single step method. It requires
66
00:07:16,389 --> 00:07:25,660
only the previous value of x in order to get
the new value of x. For linear systems, trapezoidal
67
00:07:25,660 --> 00:07:31,020
rule give says that this system is stable
if this is less than 1. Whereas, the actual
68
00:07:31,020 --> 00:07:37,430
system is stable when a less than 0 but, interestingly
this is a very interesting fact which I leave
69
00:07:37,430 --> 00:07:45,279
to you to work out; these 2 two conditions
are equivalent in the sense that, if a less
70
00:07:45,279 --> 00:07:51,719
than 0 it also means that absolute value of
1 plus A h by 2 divided by absolute value
71
00:07:51,719 --> 00:07:56,490
of 1 minus A h by 2 is less than 1.
So that is an interesting thing. Trapezoidal
72
00:07:56,490 --> 00:08:01,750
rule, whenever your simulating linear systems
will never give you a wrong result as far
73
00:08:01,750 --> 00:08:09,909
as stability is concerned. So, that is something
interesting and useful. But, of course, one
74
00:08:09,909 --> 00:08:15,770
important point is, implicit methods like
backward Euler and forward Euler are difficult
75
00:08:15,770 --> 00:08:20,439
to solve because, especially when you come
to non-linear systems, because you will need
76
00:08:20,439 --> 00:08:28,589
to use a numerical method to compute x k plus
1 from x k. It is a iterative method. At every
77
00:08:28,589 --> 00:08:34,220
time step itself, you will have to iterate
and solve a numerical solve a non-linear algebraic
78
00:08:34,220 --> 00:08:40,940
equation. When you are talking of linear systems
of course, it would mean typically a matrix
79
00:08:40,940 --> 00:08:46,140
inversion or a solution of a linear system
of equations.
80
00:08:46,140 --> 00:08:53,470
We also discussed in the previous class; not
trapezoidal, this method which I am talking
81
00:08:53,470 --> 00:09:00,279
of is Runge Kutta method. Runge Kutta method
of fourth order has the following evaluations;
82
00:09:00,279 --> 00:09:05,330
you have to evaluate k 1 at these points.
I have done this in the previous lecture,
83
00:09:05,330 --> 00:09:13,639
evaluate k 1 first, that is the function evaluation;
k 2 is a function evaluation at a slightly
84
00:09:13,639 --> 00:09:19,800
different point depending on which depends
on the first evaluation. So, f of x k plus
85
00:09:19,800 --> 00:09:28,550
h by 2 into k 1 k 1 is evaluated here. Similarly,
you have got the evaluation of k 3 which requires
86
00:09:28,550 --> 00:09:38,509
the value of k 2 and k 4 which requires the
value of k 3. So, what you see here, is sorry
87
00:09:38,509 --> 00:09:42,000
this is small error here.
88
00:09:42,000 --> 00:09:49,000
What you see here is; the final value of x
k is a combination of these. Of course, Runge
89
00:09:49,000 --> 00:09:53,690
Kutta requires a lot many more evaluations
and can be proven that this is the fourth
90
00:09:53,690 --> 00:09:58,889
order method. The 1 one I have shown you is
a fourth order method. It is of course, single
91
00:09:58,889 --> 00:10:08,360
step. The reason is that we require only the
value of x k in order to compute all these
92
00:10:08,360 --> 00:10:15,519
functions k 1 k 2 k 3 k 4 and it is an implicit
method. At least this one which I have shown
93
00:10:15,519 --> 00:10:21,120
you is an implicit sorry it is an explicit
method. I am sorry Runge Kutta fourth order
94
00:10:21,120 --> 00:10:25,959
which I have shown you here, is an explicit
method. Explicit in the sense, that all these
95
00:10:25,959 --> 00:10:30,551
functions evaluations are single shot. They
are explicit evaluations of functions. You
96
00:10:30,551 --> 00:10:39,509
do not have to solve any non-linear algebraic
equations to evaluate any any value here.
97
00:10:39,509 --> 00:10:44,680
Of course, when you are trying to numerically
integrate linear higher order systems, again
98
00:10:44,680 --> 00:10:52,769
let me point out here, that the reason why
I am studying the behavior of numerical methods
99
00:10:52,769 --> 00:10:58,389
using linear systems this is not because linear
systems cannot be sort solved. We saw in the
100
00:10:58,389 --> 00:11:03,821
lectures, previous lectures that, linear systems
in fact are amenable to exact solutions. You
101
00:11:03,821 --> 00:11:07,630
can write down the solution in terms of very
well known functions. So, the only reason
102
00:11:07,630 --> 00:11:12,889
why I am using linear systems a numerical
integration of linear systems, the study of
103
00:11:12,889 --> 00:11:18,490
that is because I would like to know how the
numerical methods behave.
104
00:11:18,490 --> 00:11:25,339
Now, if you have got a linear higher order
system like, x dot is equal to x a is a matrix
105
00:11:25,339 --> 00:11:33,870
x is vector. It is stable you know, it is
stable if real of the Eigen, any of the Eigen,
106
00:11:33,870 --> 00:11:40,640
all Eigen values should have the real part
less than 0 then we say its stable. But, if
107
00:11:40,640 --> 00:11:47,820
I numerically integrate this higher order
linear system; in that case what you have
108
00:11:47,820 --> 00:11:55,139
is, x k plus 1, this is a vector is 1 plus
A h times x k, which is again a vector. This
109
00:11:55,139 --> 00:11:59,800
is a matrix. Now, this system is stable, this
is a discrete time system. It is stable if,
110
00:11:59,800 --> 00:12:06,889
any the Eigen values of 1 plus A h is less
than 1. This not something I am explicitly
111
00:12:06,889 --> 00:12:11,990
proving but, it is not difficult to prove.
You can use the model analysis or linear systems
112
00:12:11,990 --> 00:12:17,160
analysis, Eigen values and Eigen vectors to
analyze discrete time systems as well.
113
00:12:17,160 --> 00:12:25,790
Now, if the absolute value of the Eigen value
of 1 plus A h is less than 1, then it is stable
114
00:12:25,790 --> 00:12:30,300
system. That is what Euler method says. The
original system is stable if all Eigen values
115
00:12:30,300 --> 00:12:36,930
are the real part less than 0. So, for a higher
order linear system, Euler method will give
116
00:12:36,930 --> 00:12:45,820
will give a result that, if 1 plus lambda
h the modulus of lambda h is less than 1,
117
00:12:45,820 --> 00:12:52,730
then it is stable. So, lambda is the Eigen
value of a. So, let me just repeat; if the
118
00:12:52,730 --> 00:12:58,980
Eigen values of 1 plus A h have the modulus
or absolute value less than 1 then, Euler
119
00:12:58,980 --> 00:13:04,860
method will show that the system is stable.
It may be actually a stable system and Euler
120
00:13:04,860 --> 00:13:09,480
method may show it to me a unstable system
if, this becomes greater than 1. It really
121
00:13:09,480 --> 00:13:15,569
depends on the value of h here.
You can show that mu is nothing but, 1 plus
122
00:13:15,569 --> 00:13:20,649
lambda h well where lambda are the Eigen value
of a. So, the condition, this condition can
123
00:13:20,649 --> 00:13:27,259
be rewritten as the condition for this discrete
time system to be stable is, 1 plus lambda
124
00:13:27,259 --> 00:13:34,089
h, the absolute value should be less than
1. So, again let me just clarify what I am
125
00:13:34,089 --> 00:13:38,440
trying to say; I am trying to say that, even
if the original system is stable, if this
126
00:13:38,440 --> 00:13:44,980
condition is not satisfied then, your numerical
method will give you a qualitatively wrong
127
00:13:44,980 --> 00:13:49,790
answer about stability. So, you have to be
careful when you do this simulations that
128
00:13:49,790 --> 00:13:58,310
you do not get wrong answers.
Now, moving on we will do a few examples in
129
00:13:58,310 --> 00:14:07,850
psi lab. It would be wonderful if you have
got psi lab or mat lab or even if you program
130
00:14:07,850 --> 00:14:13,149
and see your program FORTRAN and verify what
I am trying to say here.
131
00:14:13,149 --> 00:14:18,060
What I will do is, I will kind of make an
algorithm for doing the numerical integration
132
00:14:18,060 --> 00:14:28,600
first of a system x dot is equal to minus
5 x. This is a stable system. So, I will first
133
00:14:28,600 --> 00:14:32,899
numerically integrate it using Euler method.
Actually, if Euler method is used remember
134
00:14:32,899 --> 00:14:41,020
that it will you know it is likely to give
correct results. 1 One of the intuitive is
135
00:14:41,020 --> 00:14:46,620
of choosing your h is that, look at this Eigen
value this Eigen value will probably will
136
00:14:46,620 --> 00:14:56,790
give in fact a solution of this kind.
We will just to try to verify whether we get
137
00:14:56,790 --> 00:15:01,540
that solution or not. In fact you have to
chose choose your h carefully. If you look
138
00:15:01,540 --> 00:15:09,139
at, this particular you know, function it
will decay almost to 0 in 5 times the time
139
00:15:09,139 --> 00:15:13,820
constant of this system. The time constant
of this system is in this particular case
140
00:15:13,820 --> 00:15:20,769
1 up on 5 that is, point tow two. It is a
time constant of this system. So, we expect
141
00:15:20,769 --> 00:15:26,649
that, it will decay. This particular thing
will decay if I give a non equilibrium initial
142
00:15:26,649 --> 00:15:35,769
condition. You will go to 0 which is the equilibrium
in this case in about 1 about 1 one second.
143
00:15:35,769 --> 00:15:41,830
So, I can just intuitively think that, if
I want to get an nice, at least a nice picture
144
00:15:41,830 --> 00:15:47,829
of what is going to happen; I should chose
h at least you know, something like 0.0 1
145
00:15:47,829 --> 00:15:56,519
or so. We will just try out with various values
of h and see what we get. I had done this
146
00:15:56,519 --> 00:16:00,059
in the previous class also. But, it will was
somewhat in a hurried fashion. So, let us
147
00:16:00,059 --> 00:16:04,529
do it bit calmly in this particular lecture.
148
00:16:04,529 --> 00:16:15,009
So this is the basic psi lab window. I have
already invoked psi lab. What we will do is,
149
00:16:15,009 --> 00:16:22,449
first of all, let us look at the solution
of, so we will do a numerical integration
150
00:16:22,449 --> 00:16:31,170
using Euler method. First of all, we will
clear this variable x. We will chose choose
151
00:16:31,170 --> 00:16:38,910
a time constant let us say of a time step
of 0.01. Let us simulate or numerically integrate
152
00:16:38,910 --> 00:16:43,779
this system till 3 seconds because we know
that it is for lightly to decay within that
153
00:16:43,779 --> 00:16:49,660
time. So, we have got a number of discrete
time steps at which x has to be evaluated.
154
00:16:49,660 --> 00:16:55,699
So, let us assume the initial condition on
x is 1. In fact, I should have written this
155
00:16:55,699 --> 00:17:02,170
x of 0. But, psi lab does not allow 0 as an
index. So, that is why I have called it 1.
156
00:17:02,170 --> 00:17:07,819
So, this is the value in fact at time t is
equal to 0. So, I have given initial condition
157
00:17:07,819 --> 00:17:15,350
of 1. Remember x dot is equal to minus 5 x
has got in equilibrium value of 0. So, we
158
00:17:15,350 --> 00:17:19,070
are giving an initial condition which is not
the equilibrium value. So, you will see a
159
00:17:19,070 --> 00:17:28,440
transient, the value of a is 5 as I discussed
some time back. So, what I will do is, I will
160
00:17:28,440 --> 00:17:36,560
evaluate the behavior of this system using
Euler method. So, x of I is equal to 1 plus
161
00:17:36,560 --> 00:17:42,990
a into h into x of i minus 1. This is what
Euler method will do.
162
00:17:42,990 --> 00:17:56,330
So if you recall, we have got x k plus 1 minus
x k up on h is equal to minus 5 of x k. So,
163
00:17:56,330 --> 00:18:07,410
this will effectively give us x k plus 1 is
equal to 1 minus 5 h into x k. In fact, minus
164
00:18:07,410 --> 00:18:13,430
give five h is nothing but, plus A h. This
is nothing but, this is 1 plus A h, a itself
165
00:18:13,430 --> 00:18:16,690
is negative. That is why we are getting 1
minus 5 h.
166
00:18:16,690 --> 00:18:35,270
So, what I will do is, I will evaluate this
and plot it. So, the way to do it is, so we
167
00:18:35,270 --> 00:18:43,060
evaluate this. Hoops Oops! It is giving a
growing function. Now, why is this happening?
168
00:18:43,060 --> 00:18:49,060
We will just look back. Perhaps, I have done,
I have made some error here.
169
00:18:49,060 --> 00:18:57,550
H is equal to 0.01 that is oh okay. I am have
not saved it. So, I need to save it first.
170
00:18:57,550 --> 00:19:00,960
yeah I have saved it first. yeah I have saved
it. So, I have given the value of h is equal
171
00:19:00,960 --> 00:19:07,160
to 0.01 which I have saved changed from the
earlier value and now I will run it again.
172
00:19:07,160 --> 00:19:15,170
So let me just run this again. s So, I am
doing this of course, simulations of 3 seconds.
173
00:19:15,170 --> 00:19:21,950
Yeah Now, what we are getting is what we kind
of expect or exponentially decaying e raise
174
00:19:21,950 --> 00:19:29,470
to minus 5 t kind of response. The initial
value is one, in case you are not able to
175
00:19:29,470 --> 00:19:38,910
see it clearly and this time period is from
0 to 3 seconds. So, what I will do next is,
176
00:19:38,910 --> 00:19:56,830
I will change the time step to 0.1 and just
to make it a bit easier, to make see, we save
177
00:19:56,830 --> 00:20:03,760
it we plot it as a red.
178
00:20:03,760 --> 00:20:12,750
So if we look at this, you see that by changing
the time step to 0.1, we have slightly diff,
179
00:20:12,750 --> 00:20:18,390
the behavior is has slightly changed. You
see this red? If you look at this red response
180
00:20:18,390 --> 00:20:20,260
it is slightly different.
181
00:20:20,260 --> 00:20:45,780
If I go ahead and I am sorry let may let me
make this 0.5 s and look at the response.
182
00:20:45,780 --> 00:20:53,190
Well, you look what has happened. Your response
is absolutely different from what is expected.
183
00:20:53,190 --> 00:21:04,350
So, if you chose a time step like 0.5 in which
is in fact greater than 2 times 1 up on 5,
184
00:21:04,350 --> 00:21:11,180
so I had mentioned that, if h is greater than
2 by a the absolute value of 2 by a; in that
185
00:21:11,180 --> 00:21:16,140
case your response is altogether different.
You really have a situation where you are
186
00:21:16,140 --> 00:21:20,700
getting an unstable response. This is not
of course, the correct response. We say the
187
00:21:20,700 --> 00:21:26,860
numerical integration is blowing up. So, though
the actual system is stable, the numerical
188
00:21:26,860 --> 00:21:31,300
method gives you a wrong answer. So obviously,
you chose your time step carefully especially
189
00:21:31,300 --> 00:21:34,760
if you are using Euler method.
190
00:21:34,760 --> 00:21:40,630
Let us now go on and try to simulate the same
system using backward Euler method. So, let
191
00:21:40,630 --> 00:21:47,710
us just have the same thing. This is backward
Euler method. h is equal to 0.01, t final
192
00:21:47,710 --> 00:21:56,520
is 3, number of steps of course, is the depends
on the duration as well as h the initial value
193
00:21:56,520 --> 00:22:04,020
of x is 1 a is minus 5 so I will just save
this. yeah
194
00:22:04,020 --> 00:22:11,550
So as per backward Euler method, you will
have X of k or X of I is equal to i minus
195
00:22:11,550 --> 00:22:18,420
divided of 1 minus a h. You can easily derive
it from the basic formula for discretizing
196
00:22:18,420 --> 00:22:26,990
or you know using backward Euler method. Now,
this particular sentence has been commented
197
00:22:26,990 --> 00:22:32,980
so, I am not actually going to implement this
particular you know this particular command.
198
00:22:32,980 --> 00:22:37,410
What I am going to implement is this algorithm,
X of k is equal to X of k minus 1 divided
199
00:22:37,410 --> 00:22:42,990
by 1 plus a h.
So, let us just numerically integrate and
200
00:22:42,990 --> 00:22:48,990
see what we get. This is with a time step
of, I believe of, so we are getting the correct
201
00:22:48,990 --> 00:23:02,280
answer of course, the time step is 0 0.01
and a is minus 0.5. Now, I can make this also
202
00:23:02,280 --> 00:23:07,720
5. We had done that before. Let us see what
happens.
203
00:23:07,720 --> 00:23:15,670
Well, the solution which you get is this.
The interesting thing is that of course,there
204
00:23:15,670 --> 00:23:20,930
is an a kind of error you know. The earlier
with a lower time step we got this response,
205
00:23:20,930 --> 00:23:27,070
now we are getting this slightly erroneous
but, of course, it is still stable. So what
206
00:23:27,070 --> 00:23:33,950
what one one small point which we should remember
is that, backward Euler a stable system will
207
00:23:33,950 --> 00:23:35,540
in fact be shown as a stable system.
208
00:23:35,540 --> 00:23:42,010
So that is no worry. But, a very interesting
thing is what if I have got an unstable system?
209
00:23:42,010 --> 00:23:51,320
A is equal to 5 and I chose my time step more
than 2 times 1 upon 5. An interesting thing
210
00:23:51,320 --> 00:23:58,920
which you will see here is, that this system
which ought to be unstable, you know after
211
00:23:58,920 --> 00:24:04,431
all a is equal to plus 5. This is the real
system is unstable, e raise to 5 t into x
212
00:24:04,431 --> 00:24:08,900
of 0 that is the solution.
What you will find is, backward Euler, when
213
00:24:08,900 --> 00:24:13,770
you use backward Euler it is actually showing
this system to be a stable one. That is very
214
00:24:13,770 --> 00:24:19,790
interesting. So, it really this something
you should keep in mind that is backward Euler
215
00:24:19,790 --> 00:24:27,170
in fact is gives you, can also give you wrong
qualitative answers if your time steps stamp
216
00:24:27,170 --> 00:24:30,211
are larger this is. Of course, not true if
you chose a time step stamp to be small. So
217
00:24:30,211 --> 00:24:38,100
for example, if I chose 0.1 backward Euler
of course, should give you the correct answer.
218
00:24:38,100 --> 00:24:43,530
Correct in the sense that is an unstable system
so you see that the unstable system.
219
00:24:43,530 --> 00:24:49,050
So it just grows up and if your, it is not
visible on your screen, this is almost 1.2
220
00:24:49,050 --> 00:24:56,390
into e raised to nine. So 1.2 into 10 raised
to nine so, that really means that it is really
221
00:24:56,390 --> 00:25:02,830
blown out in 3 second itself. That is not
very surprising.
222
00:25:02,830 --> 00:25:10,740
So, this is about backward Euler. What if
I do the same system simulation using trapezoidal
223
00:25:10,740 --> 00:25:14,731
rule? Trapezoidal rule has an interesting
property that, whether your stable if you
224
00:25:14,731 --> 00:25:22,740
have a stable system, it will show be shown
as stable. So for example, this is 0.01 just
225
00:25:22,740 --> 00:25:28,600
a sec I will just check whether I have done
the right thing. oh I have to comment this
226
00:25:28,600 --> 00:25:35,110
statement and enable this statement which
is essentially implementing trapezoidal rule.
227
00:25:35,110 --> 00:25:41,890
This is x of k is equal to x of this is x
of I is equal to x of 5 minus 1 divided by
228
00:25:41,890 --> 00:25:47,640
1 minus h by 2 into 1 plus h by 2.
So you can actually work out this that you
229
00:25:47,640 --> 00:25:52,420
get this particular algorithm if you apply
trapezoidal rule to x dot is equal to a x.
230
00:25:52,420 --> 00:26:01,120
So, if you have got a is equal to minus 5
and h is equal to 0.1 as we saw just now.
231
00:26:01,120 --> 00:26:03,260
I save this. Run it.
232
00:26:03,260 --> 00:26:08,160
The Remember, although I called this file
backward Euler; still it is doing trapezoidal.
233
00:26:08,160 --> 00:26:15,910
It is using trapezoidal method of integration
now. So, what if I increase the time step?
234
00:26:15,910 --> 00:26:27,400
So, I am using trapezoidal rule. Just remember
that. So, for example, 5 this was unstable
235
00:26:27,400 --> 00:26:34,040
for example, Euler method when we try to integrate
it with h is equal to 5, it was unstable.
236
00:26:34,040 --> 00:26:39,100
What does trapezoidal rule, is also stable
of course, the errors in the response, so
237
00:26:39,100 --> 00:26:44,020
if I actually go on increasing h, the errors
increase. But, as we shall see, I will just
238
00:26:44,020 --> 00:26:54,900
do something more sillier. I will make h is
equal to 1 second. Even if I make h is equal
239
00:26:54,900 --> 00:27:00,500
to 1 second, although I am going to get an
error in the response it is not unstable.
240
00:27:00,500 --> 00:27:11,570
This is my new response. Similarly, if a is
equal to plus 5, if I chose h is equal to
241
00:27:11,570 --> 00:27:22,010
1, this is an unstable system.
We will just do it one on a fresh plot. This
242
00:27:22,010 --> 00:27:28,190
is an unstable system and trapezoidal rule
is showing it to be an unstable system. So,
243
00:27:28,190 --> 00:27:34,830
really trapezoidal rule you know, at least
in a linear system it kind of preserves the
244
00:27:34,830 --> 00:27:41,299
stability information. So I have of course,
if I make this time step smaller, the answer
245
00:27:41,299 --> 00:27:50,730
is of course, closer the right answer. yeah
So, you see this increasing here like this.
246
00:27:50,730 --> 00:28:01,520
So to summarize whatever I have shown so far,
is that when we are trying to numerically
247
00:28:01,520 --> 00:28:09,840
integrate a system, x dot is equal to minus
5 x; if you chose a Euler method, if you chose
248
00:28:09,840 --> 00:28:16,300
A h to be large, in fact if it is greater
than the absolute value of 2 by a, a is equal
249
00:28:16,300 --> 00:28:22,140
to minus 5, in that case you in fact get a
wrong qualitative picture about the stability.
250
00:28:22,140 --> 00:28:28,090
Of course, as you go on increasing the time
step, a numerical methods accuracy versions
251
00:28:28,090 --> 00:28:34,820
worsens. So, it is not unexpected that you
will get a better results if you chose h smaller.
252
00:28:34,820 --> 00:28:38,920
But, the important point which I can emphasize
is that, Euler method may give you a qualitatively
253
00:28:38,920 --> 00:28:43,780
wrong result also. It may show you a stable
system in fact to be unstable if your h is
254
00:28:43,780 --> 00:28:45,980
larger.
This does not happen with backward Euler or
255
00:28:45,980 --> 00:28:51,900
with trapezoidal rule. In fact, a stable system
is always shown to be a stable whether one
256
00:28:51,900 --> 00:28:57,010
uses backward Euler or trapezoidal rule. But,
interestingly with backward Euler method,
257
00:28:57,010 --> 00:29:03,260
I am essentially unstable system can be shown
as s stable system if your h is very large.
258
00:29:03,260 --> 00:29:11,330
So, if h is greater than point 2 upon a, the
modulus of it and unstable system for example,
259
00:29:11,330 --> 00:29:16,880
a is equal to plus 5 is shown to be stable
1 one whereas, trapezoidal rule preserves
260
00:29:16,880 --> 00:29:23,570
the stability information of a linear system.
So, that is one thing a very interesting point
261
00:29:23,570 --> 00:29:32,030
which we must note. Now, the basic you know,
the 3 three method which I have shown you
262
00:29:32,030 --> 00:29:37,200
to get accurate solution in fact for example,
x dot is equal to minus 5 x, it makes sense
263
00:29:37,200 --> 00:29:43,270
in fact to use methods with time steps of
0.1 0 also. You will get a reasonably accurate
264
00:29:43,270 --> 00:29:47,820
picture you know. It is another thing about
you know, if you chose h is equal to 0.5 etc
265
00:29:47,820 --> 00:29:53,450
etc. Not only you are going to inaccurate
picture but, in case of Euler method you will
266
00:29:53,450 --> 00:29:59,480
get even a qualitatively wrong picture.
So, you will use a it is a good idea if you
267
00:29:59,480 --> 00:30:04,090
have got a system with a time constant of
point 2 or as in this case we use a time step
268
00:30:04,090 --> 00:30:10,990
of 0.01 and so on. But, you know 1 you know
for example, if your situations were a system
269
00:30:10,990 --> 00:30:16,010
is marginally stable you know, you have got
a system in which the real part of the Eigen
270
00:30:16,010 --> 00:30:22,840
value of the continuous time system 1 near
0; Euler method is absolutely hopeless. You
271
00:30:22,840 --> 00:30:27,850
will find that it does not give you a good
you know, it gives you a qualitatively wrong
272
00:30:27,850 --> 00:30:34,570
picture. What I will do now is, try to simulate
1 one of the systems which we have simulated
273
00:30:34,570 --> 00:30:41,540
before, try to understand qualitatively, before
that this swing equations of a power system.
274
00:30:41,540 --> 00:30:48,590
So the swing equations of a single machine
infinite bus system or delta dot is equal
275
00:30:48,590 --> 00:30:53,280
to omega minus omega dot naught and the derivative
of omega minus omega naught which is nothing
276
00:30:53,280 --> 00:30:57,680
but, omega dot itself because, omega naught
is a constant which is the frequency of the
277
00:30:57,680 --> 00:31:03,750
infinite bus is equal to omega B by 2 H p
m minus k sin delta.
278
00:31:03,750 --> 00:31:10,470
K is in fact, we assume it to be a constant.
This is the internal voltage of the synchronous
279
00:31:10,470 --> 00:31:19,370
machine. The infinite bus voltage divided
by the reactants. If you have got a single
280
00:31:19,370 --> 00:31:32,710
machine connected to an infinite bus E 2 angle
0 and this is E 1. Well, this is a synchronous
281
00:31:32,710 --> 00:31:48,850
machine which is modeled as an internal. So,
this is our synchronous machine. So this is
282
00:31:48,850 --> 00:31:54,540
E 1 angle delta ,this is the internal voltage,
this is the rotor angle, x dash is a transient
283
00:31:54,540 --> 00:31:59,450
reactant. So, this these are the equations
swing equations. This is a non-linear system
284
00:31:59,450 --> 00:32:04,140
of equations.
So, if I am trying to solve this system of
285
00:32:04,140 --> 00:32:10,940
equations by Euler method, what is the kind
of response I will get? We saw for small disturbances,
286
00:32:10,940 --> 00:32:17,630
we have done the analysis of this system and
we saw that, in fact, it is for the equilibrium
287
00:32:17,630 --> 00:32:26,800
corresponding to which is less than ninety
degrees this one, for small disturbances around
288
00:32:26,800 --> 00:32:34,440
the equilibrium we will get essentially an
oscillated response. So, we know that, this
289
00:32:34,440 --> 00:32:39,750
if you are at this equilibrium, you will get
an oscillated response. If I give of course,
290
00:32:39,750 --> 00:32:47,530
a very large disturbance from this then, the
system may even go unstable. The synchronous
291
00:32:47,530 --> 00:32:53,980
machine may fall out of step or lose synchronism.
So, we know the kind of responses which can
292
00:32:53,980 --> 00:32:57,660
be expected at least for small disturbances
and qualitatively of course, we know that
293
00:32:57,660 --> 00:33:03,950
for large disturbances the system may go unstable.
So, with this this is the non-linear system,
294
00:33:03,950 --> 00:33:11,140
by the way this is not a linear system because
your delta appears within the sin function.
295
00:33:11,140 --> 00:33:19,580
So, what I will do is, try to integrate the
swing equations by Eulerâ€™s method.
296
00:33:19,580 --> 00:33:26,500
So, you just take, let us assume P m is 1,
h is 3 at this capital h is of course, the
297
00:33:26,500 --> 00:33:35,190
inertial constant, k let us assume is e 1
e 2 divided by 0.4, this is omega and omega
298
00:33:35,190 --> 00:33:38,200
base we assume is the same as this frequency
of the infinite bus.
299
00:33:38,200 --> 00:33:46,480
Let us take time step of 0.01, the final time
is 3, the number of steps of course, can be
300
00:33:46,480 --> 00:33:54,550
calculated from t final n h. Let us assume
that delta is equal to sin inverse, this sin
301
00:33:54,550 --> 00:34:02,510
inverse p m by k. In fact this is the equilibrium
value of delta. Omega is nothing. Omega, let
302
00:34:02,510 --> 00:34:10,119
us assume is omega b plus 1. Now, remember
one important point is that, although delta
303
00:34:10,119 --> 00:34:17,020
is at the equilibrium value of that the value
of delta is actually its equilibrium value.
304
00:34:17,020 --> 00:34:22,679
The value of omega is not its equilibrium
value. So, as a result this we are not actually
305
00:34:22,679 --> 00:34:27,040
at the equilibrium because both states have
to be at the equilibrium value if we are not
306
00:34:27,040 --> 00:34:31,629
to have any transient.
So, by adding this plus 1 we have in fact
307
00:34:31,629 --> 00:34:40,070
given an initial condition which is not the
equilibrium value of the states. So, we are
308
00:34:40,070 --> 00:34:45,590
going to get a transient. So, I have given
you initial condition which is not equal to
309
00:34:45,590 --> 00:34:49,450
the equilibrium value if I have given you
the equilibrium value, the system would simply
310
00:34:49,450 --> 00:34:55,500
not move. It would just remain where it is.
So, if I apply Euler method, you have got
311
00:34:55,500 --> 00:35:05,000
delta is equal to delta of i is equal to delta
of i minus 1 plus h. This is in fact simply
312
00:35:05,000 --> 00:35:09,880
omega minus, so we are just numerically integrating
the system by Euler method. The omega I is
313
00:35:09,880 --> 00:35:17,870
equal to omega i minus 1 plus h times omega
base divided by 2 times the inertial constant
314
00:35:17,870 --> 00:35:23,620
into p m minus k sin delta.
This is easy to evaluate if I know. What the
315
00:35:23,620 --> 00:35:28,470
initial values of delta and omega r, I can
find out the new values since we are not at
316
00:35:28,470 --> 00:35:32,570
equilibrium. Since, I have given in initial
condition which is not the equilibrium condition
317
00:35:32,570 --> 00:35:37,900
of course, if I made this 0, this 1 is 0.
We would be at equilibrium. So, what I will
318
00:35:37,900 --> 00:35:48,040
do is, I will just simulate this system. So,
what I will do is again, so I am evaluating
319
00:35:48,040 --> 00:35:52,420
this system using Euler method.
320
00:35:52,420 --> 00:36:02,630
So, if you look at the behavior of a system,
it looks as if it is an growing oscillation.
321
00:36:02,630 --> 00:36:07,790
In fact, if you look at, we have done the
small disturbance stability of this particular
322
00:36:07,790 --> 00:36:14,000
system and we know that response if in fact
not a growing oscillation but, a steady oscillation
323
00:36:14,000 --> 00:36:18,920
if those disturbance is small.
So obviously, there is a problem in the method.
324
00:36:18,920 --> 00:36:24,220
Let us just try reducing the time step. What
is the problem in reducing the time step?
325
00:36:24,220 --> 00:36:28,140
There is no problem. You will have to wait
longer for the simulation to get over in some
326
00:36:28,140 --> 00:36:33,720
time. That becomes a bit. That may become
a limiting factor for doing any kind of study
327
00:36:33,720 --> 00:36:35,800
if you if you make the time step too small.
328
00:36:35,800 --> 00:36:41,860
Now, you are getting you know, a different
answer that that shows that you have to actually
329
00:36:41,860 --> 00:36:43,350
reduce the time step.
330
00:36:43,350 --> 00:36:53,450
So, I will just redo this. In fact by using
a smaller time constant time step, we are
331
00:36:53,450 --> 00:36:58,870
in fact getting a better answer. But, you
see that this is still, it is not it is not
332
00:36:58,870 --> 00:37:03,450
you know, you this value is slightly less
than the value you have here.
333
00:37:03,450 --> 00:37:09,240
So what we see is that, if you have got you
know, marginally dam systems; Euler method
334
00:37:09,240 --> 00:37:15,020
in fact does not seem to be very nice because
whatever I reduce, if I reduce the time step
335
00:37:15,020 --> 00:37:21,200
I can do that even further. Now, it is going
to be painful because you will have to wait
336
00:37:21,200 --> 00:37:32,260
for a longer time. Now, it is still, it will
take some time for this to come up. yeah So,
337
00:37:32,260 --> 00:37:35,830
you see that there is slight difference in
the what I got with the previous time step
338
00:37:35,830 --> 00:37:41,900
and the time step you are having now. 1 One
thing, we can simulate for longer and longer
339
00:37:41,900 --> 00:37:49,500
time with this time step that is 0. 2 0 1.
But, the point you can actually try it out
340
00:37:49,500 --> 00:37:53,140
whatever time step you take, you will find
this eventually starts growing. If you try
341
00:37:53,140 --> 00:37:57,871
to simulate beyond 3 seconds, say twenty seconds,
even with this time step you will find that
342
00:37:57,871 --> 00:38:02,320
it is growing.
So, Euler, forward Euler method is not very
343
00:38:02,320 --> 00:38:10,090
nice to simulate systems with (( )) which
are marginally damped or have poor damping.
344
00:38:10,090 --> 00:38:20,640
So, let us just rerun this system using trapezoidal
rule. So I am using exactly the same system
345
00:38:20,640 --> 00:38:30,140
with a time step of 0.01 and 1 one problem
comes here, now if I am going to try to do
346
00:38:30,140 --> 00:38:37,070
integrate this method using trapezoidal rule,
what is going to be the problem? The problem
347
00:38:37,070 --> 00:38:44,090
which you see here of course, is that the
statements in the program seem to have become
348
00:38:44,090 --> 00:38:45,870
much larger large here.
349
00:38:45,870 --> 00:38:53,840
One other thing you will notice is the delta.
I will not be easily obtained from the previous
350
00:38:53,840 --> 00:39:01,540
value that is delta i minus 1. So, let us
just pay our attention on what we get if we
351
00:39:01,540 --> 00:39:05,460
try to discretize swing equation using trapezoidal
rule.
352
00:39:05,460 --> 00:39:11,970
So you have got delta dot is equal to omega
minus omega naught. So, we will have delta
353
00:39:11,970 --> 00:39:29,650
k plus 1 is equal to delta k plus h by 2 times
omega k plus 1 minus omega naught plus omega
354
00:39:29,650 --> 00:39:44,700
k minus omega naught. So of course, we require
the value of k plus 1 here. Omega k plus 1
355
00:39:44,700 --> 00:40:01,480
will be equal to omega k plus h by 2 into
omega b by 2 h this is the into p m minus
356
00:40:01,480 --> 00:40:13,830
k sin delta k plus 1 plus p m. we assuming
of course, p m is a constant sin delta k.
357
00:40:13,830 --> 00:40:18,270
So, this is the discretization which you will
have to use and that is what is really been
358
00:40:18,270 --> 00:40:22,800
programmed there. Now, the problem here of
course, is that you do not know the value
359
00:40:22,800 --> 00:40:29,080
of delta k plus 1, you do not know what omega
k plus 1 is. But, you require that in order
360
00:40:29,080 --> 00:40:34,280
to evaluate delta k plus 1. So, if I want,
if I given have given you delta k and omega
361
00:40:34,280 --> 00:40:40,330
k how do I get omega k plus 1 and delta l
plus 1? You will have to use an iterative
362
00:40:40,330 --> 00:40:46,430
method. So that is what exactly I have done
in the program. What I have done is, try to
363
00:40:46,430 --> 00:40:57,920
use within every time step within every time
step solve this using Gauss Seidel method.
364
00:40:57,920 --> 00:41:06,660
So, what we do is get? Plug in this value
of delta k and omega k. You will plug it in
365
00:41:06,660 --> 00:41:15,110
here, here, here and here. You plug it in.
Omega k plus 1 is something you do not know.
366
00:41:15,110 --> 00:41:22,770
So, let us guess that it is equal to omega
k in the first instant and this is delta k.
367
00:41:22,770 --> 00:41:29,520
So, you run this. You you evaluate this assuming
this is omega k and this is delta k. You will
368
00:41:29,520 --> 00:41:37,080
get a new value of delta k plus 1 and omega
l plus 1. Pug in this here and you plug you
369
00:41:37,080 --> 00:41:42,130
plug in this here, you plug in this here and
reevaluate delta k plus 1 and omega l plus
370
00:41:42,130 --> 00:41:46,380
1. They basically, we are following Gauss
Seidel method of solving this algebraic equations
371
00:41:46,380 --> 00:41:51,440
set of algebraic equations.
You go on doing this till, there is no difference
372
00:41:51,440 --> 00:41:58,460
between a previously evaluated del delta k
plus 1 and the newest value of delta k plus
373
00:41:58,460 --> 00:42:03,120
1 which you are getting. So, that is exactly
what I am doing this here. The program is
374
00:42:03,120 --> 00:42:07,320
a bit complicated and probably you will find
it difficult to follow this these steps but,
375
00:42:07,320 --> 00:42:10,750
I will just indicate what I am doing.
376
00:42:10,750 --> 00:42:23,490
So, for what I do here is, essentially in
these steps. I am in fact using Gauss Seidel
377
00:42:23,490 --> 00:42:30,040
method. I have taken a maximum number of iterations
per time step as 10 hope with a hope that
378
00:42:30,040 --> 00:42:35,510
of course, that will converge. So, in fact
I am checking out the convergence criterion
379
00:42:35,510 --> 00:43:12,730
as well. So, let just me re go through the
whole program slowly. This is the check on
380
00:43:12,730 --> 00:43:32,630
the convergence. So, I will just run this
of course, one small point will give the same
381
00:43:32,630 --> 00:43:41,680
disturbance as before, in the sense that the
this was one in the previous. So I have given
382
00:43:41,680 --> 00:44:03,140
a small disturbance. We are away from the
equilibrium. Save this. yeah So, this is basically
383
00:44:03,140 --> 00:44:15,180
what you get. We actually evaluated just to,
it this is not the what we get with trapezoidal
384
00:44:15,180 --> 00:44:16,180
rule.
385
00:44:16,180 --> 00:44:21,651
Yeah Now, we will get what we should. So,
this is the response and actually this is
386
00:44:21,651 --> 00:44:27,520
the sinusoid does not grow with time. So,
this is a, this is a true. In fact, reasonably
387
00:44:27,520 --> 00:44:34,460
accurate picture of how delta varies if the
system is given a small disturbance so even
388
00:44:34,460 --> 00:44:39,950
if I increase the time step remember, here
this is this will not really be affected.
389
00:44:39,950 --> 00:44:50,250
So, I just changed my time step to 0.1 10
times more. Save this. Rerun it.
390
00:44:50,250 --> 00:44:54,780
Because of course, I have chosen a larger
time step. There is slight you know, kind
391
00:44:54,780 --> 00:45:05,040
of poki poki solution here. This particular
sinusoid has got rough edges but, on the whole
392
00:45:05,040 --> 00:45:11,010
it does not give you a wrong picture of stability.
It is a undammed undamped sinusoid. Delta
393
00:45:11,010 --> 00:45:15,530
is an undammed undamped sinusoid. So, this
is basically what you get with trapezoidal
394
00:45:15,530 --> 00:45:26,180
rule. If I increase the disturbance, I increase
the disturbance I increase the disturbance
395
00:45:26,180 --> 00:45:40,610
in that case you even may go unstable. So,
I am now, I am going to give a large disturbance
396
00:45:40,610 --> 00:45:43,040
relatively wrong.
397
00:45:43,040 --> 00:45:49,530
You see that, because of this larger disturbance,
the delta oscillations are become larger.
398
00:45:49,530 --> 00:45:53,320
And one interesting thing that this is no
longer sinusoidal looking you see this top
399
00:45:53,320 --> 00:45:59,570
part looks much fatter than the bottom part.
So, obviously with the if you make the disturbance
400
00:45:59,570 --> 00:46:07,860
larger, your behavior starts moving into the
non-linear zone. And you do not get the sinusoidal
401
00:46:07,860 --> 00:46:11,380
response as predicted by linear analysis.
402
00:46:11,380 --> 00:46:24,900
So, if I actually make this disturbance even
larger, you are going to get, the system looses
403
00:46:24,900 --> 00:46:29,230
synchronism. You see that, this angle just
goes on increasing with time. So, if I go
404
00:46:29,230 --> 00:46:36,420
on increasing the, you know deviation from
the equilibrium. We go into non-linear behavior
405
00:46:36,420 --> 00:46:43,750
and this is off course captured quite well
by this numerical integration method. So,
406
00:46:43,750 --> 00:46:49,500
this is basically what I wish to tell you
about using trapezoidal rule. For small disturbances
407
00:46:49,500 --> 00:47:01,400
it gives you the correct response. If I use
backward Euler method for the same problem,
408
00:47:01,400 --> 00:47:07,190
Backward Euler method is also an implicit
method. Remember, so I have to use an iterative
409
00:47:07,190 --> 00:47:14,810
solution to get this iterative method to get
the values of the states at every time step.
410
00:47:14,810 --> 00:47:35,940
So, let us just do this again with backward
Euler. So, an interesting thing to observe
411
00:47:35,940 --> 00:47:42,420
here is that, when I use backward Euler on
a marginally stable system which is just having
412
00:47:42,420 --> 00:47:49,550
0 damping; Euler, backward Euler shows it
as if it is a stable system as a damped oscillation.
413
00:47:49,550 --> 00:47:57,030
So, that is an interesting point about backward
Euler. It really changes the picture. It s
414
00:47:57,030 --> 00:48:02,520
means especially your system is a marginally
damped, it is got very little damping. So
415
00:48:02,520 --> 00:48:07,590
it is showing you the system to be a stable
system whereas, actually it as undammed oscillation.
416
00:48:07,590 --> 00:48:11,550
So this system is a showing as if it is dying
down.
417
00:48:11,550 --> 00:48:17,720
So that is one interesting observation about
backward Euler method. Of course, if I use
418
00:48:17,720 --> 00:48:24,350
backward Euler method with a very small time
step, you are going to get reasonably accurate
419
00:48:24,350 --> 00:48:30,160
solution. So, I just numerically integrated
with a small smaller time step.
420
00:48:30,160 --> 00:48:35,900
Well, we are coming, the answer is slightly
different and we are coming closer to the
421
00:48:35,900 --> 00:48:41,940
solution we have in mind. So of course, reducing
time step is 1 one of the ways you can get
422
00:48:41,940 --> 00:48:55,710
accurate solution or a more accurate solution.
So, with this I end my discussion about the
423
00:48:55,710 --> 00:49:04,300
numerical integration of swing equations.
My discussion about numerical integration
424
00:49:04,300 --> 00:49:11,380
method would be quite incomplete if I do not
consider what are known as stiff systems.
425
00:49:11,380 --> 00:49:16,510
We have learnt about stiff systems. Stiff
systems are systems in which the Eigen values
426
00:49:16,510 --> 00:49:23,010
of the of linear system especially are far
apart. In the sense that there were small
427
00:49:23,010 --> 00:49:29,690
if there is large Eigen value which translates
to the facts that the systems has got both
428
00:49:29,690 --> 00:49:36,290
fast and slow transients. And if you are actually
interested only in the slow transients, then
429
00:49:36,290 --> 00:49:42,950
you really need not modulate in great the
fast transients in detail. In fact I showed
430
00:49:42,950 --> 00:49:48,480
to you in 1 one of the previous lectures that
when you are trying to analyze the system
431
00:49:48,480 --> 00:49:52,780
with slow and fast transients the fast transient
differential equations corresponding to the
432
00:49:52,780 --> 00:49:59,820
fast transients of the states associated with
the fast transients, can be converted to algebraic
433
00:49:59,820 --> 00:50:06,150
equations and we really do not loose out much
as far as the information on the slow transients
434
00:50:06,150 --> 00:50:12,810
is concerned.
A similar issue will kind of crops up when
435
00:50:12,810 --> 00:50:19,610
you have got even a numerically integrating
stiff system. So if you got a stiff system,
436
00:50:19,610 --> 00:50:26,830
you got large and small Eigen values and the
question comes, what is going to be the choice
437
00:50:26,830 --> 00:50:31,180
of time step? So, if you have got of a stiff
system and you have got fast transients, you
438
00:50:31,180 --> 00:50:38,230
will have to chose a time step which is compatible
with the fastest transient of the system in
439
00:50:38,230 --> 00:50:46,609
order to accurately get the response. So,
if you got a system with slow as well as fast
440
00:50:46,609 --> 00:50:53,750
transients inherent then, you have to chose
a time step to be compatible with the fastest
441
00:50:53,750 --> 00:50:57,600
transient so that, your fast transient are
accurately obtained.
442
00:50:57,600 --> 00:51:02,720
But consider situation where you are not really
interested in the fast transients. So, how
443
00:51:02,720 --> 00:51:12,710
will you chose a time step? In Euler method,
it is very clear that if 1 plus lambda h the
444
00:51:12,710 --> 00:51:18,450
absolute value is greater than 1 then, the
numerical integration method is going to blow
445
00:51:18,450 --> 00:51:24,570
up. So even if your are fast transients are
stable, if I chose my time step which is not
446
00:51:24,570 --> 00:51:31,220
compatible with this particular condition
then; you your numerical integration as I
447
00:51:31,220 --> 00:51:36,080
said will blow up. But, what if I am interested
only in the slow transients unfortunately
448
00:51:36,080 --> 00:51:42,390
with Euler method will be forced to take a
time step which is small, very very small
449
00:51:42,390 --> 00:51:47,790
which is compatible with fast transient.
So, that is one of the problem. So, even if
450
00:51:47,790 --> 00:51:52,010
you interested in the slow transients, you
have chose to choose very small time steps
451
00:51:52,010 --> 00:51:56,720
so that is 1 one of the problems which you
will encounter when you try to use Euler method
452
00:51:56,720 --> 00:51:58,790
to integrate a stiff system.
453
00:51:58,790 --> 00:52:01,480
.
So, let us again look at a stiff system. We
454
00:52:01,480 --> 00:52:27,859
have done this example before. This is i l
2, this is i l 1, this is v c. We know that
455
00:52:27,859 --> 00:52:35,440
this particular system has both fast and slow
transients. So this is 10 mille Henry, this
456
00:52:35,440 --> 00:52:40,140
is 1 Henry. In fact, the state associated
with the slow transient is this and the states
457
00:52:40,140 --> 00:52:44,140
associated with the fast transient are these
from the participation factors which we had
458
00:52:44,140 --> 00:52:51,920
evaluated in a couple of lectures before.
So, if you have got this particular system,
459
00:52:51,920 --> 00:52:58,250
this contains slower and faster transients.
In fact how do you know that it contains slow
460
00:52:58,250 --> 00:53:04,200
and fast transient? affect In fact if you
have, if you know the a matrix you have got
461
00:53:04,200 --> 00:53:16,600
X dot is equal to A X where X is equal to
i l 1 i l 2 and v c then this a is equal to
462
00:53:16,600 --> 00:53:28,750
minus 10 0 minus hundred 0 0 1. We have done
this in the previous class so, I am not reworking
463
00:53:28,750 --> 00:53:39,220
out this, minus 10000 and 0. This is a, the
Eigen values of this system are approximately,
464
00:53:39,220 --> 00:53:48,710
minus 5 plus or minus j lambda 1 and lambda
2. So they are the two 2 actually complex
465
00:53:48,710 --> 00:53:56,240
conjugate pair.
So, this is a stable system but, one of the
466
00:53:56,240 --> 00:54:01,280
Eigen values is very small and the other one
is very large. This is the example of a stiff
467
00:54:01,280 --> 00:54:07,960
system. If I try to numerically integrate
the system stiff system using Euler method,
468
00:54:07,960 --> 00:54:20,510
I will have to chose my time step so that,
1 plus lambda 1 h is less than 1. So, h should
469
00:54:20,510 --> 00:54:25,900
be such such that this is satisfied. In fact,
it should be satisfied for all the 3 1 2 and
470
00:54:25,900 --> 00:54:31,910
3, all the 3 three Eigen values. So 1 one
thing you will notice that, if you want to
471
00:54:31,910 --> 00:54:36,700
satisfy it for all the 3 three Eigen values
then, h will have to be very small. In fact,
472
00:54:36,700 --> 00:54:43,910
it becomes completely unstable in case you
chose your h say as 0.0 0 1 or so.
473
00:54:43,910 --> 00:54:49,480
So, I will just do this numerical integration
using Euler method and then we will conclude
474
00:54:49,480 --> 00:54:57,040
this particular lecture. So, let us do this
1 one integration. This is a system; this
475
00:54:57,040 --> 00:55:06,310
is a matrix I have chosen a very small time
step. I will just simulate it for a short
476
00:55:06,310 --> 00:55:10,290
while because, I know it is really not going
to work out with Euler method. Even with this
477
00:55:10,290 --> 00:55:16,410
small time step, even with this small time
step we have a problem with Euler method.
478
00:55:16,410 --> 00:55:21,240
So, the first thing I will do is integrate
it using Euler method. This is a remember,
479
00:55:21,240 --> 00:55:28,869
a higher order system. So, instead of 1 plus
small case a H you have got identity matrix
480
00:55:28,869 --> 00:55:37,240
plus A into H into X plus H into u.
So, the I am implementing Euler method for
481
00:55:37,240 --> 00:55:42,590
this system. Of course, there is an input
here. So, I have used Euler method, I have
482
00:55:42,590 --> 00:55:48,720
generalized it to a situation where you have
got an input is as well. So, if I numerically
483
00:55:48,720 --> 00:56:04,420
integrate this system; hoops Oops! There is
an error here. So, I will just check out what
484
00:56:04,420 --> 00:56:32,240
is the problem and we will redo it. The problem
probably is, there is some error here. So,
485
00:56:32,240 --> 00:56:39,590
we will what we will do is, it do is redo
this example again in the next lecture, will
486
00:56:39,590 --> 00:56:44,470
just debug the error which is there in the
program and rerun this example in the next
487
00:56:44,470 --> 00:56:47,500
class.
So, in next class we will continue with this
488
00:56:47,500 --> 00:56:52,020
particular example, see how it is in fact
a very interesting and important thing. We
489
00:56:52,020 --> 00:56:58,440
are going to learn in the next lecture that,
for stiff systems of this kind; you have to
490
00:56:58,440 --> 00:57:04,080
use specific numerical method. Even you know,
explore if you want to capture. In fact, both
491
00:57:04,080 --> 00:57:08,430
the fast and slow transients you may even
have to explore things like variable time
492
00:57:08,430 --> 00:57:15,170
steps and so on. So, let us conclude this
lecture at this point. We will redo this example
493
00:57:15,170 --> 00:57:18,720
which did not work out in this particular
lecture again in the next lecture.