Tip:
Highlight text to annotate it
X
We have been looking at partial differential equations and the methods of solving these
equations numerically in which we saw one particular case of an elliptic partial differential
equation of the form del square u by del x square plus del square u by del y square equal
to 0 and we solve these equations by writing down the second derivatives as difference
equations and we wrote del square u by del x square as u at x plus delta x minus 2u x
plus u at x minus delta x divided by delta x square and similarly for delta y square
and we found that these equations, we can solve iteratively by writing it in the following
form that is if you substitute these equations for del, the second derivatives in the difference
form in these equations.
The equation which we have if you have if you substitute these terms in to that and
this in the case where delta x and delta y are the same. So when delta x equal to delta
y we can write these equations as, so that is what we saw in the last lecture that I
can substitute this difference equation into this and then write this equation, this particular
equation in this form that is the u at xy is written as ux plus delta x, y, u x y plus
delta y u x minus delta x y and u x y minus delta y whole divided by 4 and so that is
what we have written out and we will solve this iteratively that is to start with we
divide the space into a grid.
So we have this we divide this into a square grid and then we had the values specified
at each of these grid points to begin with and then at the next step, we would evaluate
the values at each grid point by using this procedure, this equation and we will continue
this procedure till we have the left 2 sides are the same that is when u x, y obtained
in the next iteration loop is the same as what was in the previous step that was the
way of solving the elliptic equations which we did in the last class.
So now we look at slightly different form that in which now we look at equations which
are time dependent. So this particular equation we solved in the context of the distribution
of temperature that is we said that if you have a slab in which a slab a conducting slab
and which is kept at whose boundaries are kept at different temperatures for example,
this is kept at let us say hundred degrees and this is 0 and 50 and 75 and the slab is
insulated from the top and the bottom, so the whole heat flow is only inside the plane
and then we wanted to know what the temperature distribution would be and that we saw that
leads to an equation of this form where u is the temperature.
Now the same problem if you want to look at the case where it is time dependent okay that
is we want to know how does now this is steady state equation. So now you want to know how
starting from a given initial condition that is this is kept at 175, 50 and 0 and everywhere
else also is 0 to begin with and you want to know how does it go in to the steady state
and in that kind of cases we will have to solve a time dependent equation. So now in
this course we will only look at one dimensional time dependent equations that is one space
dimension and then one time dimension that we will look at only those particular cases,
okay.
So that is equations of the form we will be looking at equations of the form del u where
u is again temperature is equal to some coefficient times t is little t is the time here some
coefficient times del square u by del x square. So we will look at equations of this form
so here is one example where this arises in the case of a heat conduction through a thin
rod.
So I can assume that this to be a very thin rod, so that across the cross section we can
assume the temperature to be uniform and let us say it is insulated all over its length.
So what is exposed is only the two boundaries and the two boundaries are kept at two different
temperatures and we want to know what will be the temperature distribution along the
length of this rod as a function of time. So that problem would lead to an equation
of the form which I just wrote.
So we are not looking at steady state here we are looking at this long thin rod thin,
we assume it to be thin, so that we can we take it as thin. So that we can assume in
the cross section across the cross section the temperature is uniform and it is insulated
all along its length. So the heat enters here and goes out here so you take any element
where the heat comes in here at is qx and what heat goes out here is q at x plus delta
x. So that will be the heat flux in to the cross sectional area which is just take it
is as a.
So the heat flux let us take it as q here and then q of x and then what goes out will
be q at delta x as we saw earlier so now we know that the q is the gradient of temperature
that is del u by del x in this particular case.
So these equations would lead to something like this we consider an element dx and then
the heat enters at x and it goes out at x plus dx, now the amount of heat that flows
into this thing is qx, this qx is k a in to del u by del x. We know that that is what
we, this is the Fourier's law. So what the heat enters is qx at this point here. So that
is k times a into the gradient in the temperature, u is the temperature here. So gradient is
the temperature "a" is the area and "k" is the heat conduction coefficient, so that what
we have and then we can write and we know that this heat which goes out is qx plus delta
x.
So the net change in that what comes in and goes out would be equal to del, let me write
this in this form. So you have this rod, so we are considering this particular section,
so what comes in here is qx and what goes out here is qx plus dx for this thickness
and so the change would be qx plus dx can be written as q of x plus del q by del x into
dx. So what goes out of what is retained here the heat which is retained here in this section
is qx plus dx minus qx which is del q by del x del x del q by dx.
So qx plus dx minus q of x is del q by del x into dx that is what we have. Okay now this
q itself, we know is minus k times a times del u by del x. So substitute that here so
the heat retained here per unit length per unit area, so the heat which is retained in
this area per unit length per unit cross sectional area would be given by the heat retained would
be per unit length per unit area would be del q by del x right.
So that will be the heat which is retained here, okay per unit are per unit length would
be so 1 over into 1 over a, that will us k times del square u by del x square, now that
heat which is retained causes a change in temperature there.
So that leads to the equation that del square u by del x square into k is del u by del t.
So that is one case where you get this equation, that is take a long thin rod insulated on
the sides, okay so heat comes in here goes out there and the cross section we assume
the temperature to be uniform and then the heat comes in here if you take an element
of length dx what comes in here is qx what goes out is qx plus delta x.
So retained is del q by del x into dx, so the heat which is retained per unit area,
unit area of cross section, per unit length would be 1 over a del q by del x but we know
q is KA into del u by del x. So we get the heat retained is actually del square u by
del x square and that retained heat is actually used for increasing the temperature there.
So the rate of change of temperature should be proportional to the rate of retaining of
the heat. So that gives this equation del xu by del x that is 1 example, so we saw the
steady state case in the Laplace's equation now we are looking at this case where the
temperature changes as a function of time. So now to conclude what we want to look at
is an equation of this form and we want to know how you can solve this numerically.
So again we could write del square u by del x square as u at x plus delta x minus 2 ux
plus ux minus delta x divided by delta x square and we could write del u by del t as ux t
plus delta t, okay now this is all, so now I can write like this and then minus ux, t
divided by delta t. So that is the forward difference method of writing the derivatives
del u by del t, now the question is so now we can do this.
Okay so we can there is a choice here that is when you write substitute this in to this
equation, since this is only a derivative with respect to time I can write this as u
x t plus delta t minus u x at t divided by delta t at the same position x, now the right
hand side is only derivative with respect to space variable x then the question is whether
we should take this derivative at time t or time t plus delta t, so there is a choice
there.
So what I am saying is that you can actually evaluate this derivative at either t or t
plus delta t. So that leads to two different methods if I evaluate this at t, I would call
that as an explicit method and if I evaluate this derivative at time t plus delta t, I
call that implicit method. So these are the two methods which you should be looking at
today.
Today our aim is to solve this equation this particular equation using two different methods
that is substituting for these derivatives, their difference equations but taking the
space derivative in one method at time t that is called the explicit method and another
method in which we would take this derivative at time t plus delta t and that would be called
the implicit method we will look at these two methods, these two equations. So let me
summarize this again.
So we want to look at equation of this form which arises in the case of thermal transport
in a long thin rod and the time dependent equation for that and we want to solve this
and we will discuss two different methods and one is the implicit method in which the
second derivative here u is evaluated at time t plus delta t and there is another explicit
method in which this is evaluated at time t.
So here is a explicit method, so it is very similar to the one which we use for the steady
state case, so we replace the second derivative here the del square u by del x square by this
difference equation. So now the notation we use is the following that we descriptize space
and with the space positions are known l, so u of x we call as ul and u x at time t.
So u x at time t you would now denote by ult, so u x plus delta x at time t will be written
as ul plus 1 which is equally spaced points we are going to take we are going to take
equally spaced points in the spatial variable, so it is u l plus 1 time t.
Similarly, t plus delta t will be t plus 1here.Okay so here we are taking all second derivatives
at time t. So that will be utl plus 1 minus 2 ult plus ul minus 1 t divided by delta x
square and now there is the time derivative del u by del t as ult plus 1 minus lt by delta
t. So now these 2 equations substituting into our time dependent partial differential equation,
we get we can write that as this one that is u, now what we did was substitute these
two into this equation. So we are going to write it in time t, so we are doing explicit
method so all this is evaluated at time t, so what do we get if you do that.
So we are going to get ul at t plus 1minus u l at t divided by delta t is equal to u
l plus 1at t minus 2 ult plus ul minus 1 t divided by delta x square and there is a KA
coefficient here it is just k. So that is the equation we have, so we can rewrite this
in terms of, so what we want to find out is what is the temperature at t plus 1 at time
t plus 1, so we can write ul at time t plus 1 is equal to k times delta t by delta x square
into ul plus 1t plus 1 minus 2, minus 2 ult plus ul minus 1 t plus, so now we have so
now plus ul at t.
So that is the equation which we will be writing, so this is easy to solve because at given
time t, that is at time t equal to 0, for example or at any starting time we are given
the temperature at all the points. So we have this particular rod long thin rod, so we have
been given the conditions at the boundary and then we write we descriptize this into
equal intervals in space with interval delta x and at time t equal to 0 we have given the
temperature at all the points.
We know what the temperature at all the points at time t equal to 0 is and then at time t
plus 1 that is t plus delta t what is the temperature is given by this equation. so
the right hand side is completely at time in terms of temperature at t and the left
hand side gives us the temperature at t plus 1. So using this equation we can easily solve
this, so now this term is called the parameter lambda so that is called lambda this equation
this particular parameter is called lambda k times delta t by delta x Whole Square.
So now note that one problem with this equation is that, it is second order in delta x in
space variable and in first order in time so because of this particular reason it is
not equally its sensitivity is not the same in terms of the spatial descriptization and
time descriptization that is reflected in this particular parameter.
So that is when we change delta x what is the effect which it has on this equation is
not the same as when we change delta t the time descriptization .So we have descriptized
both time and space here but the effect of the descriptization on the stability of these
equations are not the same for time and space and that is because it is first order in time
and second order in space derivative.
So this simple equation can be used to solve the heat equation but the only problem is
that this particular equation is found to be extremely sensitive to the value of lambda
which we use that is this if you use a very large value of lambda this is unstable.
We have to use small value of lambda that is we have to use lambda typically less than
or equal to ".5" so which is less than equal to half to get good results on this equation
that is the biggest problem with this equation otherwise, it is extremely simple that we
know at all time starting time the initial time we know the temperature at all points
on this rod we can use this equation to get the temperature at the all later times using
this once we have fixed this quantity lambda.
So that is what is written summarized here. So we have this spatial derivative written
as this discrete equation and this time derivative written using the again forward derivative
forward difference equation and we substitute that into the heat equation and we obtain
a simple equation of this form, where lambda given by k times delta t by delta x square.
So now we can use given a time the temperature at t minus 1 we can compute the temperature
at time t using that equation.
So then as I said the problem with this method is that this method is convergent that is
stable only when delta t is less than or equal to delta t square by 2k or lambda has to be
less than or equal to half lambda being k times delta t by delta x square has to be
less than half and that is what we see. So now that is the problem with this explicit
method.
So if you want the implication of this is the following that if you actually change
the grid size delta x in spatial grid, we also have to change the delta t correspondingly.
So that is time we want to keep lambda a constant that is, so that we choose a lambda for which
the equation is stable and we give stable solution and we want a little more, let us
say spatial resolution. So you want to know so we choose a particular delta x and we solved
the problem using this solve this equation using this particular explicit method and
then we want to increase the resolution in space.
We wanted a temperature little more to each other that means you want to reduce the delta
x. So what happens is if you change delta x here to keep the delta, to keep the lambda
the same we also have to change the delta t so you can see that if I actually halve
the delta x so to get the same stability and to see that this equation is stable I have
to go 1 by 4 in delta t because this is delta x square this is delta t.
So if I want to increase the spatial resolution by going to half the space distance the descriptization
half my descriptization in space then I have to go 1 by fourth in time. So my computational
expenditure will go up by 8 times. So to solve that for a given time t we have descriptized
both time in terms, in intervals of delta t and space in terms of delta x if you want
to go to higher accuracy in the spatial variable and if you halve the delta x keep the same
lambda, so that the equation is stable we have to go 1 by fourth in delta t. So that
means we will increase our computational expenditure by 8 times so this is one of the drawbacks
of this particular explicit method.
We will see this method implemented in a program before we go and see a little more accurate
method which are the implicit methods. So we look at this program which uses the explicit
method for solving the heat equation. O kay so we have we are using a lambda which is
equal to ".5" and what we do is we descriptize the space into 5 points and you use the explicit
equation which is saying that s at j it is actually saying that s at j plus 1 plus s
at j minus 1 and 2 s j minus 2 s j multiplied by lambda added to s j is the next one which
I had to store as d of j here.
So d is time t plus 1 and s is time t that is what you are going to use. So before we
start that we need some boundary condition that is we have these equations and I am1
as d the temperature and then we need to implement this equation at every point on this descriptized
space descriptized space, the problem is again as we saw in the Laplace's equation that we
cannot use these equations at these 2 initial, 2 end points here, we do not know what l plus,
l minus 1 is and here we do not know what l plus 1 is.
So either we have to use we have to use some boundary conditions, the boundary condition
could be that this the two end points the temperature is specified. So we could say
that this is equal to 100 and u is 100 here and u is 0 here or we could use the derivatives
at the 2 boundaries fix the derivatives at the 2 boundaries and get the temperature corresponding
to 2 fictitious points at that side. So that the derivative can be fixed and then we would
have to solve for all the points including the boundary points.
So in the case where the boundary condition is specified in terms of the variable at the
2 boundaries that is what the program is for we have to solve only for the interior points.
So these two boundaries the value of u or the temperature is fixed all time for all
time it is fixed as 100 and 0 and we solve for all interior points using a lambda equal
to ".5". So that is what we do so we have time going from 1 to 50 and we solve this
for every time step so to start with everywhere else is put as 0, so one end is kept at 100
and every other point is kept at 0 to begin with and the last point is always maintained
at 0.
So that is the two boundary conditions which we have, so now we solve this using this this
is the time t plus 1 and this is the time t, so time t plus 1 the temperature is equal
to time t plus lambda times this particular second derivative of the function and then
after that calculation is over I transfer the time t plus 1 to time t and then goes
back and do the next time step etcetera. So every five time step I write that temperature
profile into a file named as explicit one explicit 2 etcetera now this is the method
of writing the files like that.
So I can use this function which you have not seen before may be it is juts called s
printf. So I have a character string here and into
this character string I am putting a name for a file the name of the file is explicit
and a number here and the number is given by ib the number is percentage d is given
to put the number here. So now this character string will get this name explicit and whatever
the number which ib takes.
So ib takes t by 5, so when t is 5 it is 1 and t is 10 it is 2 etcetera, now this function
can be used to automatically generate file names this generates file names explicit 1,
explicit 2, explicit 3 etcetera as we run this program. So every 5 time step it will
store the temperature into a file whose name is generated as explicit 1, explicit 2 etcetera,
so now it is storing the temperature and the l value or the i value in this case which
corresponds to discrete points on the along the lot.
So we run this program and so we compile that and we run this program and then now we have
all these data files which I have written as explicit 1, 2, 3, 4, 5 etcetera this corresponds
to temperature at time 5, 10, 15 up to 50, so 0 to 50 is what we have so now we plot
this one of the things that means we get a time t equal to 5 that is to it the temperature
goes from 100 to 0 as we go from 0 to 5 in space.
So what we see is that the temperature drops and so that is what we expect and it goes
to 0 in the continuous fashion. Remember, we have only 5 points in the interior and
now we could plot this time two also so we write the next time step that is explicit
2 or explicit 4 may be now that is after few more time steps.
So you can see that by the 4 that is 20 time steps we have actually reached the value,
the steady state value where this linearly drops that is what we the steady state solution
would be that the temperature drops linearly from one end to the other end and that is
what we see at the time 20 time steps. So this method definitely works for lambda equal
to .5 and gives us the evolution of the temperature profile as a function of time.
So now let us look at the same program, let us look at the same program when we change
the lambda to "2.7" so we increase to, we change this to ".71" that is what I have written
in the comment here then you will see that it does not actually work very well for lambda
equal to ".71".
So what we have done is we change the lambda we increase the lambda that is we either increase
delta t or we decrease delta x to increase our spatial resolution probably we increase
delta x or decrease delta x and that is increasing lambda and then we try to run this again.
So let us see what we get so we run this program and then we again plot this as a function
for the first time step so and that is what we see.
So then we see that when we increase lambda instead of getting a smooth profile as we
obtained for smaller lambda, now we have oscillations in the temperature and we are not supposed
to be getting these oscillations. We know that if I keep a rod at 1 at 100 degrees and
other end at 0, I expect the temperature to drop in a monotonic fashion not oscillatory
in space and actually these oscillations would blow up and we would get by next time step
you would see that we cannot actually get the solution correctly. So you can see that
the temperature now has gone completely crazy, that is now we have temperature going to minus
2000 plus 6000 etcetera, so which is obviously wrong.
So this particular method does not work for lambda larger than ".5" I just use a larger
which is higher to get probably a better spatial resolution and we see that the method does
not actually work.
So this is an example to show that even though very simple this method is not very efficient
in solving equations of this form. So we go back to our next description and look at a
little better method called an implicit method.
So as I said earlier, the difference is that in this method the spatial derivative in a
simple implicit method, the spatial derivative on the right hand side of this equation which
we have here, this spatial derivative is now evaluated at time t plus 1. So that is we
will evaluate all these derivatives at time t plus 1. So this derivative would be evaluated
at time t plus 1, so that means that when you write this equation. So that is the discrete
form of the heat equation and all these derivatives are taken at time t plus1. So now obviously
this is, this will not work anymore, so we have to write a different equation for that.
So when this is all evaluated at time t we have a very simple equation for the temperature
at time t plus 1, okay but now we have the right hand side evaluated at time t plus 1
and because of that now we have to when we write this equation for time, when you write
this equation separating t plus 1 and t time variables, so you will have ul minus 1 t plus
1 and plus 1 plus 2 lambda minus lambda. So we have lambda is remember k delta t by x square, so we have minus lambda
ul minus 1 t plus 1 and then you would have 1 plus 2 lambda ul t plus 1 and minus lambda
ul plus 1 at t plus 1 is equal to ul at t.
So now we have the equation in this form and we have to solve for their time at t plus
1. So we write this equation for every spatial
point l and then we will get a set of linear equations. So if you write this equation for
every spatial point, we get a set of linear equations and that set and then has to be
solved using matrix methods to get the correct answer. So again if you write this for example,
l equal to 0 we do not need, so we have again u0 let us say is 100 and we say u some n is
equal to 0 that is fixed at all times for any t this is fixed and then we write this
equation starting from l equal to 1 to n minus 1.
So for l equal to 1 is a special case because then this is 0, so then for l equal to 1we
have it as 1 plus 2 lambda u1 at t plus 1 minus lambda u2 at time t plus 1 is equal
to u0 at time t plus lambda times u, u1 at time t and u0 at time t that is l equal to
1, for all other l that is l greater than 1 and less than n minus, n less than n, n
minus 1 we have this equation that is minus lambda ul minus 1t plus 1 plus 1 plus 2 lambda
ul at t plus 1minus lambda ul plus 1 at t plus 1 equal to ul at t.
So that works for all l between 1 and n minus 1 for l equal to n minus 1 again, we have
problem on this end. So we can write that now as minus lambda un minus 2 at time t plus
1 plus 1 plus 2 lambda times un minus 1 at time t plus 1 is equal to ul un minus 1 at
time t plus lambda times un at time t. So lambda times lambda u0 lambda un are fixed
by this, so these 2 equations are slightly different from the other point equations.
So that is, so that will be this will give us a set of, so we have a set of linear equations
and the variables are the unknowns are like ul minus 1, ul and ul plus 1. So you have
a tri diagonal matrix to solve. So we can write this in a matrix form and you get a
tri diagonal matrix with elements which will be of this form. So I write that here, so
you will have a equation of this form 1 plus 2 lambda minus lambda and then you have 0
0 0 etcetera.
So let us write for 5 points so you have 1 plus 2 lambda minus lambda 0 0 0 and the next
1 would be minus lambda 1 plus 2 lambda minus lambda 0 0 and then you will have 0 minus
lambda 1 plus 2 lambda minus lambda 0 and etcetera. So you will have until you have
the last point would be okay I will write of it, so 0 0 minus lambda 1 plus 2 lambda
minus lambda and then 0 0 0 minus lambda 1 plus 2 lambda that is our matrix and this
multiplied by so this matrix will have to be multiplied by we will have then u1, u2,
u3, u4 and u5. So that is what you will have to get, so we have to get these 2, this matrix
multiplied by this column vector and the right hand side of that will be given by z, t plus
1.
So this is all at t plus 1 which are our unknowns, so we u z time t plus 1 is our unknowns and
the right hand side of this equation would then be this, the first one will be u1 at
time t lambda u0 and the rest would be u2 u3 u4 at time t and the last one would be
u4 at time t plus lambda u5 that will be right hand side of this matrix, this is the matrix
which we have to solve, so we solve this. So we have the equation in this form that
is the first one as the right hand side as, so matrix equation now is u1 at time t plus,
u0 at time t plus lambda u0 at time t, u1 at time t and then u2 at time t, u3 at time
t, u4 at time t and u5 is u4 at time t, u plus lambda times u, actually u5 at time t
and u6 at time t, now this point is solved by this and this are fixed these two values
and other and these also and these are values at time t. So we know all these the right
hand side, so we can solve for these equations at time t plus 1.
So that is the method and we will see an implementation of this method in this particular program.
So here is another program which does that we call it implicit method. So it is simple
implicit method, so here I am using a lambda which is large which is ".7" just to show
you that it works for this large lambda and other things are same as the temperature is
2, boundaries are fixed at s0 and sn are the boundary temperatures that is 100 and 0 and
all temperature at time t equal to 0, everywhere else is 0. So now I construct I have to construct
this matrix so here is the matrix which is constructed.
So I have to start with all matrix elements are put equal to 0 and then I give the matrix
elements these values as lambda 1 plus 2 lambda and minus lambda along the diagonal, that
is along the diagonal, along this is a tri diagonal. So along the diagonal I have 1plus
2 lambda minus lambda 1 plus 2 lambda and minus lambda as the matrix elements.
So that is what is been given here for each of these columns I have i minus 1, I, i minus
1 I, i and i, i plus 1 as minus lambda 1 plus 2 lambda and minus lambda and the right hand
side for the matrix which is given of this equation is given by d i equal to s i that
s remember is the temperature at s time t and d in this particular case s is the temperature
at time t.
So here they are given that and the two special cases, so the first row the first row I have
to this is for the interior points and the first row I have only 2elements that is 0
0 is 1 plus 2 lambda and 0 1 is minus lambda and the right hand side is s0 plus s of 0.
So that is 1 more it is actually for the first row and for the last row again we have only
2 points that is n minus 1, i minus 1 and n minus 1 i. Okay now this is the last part
that is actually i is n here.
So n minus 1, n minus 1 and n minus 1 n, i is n minus 1. So n minus 1, n minus 2, n minus
one n minus 1. So these are the two last points the last row, so I am just writing these 2
elements and these 2 rows especially that is here I have n minus 2 and n minus 1 and
here I have 1 and 2. So a special case and similarly this right hand side equation this
is u5 plus lambda u6 which I call u n in the program and this is u1 plus lambda u0.
So that is the equation, so then I call this function which finds the inverse of this matrix
d, so n is the dimension of this matrix. So the
idea is that I would just find the inverse of this matrix right and multiply with this
then I will get that temperature at time t plus 1. So the advantage of that is see this
is a constant matrix this does not depend on the time the matrix here does not depend
on the time.
So what depends on time is only these 2column vectors, so once I have constructed this matrix
for the given descriptization I have and I found the inverse of this matrix and then
I can just go through a loop to find the t plus 1 as a function of t and that is what
I will be doing in this program. To find the inverse of the matrix, we know we have this
is a function here which finds the inverse of this matrix and so we look at their functions.
So this is the function to which that matrix d as pass as an array of pointers which you
have learnt earlier and the method this function uses is to find, do the l u decomposition
and so we find the inverse by using the l u decomposition which again we have found
earlier. So we have the l u decomposed form of this matrix d which we have seen earlier
in our earlier lectures and then use that l u decomposed form of this matrix and supply
the right hand as 1000100 etcetera.
We give the d matrix the d column vector in this form and find the inverse we have gone
through this earlier. So we have the inverse here this inverse is written back into this
program. So once when I call this inverse function here in this program I call the inverse
function and I pass the matrix d and when it comes back the d is replaced by its inverse.
So because we do not need the d matrix what we need is the d inverse and then I do it
for 10 time steps here. So I am doing for time equal to 1 to 10 and each time step what
I do is I multiply that matrix, since this is a matrix multiplication I multiply now
the temperature at time t plus 1 is temperature at time t plus the matrix multiplied by the
d and then I replace d by that particular time t plus 1 and again treat the 2 end points
separately and goes back into this loop again.
So the method is this, I find the inverse of this matrix and multiply this vector column
vector here by the inverse of this matrix and I find at time t plus 1 what the values
are. Then I transfer the time t plus 1 values in to time t in to this vector and again multiply
it by the inverse of the matrix to get the next time step.
So that is what we are doing here. So to start with you have the temperature 0 every where
except at the 2 end points and then I multiply the matrix the right hand side, the d column
vector by the inverse of the matrix d. This is now replaced by its inverse and then I
transfer what I get into the d column vector and I treat the 2 end points specially because
they are not s m but they are they are not u at time t plus 1 but they are s0 times lambda
plus the u at s values and similarly the last point and then I write those temperature as
a function of time into 2 different files here again I use the same technique to generate
different file names and then I store them all into different file names called implicit
1, implicit 2 etcetera for every time step and remember I am going to run this for lambda
equal to ".7".
So we will now run the implicit program, so instead of explicit we will run it now implicit
program. So that is, so now we will now plot the implicit method. So okay and you see that there is no oscillations.
So the temperature this is at time 1 and this is at time 4 etcetera. So there is no oscillations
here and it slowly keeps improving. Okay so the temperature at any point in between the
2 ends is increasing as a function of time and it goes towards the straight line. So
what I wanted to show you here is that again I used only 5 points, so I used 5 points between
the 2 limits and I use lambda which is larger than ".5" that is why I use lambda ".7" for
which our explicit method did not work and you saw that the implicit method works for
this particular value of lambda.
So what is the difference, in the implicit method we wrote the right hand side evaluated
the right hand side that is the spatial derivative at time t plus 1 and while in the explicit
method the right hand side was purely a function of a temperature at time t. So that is the
difference between the implicit and the explicit method. Okay there are better schemes than
this.
So if you look at this this particular method is again an Euler scheme itself but the only
thing is that we have evaluated the right hand side at time t plus 1. So we could do
even better than this by actually taking a mixture of the type derivatives at time t
and time t plus 1, so that is also that will be what you will be looking at next.
So to summarize here, so we had evaluated the derivatives at time t plus 1 and then
we wrote down the equations for all the interior points and the 2 boundary points were specified
as boundary conditions for all time and then we write down the equations and we write down
the equations for all interior points and then solve them using matrix methods.
So we find the inverse of this matrix and multiply the right hand column vector by the
inverse of that matrix to find the temperature at time t plus 1. So the other method in which
we can use, other implicit method which uses a mixture of these two derivatives that is
we could use the derivative on the right hand side the second derivative that is del square
u by del x square evaluated at both time t and t plus 1 and take the mean of that.
So this is called the Crank Nicholson scheme. So in this method this is again an implicit
method in the sense that the right hand side of the equation del u by del t equal to k
del square u by del x square is still a function of temperature at t plus 1 at time t plus
1 but now it is a mixture of time t and t plus 1.
So again we could use the same scheme So then we would write the equation now in this form
that is ul at t plus 1 minus ul at t by delta t, so this is our del u by del t now the del
square u by del t square is written as mean of the second derivative evaluated at time
t and at time t plus 1 and then we collect all the terms here and write the equation
now for the interior points as this one.
So we have at the interior point equations as t plus 1 minus lambda ul minus 1, t plus
1, 2 times 1 plus lambda ul t plus 1 minus lambda ul plus 1, t plus 1 and the right hand
side now is now a function of l minus, a temperature at l minus 0, l and l plus 0. So this is actually
l plus 1, so that is the difference so we have instead of simply writing it as so what
we wrote in the simple implicit scheme was at ul at time t plus 1 minus ul at time t
divided by delta t. We wrote that as k times ul plus 1 t plus 1 minus 2 ul at t plus 1
plus ul minus 1 at time t plus 1 divided by delta x square.
So now instead of that now you write ul at time t plus 1 minus ul at time t as equal
to lambda which is delta k delta t by delta x square into ul plus 1 time t plus 1 minus2
ul at time t plus 1 plus ul minus 1 at time t plus 1 plus now ul plus 1 at time t plus
ul, ul plus 1 at time t plus minus 2 ul at time t plus ul minus 1 at time t then divided
by half the mean.
So we collect terms here and write this as now, so we have all the terms at t plus 1
on to the left hand side and all the terms at t on the right hand side so you have ul
plus 1 at time t, ul minus 1 at time t plus 1 into, so we have as minus lambda by 2 and
then you have ul. so we have1plus lambda into 21 plus lambda ul at time t plus 1, 1 plus
lambda.
So that is 1 term coming from here and one term here and then minus lambda by 2 ul plus
1 at time t plus 1 on the left hand side and the right hand side then we have equation
as lambda by 2 into ul minus 1 time t plus ul plus ul plus 1at time t minus 2ul at time
t then we have that that is plus ul at time t.
So that is our equation now to solve, so again this is in the matrix form. So we have all
the time t plus 1 on the left hand side all the t on the right hand side. So the only
thing this is exactly the same as simple implicit method except that the right hand side now
has to be evaluated using the temperature at l minus 1, l plus 1 and l at time t, that
is the only difference and then you could use exactly the same inversion scheme, so
once you have set up all the matrix again you could invert that matrix and then for
every time you could multiply the right hand side by that matrix and get the temperature
at time t plus 1.
So that is also when it is one spatial variable and we will look at in the next class what
happens, if you have a two spatial variables that is x and y and how do we deal with implicit
and explicit methods in that particular case that is what you would be looking at in the
next class.