 In the last lecture we started talking about methods for time integration of equation of motion and we introduced the basic terminologies and the role of finite difference schemes in developing these methods, we started talking about the forward difference method or the so called forward Euler method, so here we consider the equilibrium equation we are restricting our discussion to linear systems, MU double dot plus U dot plus KU equal to F of T with some prescribed initial conditions and the time duration of interest is between 0 to TF, a forward difference approximation to the derivatives lead to these equations for displacements velocities and I mean substituting that we will get this equation for the evolution of the displacement, we were discussing the issue of how errors grow as the time marching takes place, so we to illustrate that we considered a single degree freedom system which is driven by F of T, the idea is that even if you have a multi degree freedom system moment you do a normal mode analysis and transform the coordinates to the natural coordinates in the uncoupled state the equations of motion would be a set of single degree freedom systems, so to understand growth of errors it is enough to focus on how the errors would behave when we develop the scheme for a single degree freedom system, so this is the equation based on forward difference scheme for this oscillator and we can put it in a matrix form, capital Y N plus 1 is the vector consisting of a 2 cross 1 column consisting of U N plus 1 and U dot N plus 1 and L N is this 0 delta T into F N, so A is this matrix. Now what we do is at the nth step we assume that an error gamma N is introduced because of truncation and round off and issues like that, now we substitute this equation into the governing equation and we want to develop a scheme on understanding how the errors grow, so from this equation we see that the governing equation for evolution of error is given by gamma N plus 1 is A into gamma N. Now to be able to understand how the solutions to this type of equations behave we considered an equation of the form X N plus 1 is A X N and we introduce the transformation on X which is Phi Z, we introduce the new coordinate system Z by introducing the transformation X equal to Phi Z and Phi was selected to be the matrix of eigenvectors of A matrix and we showed in the last class that the condition for the X, the modulus of X, the jth element of that to go to 0 as N tends to infinity is that the maximum eigenvalue, the absolute value of maximum eigenvalue of matrix A must be less than 1, so this quantity is known as a spectral radius of matrix A, so we want that the spectral radius must be less than 1. An alternative way of looking at it is if you plot the eigenvalues on the complex plane this is real part of the eigenvalue and imaginary part of the eigenvalue and if this is a unit circle we want that for stability the roots must lie within the unit circle that is for asymptotic stability for errors to go to 0, if errors not to grow it is sufficient if the roots lie on the unit circle and if they are outside that the solution would be unstable in the sense the errors will grow as we march in time. Now let us return to the analysis of forward difference scheme and we had this equation YN plus 1 is AYN and from this I developed this equation gamma N plus 1 is equal to A gamma N, now let us rewrite this in terms of gamma naught, so we get that gamma N plus 1 will be A to the power of N into gamma naught, how do we see that? If you write now gamma 1 this will be A into gamma naught, if I write gamma 2 it will be A into gamma 1 which is A square into gamma naught, so continuing in this manner we get this equation. Now this equation we can see that it describes how error introduced at N equal to 0 grows in time, now what we do is we seek the solution of this equation in the form gamma N is equal to alpha N gamma naught, where alpha is a scalar quantity which could be complex valued, so upon substituting into this equation I get alpha N plus 1 gamma naught is A alpha N gamma naught, or in other words for non-trivial solutions we see that this leads to an eigenvalue problem as stated here, and for non-trivial solution the determinant of A minus lambda I must be 0, so the roots of this equation if I write it as alpha naught e raised to I theta because we expect alpha to be typically complex valued so I can write the complex number in this form where alpha naught is amplitude and theta is a phase, so gamma N will be alpha naught N exponential I theta into gamma naught. So now if we are interested in the asymptotic behavior of gamma N as N tends to infinity we want that if this amplitude of these eigenvalues satisfy this requirement that the maximum value of this absolute value of this eigenvalue is less than 1 then it guarantees that the solution decays to 0, so the behavior of the growth of errors in forward difference scheme is therefore governed by eigenvalues of this matrix A, what are the parameters that enter this matrix A? Delta T, omega and eta, now eta and omega are parameters dictated by the problem on hand, on that we don't have any control when we are developing the numerical scheme, so the only algorithmic parameter that we have is Delta T, so now clearly for this condition to be met it will place the restrictions on Delta T, now if we to see that we will actually find the roots, so if I now use this equation determinant of A-alpha is 0 I get this characteristic equation and this is a quadratic equation in this case we can quickly solve that and if we solve this we see that these two, the two roots are given by this. So upon simplification we can show that the roots are given by 1-eeta omega T plus-i omega D Delta T where omega D is the damped natural frequency omega into 1-eeta square, the beta of algebra will lead to this, now therefore what is the maximum value of the amplitude of the eigenvalues, we can compute that and we get that to be this, 1 plus omega square Delta T square minus 2 eta omega Delta T, for stability asymptotic stability we want that this maximum value should be less than 1, so if we impose this condition we get the requirement that Delta T should be less than 2 eta by omega, for the condition of stability to be satisfied, asymptotic stability to be satisfied, so therefore we conclude that the forward difference scheme is conditionally stable for damped systems, this critical step size that is 2 eta by omega this one as we see here is a function of eta and omega, so if system has damping then Delta, the step size that we use should be less than 2 eta by omega to ensure that the errors don't grow in time, but if system is undamped we see that there is no step size which guarantees stability, so this leads to the conclusion that the scheme is unstable for all choices of step size for undamped systems, so in an undamped system this would mean that the scheme introduces an artificial negative damping into the discretized model which alters the basic physics of the problem, so that alteration in the behavior of the system is essentially due to the way we have discretized the governing differential equation, so this you would indicate that forward difference scheme has many pitfalls and we should be cautious if we want to use this. Now let's look at another first order method that is the backward Euler method, here instead of using the forward difference scheme for approximating derivatives we use the backward difference scheme, as you would see this switch results in a substantial change in the qualitative behavior of the integration scheme, so again to see that let's again start with our multi-degree damped forced oscillation problem MU double dot plus U dot plus KU equal to F with prescribed initial conditions and time range of interest from 0 to TF, now let's consider the time instant T plus delta T, and we write now the derivative at that point it will involve the displacement at T plus delta T and at T, so this is the approximation according to the backward difference scheme. Similarly acceleration can be written in terms of velocities and suppose if we use now the expression for U dot at T plus delta T we see that if I solve for U N plus 1 it will be delta T into U N dot N plus 1 plus U N, similarly I get from the backward difference approximation to the acceleration I get U dot N plus 1 is delta T U double dot N plus 1 plus U N dot, now clearly you see that this is an implicit scheme because the right hand side contains unknowns U dot N plus 1 and U double dot N plus 1, so it's an implicit scheme, now let's return to the governing equation what we did just now was to obtain forward difference approximation to the displacement and velocity, now let's return to the governing equation of motion and substitute those approximations, so for U N dot I get this, for N U N I get this, that is I am considering equilibrium equation at time T N plus 1, now this is the equilibrium equation as per the forward difference scheme and by rearranging these terms and collecting coefficients of U double dot N plus 1 I get the time evolution of U double dot N plus 1 is given by this equation, okay. So in this let me again emphasize we have got, we have used backward difference scheme for approximating velocity and acceleration terms, so how do we implement this? A new version of this implementation is we start with T equal to 0 and set N equal to 0, read U naught and U naught dot that is the information we have at T equal to 0 and use this equation in this sequence, first we write the expression for acceleration that is given by this, so when N equal to 0 I am writing U double dot 1, so on the right hand side I have U naught dot which is known, U naught which is known, F N plus 1 is known because F of T is given to us. Now once I, then I go to the expression for U dot of N plus 1 I would have computed this U double dot N plus 1 in the previous step so I can use that, so similarly when I come to U N plus 1 I would have computed U dot N plus 1 in the previous step I will use that, then I increment N and I will stop this process if I cross T F, otherwise I return to 2 and go through this. Now a more useful scheme of for implementation is if we go back here we notice that at every time step we need not invert these matrices especially if delta T is constant for all time steps, then this can be inverted only once and that can be stored outside this loop. Similarly the product of this matrix with this C plus delta TK and K can be evaluated only once and kept outside the loop, this is of course valid only if delta T is independent of steps N, for all N I am using constant delta T, so in this case I will first store all these matrices and initial conditions and the forcing functions and then I will evaluate this matrix M plus delta T C plus delta T square K inverse and product of that into C plus delta TK and K I will evaluate outside the loop, then I will start my count N equal to 0 and then use these three relations, so here I am using now B, D and A which have been computed outside the loop. Now I increment N and I will stop if I reach the final time instant otherwise I return to this, so that means as I am running through these steps 4, 5, 6 I am not evaluating A, B and D repeatedly, so this will obviously enhance the speed of computation. Now let's see what are the stability characteristics of the Euler backward Euler method. Now in the previous study on stability of forward difference scheme we notice that the presence of forcing function doesn't affect the matter related to stability, so we need not consider a forced oscillator but we can consider now an single degree freedom system which is unforced, so the final difference schemes that we are using are listed here and by reorganizing these terms I can write the evolution equation as some matrix which is A in our earlier case and UN plus 1, U dot N plus 1 is a previous step, so the next step if you want you have to invert this matrix and get this, okay. Now this matrix we have to now find the eigenvalues of this matrix, that means this 1, minus delta T this inverse, so this amplification matrix here is inverse of this, this is A, now this matrix we have already encountered in the application of forward difference scheme, so we have found eigenvalues of A there, now we want the eigenvalues of inverse of that, so we can see that eigenvalues of A can be obtained as reciprocals of eigenvalues of A inverse, so that you can verify why that is true. Then the characteristic equation we have already solved and these were the roots that we obtained for A matrix in the earlier case, so this is this and the reciprocals of that would be this now. So now if you separate real and imaginary parts I get this expression for the two eigenvalues, so the spectral radius which is maximum value of lambda bar is given by this expression, now we want that this spectral radius should be less than 1 for stability, now you observe this there is one here and there are other quantities in the denominator which are all positive, so that would mean the denominator here is always greater than 1 for any choice of delta T, so that also means that the spectral radius is less than 1 for any choice of delta T, so this means the integration scheme is unconditionally stable, so you can select delta T in whichever manner you want the errors won't grow. Now let me show some numerical results, so let us consider an undamped single degree freedom system with frequency, a natural frequency 2 pi the period will be 1 second and we are starting from zero displacement and unit velocity that means we are basically computing the impulse response function of the system, so let us start with a forward difference scheme and I will take about 100 points in one cycle of oscillation, so not withstanding that you see that the blue line that I am showing here is the result obtained from numerical simulations and we have already seen that forward difference scheme is always unstable for undamped system and that you see here the error is growing, okay, so the system the discretized system has a negative damping, now let us introduce the damping suppose about 3% damping is introduced and we know the delta T critical we can find out 2 eta by omega is the critical step size, and if I now take for illustration delta T to be 2 times the delta T critical I expect solutions to be, solutions to grow as time passes and that's what happens here because forward difference scheme is unstable if delta T exceeds delta T critical. On the other hand if it is exactly equal to delta T critical you see that the blue line amplitude remains constant but still it's not a good approximation for the red line which is the exact solution, okay. Now the physical system here has a natural damping but my numerical solution appears as if the system is undamped, so this is again a influence of discretization which is undesirable, now how do you remedy this if you are using the same method you reduce step size, so I have half the step size delta T critical by 0.5, so again red is exact and blue is the approximate solution, now the blue lines start showing the decay that we expect but still the results are not accurate enough. Now I take delta T critical by 20, now I see that the two solutions start comparing quite well, so the forward difference scheme in this case with the right choice of step size produces stable and acceptable solutions. Now let me take the backward difference scheme, so it's a damped system, this system the backward difference scheme is stable for any delta T, so now if I take 25 points in a cycle it looks adequate from an engineering point of view but the numerical solution show that that's not quite adequate because the red line is the exact solution and blue is this, so it's not showing, it's showing the comparison is pretty bad. Now instead of 25 points suppose if I take 50 points in a second there is some improvement but still not impressive, 100 points still not good enough, 200 points there is improvement but still not acceptable, 400 points we are nearing acceptance but still not quite there, 800 points seems okay but still there are differences, 1600 points the comparison is pretty good, so that means if you want to use backward difference scheme you have to take about 1600 points in a cycle for this case to get a good approximation to the impulse response of the system, this also points towards the fact that these methods are first order method so you need very small step size to get acceptable solutions. Now I have plotted the norm of the error for the backward difference scheme and the norm that I have used is defined here since we know the exact solution it is X exact minus X numerical whole square plus X dot exact minus X dot numerical and there is a, to make units consistent I have divided by omega, so this is a norm I have shown and this is with 200 points in a cycle the value of this error is about 0.035 and with 1600 points the error comes to about you know 5 into 10 to the power of minus 3, so this is how the backward and forward difference schemes perform with respect to a single degree reference system, so but our objective is really to apply these methods to multi-degree forced oscillation problems. Now we can make few remarks now, what are the modes of unsatisfactory performance of a given scheme? Lack of stability that means errors grow, the numerical dissipation and amplitude distortion even if errors don't grow there may be undesirable damping induced by the integration scheme which may distort the amplitude of the response and it also could lead to distortion of the frequency, then low order of accuracy resulting in requiring smaller step size to get acceptable solution, so these are some unsatisfactory features of the integration scheme that we have seen so far, and one more thing that we have to notice is stability of the scheme does not guarantee that solutions obtained would possess acceptable accuracy, okay, we have seen now we have used step sizes which are quite satisfactory from the point of view of stability requirement for both forward and backward cases for difference schemes, but accuracy was not always acceptable. Now when we apply the integration schemes to multi degree freedom system we need to make you know take into account a few facts, so I will try to start discussing that, suppose we have equation MU double dot plus CU dot plus KU equal to F of T with specified initial conditions and if we are integrating this equation numerically that is directly that is without introducing transformation on the dependent variable that is known as direct integration, so here we would be integrating a set of N coupled ordinary differential equation, that is direct integration, an alternative would be we will consider this equation but we will make this transformation, U equal to Phi Z where Phi is the matrix of eigenvectors of the underlying stiffness and mass matrices, and Phi transpose M Phi is identity, Phi transpose C Phi is diagonal we assume damping to be classical, and Phi transpose K Phi is diagonal with elements being the squares of the natural frequencies. With this normalization scheme the governing equations for Z will result in a set of uncoupled second order ordinary differential equation this we have seen this, so now P of T is the generalized force obtained by using the relation Phi transpose, P of T is equal to Phi transpose F, and similarly we can derive the initial conditions on Z using the orthogonality relations and the mass matrix of the system. Now what we could do is we can integrate each one of these equations numerically find out ZN of T as time history and our interest would be to find elements of U of T, so we have to go back to this equation and find elements of U of T, so in a more superposition method that the numerical solution of a differential equation is carried out in the Z domain, and to get back to U domain you have to again go back to U equal to Phi Z, the relation U equal to Phi Z. Now suppose we are doing modal superposition and we consider this set of equations N equal to 1 to N, in N degree capital N degree of freedom system there will be N such oscillators, and suppose if we integrate all these equations with a common step size of delta T, delta T is common not only for each mode for at all times but also the same delta T is used to integrate Z1, Z2, Z3 and Z capital N, this would be equivalent to evaluating U of T by integrating directly the governing equation without doing mode superposition with the same step size delta T, so this is something that is very important to take cognizance of. Suppose we are using forward difference scheme which is conditionally stable to integrate these equations MU double dot plus U dot KU equal to F of T, in order that we get stable solutions we need to use a step size of delta T which is minimum of 2 eta N omega N, where this N runs from the first mode to the last mode. Now suppose if eta N is constant for all modes then the critical step size is given by 2 eta omega capital N, that means the highest natural frequency plays a crucial role in determining the critical step size if you are using the forward difference scheme, so this places severe demand on the step size, because as you go higher up in the frequency the time periods become small and our requirements on step size become more and more stringent. On the other hand if you use backward difference scheme which is unconditionally stable to integrate this equation, the scheme would provide stable solution for any step size, so there is no requirement on step size from point of view of stability, so the choice of step size now step size to be used in this case is governed by the considerations on accuracy, so here the analyst has to exercise his or her judgment in selecting delta T, so this calls for engineering judgment, that means we should have an idea on how to select step size, because any step size that we use will result in stable solutions. Now from the numerical schemes and the illustrations that we have shown so far we see that about 800 to 1600 steps within a cycle of oscillations at the highest frequency expected in the response has to be used to get acceptable accuracy. Now when an unconditionally stable explicit scheme is used for analyzing a multi degree freedom systems, the requirement on step size to ensure stable solution typically results in stringer requirements than requirements on step size imposed by consideration of accuracy, that means if you ensure that the solution scheme is stable in a multi degree freedom system that means you have to go to the highest natural frequency and find out the critical step size and use the step size which is less than the delta T critical, then by and large the requirement on accuracy would be met because in the response of multi degree freedom systems if we use modal superposition often we see that the majority of the contribution to the response would be made by first few modes, so if you use now the highest mode in determining step size then you will have adequate number of points within a cycle of lower order modes, so that is likely to lead to acceptable solution, so the unconditionally stable schemes therefore accuracy is likely to be guaranteed if you use step size which is less than critical step size, but this is not a universally valid statement but one could typically expect this to happen. Now in the methods that we discussed so far we have used first order approximations and that led to you know for getting acceptable solution from point of accuracy we saw that we need very fine step size to be used in the scheme, now that requirement can be mitigated if you use now second order finite difference schemes to approximate the derivatives, so we start discussing about what is known as central difference method, so first let us review quickly what is the central difference approximation to a derivative, suppose I consider U dot of T I can write this in terms of value of displacement at T plus delta T by 2 and T minus delta T by 2, so delta T is the time difference between these two time instance and this will have accuracy of order delta T square, similarly U double dot of T I can write in terms of velocities as this, now for each of these U dot and U dot T plus delta T by 2 and T minus delta T by 2 I can further use the finite difference approximation, so that means what we are doing here is we have time instant T here and we have T plus delta T by 2 and we have T minus delta T by 2, so at this point we are writing U dot so that will be the displacement here at T plus delta T by 2 minus T minus delta T by 2 by the difference delta T, similarly acceleration can be written, but now when it comes to approximating velocities I need velocities at this point as well as at this point, what I do? I consider other time instant T plus delta T and one more time instant T minus delta T, so they will start appearing in our approximation, because we need to get velocity here again central difference method requires the displacement here minus displacement here, so that's what we are doing here, so this is the approximation to acceleration this is approximation to velocity, now if I now return to the governing equation of motion this is MU double dot plus U dot plus K equal to F of T, now we will consider equilibrium at time T and we will use for velocity this approximation and for acceleration this approximation, this is what we have derived just now, substitute for U dot and U double dot into this equation we get this equation. Now if I rearrange these terms and collect terms which multiply U of T plus delta T and so on and so forth I get this equation, now the equation for U of T plus delta T can be reduced from this and we get this equation, therefore U N plus 1 in this approximation will be this matrix inverse into this column vector, so to implement this scheme what we should do we'll see that, now if we now there is one point that we should notice so this is the expression for U N plus 1, now if N equal to 0, U 1 will have if you see carefully U naught and U of minus 1, okay this will be U of minus 1, now what is U of minus 1, how do we get U of minus 1, U of minus 1 is a hypothetical quantity, so this also means that this scheme is not a self-starting scheme, what we do is we manipulate these expressions and get an approximation to U of minus 1 as follows, so we have U naught dot is U of 1 minus U of minus 1 by 2 delta T, from this I get U 1 as this, similarly U naught double dot will be this, from this I get this expression, so now by solving for U of minus 1 I get this, so this is what I will be doing and we should notice that U double dot of 0 is reducible from the governing equation, so this will be simply M inverse F of 0 minus C U naught dot minus K U naught, so U double dot of 0 is known from the definition of the problem, therefore the term that lie on the right hand side of this equation are also known. Now how do we implement this method, so we can write now the steps, so we set N equal to 0 and obtain U double dot of 0 from the governing equation, because we know F and F of T for all T, then from this I deduce U of minus 1 which is this, then I evaluate these matrices I am assuming that delta T is constant for all times, that would mean the inversion of this matrix need not be done within the time loop, within the loops that through which the time advances, so that can be done outside, so I evaluate all these matrices which I need subsequently outside the loop. Now in this step I evaluate U 1 plus 1, U N dot and U N double dot, for which we already deduce the equation, now we advance N and we'll stop if it crosses the final time instant otherwise we go back and continue advancing the time, so this is the scheme, so again you should notice that the matrices A, B and D are evaluated only once outside this looping, of course this is possible only when if we are dealing with linear systems, if the system is nonlinear then at every time T we need to solve an algebraic nonlinear equation and we'll deal with that sometime later. Now before we consider the questions about stability of the central difference scheme we'll digress a bit and we discuss what is known as Jury's criterion, in our stability analysis we have already seen if A is the amplification matrix the question that we need to answer is whether the spectral radius of A is less than 1 or not, we don't really need to know the numerical values of Eigen values of A, we need simply to check an inequality, so that means the question we are asking is very simple, far simpler than asking what are the Eigen values, we are simply asking whether spectral radius is magnitude of spectral radius, the magnitude of the highest Eigen value is less than 1 or not, to answer that question there is a simpler strategy which avoids the actual determination of Eigen values and that strategy uses what is known as Jury's criterion, so to illustrate that we'll consider the polynomial P of Z as A naught Z to the power of N A1 ZN minus 1 and this is the polynomial where A naught is greater than 0. Now Jury's criterion provides a necessary and sufficient conditions for the roots of the polynomial to lie within the unit circle in the complex plane, so the necessary conditions are P of 1 must be greater than 0 and minus 1 to the power of N P of minus 1 should be greater than 0, now to construct sufficient conditions we construct what is known as Jury's table, that table has this form, so from the given polynomial we start constructing these rows and the number of rows to be constructed is 2N minus 3, suppose this is a second order polynomial there will be only one row, so the rows are shown here, the first row clearly provides the coefficients of Z to the power of 0, Z to the power of 1, etc. The second row is also derived from that and subsequent rows have B's and C's etc, where the B, C, K up to Q this continues, so this is a sequence continues is given by this set of relations, so from the given polynomial the first step is to construct this table, then the necessary conditions we have already seen P of 1 greater than 0, P of minus 1 to be greater than 0 of N E 1 otherwise less than 0 for N equal to R, the set of sufficient conditions is given by terms contained in this table and they are listed here, so the number of sufficient conditions is actually N minus 1, now the criteria do not provide information of the value of the roots, that means we are not finding the eigenvalues but we are only checking whether the roots lie within unit circle, now. So let's try to now return to the question of stability analysis of the central difference method, now let's consider this equation X double dot plus 2 eta omega X dot plus omega square equal to R of T and if we apply now the central difference method I get for X N double dot this relation and for X N dot this, so substituting that this we have done we get this equation, for reasons that will become clear we will combine this equation with the identity X N equal to X N and we get this equation, so this is a vector of system states at N and N plus 1 and this is the amplification matrix given by this, so the stability of the scheme depends on eigenvalues of A and eigenvalues of A there in turn depends on the natural frequency damping and delta T, natural frequency and damping are parts of the problem specification, so it is not a, they are not the parameters that we can choose, so the only parameter that we can choose is delta T and by examining eigenvalues of this we will be able to tell what is the critical delta T, so we get by formulating A minus lambda I determinant of that equal to 0 this polynomial equation, this is a characteristic equation. Now of course I can solve for this lambda and find out what the eigenvalues are but let us now use Jury's criteria and see what we get, now P of lambda therefore is lambda square minus lambda this, this is the polynomial associated with the characteristic equation, the necessary condition is P of 1 must be greater than 0 so that leads to the requirement omega square delta T square divided by 1 plus eta omega delta T must be greater than 0, numerator is positive, eta is positive, omega is positive therefore this requirement is met, similarly minus 1 whole square P of minus 1 if we manipulate we get this expression and this has to be greater than 0, so for this to happen we get a condition on delta T and that turns out to be delta T less than 2 by omega which is this 2 by omega is nothing but Tn by pi where Tn is the period 2 pi by omega n. Now the sufficient conditions how many are there 2 minus 1 which is 1, so that is modulus of A naught must be less than modulus of A2 that means 1 should be less than 1 minus eta omega delta T divided by 1 plus eta omega delta, so this is always satisfied we can see this, examine this and we can verify that it is always satisfied. So now delta T critical is obtained as 2 by omega, so the step size that we select must satisfy the requirement that delta T should be less than 2 by omega and we can immediately notice that this critical step size which is 2 by omega is independent of damping in the system, so in the forward difference scheme if you recall we got 2 eta by omega as the critical step size, backward difference scheme of course was unconditionally stable but here we get the requirement that delta T should be less than 2 by omega. So therefore we can observe few properties of this method, the method is an explicit method with second order accuracy, the method is conditionally stable with critical step size of 2 by omega, the critical step size is independent of damping, now the method requires a special starting procedure it is not a self-starting scheme and we have assumed in our discussion that the step size is constant, now if M is diagonal and C is proportional to M then this matrix A, if M is diagonal and C is proportional to M this matrix will be a diagonal matrix, therefore the implementation of the method does not require inversion of any matrix, so that is when you are dealing with large system that this is something that we should take note of and it is also important to note that this requirement may, that M is diagonal and C is proportional, the observation that the method does not require inversion of the matrix is not particularly significant if you are dealing with linear systems, if you are dealing with either multi-step methods or non-linear systems this becomes an important observation. If this condition is not satisfied then we need to invert this at least once if you are using same delta T, otherwise you have to invert this at every delta T, that is if delta T is not constant the matrix needs to be inverted at every step, so this can place severe demands on computation if you are dealing with large scale problems. The method is easily implementable even if the system is non-linear, because it is an explicit method there is no solution of non-linear algebraic equations as time progresses. Now delta T is the only algorithmic parameter to be selected and that has to be less than 2 by omega, so some numerical illustrations, so let's start with an undamped free vibration of a system with natural frequency 2 pi eta is 0 and we will select a step size which is about 5% more than the critical step size, so we see that the red line is the finite difference approximation, so the error is growing and it swamps the true solution, so on this scale the true solution appears like a flat line. Now on the other hand if delta T is about 95% of the delta T critical, so things look at least the errors don't grow but the solution obtained is quite unsatisfactory, okay, blue is the exact solution, red is the solution from central difference scheme. Now if I now take step size to be delta T critical by 5 the solution starts looking much better, but still as you see as you progress in time there is still a distortion in phase and there is some observable difference between blue and the red lines. If you take now delta T to be delta T critical by 10 then the solution looks quite good, the red and blue look almost identical, of course one can plot the error norms and examine what exactly are these differences but through these graphs we can see at least qualitatively that the solutions are agreeing quite well. Now let's consider damped system, this critical step size is independent of damping that remains constant, suppose if I take now delta T to be 1.05 times delta T critical which is not a right choice to make the solutions will become unstable that is what we see here. On the other hand if I take 95% of delta T critical the errors there are no growth of errors but still the solution is, numerical solution is very poor, so with factor 5 that means delta T is delta T critical by 5 we start getting good solutions and with delta T critical by 10 I get quite pretty good solution. Now we want to now study the stability of other integration schemes that I will be introducing shortly, central difference method is an explicit second order method it is conditionally stable, so the next question that we need to ask is are there any schemes which are implicit and unconditionally stable and have second order accuracy, so that discussion we will take up in the next class, as a prelude to that we need to have some idea on what are known as Z transforms that we will be using in analyzing those other schemes and we can also as well analyze what we have already done using this approach, so let's quickly recall some preliminaries about that. So let F of T be a function of time such that F of T is 0 for T less than 0 and let S be a complex number, we define the Laplace transform of F of T as F of S 0 to infinity F of T e raised to minus S T dt, now if F of T is a sample that means it is multiplied by a comb like function a series of direct delta functions, so how does this function look like? F of T is a series of direct delta functions placed at a common interval capital T, now if F of T is written like this that means it is the F of T is sampled at all these time instance the sampled waveform can be written in this form, now let us consider Laplace transform of that, so I substitute this into this, so upon carrying out the integration I get Laplace transform F of S as N equal to 0 F of N T e raised to minus S N capital T, now this can be written as F of N T exponential S T to the power of minus N, suppose if I call this quantity e raised to S T as Z, I can write this as F of N T Z to the power of minus N, so this quantity is known as Z transform of F of T and it is written as F of Z as this, so this is a discrete version of Laplace transform, so this is the definition, now we need to recall some result from geometric progression a series like 1 plus R plus R square up to R N is called a geometric progression, and if I multiply now this by R I get R plus R plus square so on and so forth plus R to the power of N plus 1, suppose if I subtract this I get 1 minus R into S N which is equal to this from which I get a sum of first N terms in a geometric progression to be given by this, so if modulus of R is less than 1 then as N tends to infinity S N goes to S which is 1 by 1 minus R, now we can quickly see few examples, if F of T is a step function, F of N T will be 1 for all N otherwise it is 0, now using that and using this result we can show that the Laplace, the Z transform of a unit step function is given by Z divided by Z minus 1, similarly if I have F of T is KT that means a straight line then F of N T will be N capital T and the Laplace transform of that we can show through some manipulation that it is given by T Z divided by Z minus 1 whole square, this you can verify, there are few more examples, for example we having exponential if you are having a F of T as E raise to minus A T then F of N T will be exponential minus A N capital T, now the Laplace transform of this if I do I will get this as Z divided by Z minus exponential minus A T, similarly we can verify that the Z transform of sine omega T and cos omega T are obtained through these expressions, so I leave this as an exercise for you to verify. Now a few properties of Z transform can be quickly recalled, the Z transform of A of F of N T is A into F of Z, similarly the Z transform of A to the power of N, F of N T is F A minus 1 Z, Z transform of F of T minus N T is Z minus N F of Z, Z transform of F of T plus N T is given by this, okay. Now if we now I mean this can be used for K plus 1, K plus 2, K plus N so on and so forth, suppose if I now consider Y of K to be summation N equal to 0 to K X of N, if we have this Y of K I can write Y of K as Y of K minus 1 plus X of K. Now if I take Z transform, Z transform of Y of K minus Y of K minus 1 is Z transform of X of K which is X of Z, so from this I get Y of Z to be X of Z divided by 1 minus Z, so if there is a dynamical system which evolves as YK is equal to YK minus 1 plus XK then the X of Z divided by 1 minus Z would be its transfer function. Now to make that clear we will consider a dynamical system with input X of T and output Y of T, let the input-output relation for the sample data be given by this equation that may Y of N plus A1 YN minus 1 plus so few terms up to AK YN minus K, this is equal to B naught X of N plus B1 XN minus 1 and so on and so on up to N minus M. Now if I take now Z transforms on both sides and use the identities that I briefly mentioned we can show that the frequency response function in the Z domain is given by this ratio and we get this as the frequency response function. Now the idea of introducing this Z transform is basically to visualize the finite difference schemes that we are going to develop as a kind of input-output relation of this form, so consequently associated with the finite difference scheme that we use we can define a frequency response function. Now this frequency response function will be in the form of a ratio like this and the zeros of the denominator are known as zeros and the zeros of the numerator are known as poles, so the zeros of the frequency response function must lie within unit circle is a requirement that we need to verify. So here if you are following this approach we need not formulate the A matrix, the amplification matrix instead we can formulate the frequency response function in the Z domain and look at the zeros and poles of the frequency response function. Now to quickly see that let's revisit the central difference scheme and we got this as a evolution equation which is this, now if I write this for N I get this equation on the right hand side there will be N minus 1, now if I take Z transform I get this equation. So for non-trivial solutions you know the this bracket terms inside this bracket must be equal to 0 and this leads to the polynomial equation that we are looking for, and this is identical to the characteristic equation we obtained by considering eigenvalues of the amplification matrix. So now what we can do is instead of applying Jury's criteria on this, criterion on this we can apply the Jury's criterion on this equation, so if in this case we will get the same solution as we got before because we applied Jury's criterion on this equation to derive the critical step size, and if we had done this using this logic would have got the same answer. So now let's conclude this lecture at this stage, in the next class I will ask the question are there any second order schemes which are implicit and which are unconditionally stable, so that takes us to discussion of what are known as a Newmark family of methods and we'll take up the discussion in the next lecture.