 We have been discussing solution of equilibrium equations in dynamics using time-stepping methods, so in today's lecture we will talk about second order implicit methods. So far we have discussed this forward Euler scheme which we have seen to be an explicit scheme it is conditionally stable for damped systems and it is unstable for undamped system for any step size, when stable it has first order accuracy, the backward Euler scheme is implicit it is unconditionally stable and it possesses first order accuracy. The central difference scheme is explicit it is conditionally stable it has second order accuracy it is not a self-starting scheme whereas these two schemes are self-starting. Now we will continue with this discussion we are interested in developing methods that possess second order accuracy and let's consider if we can develop methods which are implicit in nature. Now before we take up that topic we will revisit some issues about choice of critical step size in conditionally stable schemes. See we had for forward difference scheme the critical step size was 2 eta by omega, so in terms of period it is eta T by pi, for central difference method which is again conditionally stable the critical step size was T by pi where T is the period. Now if we apply any of these schemes to multi degree freedom systems it means that delta T must be selected such that these conditions are satisfied for all the participating modes, that means all the modes in the system in a n degree freedom system there will be n modes and this condition should be satisfied by all the n modes. If you are doing direct integration the question arises how do you estimate the maximum frequency? So delta T critical is actually for forward difference scheme the minimum of this taken over all the possible modes, so if damping is constant for all modes it is 2 eta by omega n, if eta n is the same for all modes, for central difference scheme it is 2 by omega n which is 2 by omega capital N if central, if we use central difference scheme. For mode superposition method each mode can be integrated at least in principle with different step sizes with higher modes requiring finer step size, but when we construct back the solution in the original coordinate space you need to take care of adding the modal contributions at appropriate times. Now given a finite element model for a multi degree finite element model how does one estimate what is the highest natural frequency? Typically in professional softwares you will not be able to compute the highest natural frequency if you do an eigenvalue analysis. Suppose omega n is the highest natural frequency of an assembled FE mesh with NE elements, let KE and ME be the element matrices for E equal to 1 to NE, now what we do is we consider the eigenvalue problem for the Eth element, so KE phi E is omega E, omega E square ME phi E, now let omega NE be the largest frequency of element E, now we can show that the highest natural frequency of the model is actually the maximum of the highest natural frequency of the all the elements, thus for central difference method and a bound on delta T, the delta T critical can be obtained as shown here, so that means for each element if we were to do an eigenvalue analysis and find the highest natural frequency for every element the maximum natural frequency of the assembled system will be less than or equal to the highest natural frequency of the elements taken over all the elements, so how do we show that? So let's start with the eigenvalue problem, so omega n square that is the highest natural frequency is given by this, this is exact there is no problem, now for K I will write it as KE summed over E from 1 to NE, now similarly for ME I will write this, so I can take now this phi n transpose inside I will get E equal to 1 to NE phi n transpose KE phi n, now let's call this phi n transpose KE phi n as UE, and similarly phi n transpose ME phi n as IE, now for the Eth element we have by property of Rayleigh's quotient, phi n transpose KE phi n phi n transpose ME phi n which is UE by IE will be less than or equal to the highest natural frequency of that element, so this would mean UE will be less than or equal to this quantity, now let's return to this expression for omega n square and I can write this now in the denominator since UE is less than or equal to this for UE I will write this, so instead of an equality I will have now the inequality, now if I take this term outside the summation it will be maximum of omega NE whole square from 1 to E equal to 1 to NE that goes out and I am left with this ratio which gets cancelled which is equal to this basically, so from this I get the required relation that the square of the highest natural frequency of the assembled system is less than or equal to the square of the highest natural frequency of taken across all the elements. Now here you should understand that KE and ME the sizes of KE and ME are taken to be that of K and M, that means if you recall that AS transpose MEAS is KE it's because for example if a structure has say 20 degrees of freedom and beam element if you are using it may have 4 degrees of freedom but this KE will have 20 degrees of freedom not 4 degrees of freedom, so many eigenvalues the initial first few eigenvalues will be 0 but the highest eigenvalue there won't be any problem, I mean you will get the required result. In solving element level eigenvalue problem we have removed the boundary conditions and interactions with other elements that also you should remember. Now so now we can now look at say for example actually vibrating rod this is the stiffness matrix and this is a consistent mass matrix, so delta T should be less than if you find the eigenvalues of this 2 by 2 stiffness and mass matrix you will get delta T should be less than this square root E by rho is nothing but the speed of the wave velocity in the rod. So what this requirement basically means is that the step size should not be less than the time taken by the wave to pass through the element, so if you are using lumped mass model instead of consistent mass matrix model the delta T will be actually L by C where C is square root E by rho, this is a classical result available in mathematics literature attributed to Courant. Now similarly for Euler-Bernoulli beam with consistent mass matrix we can show that the delta T must be less than this, so for a given finite element model we can find out these critical step size for all the elements and pick the one which is smallest, so we need not really estimate the highest natural frequency this will provide a bound, okay. Now let's take up the problem of developing second order implicit schemes, so let's start with what are known as average acceleration, linear acceleration and Newmark beta methods, now in this diagram X is the time axis and Y is the acceleration axis and I am considering two time instant TN and TN plus 1, so at TN acceleration is UN double dot, and at TN plus 1 the acceleration is UN plus 1 double dot, the question is at a time instant tau lying between TN and TN plus 1 how do we approximate U double dot of tau? According to average acceleration method U double dot of tau is taken to be the average of these two values, so that is half of UN plus 1 double dot plus UN dot, UN double dot, now to get velocity from this I integrate this with respect to tau so it will be UN dot initial velocity plus integral of this is tau by 2 into this is a constant it will remain as it is, now from this it follows U dot of N plus 1 will be this, now if I now integrate this further to get displacement I have to integrate this now, so I will get UN plus tau UN dot I am integrating this and similarly I have to integrate this now tau square by 2 that will be tau square by 4 into this, now for tau if you put you know you reach this point you will get this as UN plus delta T UN dot plus delta T square by 4 and U double dot N plus 1 plus U double dot N, so we have now the expression for the velocity and the displacement in which we are assuming acceleration at any point tau is average of these two, a slightly different version of this is instead of assuming that acceleration here is average of these two we could as well take that to be linear, assume it is a linear variation, so in which case U double dot of tau in terms of these two values I will fit a straight line and the equation for straight line is this, so this is known as linear acceleration approximation, now from this again I can get velocity by integration so the initial velocity is UN dot, so if I integrate this I will get, first term will lead to UN double dot of tau this is tau square by 2 delta T into this constant, from this if I now evaluate velocity at N plus 1 this will be the expression, and now you return to this and integrate with respect to tau I will get U of tau so I get this expression, so again with tau if I put now if I reach T N plus 1 I will get this as an expression for displacement, okay. Now what Numark beta method, Numark method is also known as Numark beta method because a parameter beta appears in its description I will simply call it as Numark method, now if we consider U dot of T we can write U dot of T plus delta T minus U dot of T divided by delta T as U double dot of T plus delta T for some delta belonging 0 to 1 this is by mean value theorem, so from this I get U dot of T plus delta T is U dot of T plus delta T U double dot T plus delta delta T for some value of delta lying between 0 and 1. Now suppose we now approximate U double dot of T plus delta T as a weighted sum right as shown here I will get now U dot of T plus delta T will be given by this, okay, so what I am doing is for this term I am writing an approximation, so delta is now an algorithmic parameter, okay which has to be selected, so now this means in a discrete version it is U dot N plus 1 is given by this, now if we now consider we have to now make a model for displacement, now if you consider Taylor's expansion for displacement I get U of T plus delta T is U of T plus delta T U dot of T plus delta T square by 2 into U double dot T plus 2 alpha delta T for some alpha lying between 0 to 0.5, again by mean value theorem. Now this second term now I approximate as a weighted sum with weights 2 alpha and 1 minus 2 alpha as shown here, and now if I substitute back into this I will get a model for displacement. So here alpha is another algorithmic parameter which is taken to be independent of delta, so in the Newmark's method what we have done is we can take a look at the three methods, average acceleration method, linear acceleration method and Newmark method, so in Newmark method what we have done is we have introduced two algorithmic parameters delta and alpha, in addition to delta T which is a parameter that we already have. Now by adjusting this value of alpha and delta for example for alpha equal to 1 by 4 and delta equal to half Newmark's method becomes average acceleration method, similarly sorry yes correct, then if alpha equal to 1 by 6 and delta equal to half Newmark's method become linear acceleration method, so you can see these expressions and you can see that Newmark's method is an implicit method because U dot N plus 1 will have U double dot of N plus 1, and if I write this in terms of the velocity it will involve U dot N plus 1 terms, okay, so this is implicit, and it's a self-starting scheme, it's a single step scheme. Now we can study Newmark's method and reduce properties of the other two method we need not have to consider these two methods separately. So how do we proceed to implement Newmark's method? These two expressions we have already got for U dot N plus 1 and U N plus 1, now I write the equation of motion at T equal to N plus 1 I get, so I get this, so I am considering basically linear systems in this discussion. So now I will solve for U double dot N plus 1 from equation 2 and I get this, and what I get in equation 4 if I substitute it in equation 1 I get this expression, so this is some simple manipulations you need to sit with pen and paper and work this through. Now if you look at equations 4 and 5 they contain U N plus 1 which is unknown, now if I substitute 4 into 5, 4 and 5 into 3, 3 is what? 3 is this equilibrium equation at T equal to T N plus 1, so I have now U dot N plus 1 equation 4 is U double dot N plus 1 and U dot N plus 1 I have, so if I substitute that I will get an equation for U N plus 1, so after rearranging the terms we get the equation for U N plus 1 as this, so this coefficient matrix consists of MCK matrices and the algorithmic parameters of the Newmark method, and on the right hand side we have A0, A2, A3 and A1, A4, A5 their meaning is displayed here, okay, right. So now how do we implement this? We can look at this expression and a simple minded implementation would be start with N equal to 0, read the initial conditions and derive the initial acceleration from the equation of motion, then use the expression for displacement which is given by this, which is basically this equation, and then you compute U double dot N plus 1 and U dot N plus 1 moment U N plus 1 is known, these terms come on the right hand side here, so you are able to get displacement, acceleration and velocity, then increase N to N plus 1, stop if you have crossed the final time limit otherwise go back to 2, now this can be refined a bit because if delta D is constant, I need not have to invert this matrix and multiply this matrix with MC etc at every time step, so that can be done offline outside this loop, so a more efficient implementation would be you start with N equal to 0, read the initial conditions, compute the initial acceleration, now compute these matrices A, right, then B is A into M, D is A into C, so U N plus 1 will be now AFN plus 1, B into this, plus D into this, so what I am basically doing here is that I am avoiding this inversion at every time step and this matrix multiplication at every time step, so next I get acceleration and then velocity, I increment time and if we reach the final time instant I will exit otherwise I go here, so 3, 4, 5 are the steps that will be repeatedly implemented as we advance in time and step 2 is done outside that loop, so this is a more efficient implementation. Now one of the concerns that we have about integration schemes is the questions about the growth of errors, so for studying that as we have been doing we consider a test system which is basically an undamped, sorry damped single degree freedom system, so let's consider X table dot 2 eta omega X dot plus omega square X equal to 0, now let's implement the Newmark method on this system, so acceleration is given by this, velocity is given by this and we can manipulate this, this is the equation for displacement, we can do a bit of manipulation and we can recast this equation in a matrix form if I define a velocity and displacement as a state vector the Newmark's method leads to a matrix equation of this form, there is a coefficient matrix here let's call it as A, A into suppose the state is called as capital XN this will be A into X dot XN plus 1 and B is this matrix which is on the right hand side B into XN, this is XN, okay, so this is how the solution advances. Now the amplification matrix is A inverse B, now we have to find eigenvalues of A inverse B and look at the spectral radius and see whether that is less than 1, condition for stability is the spectral radius must be less than 1, we can as well consider the characteristic equation B minus lambda A equal to 0, so this is a 2 by 2 matrix, this characteristic equation will be a quadratic in lambda, so I can write with lambda as Z we can use the Z transform method and I can write the polynomial equation as A naught Z square plus A1Z plus A2, now by expanding this determinant I can identify these terms A naught A1 A2 and this is the polynomial, now according to the Jury's criterion what we should do is we have to find P of 1, so P of 1 will be omega square delta T square, then minus 1 to the power of N P of minus 1 is given by this, so the necessary condition for stability is P of 1 must be greater than 0 that means all eigenvalues should lie within unit circle in the complex plane, so P of 1 is omega square delta T square which must be greater than 0, so this is satisfied because omega is positive, delta T is positive, omega square delta T square is positive. Now the second condition if you look at this is the second condition, now this will be satisfied if 2 delta minus 1 is greater than 0 and 4 alpha minus 2 delta is greater than 0, so that leads to the condition delta greater than 1, and alpha greater than delta by 2, okay. Now sufficient conditions this is a second order polynomial therefore there will be one sufficient condition that we need to check that is A2 is greater than A naught, so for that to happen this is the equation that we get. Now this again we can see that if delta is greater than 1, this will be satisfied because all these quantities are positive, so we get the now the conclusion that the method will be unconditionally stable if delta is greater than 1, and alpha is greater than delta by 2. Now these are algorithmic parameters alpha and delta, they are nothing to do with system natural frequency and damping, therefore these requirements if they are met by choosing the algorithmic parameter so that these conditions are met then the method becomes unconditionally stable for all damping and natural frequency of the system, otherwise this is not the only requirement we can derive the conditions for stability as given by this, so this is a stronger requirement we can relax this and still be able to get stable algorithms provided these conditions are met, these are conditional in which case the step size is function of, critical step size is function of natural frequency and damping and the method becomes conditionally stable. Now we have already seen that the Newmark method for alpha equal to 0.25 and delta equal to 0.5 is average acceleration method and we can verify that this is unconditionally stable, then linear acceleration method becomes conditionally stable, for delta greater than 1.5 and alpha greater than delta by 2 the method actually displays what is known as high frequency dissipation characteristics, that is the spectral radius drops below value of 1 for larger omega if we keep delta T constant, that means you look at the amplitude of the highest eigenvalue as a function of omega by keeping delta T fixed, so the spectral radius for higher omega drops below 1, so that means if you are using this scheme for direct integration for a fixed value of delta T the contribution from higher normal modes which are spurious actually will be dissipated out, whereas the lower modes which are important for calculating the response correctly hopefully are not dissipated, so this is one of the requirements of a integration scheme that the higher order mode should be numerically dissipated, especially if system response is governed by contribution from lower modes. Now if delta equal to 0.5 the high frequency dissipation characteristics is not witnessed, you can see that the spectral radius becomes 1 if delta is 0.5 and we don't get any high frequency dissipation characteristics, therefore that is not a very desirable feature. Now if damping is 0 we can work through the requirements on step size and we get this as a requirement. Now one important observation which we need to make, we can make at this stage and later on substantiate is we have seen that for delta not equal to 0.5 there will be higher high frequency dissipation possible, but on the other hand if delta is not equal to 0.5 the order of accuracy becomes order delta T, it won't be order delta T square as we wished it would be, it will be order delta T, so there is a tradeoff between achieving second order accuracy and achieving numerical dissipation for higher modes, so numeric beta method cannot achieve both, you have to sacrifice one of them, so I have plotted a few results for spectral radius, on X axis we have omega into delta T and here I have the spectral radius, so ideally the spectral radius should be less than 1, so I have selected delta equal to 0.5 and I have varied alpha, alpha should be less than 0.25 for stability and we can see here so for example the green line alpha is 0.2, alpha should be greater than 0.25 for stability, if alpha is 0.2 this becomes unstable and similarly alpha is 0.4 which is okay the spectral radius is less than 1, red one alpha is 0.24 it crosses this, so the higher modes become unstable so the errors would grow, okay, so that is the problem, so on the other hand 0.25 actually that color scheme shown here that graph is not clearly seen it will be 1, now here what I have done is I have fixed alpha and varied delta, so this is a spectral radius on X axis it is omega delta T, suppose I fix delta T and vary omega, so for low frequency up to about say 10 to, that is 1 hertz the spectral radius is constant, as you go higher in the frequency the spectral radius drops and the contribution from higher modes get numerically dissipated, so this is a desirable feature. Now we talked about order of convergence, the order of approximation and we can work through some details for a couple of schemes, I will illustrate this with convergence of the Euler's method and the reference for this is a book on numerical analysis by Konte and De Boer which I have mentioned here, it is instructive to go through this is a relatively a simple derivation we can understand quite a bit by understanding what this tells us, so let's consider a differential equation X dot of T is A of X of T, T with initial condition X naught, so X of T is 2 N cross 1, so it's a vector differential equation. Now discrete time approximation let us take TN to be N delta T plus T naught, according to Euler's method we approximate YN plus 1 minus YN by delta T as A of YN TN for N equal to 0, 1, 2, 3, etc., where see the notation we have used is X of TN is a true solution of this equation at T equal to TN, whereas YN is the solution obtained through Euler's method, suppose YN is approximate solution to this equation at T equal to TN, therefore the error would be X of TN minus YN, that's a discretization error, since we expect to specify Y naught exactly that is initial condition there is no truncation there E naught would be 0. Now the objective of this discussion is can we derive an equation for growth of E N, that is first question, second question is can we find a bound on E N by solving that equation for evolution of E N. Let's consider X of TN plus 1 which is X of TN plus delta T, now Taylor's expansion tells us that is X of TN plus delta T X dot TN plus delta T square by 2 X double dot of some XIN where XIN lies between TN to TN plus 1, this is called delta T square X double dot XIN is called the local discretization error, because this is where you will truncate the Taylor series, so this is a discretization error, that is it's the error committed in the single step from TN to TN plus 1 assuming that X of T and X dot of T are known exactly at T equal to TN. Now we are assuming that in this discussion the round off errors that occurs due to finite precision calculation on computer is ignored, you may use double precision or single precision, so there will be round off error at every step, I am not talking about that, I am only talking about the errors that result from truncating the Taylor's expansion at a given term. Now therefore now let's consider the Euler approximation, this Euler approximation, now from 1 and 2 if I now subtract, if I now write an equation for YN plus 1 and subtract from this I will get X of TN plus 1 minus YN plus 1 would be X of TN minus YN plus this will be delta TA and there will be this term, so I will get delta T into X dot TN minus A of this plus this term, okay, fine. Now this first term is nothing but EN plus 1, the left hand side, this is EN, what multiplies delta T is X dot which is nothing but A of X of TN, TN minus A of YN, TN plus this dot, so this is the equation we have got. Now let us consider this A of X of TN, TN and this I can write it as YN plus X of TN minus YN, okay, so then I will do a Taylor's expansion on this I will get A of YN TN plus X of TN minus YN into A prime of this for some XN bar between X of TN and YN, again I am using mean value theorem here, okay, so I can thus write the difference A of X of TN comma TN minus A YN TN can be written in this form. So consequently the equation for evolution of error is given by EN plus 1 is EN plus delta T A prime into this, into EN plus this with E naught equal to 0. Now if we assume that our A is such that I can place an upper bound on modulus of A prime as L and the second derivative as capital Y I can write that this error would now be less than or equal to I am writing absolute values into this, where this L and capital Y are bounds on this A prime and X double dot term present here. The idea here is not so much to evaluate L and Y but to understand how the errors behave, okay. So now the equation for the growth of error is an inequality that I get like this. Now let us consider the difference equation where this inequality is replaced by equality I get XIN plus 1 is XIN 1 plus delta T EL plus delta T square by 2 Y with XIN naught equal to 0, let's consider this difference equation. Now at N equal to 0 XIN naught equal to 0 therefore E naught equal to and E naught is already 0, therefore at least at N equal to 0 this XIN is greater than or equal to modulus EN, then in fact equality holds good here. Now let's assume that this is true for some N greater than 0, I am using induction here, so then you can see that if this is true I can write for EN I will write this and consequently I get this, so it follows that this error is bound by this XIN plus 1. Now this is a difference equation which can be solved, okay, so thus to obtain a bound on the error we need to solve the difference equation which is equation 5, how do I solve that? You consider the difference equation, first you consider you remove this constant term, if you consider only XIN plus 1 is XIN into some constant we can assume solution to be in the form of C beta to the power of N and from this I get C beta N plus 1 is C beta N 1 plus delta T L therefore beta is this. Now what I will do is for this equation I will assume solution in this form, XIN is C into 1 plus delta T L to the power of N plus B for some B, now if I substitute this back into this and find a B which satisfy that equation then I am done, indeed it is possible, so I will write for XIN plus 1 this term with N plus 1 here and on the left hand side XIN I will write this I get this equation, now after simplification I will get that this equation will be satisfied if B is given by this, so consequently I get XIN is C into 1 plus delta T L to the power of N minus this, now C is a constant which is yet to be determined and since we know XINOT is 0 by substituting this here I get C also, so the final solution would be given by this, okay. Now you recall E raise to X can be written as 1 plus X plus X square by 2 into some exponential of Zeta where Zeta lies between, Zeta is a you know again mean value theorem is applied here 0 to infinity or range of X, therefore from this I can conclude that exponential of X is greater than 1 plus X for all X, so therefore for this 1 plus delta T L I can write it as 1 plus delta T L is less than exponential delta T L and therefore Nth power of that will be less than exponential N delta T L, so using these facts I will be able to get a bound on XIN which is given by this, consequently I get a bound on the way the error grows, okay. So this is one instance where this type of analysis is possible we are getting this is the error and we see that this is the global error is of order delta T here, that is the error tends to 0 as delta T tends to 0 like C delta T for some constant C if T equal to TN is held fixed. In the expansion however if you look at the local error is order delta T square before integrating before finding EN, so this is a typical situation where the local error is if it is order delta T square global error will be of the order delta T, now this function this inequality provides an upper bound, it may not be a realistic upper bound, the objective is never to bound this error and use it as a prior estimate of error or things like that, the objective here is to understand how the error grows and what is the nature of the error, so the main conclusion from this is contained here that the global error is of order delta T, now how do we study make a similar study for the numerics method? Now let's consider the exact solution to a governing equation say MU double dot plus U dot plus K equal to F of T, they let the exact solution be U of T, U dot of T and U double dot of T, displacement velocity acceleration, let the numerics approximation to these quantities be respectively U1, U1 dot, U1 double dot, now according to the numerics approximation U dot N plus 1 is given by this, now let us consider U dot, U double dot of T plus delta T, that's what I need here, let us see what we get, this is U double dot of T plus delta T into U triple dot of some Zeta for Zeta lying between T and T plus delta T, again use of mean value theorem, so now you substitute that result into this I will get U dot N plus 1 is U1 dot plus delta T, 1 minus delta U1 double dot plus delta into for this term I will write this that will be delta U1 double dot plus delta delta T U triple dot Zeta, similarly we have for U dot T plus delta T, U dot of T plus delta T U double dot of T, delta T square by 2 for U triple dot Zeta for Zeta lying, now if I subtract these two I will get now the error term this is EN dot, EN plus 1 dot that is U dot of T plus delta T minus U dot N plus 1, this is the exact solution this is the numeric beta solution, numeric solution, this is again error and this is the you know terms that helps us to understand the error, so the equation I get is E dot N plus 1 is E dot N plus delta T E double dot N plus this term. Now similarly I can do the analysis for U1 plus 1, you can go through this and we get for U that a displacement U of T plus delta T minus U1 plus 1 is given by this, so in terms of error notation E I get from the first equation E dot N plus 1 is this and from the equation for displacement I get this, now if you see now if delta is half the velocity the error term will be of delta T cube, because this is half the remaining term will be of the order delta T cube, similarly if alpha equal to 6, 1 by 6 displacement is order delta T to the power of 4, these are local errors, and this is for linear acceleration method. Now if delta is greater than half the method displays high frequency dissipation characteristics which is desirable, but the global error will be of the order delta T, if suppose this is not 0 then the local error is of the order delta T square, global error will be of the order delta T, okay, so for the method to display numerical dissipation I need delta to be greater than half, in which case the term inside the parenthesis would not be 0, and the error the local error becomes order delta T square, and the global error will become order delta T, okay, so that is precisely what I was telling is that tradeoff between the desirable feature that the method should have high frequency dissipation characteristic, and also it should have a higher order accuracy, so both cannot be met simultaneously. So that takes us to the development of another method in which this weakness of Numer beta method is overcome, this is called HST alpha method, the HH, the alphabets HHT refer to the names of the 3 authors Hilbert, Hughes and Taylor, this is the reference that I have provided, and this is another related reference on the method, so the motivation of developing this method is to develop an unconditionally stable second order accurate method with acceptable high frequency dissipation characteristics, okay. So what did we see in Numer beta method? It was unconditionally stable alright, but there was a question of either achieving second order accuracy or achieving a high frequency dissipation, not both, now can we achieve both is a question, so indeed that is possible I will not get into all the details, I will just mention what exactly is done here, so in this method or in this class of method there are generalizations, what is done is in Numer beta method we used this representation for velocity and displacement and we introduced algorithmic parameter delta and alpha as we shown earlier. Now in HST alpha method a similar scheme is introduced in writing the equation of motion, so the MU double dot is written as a weighted sum of 1-alpha M into MU double dot N plus 1 and alpha M, MU and double dot, so this is an approximation to MU double dot at N plus 1, so alpha M 1-alpha add to 1, so it's a weighted sum, similarly the damping term is written in this form, and stiffness term is written in this form, and the external force is written in this form, so here we see now that we have now two more algorithmic parameters namely alpha M and alpha F, in addition to this delta and alpha, so that helps us to now examine make choices by suitably making choices of these parameters we can hope to achieve the objective of achieving unconditionally stable algorithm with second order accuracy and desirable high frequency dissipation characteristics. So we can run through this I will not provide all the details, so these three equations can be recast in a matrix form as shown here, and we can write this as AXN plus 1 equal to BXN plus FN, and we can write out the steps for implementation, it's not very difficult it can be done, if you follow, if you rearrange these three equation in this form we get you know these are the steps for implementation, we read algorithmic parameters I have to come to issue of how to select them I'll come there shortly, and read the initial conditions and evaluate the initial acceleration, and we form these matrices if delta T is constant we need to do this only once, and we evaluate this A tilde which is A inverse B which again we need to do only once if delta T is constant, now then set N equal to 0 and XN plus 1 is given by this, and increment N stop if you reach the final time instant otherwise simply go to step 6, that means the looping is only over step 6 and 7, so all other calculations are done outside this loop, okay, so this can be done later on we will consider some examples, but now we are focusing on qualitative features of the integration schemes, so let's now ask the question how about the stability characteristics of this method, so again the test system is a undamped single degree freedom system, so I write now the HST alpha you know equations for acceleration velocity and displacement, and we can assemble them in this form as we obtained earlier I will get AXN plus 1 is BXN, so the characteristic equation will be some AZ minus B determinant equal to 0, I use Eigen value as Z instead of lambda which I was doing earlier it doesn't matter, so the question is do the roots lie within the unit circle of the complex plane, that is a requirement for stability analysis. Now since this is a 3 by 3 equation we will have a cubic polynomial and that is given by this A naught ZQ plus A1Z square plus A2Z plus A3, now again we use the Jury's criterion so I have A naught A1 A2 A3 from the expansion of the determinant we can get these terms you can use symbolic maths you know, symbolic you know mathematics softwares to do this this is precisely how I have done this, the necessary condition will be P of 1 is P of 1 turns out to be omega square delta T square it must be greater than 0 that is satisfied and minus 1 to the power of 3 P of minus 1 which is given by this this must be greater than 0, so again we can examine this equation and see that one way to satisfy this is to take the parameters to satisfy these inequalities, there will be two sufficient conditions now modulus of A3 should be greater than modulus of A naught and modulus of B naught must be greater than B2, so we can examine all this and we can place now restrictions on the range of parameter values to achieve unconditional stability and that is summarized here. So the idea the important the details can be understood but the point is that we are now adjusting the values of algorithmic parameters to achieve unconditional stability, now I will illustrate some of these algorithms with simple single degree freedom system, so I start with a system with period 1 second that is omega is 2 pi and I am using Newmark's method in an acceptable region with 10 points in a cycle I get this solution, blue line is the exact solution and red is a Newmark's solution this is for undamped system and this is for damp system, so 10 points in a cycle doesn't seem to give good results but if you take 100 steps both give reasonably good solution, we will come to guidelines on selecting step size etc sometime later but we will just observe now outcomes of these simple numerical experiments. This is now HST alpha method again this is a parameter that I have taken which guarantees unconditional stability, 10 points in a cycle this is the result for undamped system and this is for damp system, now with 100 points in a second both the methods I mean the HST alpha method gives the result from HST alpha method agrees with the exact solution as shown here. Now this line of development of integration schemes to achieve you know desirable features in the schemes has been pursued by several authors and I just want to give an example of one such study, so in the HST alpha method if you go back we had a parameter alpha M to approximate the inertial forces and the same alpha F was used for damping stiffness and external forces, so one can relax that introduce more algorithmic parameters, so a study to illustrate that has been carried out by this person in 2008, so what he does is introduce alpha M here and alpha F here alpha K here and that alpha K is used here again, but even in writing the expressions for velocity and displacement the new parameters are introduced, now these become algorithmic parameters, seven algorithmic parameters, so now how do we select them, so that takes us into the question of what are the desirable properties of a integration scheme that is what we have been talking about, so typically we want at least second order accuracy that is at the global level and when applied to linear time variant systems we want unconditional stability and then controllable algorithmic damping in higher modes, so what we do is we examine the spectral radius as frequency becomes large, what we want is the spectral radius should drop, so that higher contributions from spurious higher modes are dissipated, this is desirable as I already said if response is dominated by low frequency mode, then for however for low frequency reach and the spectral radius should be 1, otherwise the contribution from lower modes also get distorted which is not desirable, we don't want numerical dissipation to affect the contributing modes, the significantly contributing modes should not be disturbed by the numerical dissipation, but the higher modes it helps if they are quickly dissipated numerical, there is another phenomenon that is observed is known as overshoot that is this is excessive oscillations during the first few steps, so HST alpha method has this deficiency and the further improvements on that has overcome some of these deficiencies, and we want the algorithm to be self-starting and we have not talked about nonlinear system but when we talk about nonlinear systems we don't want to solve more than one set of implicit equations at every time step, so the idea here is that by formulating the integration scheme in different ways and introducing different algorithmic parameters we can find regions in the parameter space where if these parameters are selected from those regions these objectives are met, so this has been the approach to developing these type of schemes, so in the next lecture what I will do is we will make some additional comments on these methods and we will try to apply these methods to some multi degree freedom systems just for illustration and see how they perform, so we will conclude this lecture at this stage.