Tip:
Highlight text to annotate it
X
Hello and welcome to lecture 6 of module 7, where we were discussing numerical methods
to solve ordinary differential equations - initial value problems. What we started off with was
Euler's equation, and in yesterday's lecture, we considered one more example with Euler's
equation and saw that Euler's equation under certain condition becomes unstable, that means,
we go to minus infinity or plus infinity or we have highly oscillatory behavior.
So, what I am going to do in today's lecture is going to talk about the stability issues
with respect to the Euler's method, the implicit method and the explicit method. Before actually
going to that, what I will do is just spend a few minutes recapping what we have done
so far; I think, every lecture of this module, we are just starting with recapping all of
these methods simply to put all these methods, because we are covering a fair number of methods
in this particular module compare to any of previous modules.
So, I want to put all these various methods that we are considering into perspective.
For example, when it comes to some of my own work, my own research work, I tend to use
not just one method, but I tend to try a few of these methods before figuring out that,
this method is going to be work the best for me. There are cases, where I have ended up
using the second order RK method; there are cases, where have ended up using the forth
order RK method and so on. As far as possible, we start with RK 4 method
and wherever we find certain problems using RK 4 methods, we will go to some of the less
accurate methods, but some of them which are actually more stable, for example, through
implicit methods. So, I will talk about the implicit methods and the explicit methods
compare their stability, we will just cover it with respect to Euler's method and then
make general conclusions about how this is going to work.
So, that is going to be our strategy for today's lecture. So, let us consider Euler's method.
And - Euler's method - explicit method was, y i plus 1 equal to y i plus h multiplied
by f (y i, t i); so, this was Euler's explicit method. And Euler's implicit method, y i plus
1 equal to y i plus h times f of y i plus 1, p i plus 1. Then we went on to the Runge-Kutta
family of method; and in that case, I will talk about Heun's method - RK 2 version.
So, in Heun's method, the RK 2 version of Heun's method, what we do is, we write our
k 1 equal to f( y i, t i); k 2 is f of y i plus h times k 1, t i plus 1; so, at this
stage, what we can do is, we could perhaps not we could perhaps, we could write y i plus
1 bar 0 equal to f of i plus i y plus 1 bar 0, we could write this equal to y i plus k
1; and then this particular function, we will be able to write this as f of y bar 0 i plus
1. And our y i plus 1 is y i plus h times s and
s is h by 2 multiplied by k 1 plus k 2; so, y i plus h by two multiplied by k 1 plus k
2; k 1 is nothing but f( y i, t i); and k 2 is nothing but f of y i plus 1 bar 0, t
i plus 1; so, this is the RK 2 version of Heun's method. The predictor-corrector version
of Heun's method is, this becomes a predictor equation and recursive of use of this equation
becomes a corrector equation. So, we have y i plus 1 bar 0 equal to sorry
y i plus h times - I miss then h - I miss the h over here, so I will just correct this
thing plus h times f( y i, t i); and this is the predictor; and the corrector equation
is, y i plus 1 bar m plus 1 equal to y i plus h times f of y i plus 1 bar m, t i. So, that
is the corrector equation, sorry this going to be h by 2 this plus f (y i, t i); so, this
is what we are going to use as the corrector equation for the predictor-corrector Heun's
method; and if we keep using this particular equation until convergence, we will get the
Crank-Nicholson method.
So, what we mean by the Crank-Nicholson method is, when this particular quantity becomes
equal to this particular quantity, we call that particular iteration has converged and
we denote this as y i plus 1; so, we replace y i plus 1 and y i plus 1 inside of y bar
m and y bar m plus 1 and we will get this equation; and it is a non-linear equation
that needs to be solved in this Crank-Nicholson method, which is a semi implicit method; so,
that is the RK 2 version of Heun's method predictor-corrector version and the semi implicit
method, which is essentially the Crank-Nicholson method for solving this particular problem;
what we realize over here, when we do this predictor-corrector method is, this is a trapezoidal
rule, this being a trapezoidal rule has an order of h cube accuracy; but this is just
an explicit Euler's method; being an explicit Euler's method this has an order of x square
accuracy. So, the question that we are trying to ask
now is, can we actually do better in the predictor-corrector Heun's method then using an x square accuracy
for the predictor; and the answer to that, of course, is yes and we will have a more
accurate version of the Heun's method.
So, modification to Heun's method and this particular modification will be implemented
in order to improve the accuracy of the predictor equation. In this case, we will recall our
result that the central difference formula during numerical differentiation were more
accurate then the forward difference formula using basically parallel arguments; to those
arguments we will convert this particular formula from y i 0 equal to y i, h of h multiplied
by f ( y i, t i), we will convert it into a slightly better formula.
And what that better formula is, is the numerical difference dy by dt we can write that equal
to y i plus 1 minus y i minus 1 divided by 2 h, now dy by dt was equal to f (y i, t i);
as a result, our improved - our - method with an improved accuracy is going to be, y i plus
1 bar 0 is going to be equal to y i minus 1 plus 2 h multiplied by f ( y i, t i). So,
that is the predictor equation. And the corrector equation remains the same
y i plus 1 bar m plus 1 equal to y i plus h by 2 multiplied by f of i plus f of i plus
1 bar m; I am just using shorthand notations over here, because I do not want to write
this entire expression all over again. So, we have the equation dy by dt equal to f subject
to y at 0 equal to sum y 0. So, when i equal to 0, we are going to get
y 1 bar is going to be equal to y minus 1 plus 2 h times f(y 1, t 1); but now the problem
over here is, y minus 1 is not known, we have started with - sorry not y 1 comma t 1 - y
0 comma t 0; we have started with an initial condition for y 0 y 1 represents the value
at one time prior to this y 0 and we do not know this particular value and this particular
value is not given to us; this particular value is unknown; as a result, this predictor
equation cannot be used at i equal to 0. So, the predictor equation cannot be used
at i equal to 0; as a result, this improved Heun's method is non-self-starting method;
in other words, what happens is, we need a method - a different method - to be implemented
at i equal to 0, after that from I equal to 1 2 3 and so on, we can actually start implementing
then all itself starting method. So, how to get out of this problem? Well,
to get out of this problem is straight forward at i equal to 0, that means, the first time
that you are going to take the first step in this integration problem; we do not implement
the modify Heun's method, but instead we implement the original Heun's method; so, at i equal
to 0, the idea is implement the original Heun's method; at I equal to 1 2 3 and so on, implement
the modified Heun's method. So, the ad-hoc solution we can say...
So, the way to handle this particular problem is, at i equal to 0 we will use the original
Heun's method, which means the predictor will be of the form y i plus 1 bar equal to y i
plus h times k 1; so, in other words, y 1 bar 0 equal to y 1 plus h time k 1, that is
what we use for the predictor. We will use this particular equation for the corrector,
keep in mind that, the corrector equation has not change from Heun's method to the improved
Heun's method. Now, when after as we solve this particular
problem at i equal to 0, we move onto i equal to 1; and i equal to 1 we can actually use
this predictor, because y bar 2 equal to y bar 0, which is y 0, which is a known quantity
plus 2 h times f(y 1, t 1); y 1 is known; t 1 is known; as a result, we can now use
this method. So, a non-self-starting method and Heun's
method is just one example of non-self-starting method. In the next lecture, we are going
to consider a few more examples, specifically, we will consider the Adam-Bashforth and Adam-Moulton
family of methods, which are actually non-self-starting multistep methods and over there as well we
are going to use ideas similar to this is, for i equal to 0 or i equal to 0 and 1 so
on and so forth. We will use a less accurate method to get
our method of choice started; and once the method of choice get started we will go on
and keep using our method of choice in order to improve the accuracy.
So, now, what we have done is, we have talked about Heun's method - RK 2 method, we have
introduced ourselves to the predictor-corrector methods. Now, we will go back and analyze
what happened in the previous lecture of this particular module, what happened when we saw
that the method using the Euler's explicit method did not converge, but instead we saw
oscillatory behavior.
So, what we saw in previous lecture? We changed our p f r equation to a first order equation
of the form dc by dv equal to minus 1 by 2 c to the power 1, instead of 1.23 we changed
to c to the power 1. When we choose again starting with c at 0 equal to c 0 equal to
1 unit moles per liter - moles for meter cube, whatever the appropriate unit might be, we
started at that. We solved this particular problem, you at
h equal to 1 and we got a stable solution; the solution was not very accurate using the
Euler's method, but at least the solution was stable; we got concentration c at 5 as
some finite value. Next, what we did was, we increase the h to h equal to 2, when we
increase the h 2 h equal to 2, c 0 was equal to 1, c 1 - me immediately saw was became
- equal to 0 and after that it remained at 0.
C 2 c 4 and so on were all equal to 0 for h equal to 2; so, what we saw happening with
h equal to 1 was, something of this sort; what we saw happening at h equal to 2 was,
in one step itself the solution came to 0 and then it stabilize at 0.
Next, we considered a higher, h equal to 2.1, at h equal to 2.1 what happened is c 2.1 was
negative, c 4.2 was positive, and so on and so forth; so, what we saw is, we saw some
amount of oscillations taking place in this particular system. Specifically, when we plot
concentration versus volume, what we observed was that the concentration drops and then finally reaches 0. So, this is what
we saw and these oscillations for h equal to 2.1 were relatively small, they got damped
and very quickly the final solution stabilizes. Then we increase the h further and we increase
the h to h equal to 4, at h equal to 4, we saw c 4 equal to minus 1, c 8 that we observed
was equal to 1, c 12 was minus 1, c 16 was 1, so on and so forth; so, what we had observed
in that particular case was, behavior of this type; so, each alternate values was negative
or positive so on and so forth. That is the behavior that we saw for at h
equal to 4; and finally, when we increased from h equal to 4 to h equal to 5, I think
we increase it to h equal to 5, but any value greater than 4 will suffice, what happened
was, that c 5 c 10 c 20 on so on that is increasing and finally the concentration blew up; so,
if you were to plot this particular guy for h equal to 5, we will actually see a plot
of this type. So, what we see is that, this guy keeps expanding
and we will either reach plus infinity or minus infinity alternatively as volume v tends
to infinity; so, this is what we saw with Euler's explicit method; so, what happens
over here in this range is, that the solution converges monotonically; keep in mind that,
this is a plug flow reactor, this is in physically in a plug flow reactor, the concentration
cannot go in a fluctuating manner, it has to actually smoothly go to the final steady
state, that is what we will physically observe in a plug flow rector in which a first order
in which actually any reaction any single reaction is taking place.
The next case was, when h was between 2 and 4, between 2 and 4 what we saw was, there
were oscillations; and when h was increased beyond 4 we saw that the system was unstable;
we will now go ahead and analyze this behavior for the Euler's explicit method, try to find
out why we get this particular type of a behavior and then we will extend the same analysis
to the Euler's implicit method.
So, I will go back over here, where I had initially written down the expression for
Euler's implicit and Euler's explicit methods. We will first we will now look at the Euler's
implicit method. The equation that we started off with was of the form dy by dt equal to
minus lambda y, where lambda in this particular example was half.
So, now, if you want to solve this particular equation analytically, we divide throughout
by y and we take dt on the other side and we integrate and we will get ln of y equal
to minus lambda t plus constant c at t equal to 0, we have y equal to y 0, so the constant
c is going to be equal to ln of y 0; and we just rearrange this and we will get y at any
time t equal to y 0 e to the power minus lambda t. Now, exponential function is always going
to be stable limit as t tends to infinity, e to the power minus lambda t is going to
tend to 0 is for positive values of lambda. So, the analytical solution is stable for
lambda greater than 0, plus the behavior of the exponential function e to the power minus
lambda t, with t is going to look somewhat like this, which basically means and finally
it should settle down at 0, it would not go below 0.
This is the function y against t starting at value y 0; so, what we will see, so the
two physical features of the analytical solution of this particular equation is 1, that the
equation settles down to t 2 y equal to 0 and it is second thing is it is settles down
monotonically provided lambda is greater than 0; if lambda is less than 0, of course, we
are going to be at unstable solution. Now, what we want from over numerical method
is, the first thing is that the numerical method of our choice should converge to a
value, what did it means by, should converge to value is that limit as i tends to infinity
or limit as t tends to infinity, y should converge to a finite value, this is of course
true, only if we have the original analytical equation being stable, we are just considering
those cases for now. So, we want the value of y i to be equal to
a finite value as i tends to infinity. So, let us look at the explicit - methods - method
first. So, we will have y i plus 1 equal to y i plus h times f(y i, t i) and the function
f in this particular case is minus lambda y, so f (y i, t i) when we substitute over
here is going to be minus lambda multiplied by y i.
So, we can write this as, 1 minus lambda h multiplied by y i; so, y i plus 1 is 1 minus
lambda h multiplied by y I, that is what we get; now, this we can write it as, 1 minus
lambda h multiplied by y i multiplied by y i minus 1. So, we will have this as 1 minus
lambda h squared y i minus 1, which we can write it as 1 minus lambda h cubed y i minus
2 and so on equal to 1 minus lambda h to the power i plus 1 multiplied by y 0.
Now, this is our value of y i. Now, let us consider the three cases, the case where y
i decrease monotonically, the case where y i had oscillations, and the third case where
y i did not have any oscillations, but sorry y i had oscillations, but it was unstable.
So, in the case of monotonically decrease is when we get y i behavior something like
this and with y i ending at y equal to 0 as i increases, for this to happen, what we need
is that, 1 minus lambda h should actually lie between minus 1 and plus 1; that is required
for stability; and this particular guy should be positive, that is required for monotonicity;
so, for this particular thing is that, let us call the 1 minus lambda h as alpha for
y i plus 1 or rather y i equal to alpha to the power i multiplied by y 0; for this particular
equation to be monotonic, we need alpha to be less than 1 and greater than 0.
This is what the alpha should satisfy in order for this particular equation to stabilize
monotonically for us to observe some oscillations, but not whole out of oscillations; if let
say, if alpha was equal to minus half, what we will get is y 1 equal to minus half multiplied
by y 0; so, y 0 starts at 1, this becomes minus 0.5, next y 2 is minus half multiplied
by y 1 and y 1 is minus 0.5, so y 2 becomes plus 0.25, y 3 becomes minus 0.125, y 4 becomes
plus 0. 0625 so on and so forth. So, for this behavior, we need alpha to be
between minus 1 less than alpha, less than less than 0, or actually I should be writing
as equal to... when the alpha is equal to 0, at that time, we immediately will hit this
particular value and then flat flatten out at y equal to 0; and for the case where we
get instability, we need absolute value of alpha to be greater than 1. That means, alpha should be less than
minus 1 or alpha should be greater than minus 1. So, these are the various three conditions.
And let us evaluate the first two conditions now. Let us, let us, evaluate the first condition,
where 0 should be less than or equal to alpha should be less than 1 and we substitute alpha
equal to 1 minus lambda h, so we will have 0 less than equal to 1 minus lambda h less
than 1, or we can take 1 on this side and we will get this as minus 1 and this will
become 0. So, we will have minus 1 less than equal to
minus lambda h less than 0 or 0 less than lambda h less than equal to 1. So, as long
as h... now, h is a positive quantity, so this becomes a redundant thing, h we always
know is going to be positive; so, as long as the value of h is less than equal to 1
by lambda, we will get monotonic behavior. So, that is the first result, that the overall
numerical solution of the Euler's explicit method is going to be monotonic and monotonically
decreasing for a linear system as long as h is less than or equal to 1 by lambda.
Now, let us then consider the next question, when do we get oscillations? We get oscillations,
we will get oscillatory but stable behavior, if alpha is less than 0 and less than minus
1; that means, minus 1 less than 1 minus lambda h less than 0 or minus 2 less than negative
lambda h less than minus 1, or 1 less than lambda h less than 2 or h should lie between
1 by lambda and 2 by lambda. So, if h is between 1 over lambda and 2 over lambda, we will get
the system to be stable but oscillatory. If h is negative, we will get the system to
be unstable; if h is greater than 2 by lambda, we will get the system to be unstable again.
If h is greater than 2 by lambda, will not only get the system to be unstable, in addition
to being the system being unstable, we will also see oscillations; so, system is both
oscillatory and unstable. So, to summarize this results...
So, summary of stability results for Euler's explicit method is that, if h lies between
0 and 1 by lambda, we will get monotonic and stable behavior; if h lies between 1 by lambda
and 2 by lambda, we will get a stable though oscillatory behavior; if h is greater than
2 by lambda, we will get an unstable behavior. So, I will start with 0 over here, I will
put a point location 1 by lambda put another location 2 by lambda and then these are the
values greater than that; if h lies in the red region, we have stable and monotonic behavior;
if h lies in this yellow region, we have stable but oscillatory behavior; and if h lies beyond
this, we have unstable behavior.
So, let us go to the result of the previous lecture. In this particular case, lambda was
equal to half, so 1 by lambda is 2 and 2 by lambda is 4, between 0 and 2 we get stable
and monotonic behavior, between the 2 and 4 we get stable but oscillatory behavior,
at 4 we get exactly oscillatory behavior; the oscillations do not increase, they do
not dampen, they keep going like this, because when h equal to 4, we get our alpha value,
that is 1 minus lambda h is going to be exactly equal to minus 1.
As a result, y i is going to be y minus 1 to the power i multiplied by y 0; so, it is
going to be y 0, negative y 0, y 0, negative y 0, it will keep switching between those
values that is exactly, what we had observed through our numerical simulations in the previous
lecture; and when h when above the value of 4, for example, when we took h equal to 5,
what we observed was that the Euler's explicit method was unstable.
So, that summarizes the stability results for Euler's explicit method; similar stability
results are applicable for most other explicit methods - all for all other explicit methods
- also; explicit methods have a limit on the step the size h that you can use before the
system becomes unstable; if you take a very large step size using the explicit method,
we will have the system becoming unstable. In case of RK 2 method, the stability envelope
for explicit RK 2 methods is similar to that we see for the Euler's method. When we go
to RK 3 method, the stability envelope is much larger compare to the Euler's explicit
method. In RK 4 method, the stability envelope is still larger.
However, in all these cases, there is a maximum value of h, which depends on the value of
lambda; if you go to a value of h, which goes beyond that particular value, we will get
this behavior of the numerical solution to be unstable; unstable behavior of numerical
solution has to be avoided at all cost for a system that we know to be stable.
Now, this is where the implicit methods become useful. So, let us go back to our derivation
that we had done for the Euler's explicit method; so, this is the derivation for the
Euler's explicit method; I will repeat the derivation for Euler's implicit method. For
Euler's implicit method, we had y i plus 1 equal to y i plus h times f of y i plus 1
and f of y i plus 1 is nothing but negative lambda multiplied by y i plus 1.
So, we will go over here and write that down and we will take this term on to the right
hand side and we will get one plus lambda h multiplied by y i plus 1 equal to y i, which
is going to give us y i plus 1 equal to 1 divided by 1 plus lambda h multiplied by y
i. So, this is the value of alpha. Now, we know
that, for positive values of lambda, the analytical solution is stable. So, we are only going
to consider positive values of lambda, the step sizes always positive; so, no matter
what value of h we choose, any non-zero value of h that we choose over here, alpha value
is definitely going to be 1 divided by 1 plus some positive number, which means that alpha
is always going to lie between 0 and 1. In this particular case, for lambda greater than
0, alpha is guaranteed to lie between 0 and 1; no matter what value of h we choose, we
know for sure that alpha will lie between 0 and 1, because alpha will lie between 0
and 1; this method, an implicit Euler's method is, what we call as globally stable.
So, no matter what value of h you choose, no matter how large a value of h we choose,
we will always get an the implicit Euler's method to be stable; stability does not say
anything about accuracy, so this stability property is the reason y the implicit methods
are actually used. Now, what I have done so far is, I have talked
about the various Euler's method, Runge-Kutta method and so on. We have discussed about
accuracy of those methods; we have discuss about stability of those methods; we have
discussed then predictor-corrector methods, what we mean by predictor corrector method;
we have also look at non-self starting methods.
So, this is what we have done so far. All of this discussion involved single variable
problems, we had dy by dt equal to function f and that particular function was a single
value function.
Now, what we are going to do is to look at the extension to the multi-variable case.
So, consider that, we have the equations of the form dx by dt equal to f 1 x dy by dt
equal to f 2 x and so on and so forth. Let us just consider three variables x y z
as the three dependent variable and t; and we require the initial conditions, an initial
conditions means x y and z have to be defined at the same value of t 0; so, let say, x at
t 0 equal to x 0, y at t 0 y 0, and z at t 0 equal to some z 0.
Now, let us consider Euler's explicit method; and in Euler's explicit method what we are
going to do is, x at i plus 1, so x i plus 1 minus x i, so dx by dt we can write it as
x i plus 1 minus x i divided by h that is going to be equal to f 1 of x i y i z i t
i; and then we multiplied by h throughout and then we take x i on to the other side
and the equation is simply going to be x i plus 1 equal to x i plus h times f 1 of x
i y i z i t i; likewise, we can write for y i plus 1 as well, y i plus 1 equal to sorry
y i plus h times f 2 of x i y i z i t I; and z i plus 1 equal to z i plus h times f 3 of
x i y i z i t i. So, this is the overall expression that we
get. Now, let us writes this in a vector notation. In vector notation, this we will choose our
vector capital x as nothing but x y and z and our vector capital f as nothing but f 1 f 2 and f 3.
So, this we can be written as, x i plus 1 y i plus 1 z i plus 1 is going to be equal
to x i y i z i plus h times f 1 f 2 f 3 computed at time i.
So, this guy is nothing but our capital x i plus 1; this guy is nothing but our capital
x i; and this guy is nothing but our capital f computed at x i, t i; and this is our multivariable
Euler's equation.
So, in the vector notation, the equation that we need to solve is going to be, d capital
X by dt equal to some capital F of capital x, t, where capital X is in general is a vector,
capital F is also a vector. The Euler's explicit expression is going to be, X i plus 1 equal to X i plus
h times F of X i comma t i; this is going to be Euler's explicit. Euler's implicit will be, X i plus 1 equal to X i plus h multiplied
by F of X i plus 1, t i plus 1; this becomes the Euler's implicit method.
The same idea we will extended to the Runge-Kutta second order method. Let us talk about, since
we have been doing Heun's method in this particular lecture, we will talk about Heun's method;
in Heun's method our K 1 capital K 1 of x, or we just write this as capital K 1 is going
to be equal to f 1 of x i, t i f 2 of x i, t i and so on.
So, capital K 1 we can write as equal to capital F of X i t i; then we can write our K 2 or
before writing K 2 we will write our Y bar or X bar sorry X bar 0 i plus 1, we will write
that as nothing but X bar sorry X i plus h multiplied by K 1, that becomes our next expression;
our K 2 is going to be nothing but function f computed at X bar i plus 1 0, t i plus 1;
and our final expression for Heun's method is going to be X i plus 1 is equal to X i
plus h by 2 multiplied by F of X i, t i plus F of X or I will just write this as equal
to K 1 and k 2, is going to be just h by 2 multiplied by K 1 plus K 2.
So, this is our equation for the Heun's method; where K 1 is computed as F of X i, t i; K
2 is computed as f of X i plus 1 bar t i plus 1; and X i plus 1 is X i plus h by 2 multiplied
by sum of K 1 plus K 2, where each of these capital letters are not scalars, but they
are vectors of the same dimension. So, that is where we finish our this lecture
of this module - of the 7th module. In the next lecture, what we are going to do is,
we will re visit the RK method and we will talk about what is known as the adoptive step
sizing. And finally, we will just introduce our self to the Adam-Moulton and Adam-Bashforth
family of methods, that is what we will do in the next lecture; and the final lecture
I will finish off with the Adam-Moulton and Adam-Bashforth family of methods and do on
overall recap of this particular module. So, I will see you in the next lecture. Thanks.