Tip:
Highlight text to annotate it
X
We begin lecture 30 with a brief review of what we have talked in the last class, namely,
the sixth order scheme of Lele that was obtained by Taylor series analysis. However, we also
pointed out that there must be some basis for choosing compact scheme. One way of doing
it was an optimization of the scheme by looking at the truncation error; as we talked about
in the last class about Haras and Ta'asans method of looking at it in the k space instead
of worrying about the physical plane matching term by term.
In the process, we can develop very accurate schemes but, whose order could be very less.
That naturally brings us to this discussion on the distinction between orders of approximation
versus resolution of the scheme. We will surely highlight what is more important is the resolution
of the scheme rather than formal order of approximation.
In the process, we also figured out why compact schemes are so accurate as compared to other
methods because of its global nature. That global nature can be found out very clearly
by an analysis that we have developed - spectral analysis tool. We can figure out that this
implicit method when written as an equivalent explicit scheme that is equivalent to 20th
or higher order accurate scheme. So, this is one of the plus points of compact schemes
as a global approximation. We have further developed the spectral analysis
tool by finding out the resolution at any node for a non-periodic problem in terms of
a projection operator that is what we are going to talk about. Once, we do that taking
the derivative is equivalent to converting the k, a wave number into an equivalent k
which we have called as k eq, which will have a real and imaginary part. The real part actually
represents the phase while the imaginary part represents dissipation or anti-diffusion.
This is rather important, because we will just show using one of the older schemes due
to Adams - near the boundary; this scheme actually shows very large dissipation near
the outflow and instability near the inflow. Then, we will also talk about carpenter scheme
of fourth order accuracy. We noticed that this method also has a problem with the boundary
closure; this is essentially due to the asymmetry of all compact schemes.
So, we have to be extremely careful in how we get this boundary stencil in the compact
scheme. This has led to opponent compact scheme which we will be talking about, but we will
in the meanwhile talk about evaluating second derivative as a sequence of two first derivative
operations and with this we will conclude our lecture 30.
We were looking at high accuracy methods of computing. In that context, we were talking
about compact schemes. Compact schemes typical example is given here in equation 9, where
the derivatives determined by these prime quantities on the left hand side are related
to the function on the right hand side. You can see that these are implicit methods
because the derivatives that you are looking for they are coupled in this equation, as
opposed to explicit methods that we have done before. So, there are some constants alpha
on the left hand side, a and b on the right hand side, h is the spacing between the nodes,
then what you can do? You can equate the coefficients of Taylor series on both sides.
When you equate the coefficients of u prime on the left hand side of course, you note,
it gives you 1 plus 2 alpha and on the right hand side, it is a plus b. So, next coefficient
is the third derivative term coefficient. In the left hand side, you get alpha and on
the right hand side you get a plus 4b by 6 and so on so forth. What happens is, if 3
coefficients are found out by solving let us say, these 3 equations 10 to 12 and that
is how you get 3 values of these constants.. This is what we usually would get; these are
unique set of values of constants for these particular case. Since, you have satisfied
this Taylor series expansion 9 up to derivative so, next drop term is the seventh derivative
that is why we called it as a sixth order scheme. I emphasize the fact that you would
always be satisfying equation 10 and that is the coefficient of u prime, because that
is what you are trying to find out. Since you are looking for u prime, there is no way
you can forget equation 10.
The satisfaction of 10 is a must and this is what we called as the consistency condition.
This was a unique scheme given here alpha a and b. We talked about the leading truncation
error; it is the seventh derivative term; that is why it is called sixth order scheme.
We talked about Lele's effort in trying to develop a fourth order scheme. Fourth order
scheme means, you would have satisfied 10 and 11 and omitted 12; that would give you
2 equations and 3 unknowns. So, you get one-parameter family scheme; let us take alpha be that parameter,
then you can get a and b in terms of alpha. We talk about this and we said to look these
coefficients that we are seeing; they are all real numbers however, there are to be
a sort of a scheme.
In the sense, they are not like direct functions that only at the chosen value you have defined
and everywhere the error is degraded. That is the question that posed here, whether they
are discontinuous function or continuous function. We did come later and discuss further more
saying that they are indeed continuous functions and that opens up possibility of obtaining
optimized schemes, you can optimize it. That is what Haras and Ta'asan did. They did optimize
scheme for evaluating the first derivative by minimizing the error committed by this
approach with respect to the spectral method. They tried to look at the full range of k
for the wave equation only and that was somewhat different than what Lele tried to do earlier.
He also tried his hand at optimization because they were 3 equations, 3 unknowns; he looked
at the departure at 3 values of k and from there you obtain the values.
Lele observed which the core of this subject is, if you take a look at these optimized
or non-optimized implicit schemes; they are far more accurate than their counter part
in the explicit schemes, this could be quiet significant. For example, a fourth order scheme
that Lele obtained was out performing a tenth order explicit scheme. So, basically those
made us observe that we really should not be bothered too much about looking at the
truncation error. Instead, we should try to find out in the k space what is happening
to error that is the whole rule of this game.
Basically, we define what we call as the resolution. Resolution is of course, what we are of, we
are not worried so much about order of approximation; resolution means, how well you are resolving
all length scale in the problems. You may have noticed that we have already talked about
the Nyquist criteria the chosen grade that is one criterion.
But, the resolution is not the Nyquist criteria, it is something looked in that and that is
what we need to figure out. We basically, would do that by looking at spectral analysis
in the k space and that is what we did subsequently. This is where you begin by expressing the
unknown let us say, at the jth node in terms of it is Fourier Laplace transform; please
understand that this is an integral, this is not a series. So, this implies that you
are able to take care of non-periodic problem; you are necessarily solving a periodic problem.
If you recall that in general, what we are doing? We are looking at a linear algebraic
equation where the derivatives let us say, the first derivative is obtained by solving
this linear algebraic equation . Most of the time you would note that, this A and B matrices
are going to be constant value matrices; they do not depend on at what level of computation
stage you are in, you can do that. So, what you can do is you can write it down instead
u prime should be equal to C u, where C is A inverse B, but please pay attention that
in your actual computation you would never do this, this A inverse operation.
What you just now noticed that A matrix - as you can see - from this say equation on 9
would be the tri-diagonal matrix. Even then, you would not be doing A inverse of B; that
is just for the sake of analysis we are writing A inverse B.
In the actual case, you will always solving this equation because this will involve, so
this is a tri-diagonal matrix here; this could be at the most a penta diagonal matrix. What
you would notice? That when you write down like this u prime equal to A inverse B into
u or what will be the bandwidth of the matrix?
This is interesting, because you see here we are solving an equation of course, on the
right hand side, you do not have to worry about whether it is penta diagonal or hepta
diagonal, because your fraction values will be given.
So, you will make this product and it will be just a vector. This will be a vector that
would be computing and then you will be solving a tri-diagonal matrix function; that is the
whole approach that you would be taking. If I know the dimension or the rank of this matrix
is N then we do 5N to 7N calculations in solving this equation. If you do A inverse of B what
do you think the bandwidth scheme matrix would be?
Well, if you go ahead and perform that you would be surprised to see that C matrix will
have non-zero entries about the diagonal, which spans about 10 or 11 points on either
side. What does that mean? It means that C matrix is something like a banded matrix with
22 non-zero elements. That would be something like your twenty second order or twenty third
order of equivalent explicit scheme. We discussed it in the last class also. Let
us say, for example, a second order scheme we require 3 points; the middle point and
that neighbors, the fourth order we need 5 points. If I have 22 points, the band-width
of C matrix is 22 then, I am actually equivalently writing out about let us say a twenty third
order scheme. So, this is the clinching issue of compact scheme that although you are solving
a tri-diagonal system, you are developing an equivalent exclusive scheme, which is of
the order of 20 plus which is a significant advantage that you should have. This is the
prime mover for this whole activity. Now, you also notice that if we are looking
at explicit scheme, what is A? A is going to be an identity matrix, because we are explicitly
finding the derivative at 1 point in terms of the function values. So, A will be just
simply nothing but, an identity matrix that is what we said. Basically, what we are trying
to do? We are trying to develop a framework, where we can analyze a scheme.
It is not necessarily has to be finite difference, finite volume or finite element all you need
to do is come to up to this level A and B time permitting; I will just tell you how
to construct finite element methods in this course or whatever, little time we would have.
There we would note that this A and B matrix once again are going to be some constant valued
matrices. This framework that we are setting up is going to be rather numeric and you could
use it for any discretization scheme that you would like to.
Now, one thing that we did not talk about is this the band-width, which we are now familiar
that C is not necessarily band limited by the band-width of A and B. So, this is a major
issue that I emphasize again. That compact scheme appears local in nature because if
you have looking at j th point on the left hand side, you have at the most j plus 1 and
j minus 1 that makes it local. However, because of this A inverse B operation makes the C
matrix band-width rather large. So that makes it basically a global approximation
and this is what we call as the global near spectral approximation it is so . Now, let
me tell you the compact schemes are so much more handy and accurate that so called spectral
element method an off-shoot of finite element method that often used by many people. Although
they call it a spectral element but, they do use Lagrange interpolation functions of
the order of 3, 4 or 5 at the most. There the accuracy is that much whereas, this provides
you with a much greater accuracy and much more flexibility in solving problems.
So, what we are essentially said in the last bullet is that if we are looking at any generic
compact scheme, we can construct an equivalent explicit scheme that is the connection between
an implicit and an explicit scheme.
For example, say you are trying to obtain the first derivative as we have seen here
on . First derivative is equivalent to multiplying the Fourier Laplace amplitude u of k by simply
i k. Then let us say if we are trying to evaluate I think there is a mistake this thing should
be u l prime.
What we are trying to do? We are trying to find out the derivative at the l th node,
if I want to write this what I should be writing here i k u of k e to the power i k x l d k.
This is the spectral representation of this . However, what we are doing here? We are
writing it like this so, what we are doing? In a sense apriori, we do not know how many
of those functions are involved. We would assume that B is full that means,
we are taking the contributions from each and every node. So, what happens is we could
project this; let us say, I could project it any other point say I could write it like
e to the power i k and then e to the power i k, I will write x l minus x j.
So, this quantity if I call this as some projection of the j th node on to the l th node that
is what we are written here in equation 17. Although by derivative is on the left hand
side it relates to the l th point whereas, I am evaluating the phase at the j th node,
all the j th node. Then what will happen is if I just look at
one particular point x j then I need to multiply by this projection operator that is what we
have written here. So, you can see with a uniform point that will be simply nothing
but, i k h times l minus j. That basically is further split into a real part and the
imaginary part here I call it as r l j and plus i, I would tend to i l j.
If I like to write the derivative at say jth node, I can write it here because on the right
hand side we have contribution coming from all the other nodes, nonetheless we try to
write an equivalent form where the phase is essentially written at the j th node. You
have to understand that whenever we do this derivative operation, these are local operation
means I am trying to find out derivative at the l th node, my phase is dictated upon by
this quantity evaluated at the l th point itself whereas, what we try to do numerically;
we try to add up all the points contribution in the domain.
That is what? These right hand sides apply that if I am trying to find out the j th point
this is not necessarily at the j th point but, this for u 1 to u n all the points are
coming into picture.
That is what we are trying to state here, I have a sort of a global operation in numerical
sense, but if that is to find out what is happening in the corresponding to theoretical
spectral auxiliary mode, then I should look at this because now, what has happened here?
This is like u j prime on the left hand side and this is my C matrix. Let me write this,
so what is this C matrix? Suppose, I want to find out u j prime what I would do? I would
take the j th row and multiply by the whole column.
So that is precisely what we have done here; that is what you are seeing here. I have taken
the j th row and I have multiplied I have done that phase shift business I project each
and every point to the j th point, u 1 has been projected to u j, why? having this and
so on so forth. That is why this quantity is within brace tells you how the projection
has been take with the help of this C matrix. C matrix is also again a constant coefficient
matrix, you do not have to worry about anything these are easily obtainable. So, having obtained
this now what happens if this is my spectral representation. What we have done in a computational
set? We would have done it like this instead of i k, we would have obtained some i k equivalent
we know now how to calculate i k equivalent, I will show you here also.
Then I would have multiplied by u of k and then ideally I would have like to be getting
at the l th node itself, but I am doing it like i k h of I am doing it l minus j; so,
this is what we are getting, this is what we want. What happens is if I look at this,
every u j or u l has been projected to the j th point and that projection operation is
given here . This is what we are getting this is the projection
and then of course, you brought the phase back to j th node that is what we are looking
at. This is essentially than the quantity within that brace is nothing but, i k equivalent,
that is your i k equivalent, so i k equivalent is nothing but, your C j l into p l j and
that you have to sum over all the nodes l running from 1 to n.
So that is what I said that you could i k equivalent at the j th node in terms of the
j th row of the C matrix multiplying by that projection operator p l j. This is the story
that you can now calculate k equivalent and remember this is something that we have achieved
now. We have been able to find out this k equivalent at all the nodes simultaneously.
It is not a local analysis; it is not a global analysis, because all the nodes l equal to
1 to n are involved in evaluating k equivalent at the j th node. This is something that we
did and tell you what this was not very old story, this we did over here it only few years
ago, one of our B tech student Gourav, he joined me in doing this. Before that this
did not exists; so, in a sense I think we did develop a spectral matrix stability analysis
over the whole domain. This is important because local analysis is
good but, in a long run what you need to know is how you are solving the problem in the
full domain. You would know how each alteration at one point can affect the other points;
this is what we need to do. Now the second point that you need to notice that in this
k equivalent that you are obtaining is essentially a complex number with a real part, the real
part will be giving you what? Here also I would be writing some k equivalent.
If I now write that equivalent to consist of a real part then, what happens is of course,
if I am looking at i k equivalent so this part is going to give a phase part that you
can easily skip, the real part actually contributes to phase . However, what does the imaginary
part do? Look at i, i was sitting outside and there is 1 inside so, this will give us
something like e to the power minus k imaginary into x.
Whenever you would estimate this k equivalent the real part would tell you how well the
phase has been represented, while the imaginary part will tell you if we have added numerical
dissipation. See that is why we have been discussing so far about central schemes. We
do not want to needlessly bring in numerical dissipation because chances are there that
will cause the solution to attenuate un-physically. So to avoid that, what we are doing? We are
looking at central scheme, but suppose we have although seen that even though we are
dealing with the central scheme that near the boundary we still have to close the system.
We said that we will have to bring in near boundary stencil, which will only get information
from inside the domain. So, on the left hand side it will come from inside that will go
from right to left and on the other side you would get it from inside towards outside;
so, we have talked about it. Those boundary closures, closing the boundary points would
still bring us to one sided schemes. As you know, you have seen that if I have
this A and B matrix, which are not perfectly central because the first row and last row
we will have to change because of lack of information from outside the domain, we will
have to keep one side a domain analysis. Then what happens? You will get a C matrix, which
will not be symmetric which will be non-symmetric and that in turn will give you a k equivalent
which will have an imaginary part. .
Why because you have seen that k equivalent comes from this C matrix here. If the C matrix
is non-symmetric, then you can be rest assured that it will be non-symmetric and you will
have an imaginary part and that can add to numerical dissipation. If you are not careful
which happened to most of the people before we we found that many of the schemes those
are being popular and used instead of numerical dissipation. They are having just the opposite
effect they were actually making the flow unstable.
The problem where the signals were propagating with the amplitude growing that mean k imaginary
instead of being negative, where actually becoming positive; I will explain again, please
do understand that well. Let me explain with the proper example to understand this role
of real and imaginary part this would be useful. Suppose, let us take our simple example here.
If I am trying to solve this equation, if I add numerical dissipation then, what do
I do? I would add some, let me call this nu . See if nu is going to be a positive quantity
than of dissipating the solution, this is like adding a physical dissipation.
If I have a method where instead of solving this equation I have equivalent equation of
this kind, what am I doing? I am actually adding a dissipative term only when nu is
positive. If nu is negative the way we are doing here, then will be actually, instead
of adding dissipation we will be pumping in energy to the system.
So, this is what? I would call a dissipative method. Let me put this instead if I get like
this a strictly negative quantity then what do I get? This will be not dissipative just
the opposite of that so, people have called it as anti-diffusive method, we are adding
anti-diffusion . If I now go back and look at the methodologies we are talking about
here.
Look at equation 22, if I have k equivalent like this then what am I getting; k equivalent
u of k then e to the power i k equivalent x of j d k, this is what? This is let say
we are finding out del u del x at x j location, that is what we are finding out. Now say k
equivalent has a real and imaginary part itself well, let me write it. This part I am just
leaving it as it is now what happens? This I could write it as i k real u of k e to the
power i k equivalent x j d k and I will get minus k imaginary u of k .
Now you see this equation, if you look at this equation, what we are getting here is
this; del u del t plus this is this part so, this is this part that is C times this and
this is additionally that has been brought in here . So I could take it on the other
side then you can see the role of k imaginary. Now you see there is no i term here like here
there is i term so this is like your second derivative brought in there.
Now, you can see what we are talking about if k imaginary is negative then what will
happen here? If it is positive it will work as a diffusive term, but if it is negative
it will work as anti-diffusive term.
This is something that you have to realize that how important it is for one to obtain
this k equivalent. Then you would like to probably plot it in a non-dimensional form
by scaling it with respect to the corresponding theoretical estimate that is why we will be
plotting k equivalent by k. Of course, by now you are comfortable with the idea that
instead of plotting it across as a function of k we will be plotting it as a function
of k h. That will make your job easier because all the time we will plot between 0 and 5
that is the universal limit said by Nyquist Limit.
So that is what I just restate here: the ratio reflects efficiency of a method with respect
to the best theoretical estimate that you can get through the spectral method. I emphasized
also this path that if you would choose a symmetric stencil for the interior node that
would give raise to k equivalent, which will not have the imaginary part. That would actually
mean that you are developing a non-dissipative scheme that is what the second bullet states.
However, we have also seen that we are forced to take asymmetric stencil for the boundary
and near-boundary points that would make the method either dissipative or non-dissipative
near the boundary. Then this is an implicit threat of compact schemes so, we have to be
varying of it and we have to guard against it.
We will have to figure it out and we have also seen that these are all matrix operation,
even if I do something or say j equal to 1 or 2 it is effect will percolate much more
deep inside the domain; so, this is what it is. Let us look at some of the things, which
people have being doing; this is taken from nick Adams thesis work, he suggested to use
a interior stencil like this. As you can see, this is the on the right hand side, we have
a 5 point molecule like what we have already stated.
This you can do it only from 3 to n minus 2 so, you need to have closures at j equal
to 1 and j equal to 2 and as we kept saying time and again this will have to be one sided.
So, at j equal to 1 it would be add mixture of u 1 prime and u 2 prime on the left hand
side and on the right hand side also it has to be one sided, right; that is what is done,
here you can convince yourself. You must do it yourself, do a Taylor series
expansion and check what is the sort of order of this team? You would be noticing that you
are matching the first and third derivative; this is essentially a fourth derivative stencil
that fourth order stencil that you see in 23. Same about this also for j equal to 2
you take a stencil, well this is somewhat slightly better on the left hand side; again
you have a symmetry though 1 4 1 and same thing about the right hand side also.
You can see there is an anti-symmetry that is what you require 3 and 1 both have amplitude
equal to 1. If the point j equal to 1 that may give rise to some problem whereas, rest
of the point we have a perfectly symmetric stencil on the left hand side that you see
1 3 1. On the right hand side you again get what you want in the anti-symmetric coefficients;
j minus 2 is minus 1, j plus 2 is plus 1, j minus 1 is minus 28, and j plus 1 is plus
28 all are divided by 12 h. So you can see that again we can do a Taylor
series expansion and you would probably note that these are sixth order. Basically, Adams
started looking at these problems using a sixth order interior stencil and fourth order
boundary stencil. What do you do? You have to do a similar thing at j equal to n minus
1 and j equal to n; j equal to n minus 1 would be nothing but, a mirror image of 24, this
will be u n, this will be u n minus 1 and this will be u n minus 2 , whereas what will
happen here? Please note this; on the left hand side the
coefficients remain as they are, on the right hand side you are going to see anti-symmetry.
So, you will write it here for 24, instead of 3 we would be writing u n minus 2 and this
will be u n but, in addition we must have a minus sign. You can convince yourself that
this would be just the anti-symmetric representation of the right hand side of 24.
So, same thing holds for j equal to n that would be a mirror image of j equal to 1; left
hand side coefficients remain as they are, replace u m prime by u n prime and u 2 prime
by u n minus 1 prime. On the right hand side, you will be writing instead of minus 5 u 1
you will be writing plus 5 u n, then instead of writing 4 u 2 you will be writing minus
4 u n minus 1 and then minus u n minus 2. So, this is the way you have the full A and
B matrix set up for you. All the time whenever you need to calculate these derivatives which
I have given on the left hand side, you will be solving this equation by tri-diagonal matrix
algorithm.
So that is the practicality of this issue. I think that is what we explained just there;
that we have near boundary point stencil which are fourth order accurate and the interior
points are sixth order accurate, formal accuracy. This is one of the schemes that were used
by Adams; there is this group from nasal army carpenter. They used this type of scheme and
what I just now said, I just am showing you here. So, what about this scheme 27? 27 is
basically something like what Adams took actually for the interior, sorry; for 24, you see Adams
took this scheme, which is a fourth order accurate scheme but, it is symmetric that
is what I mentioned. In the carpenter scheme the interior nodes
are actually this; the same equations that Adams did, this is for j equal to 2 to n minus
1 you can see they can be applied all the way up to 2 to n minus 1. You need to only
close the system by having a scheme for j equal to 1 and j equal to n. So, this 2 are
written here and what I said about how to write out the anti-symmetric scheme 28 is
just the anti-symmetric of 26 on the right hand side that you can see.
To compare with this fourth order accurate compact scheme, carpenter also looked at explicit
fourth order schemes and this is what is given here. You see when you go to explicit fourth
order scheme here on the right hand side in 31, you can see you need j minus 2 to j plus
2 so, you can write it for 3 2 n minus 2. Even though we are doing explicit method your
work actually increases because you have to add two more additional closure schemes for
j equal to 1 and j equal to 2, they are also written like this.
I suppose this 29 and 30 they are also formally fourth order accurate. The whole idea in this
exercise by carpenter and his colleague were to compare a fourth order explicit versus
fourth order implicit scheme, I will show you the result now.
I think maybe I should first show you the result, and then we will come back here again.
If I look at fourth order implicit scheme what we have here is on this side we have
plotted k equivalent by k, on this side we have plotted k h. You can see we have taken
a domain from j equal to 1 to let us say 31 and we are showing half the number of points,
because it is perfectly symmetric. If you see only one half you can imagine what is
happening on the other half because there will be coincident on the real part.
The left hand side is the real part of k equivalent by k and ideally you would want it to be equal
to 1 to as large a range of k h as possible. That is what we want because in a spectral
method that remains one although we have 2 k h equal to 5 and then falls off at the Nyquist
Limits. We would like to do that but, in this compact schemes you would notice they do remain
one to a pretty large extent and then this starts falling off and this is what we call
a filtering of the solution, you see this is what you may have done.
I do not know if any of you have taken this course on electrical engineering, they talk
about filtering low pass, high pass filter. Basically, we have seen how this filter is
work. You can tune the filter and selectively attenuate or amplify different frequency range.
It is also the same way here, instead of time domain we are looking at in space and instead
of looking at frequency, we are looking at wave number; that is the only essential difference
but, other ideas are exactly same. So what you see that compact schemes or any
discrete scheme as this attribute as we have plotted. We have seen that they have this
attitude approach to make this k equivalent k go to 0 at 5 we have done it. If we have
looked at c d 2 and c d 4 expansions in the previous part this is what happens. However,
you notice one very interesting thing about this at j equal to 1 that is the inflow part
of the domain look at this. It has a very exaggerated overshoot, what
does it mean that some of those intermediate wave numbers are weighted more than they would
have been even for a spectral method. Whether this over emphasizing from k component is
good or not we will talk about it later, we will not appreciate it now. But, we do agree
that this is not physical, physically it should have been 1. Then it should have fallen to
0 and this departure overshoot is significantly more than this filtering attribute of any
discrete scheme. You also look at say imaginary part that is
shown here . We have plotted k equivalent by k imaginary part versus k h; what you notice
that j equal to 1 line actually goes off to exceedingly high value. These positive values
I think here we have taken a sign inbuilt into it that is why all this positive values
imply numerical instability. So, what it actually tells you that in the
process of spatial discretization itself if you are not careful instead of dissipating
the solution, you actually are going to pump in energy. That would be quite catastrophic
as you can see you try to use this and try to solve even this equation, which is what
you are going to do as your next assignment. I suppose, you will get hands on experience
in analyzing schemes and calibrating it with respect to a standard problem like this. You
would be able to see what this is doing, so the first point is very badly behaved as we
can see here, this goes off to a value of k equivalent by k is 5 times.
This is exceedingly high value so, the same thing happens with other nodes also; j equal
to 2 3 and so on so forth. What you are seeing that is what I mentioned a moment ago that
although you are making a mistake at j equal to 1 taking a wrongly directional stencil
that effect is percolating inside. You are seeing that it is effect is felt at j equal
to 2 3 all the way and only when you are substantially inside, then you see it is becoming what you
wanted to be it is a symmetric scheme. You want the imaginary part to be 0 and that
is what you are seeing only near the Well, I would say j equal to may be 13 14 that is
what we have getting, even this value here I think 5 or 6 they are unstable very clearly.
You would also see that some of these things are very interesting that it goes up and down;
so, it shows you that instabilities are very selective in across a small range of k, you
have to do that. This was what you see for the implicit scheme and you can see this k
equivalent by k here is equal to 1 all they up to about 1, above it is even more than
that it is about 1.2 or something.
Now, if you look at the corresponding explicit scheme you can see what if you are getting
the real path is on top and you can see it falls off from a value of about point 6 or
4. Although you are having both the methods at fourth order scheme but, explicit scheme
allows you to maintain accuracy only up to point 5 point 6 whereas, the implicit scheme
allows you to go to twice the range. So, what does it mean? That if you are solving
a problem you should be able to do accurate calculation by taking half the number of point,
if you are adapting the implicit scheme as compared to the explicit scheme. If you have
let us say a 3 dimensional problem and each direction you get a benefit of this kind the
factor of 2 reductions. So, you are actually going to get by with may be about half the
number of points in each direction that amounts to about 8 to 10 times benefit.
In fact we will come to some schemes, which we have developed here itself. We get such
benefit in each direction by a factor of 8 to 10 and you can understand that can translate
into benefit of the order of 500 to 1000 times. That is what we have been doing for last 10
years using pc to compute problems, which others use super computers to do, so it is
possible. So, we have exploited the benefit of implicit schemes; let us go back, this
is what we talked about.
Now, let us talk about another aspect; we have been talking about the first derivative.
First derivative we are evaluating by solving this equation or we are doing this u prime
equal to C u. Suppose, I want to evaluate the second derivative then what I could do?
I could equate the second derivative with the first derivative this equation is that
but, if I write this so if I use this in here I am going to get.
Basically, we have done it twice to get the second derivative. Now, what does it do? That
is what we tried to understand and we tried to do in the initial stages when you are at
it. We tried to evaluate the second derivative by this equation 32 and what you had noticing
here is that this has been done in 2 step, first you solve a prime A u prime equal to
B u and then you again solve it again A u double prime equal to B u prime. So that is
what is implied by equation 32 here. Then let us say simplify the problem very much
we take that say some sort of a stencil, which we had shown here, given here, this is Adam
scheme.
We take n equal to 5 so, what do we do is j equal to 1, we use this equation j equal
2, we use this j equal to 3, we use 25 and for j equal to 4 we write out an expression
which is similar to 24.
For j equal to 5, we write out an equivalent expression of 23. So, if we do this now it
is a nice C matrix is going to be nice 5 by 5 matrix and then I suppose some of you may
be familiar. You can use this symbolic manipulator MATLAB tool box also as if maximize another
one, you can use the symbolic manipulator and you can work out this C square.
See all those schemes that I talked about the Adam scheme; the interior stencil was
sixth order accurate, the boundary closure schemes were fourth order accurate. We use
that combination of a fourth and sixth order and work out the second derivative analytically
using the symbolic manipulator. You can see what this is; this is your c d 2 expression
so, what happens to our obsession with order of accuracy? We used a fourth order and sixth
order and this is an exact operation symbolic manipulator comes back and tells you that
in effect you are doing actually second order accuracy.
So, this is clear evidence that worrying too much about the order of accuracy is not really
something that we like to do. So that is the last comment that I made there.
Now, if I go back to that Adam scheme, the previous example also tells you that sometimes
doing this analysis over unrealistically small number of points can give you such surprises,
see here this could be an attribute of taking only 5 points. Suppose, I would have taken
8 points or 10 points I would probably get something different. So, please do not over
generalize. I think this comes as a warning for your exam question; you always saw a quadratic
and you saw the product term is minus 1 all of you jumped the gun and said.
Product is 1 so, one is more than 1 and another is less than 1. That is not true; we have
all forgotten that both of them should be equal to 1 that is exactly what happens method.
If you would have just done it, you have seen it so, please do not generalize. This is what
we also seen that we try to see the sensitivity with n and done this exercise with the Adams
scheme with 20 points and 30 points and when we saw there is absolutely no difference between
the 2 sets of results. Then go back and adopt one of them; so, we have here on the conservative
side, all the results that you see are going to be with 30 points at least.
If you take any more number of points these are not going to be different. This global
analysis tells you how this k equivalent by k the real path varies with k h for different
nodes. Once again, you can see this first node is always will be here but, this is somewhat
much better than the carpenter's method. You saw the overshoot their reaching up to 3.5
and 4 here it has gone up to less than 2.5. We will talk about the role of j equal to
1 and j equal to n somewhat later but, let us first look at this results so, they do
not look to bad. You see in the interior, you have a sixth order scheme and that is
why you will see that k equivalent by k remains flat although up to 1.5, value is 1. It shows
you that out of the full Nyquist Limit you are able to get almost half the range. I think
this is instructive, because what we are going to do; the spectral methods are always using
uniform spacing. You cannot use spectral method using non-uniform
spacing whereas, this compact scheme we can use it after mapping your governing equation
to a transform plane and in the transform plane you have equi-spaced nodes there you
can use this. So this makes it somewhat of a utility that we will see as we go along.