Mod-01 Lec-03 Foundation of Scientific Computing-03

Mod-01 Lec-03 Foundation of Scientific Computing-03


On this third lecture of our meeting today,
we are going to discuss what constitutes physically a boundary layer and we will talk about the
various methods that have been used in studying this area. One of it is due to Wentzel-Kramer-Brillouin
and Jeffrys method which is again applicable only for linear systems.
So, leaving that aside we will basically talk about similarity transformation, that will
allow us to take the PDE directly to an ordinary differential equation and to solve this ordinary
differential equation we will require auxiliary conditions; the auxiliary conditions consists
of the initial and boundary condition. Having talked about this, we will briefly discuss
how we actually use spatial discretization method in solving this ordinary differential
equation. We also touchup on the time discretization
and show how it is different from space discretization, because of the requirement of causality, which
is distinct from the concept of upwinding that can be used in spatial discretization.
When it comes to this solution of ODEs, we classify the solution methods in terms of
whether we take one-step at a time or we take multiple steps and in doing
so we actually numerically introduce what is called the amplification factors.
And having set the goal here, we try to distinguish between multi step and single step method:
while single step method only gives rise to physical modes, multi step methods have this
additional problem of having spurious modes. We will be talking about that.
Having described the requirement of a single step method, we will talk about solution methods
for initial value problems and in this particular activity, we distinguish between explicit
and implicit methods. In resolving this flow or the problems, we talk about truncation
errors, the orders of method which defines these methods and as a follow up we will specifically
focus our attention on Euler method to describe what it is; that will be followed by two stage
or multi stage two step Runge Kutta methods. We did talk about boundary layers, right.
I think one of you after the class came and asked me that look I have not done any fluid
mechanics, is that going to be an impediment? My answer to you would be not at all; there
are various ways of learning subjects. What I prefer is associative learning. If
you have an anchor on a particular subject, you have understanding on some topic, what
I wish you to also probably look at the positive aspect of that way. Suppose, I know something
then based on that knowledge I go to something which I can build upon.
So, that is why I started with two problems: one was from mechanics and one was from fluid
mechanics to give you a glimpse of what boundary layer does. If you want to have a sort of an impersonal
mathematical way of stating the same thing, I would say boundary layer is a study of differential
equations for say, some variable y of x where the highest derivative is multiplied by a
small parameter; that in a sense we noted that a thin layer of pumps in the solutions
space. Suppose, if I am plotting say y on this axis and x on this axis like the way
we saw the velocity profile. What we found was that we have a thin layer over which you
see a rapid growth and a good computing would require you to be able to resolve such rapid
variations. So, this is one of the issues that come about in talking about boundary
layers. As I told you that it got started in the early
part of 20th century by the group at Gottingen, but people have been also looking at associated
problems notably from the point of view of the same issue where a differential equation
that has derivative term is always multiplied by small parameter. So, in one case you can get boundary layer’s
theory. There is the other method which has been looked up for a long time which is now
known by this name Wentzel-Kramers-Brillouin. Some people also add Jeffrey’s name to it;
it is called WKBJ method. This also comes under the same class of problem
where the differential equation has this kind of attribute.
The only difference here between these two cases is that in this case delta goes to 0.
I call this the thickness of the boundary layer delta goes to 0 as epsilon goes to 0
where as in WKBJ method delta remains finite as epsilon goes to 0. This is the first difference
between these two methods; both are global analysis methods. The second difference of course, is the very
relevant one that this only works for linear equations. You can study linear equations.
Well, but here boundary layer theory is little more versatile; it can also tackle non-linear
problems. In fact, you noted that the boundary layer equation that we solved was a non-linear
PDE and we could tackle it with quite ease. Despite that, you may note that WKBJ method
gets a little more sort of exposure in the literature because those three physicists
used it for studying quantum mechanical systems. So, I do not wish to or sadly with equations
originating from Schrodinger equation and talk about boundary layers and WKBJ method.
That is why, that was the reason that I took up something which engineers can visualize
a boundary layer growing in a fluid flow or you could see that you have early transient
problem in that mechanics problem, the spring mass system we saw that to be able to resolve
those early transient you have to be better equipped. Let us get back to where we stopped. We saw
that we could tackle a non-linear problem and we saw this particular attribute that
as the parameter nu, kinematic viscosity in that problem played the role of epsilon. So,
if nu becomes smaller and smaller or what some of you would know that corresponding
non dimensional parameter called the Reynolds numbers goes higher and higher; you get to
see this boundary layer shrinking – the thickness shrinking.
We saw that based on that observation, Prandlt and his group including Blasius, they exploited
the nature of the solution and Blasius went one step ahead and introduced similarity transformation.
Similarity transformation as I pointed out to you is a very versatile technique for converting
PDEs into ODEs. That is where we had landed up with this equation 16. So, that was the
converted boundary layer equation in terms of that non dimensional stream function f. psi was written down like nu U infinity x
f of eta and this is your governing equation that you had obtained. f d 2 f by d eta square
plus 2 times the third derivative of f with respect to eta that is equal to 0.
We started talking about auxiliary conditions. Auxiliary conditions refer to both initial
conditions. Well, it could be a single condition or could be multiple and then you could have
boundary conditions. The appellation really tells you that initial condition relates to
time variation. We talk about a problem; what are the conditions at t equal to 0. That is
why we relate it to the initial conditions. However, you could also use with equal felicity,
if you have the ability to solve the initial value problems using multiple initial conditions.
For example, this equation that we are noticing here has the order very clearly visible here.
It is a third order equation; so, we will require three conditions and those are the
ones we are calling as auxiliary conditions. If the stream function is given by psi here
then of course, we need to know what we are studying. If you recall that we were looking
at this problem. We have a simple geometry sharp leading edge flat plate exposed to a
uniform flow which is indicated like this and that forms the boundary layer like this
at a later stations. If this is U infinity, on the outside you have U infinity and on
the surface you have no slip condition 0 velocity. That is what you are noticing that u and v
are the velocity components given here. and these velocity components If I call this axis
as y or eta and this as the x, the similarity solution allows us to do away with the x and
y variation separately. Instead we look at it in terms of eta and what we notice that
y equal to 0 essentially corresponds to eta equal to 0 and at y equal to 0, we require
velocity to be 0. This in turn will tell you that if u is 0, the first derivative of f
that should be equal to 0. What about the other condition? If v is 0 of course, you
can see that f also has to be equal to 0. So, this is the usual easier way to visualize
that two conditions are fixed here. Unfortunately, you cannot have more conditions at this point.
If we would have been able to generate all the three conditions here then we would have
obtained what we call as an initial value problem. You can now distinguish how ODEs are viewed
in terms of initial value problem. We will call them as IVP and we have complementary
to it what we will call as boundary value problem. In essence, I just now told you that
this Blasius’s equation that we have written here will not be classified as IVP because
you do not have the all the three conditions available at one of the end of the domain. So, the domain that we will have to be looking
at as I told you the shear layer forms like this. What we do is this delta, the boundary
layer thickness that we talk about has been stretched out by the eta coordinate. That
was the whole thing that we telescope it to a larger distance. Then what happens is we
have two conditions here; we do not have any idea of any other condition that we could
generate from this. However, we also note one of the properties
is that the solution actually smoothly blends to the outer velocity.
So, u starts of from 0. We do not know what is inside it; that is what we are trying to
obtain, but we surely know that it smoothly blends into the outside velocity. That is
what we say here that the third requisite boundary conditions for solving the Blasius
equation is readily obtained from the observation that the boundary layer smoothly merges with
the outer flow. Recall, we talked about inner and the outer
variation – inner solution and outer solution. The boundary layer is the inner part and anything
that is outside of course, we call it the outer flow. That could be obtained by some
other equation by other means which we will not talk about in this course.
Since we noted also the u velocity was obtained as this at the edge of shear layer. So, we
need to fix it arbitrarily at some large value of eta, the edge of shear layer while we may
actually solve this equation not necessarily at the edge of the shear layer delta, but
may be 3 delta 5 delta 10 delta. That is what we would be doing. We will be solving this
equation over a larger value of delta and that edge of the solution domain we are calling
as eta infinity. From here you can see that at eta equal to
eta infinity, we have u is equal to u infinity and that gives you this condition – d f by
d eta equal to 1 as shown in the slide. So, basically this here is the third condition.
Here is an example of a boundary value problem – two conditions given at eta equal to 0 and
one condition given at eta equal to eta infinity. This kind of through an example demonstrates
the distinction between the initial value problem and the boundary value problem stated
like this: that all the auxiliary conditions are given at one end then we have initial
value problem; if the conditions are distributed between the two boundaries then we have a
boundary value problem. All time dependent problems are inherently initial value problems
because we have to give some initial conditions for time dependence. Now having obtained an equation like this,
next thing of interest would be for us to be able to solve this equation. To be able
to solve a differential equation, we have to be able to replace these derivatives by
some simpler information of the dependent variable and one of the easiest ways to do
is use Taylor series expansion. Say for example, what we are showing here
is if we have the knowledge of the function at y equal to x and its derivatives then of
course, we can obtain a solution at a neighbouring point which is h apart. If it is in the plus
direction then of course, we have a simple look of it like this and if it is on the other
side, we will get the odd derivatives with a minus sign; this is quite straight forward.
You all realize that computer does not allow you to directly use these differential or
integral operators. So, we will have to do something and that is what we do – discretize
the problem. Instead of having the problem in terms of the derivative we convert it in
terms of functions at different nodes and that is where Taylor series expansion comes
into picture. What we have shown here for a point at plus
h or minus h, we can do it at points which are at different distance from the point of
interest. Basically, what we are talking about is if we are looking at y, this is x and let
us say this is our domain of interest. We discretize the domain into a network of points
like this and then if I am calling this point as x then this could be say x plus h. We are
talking about uniform spacing, let us say and this could be x minus h and we can similarly
get information at x minus 2 h, x plus 2 h etcetera. Then with the help of those series representation,
we can write the first derivative simply like this. I am sure most of you are familiar with
it I am just giving a quick recap of the same. The first derivative could be written in terms
of the function at x and its right neighbour or the function with its left neighbour or
we could write it like this. We could write it in terms of the right neighbour and left
neighbour divided by twice the gap between these two successive points. Now of course,
you can visualize there would be many ways of evaluating derivatives – even these first
derivatives. They would be far too many. How do you choose?
Of course, we need to come out with some kind of rational, some kind of a guideline as to
how to choose that. One thing you can clearly see in either of these relations if you look
up here, I have gotten this first derivative relation by truncating the series beyond this
term, second derivative onwards and then of course, pulling this on this side and dividing
this by h, I got the expression for dy by dx. This tells you that the Taylor series
has been truncated at the first order itself. So, this is called a first order method.
Same thing if I pull this to the right and get this h dy by dx to the left and then of
course, ignore these set of terms I get the other expression that also happens to be a
first order accurately. Whereas, this expression that you are seeing here taken from these
neighbours on either side I subtract the second from the first then of course, y x cancels
out and I get 2 h dy by dx. Then you also notice one interesting thing
happens that these even derivatives all cancel out. In the process, what you get is this
expression in which you have been able to raise your order of approximation from first
to second order. So, it is quite obvious that if you get a higher order representation,
you have a higher accuracy; this is a legitimate expectation. So, that could be one of the rational that
let me also tell you people have overplayed this obsession with order. We will see that
there are better ways of characterizing various operations and we will learn as we go along.
That was about spatial derivative discretization. When we have time derivatives then we have
to think afresh. The reason is here – causality. Well, this is something which we all accept
at the back of our mind; we assume it to be valid. As I keep telling whenever I get an
opportunity that this is more of a philosophic view point than a mathematical statement.
Nobody has proven causality, but nobody has seen it violated too. There you have it. It
basically tells you that if you are trying to look at events at t at advance level which
we did by index n plus 1 from the knowledge of the system, from the values that we have
up to t n that would be legitimate, but we should not be hoping that somehow future can
influence present that only happens in Hollywood films anyway, but this is a rock solid statement
we have never seen it to be violated. Basically, we need to be aware of this and
let me also tell you that for spatial discretization we have seen that if I am marching from left
to right then I could take a point forward to the right. Fine, that is some kind of an
analogy of going beyond than where we are now. So, that is routinely done for spatial
discretization, but for time discretization you are not allowed to.
However, what you could do is you could write down the expression for the variable associated
with the events at t equal to n plus 1th level to be dependent on what happened in the past
at the nth level, n minus 1th level, n minus second level and so on and so forth.
What happens then that is a very interesting aspect and this is often been overlooked and
this is one of the issue that we are going to talk about in greater detail in future. Say, give a simple example, if I have a simple
equation like this with a first order time derivative written
here. What we could do is we could approximate this first derivative in the left end side
by its approximation at n plus 1th level. As I have told you I could do it in the simplest
possible way would be to relate it to the events that had occurred at just the immediate
predecessor time. If I do that then I could define an amplification
factor which I call it say g of x because I could be looking at a space time dependent
problem. This u could be simultaneously a function of space and time. I could evaluate
it at some point. We will not worry about this. What we will do is we will say that
g of t I will write it like t of n, I will write as u n plus 1 by u n. So, how much the
function has gained over one time step; that is the measure of this amplification factor.
Now, you can very clearly see that this equation if I look at it and divide the equation by
u of n then what I am getting, u n plus 1 by u n that will be g and 1 by delta t will
be L of u n by u n. So, what I find that for a first order equation like this du dt equal
to L u, g happens to be 1 plus delta t by u of n L of u n, that is what I get. Now, for a moment think of if we would have
done it slightly differently. We do not do this; instead, let us do it like this – u
n plus 1 minus u n minus 1 by 2 delta t. Remember the central difference formula that we just
now saw for the spatial operator. We can do the same thing and this does not violate causality
because we are drawing information from n minus 1 and we are going from n to n plus
1. So, this is legitimate; there is not any problem.
Now what happens? If I go through the same exercise, divide this equation on both side
by u of n then what am I going to get? Help me. What I am going to get is u n plus 1 by
u n will be g. So, I will get g then I have u n minus 1 by u n that is 1 upon g. So, what
happens? Yes, tell me.
Sir, the two will not would not be the same gs.
Why not? One is a t n and one is t n minus one
Depends on how the L u n is. You will see in most of the cases if you are not seeing
what we. We will come to that – distinguish system which are called autonomous which are
non autonomous. If you are looking at autonomous system, it will not matter.
We will go in much greater detail. This is just an introductory stage. So, please bear
with me. Consider that all gs are same; these are called stationary processes where events
are same at all time step. They are called stationary processes and you will be seeing
in most of the mathematical physics problem you will come across stationary processes.
It is very rarely you talk about non stationary processes. I might even actually give you
an assignment since you provoked me. We will find out that even for a stationary system,
careless handling of that system numerically can actually get you to what you are dreading.
It will happen and it does happen. So, what we are essentially getting here is
a kind of a quadratic equation. I will just simply write this. What it actually means
that here I just simply add g of 1; only one amplification factor, here I am going to get
1 and 2. I could write it down. I am not bothered about this for a moment.
This is what we have noted here that trying to relate events with multiple time steps
not necessarily at one time step, but more than two levels are involved. See in this
case of course, on the right hand side we have introduced only two time step; here we
are talking about three time steps method. These go on from n minus 1 to n plus 1 and
as a consequence instead of getting 1 g, we are getting 2. Now, which is correct?
Should we have 1 or should we have 2? Suppose, I would have gone on; instead of three time
levels if I would have incorporated four time levels I would have a cubic for g now. So,
I would end up with 3 values of g. Yes.
Sir, I could not see that we are deliberate discretization in terms of time, sir because
u n minus 1 and u n plus 1 are the discretization in terms of the space.
No, it is not space; it is time. We are talking about time; that is a time derivative du dt.
So, I have not written; for a moment forgot about space. Let us say this is a node e and
time is the independent variable and n is the time index. I posed a question to you
that more than two time level method gives rise to multiple values of g and this seems
to be somewhat dependent upon what is the number of levels that we invoke.
So, can someone tell me I mean what is the correct thing? I have already written it down;
you have read it. All that I am saying is that this two time level method is the correct
one, the real one. Why? It is because whatever you are doing your differential equations
are limit processes – delta t going to 0. So, you should involve as a kind of a local
or instantaneous representation of the derivative. More defused your information is, you start
involving more time step from farther and farther backward. That is the artifice of
your numerical method; that cannot be physics. That is why this one, g 1, I will call it
as the physical mode. Here I may have one physical mode, we will see we will go much
deeper into it as we go along and in this case mode will be basically numerical mode
and no doubt about it that this is spurious. So, please be attentive to whenever you are
doing temporal discretization that when you involve more than two time levels, those methods
are called multi step methods. They come with the baggage of generating spurious modes.
If you are not careful, these spurious modes can play havoc. In fact, the subject of chaos
dynamics is probably built around, fooling around wrong numerical method.
In many of the cases, you will see that if you are not careful about chaotic dynamical
system and if you do discretization mindlessly then you might see totally a different dynamics
than what you would like to see. Whereas, whenever you use only two time level
methods, we will call them as single step methods. Let us keep this in mind and then
we will go along. Let us get to a very simple task of solving initial value problems.
The initial value problem is of course, given by du dt is equal to some forcing term f which
could be a function of t as well as the dependent variable. Now, this we would require. This
is a first order in time equation; so, we require an initial condition. That is the
whole idea. That initial value problem, the first derivative tells you that you need only
one condition and that is provided to you. Well, this method that we are talking about
you can also use it for spatial initial value problems. There is not much to worry about
as I told you that multi step methods are dangerous. We will keep away for the time
being. We will focus simply on single step methods; they help you in avoiding trouble.
So, if I am trying to solve that equation du dt equal to f. This is the generic equation
by which we will be solving it. u at n plus one is equal to the predecessor at u of n
plus h, h is the time step. Time is a function. This function we will have to work upon. That
function will depend on the time levels which are involved here for the single step method.
They are t n and t n plus 1, the function at the n plus 1th level, the function at the
predecessor level and the time step level itself.
This phi or variable phi – var phi, what you call this? It is called the increment function.
So, all numerical methods revolves around figuring out this increment function and we
can see that this single step method which has a generic appearance as given here can
be further sub divided into two. If the right hand side is independent of u
n plus 1 that means we have a what we call as a explicit method because you can explicitly
calculate the right hand side and then you can march from u n to u n plus 1. In contrast, if the right hand side is also
a function of u n plus 1 then you have to worry about solving an implicit coupled equation
that involves both the left and right hand side. What happens is explicit methods are
simple to execute computationally; implicit methods are difficult. Why should we bargain
for difficulty? The reason is the explicit method comes with severe restrictions whereas,
implicit method allows you to circumvent these restrictions.
For example, one of the restrictions for explicit method is the time step h; you have to take
it very small. We will see whereas, implicit method would allows you to take much larger
time step and when we are talking much larger means we are talking about order of magnitude
improvement. So, that is the kind of payback that you can have that should give you enough
incentive to really probe into implicit methods. However, let us keep our attention focused
to explicit time integration methods because we are also going to show as we go along that
while numerical restrictions are released for implicit method, when it comes to the
genuine accuracy of the solution explicit methods are better.
So, you pay for it by large number of time steps, but you get better quality result.
That is why I would tell you about explicit method first, later on we will go to implicit
method. So, this is the way that we write numerically.
If I represent the exact value with time indicated inside the bracket with indices like t n and
t n plus 1, this is representation of the exact value of u.
What we have written is that we could evaluate that at t n plus 1 in terms of t n and then
with this increment function and of course, we have in the process left some terms behind
which we will be calling as the truncation error.
So, truncation error is nothing, but pulling this first term on the left hand side, that
is u t n. What happens is we talked about order of the method when we were talking about
spatial derivative. We saw first order method; we saw second order method. So, generic expression
for the order of the method is you look at the truncation error term. For this first
derivative, you divide it by 1 upon h and then you will see what is the order of the
polynomial that you are getting – the terms that you have left behind.
If that is the leading term h to the power p then I will call that method as pth order
method. This is what we do. Now, we wrote the numerical equivalent and
in terms of Taylor series this is the story. u of t n plus 1 written in terms u of t n
plus h u prime at k. Everything is evaluated at t n and this is the Taylor series. We are
stopping at the pth time. So, it is a pth order method.
Anything that is left behind, we can use mean value theorem and club it all in the p plus
1th order time. Please note the argument. This is somewhere in between t n and t n plus
1. That is what that theta is doing; theta is restricting you to remain in between the
current time and the next time level t n plus 1.
If I now compare with the previous page then I can simply see for a pth order method, this
increment function is given in terms of this set of terms. It is almost like a Taylor series,
but note that u of t n is missing. Anything beyond that, u of t n comes under that increment
function term. So, we have now an estimate for writing the
exact solution at n plus 1th level in terms of the value at nth time level and increment
function and plus the truncation error term. The truncation error term is readily identified
by this. When I say I have a pth order term, the leading
truncation error term will be p plus 1th term; that is h to the power. See this is where
the order actually fools us because we are simply obsessed with the order of h, but who
gives you a guarantee about this derivative of the dependent variable. See this is the
p plus 1th derivative of u. There is no way we could guarantee that successive
derivatives are going to become smaller and smaller. In fact in many physical systems
you will see that it just happens the other way. I don’t know about Quite a few of you
must have taken this other course on numerical methods in the second year level.
There you must have come across a topic called interpolation and you must have also been
told that you are trying to interpolate a function. It is always better to restrict
yourself to a lower order than to higher order. Am I making sense? No, then we will come back
to it in future. Well, you have a noisy data and you are trying
to fit it with some polynomial. It is always better to stay as low an ordered system as
possible because any small difference between the physical state and the measurement actually
amplifies more, if you try to fit it with the higher order polynomial.
Because This error if it is localized, actually excites the higher wave number or higher frequency
and that is the reason that you should be well advised to not only pay attention to
what this order of h is, but what the associated derivative is.
In fact, we will figure out a method; we have developed a theory over here over the last
5-6 years. We look at in the case space, the wave number circular frequency space. We will
talk about it as we go along and then that would all make sense. We will come back to
it again and again. Let us look at few single step methods. The
easiest one is in fact, what I have written it out here p equal to 1 and then what you
get u at t n plus 1. I am just simply relating with u of t n and the right hand side function
has been evaluated here like this. So, that is what we call as Euler method.
Of course, suggested by Euler and you would note that Euler method looks easy, but very
nasty method. Many times it will lead to numerical problems and you practically will be advised
to stay away from it, if you are just simply looking at solving ordinary differential equation.
So, what happens? You can raise the order of your method. You notice that here we needed
the first derivative in Euler method. If I want to go to higher order method, say second
order method then I will recall u double prime and since u prime is f, I could differentiate
it and I could get u double prime which is f of t plus f into f u. If I want to go to
a third order method, I will have to do a little extra work and if I want to go fourth
order method, I will have to be looking at more details.
So, this is what you have to do at the coding level. You will have to be evaluating them
and you will have to be plugging them in, write your code and you have high order methods. Now, there could be a practical method to
fix the order of the method. If I look at the truncation error term and I say, look
I am stopping at the pth order and I do not wish to cross this threshold, the error tolerance
epsilon then what happens? This is your leading truncation error term that should be bounded
by epsilon. If I do, I have to estimate this quantity
within this modular sign – mod sign. What happens? After all, I do not know what this
theta is; so, it becomes a difficult task. Instead, what you would be doing? You would
be looking at the right hand side function, its pth derivative and you find out in that
range from t n to t n plus 1. What is the maximum of that function? Well easily said
than done, but you would be happy to note that solving ODE is quite common place, quite
routine. So, if you go to CCE and look at any of You
may have already used it. Say, go to any of these packaged routines available – libraries
they would all have that and they would without asking you, do this exercise for yourself,
but we are not talking about how to use packages; we try to find out why the packages work the
way they are designed. I make an observation here that Euler method
is a restrictive method all though not completely prohibitive. Will again spend some time talking
about what exactly you mean by it. Every sentence is loaded; it would come with lot of explanation
as we go along. We will be talking about it when we come to space time dependent problem
later. Now, we noted our restriction that we would
not like to go beyond single step method, but then we find that the rudimentary method
– the Euler method, is difficult to handle. Higher order methods have the disadvantage
that you will have to evaluate all these partial derivatives not only in terms of coding difficulty,
but computing time also will increase because you will have to be evaluating all those derivatives
– high derivative terms, numerically and they would involve lot of extra work.
What happens is these two mathematicians Runge and Kutta developed simultaneously what is
called as single step method, but within that single step you break the domain further into
smaller stages. That is why the word multi stage comes in
here. So, we are talking about single step, but within that single step we will be stopping
by at various stations and different stages. This was a very brilliant way of conceiving
how we can address this issue of accuracy and not bringing in spurious modes. You can
see that right. We are still keeping ourselves routed to single step method, but by adopting
multi staging we should be able to improve the accuracy; that is the whole idea about
this method. If I am trying to evaluate this equation advanced
by one time step, this is what we do applying mean value theorem we write the expression
like this. This f, I need to know at some intermediate
state which I identify with theta. There are various possibilities of choosing the value
of theta. After having chosen the value of theta, I could also evaluate this quantity
f by multiple options and that is where you have the whole basket of different methods.
So, it is not that you would be talking about a unique method here. You would have all kinds
of freedom to look around and shop around. For example, we have already seen that if
we choose theta equal to 0 that is your Euler method. Basically, what we have done. Let
us try to look at the solution and try to see what we are aiming actually. Let us say I am plotting here u versus t and
I have a solution at t n and I am trying to go t n plus 1 and this is say the exact solution.
What we did if you can see here this increment function that we have evaluated is nothing
but, the slope of the curve, is not it? Where have we evaluated the slope? At the starting
gate. At this point if I draw a slope, Euler method basically gives me the solution here
and that is where we have so much of problem, so much of error if the function is rapidly
wearing and the slope is not truly representative of what the function is doing. In the single
time step, you will always end up with a distant point from the exact solution and that is
one of the issues. Now, Euler also came out with this method
that instead of calculating the slope at t n, you calculate the slope at t n plus 1.
What happens? I could draw a slope here. I will start from here; so, what I am doing
is I am drawing a parallel to this line which is the basically the tangent here. What I
am going to get is something like this. So, this line is parallel to this. Well, my
drawing is not at all good, but I should be mindful of that and let us say we do that.
That is your solution. So, if this is your basic Euler method, backward Euler method
would get you there, but this just happens to be the way we have drawn the variation.
In another case, this may not be as good as what it has been seen here. So, you can see
that these two possibilities are not really at all that good that you also notice that
the right hand side involves t n plus 1. That means, you will have to know the solution
to evaluate the slope. If the function descriptor is tough then you will have a tough implicit
method at hand because u n plus 1 is involved on both sides. So, the other possibility is suppose a query
we fix the mid stage at half way. I choose theta equal to half then you evaluate at t
n plus 1 and it would be obtained like this. Now, you notice that this is not a node and
we could again have the possibility of choosing our self.
For example, for this we can look at the followings of subcase. If this is evaluated by the Euler
method, that means what I am saying is I have to work out what this u at t n plus h by 2
is. What I would do? I will evaluate it as u of t n plus h by 2 into f evaluated at t
n. Then what happens is wherever I have this
function, I will write down as u of t n plus h by 2 f of t n.
In an algorithmic sense, what I would do? I would first evaluate K 1; K 1 is nothing,
but h times f evaluated at t n. Then having obtained K 1, I will evaluate the right hand
side with the help of the knowledge that the mid station, the solution is obtained as u
of t n plus h by 2 K 1. Once I have the K 2 evaluated in this fashion, I will just simply
add that K 2 to the previous time level to get the new solution.
So, this is a very simple demonstration of how a multi stage method works, if we decide
to clean our information from midpoint of the stage. Now, we could also do various other things.
For example, I could have suggested that look this quantity that we have, we could wherever
we need to have obtained this average, we will just take it as a half of what we have
at t n and what we have at the end of the… Basically, then this will be written in terms
of f of t n plus this and eventually plug that information in there and there you have
it. Then we will be evaluating K 1 in this case again the same way, but K 2 has been
evaluated in this fashion and the difference comes from the previous method is the final
step here. I just simply take the average of K 1 and K 2.
In the previous case, I just simply obtained it in terms of K 2 remember u of t n plus
K 2, but in this case what you are doing you are taking an average of K 1 and K 2. That
is a method attributed to Euler and Cauchy. So, call that as Euler-Cauchy method.
Yes. what is the second order equation
Yes. so I have Starting point was a first order method we are trying to get it better.
So, second order method is a better method than first order method; that is the whole
approach. Let us try to recap what we just now looked
at. If we have to evaluate then we have to see the suitability, applicability of the
method from many possibilities. A natural way of course, would be to develop some analysis
tools with the help of some standard benchmark problems.
That really mimics the most physical processes involved in the differential equation. We
hope to we will of course, we will address this issues as we go along.
Let us now look at solution methods as given by Runge and Kutta and that is where we will
begin from next Monday. Thank you.

One thought to “Mod-01 Lec-03 Foundation of Scientific Computing-03”

Leave a Reply

Your email address will not be published. Required fields are marked *