 We have been discussing approximate methods for vibration analysis as a buildup to developing the finite element method, so in today's lecture we will be talking about what is known as Rayleigh-Ritz method and method of weighted residuals. Before we get into the details we will quickly recall what we have been doing, exact solutions to partial differential equations with inhomogeneous coefficients is in general not possible even for linear systems, so approximate solutions are invariably needed, right now we are discussing approximate solution strategies for finding natural frequencies and mode shapes and we will now move to force response characteristics as well and this discussion is a prelude to the development of finite element method. So we introduced during the previous lectures a quantity known as Rayleigh's quotient and for a discrete multi-degree freedom system we showed that Rayleigh's quotient U, U is a N cross 1 vector where N is a degree of freedom is given by U transpose KU divided by U transpose MU, this quantity is known as Rayleigh's quotient, this has units of radian per second whole square and if this vector U coincides with the nth eigenvector of the system then Rayleigh's quotient becomes exactly equal to the nth natural frequency of the system or the nth eigenvalue which is square of the nth natural frequency. We can also show that the Rayleigh's quotient is bounded between the first eigenvalue and the last eigenvalue, moreover we have also shown that if the vector U is in the neighborhood of the true eigenvector phi N that is phi N plus epsilon Y then we can show that the Rayleigh's quotient will be in the neighborhood of the natural frequency with an error which is order epsilon square. So we can also show that as we vary U, R of U reaches stationary values in the neighborhood of true eigenvectors of the system and it reaches its minimum value when U coincides with the first eigenvector, this is known as Rayleigh's principle. The main application of Rayleigh's principle is to estimate with simple approximations to the mode shape, the first mode shape, the first natural frequency that is the main application of Rayleigh's quotient. Now this is for discrete multi-degree freedom systems for a continuous system like actually vibrating bar or Euler Bernoulli beam we have shown that Rayleigh's quotient is now a function of the trial function phi of X defined over 0 to L and for a actually vibrating rod it is given by this and for Euler Bernoulli beam it is given by this. For distributed parameter systems the Rayleigh's quotient provides a bound as shown here there is no upper bound because for distributed parameter system the highest natural frequency tends to infinity therefore we can place a bound on Rayleigh's quotient as well. So again if as we vary phi of X whenever phi of X coincides with the true eigenfunction of the system then Rayleigh's quotient corresponds to the true eigenvalue of the system which is square of the natural frequency and again this function reaches stationary value whenever phi of X is in the neighborhood of the true eigenfunction and it takes a minimum value when phi of X is equal to the first eigenfunction. Now some illustrative examples, suppose we consider a 3-degree freedom system with K as shown here, M as shown here and let's assume that the kind of problem that we wish to solve is as follows for this system shown here the first eigenvector is obtained as shown here 1.13, 1.63 and 1.78 and the first natural frequency as 4.15 radian per second. Now due to a design modification of the modification the third mass is changed by about 15% that means this mass I increase by 15%. Now the question is how to estimate the change in the natural frequency, exact analysis but that is not what is intended here we wish to find an approximation to the change in the natural frequency. So for this system for the modified system the mass matrix is shown here so this M3 is increased by 15% that means this will go to this and K remains the same as it was earlier and what we will do now is since the first eigenvector is already given to us we will use this as the trial vector and compute the Rayleigh's quotient. So for this eigenvector and with the modified mass matrix the Rayleigh's quotient turns out to be 16.10 from which I get the modified natural frequency as this. Now so therefore the estimate of the change in natural frequency is about 3.3% so a 15% change in mass results in about 3.3% change in natural frequency as per this formulation. For sake of reference we can note that if you were to do an exact analysis of the modified system you will get the first natural frequency as 4.01 to 1 radian per second, so this is a one simple situation where Rayleigh's quotient can serve with some useful purpose. Now an example for a distributed parameter system so here I am considering a simply supported beam which carries a point mass M and also it is mounted on a spring K, if this K is 0 and M is 0 we know the exact solution to the natural frequencies of this system this is amenable for exact solution that is well known. Now the objective is to find characterize the first natural frequency of the system using Rayleigh's quotient. So we write the expression for the potential energy and the kinetic energy this is the strain energy stored in the beam and this energy stored in this spring this is at X equal to B therefore this is half K V square B, T and kinetic energy will be the kinetic energy stored in the beam which is this plus the kinetic energy stored in this point mass which is half M V V dot square at A X equal to A. Now if the system undergoes normal mode oscillations it is understood that the form of the oscillation is as shown here that means all points on the structure vibrate harmonically at the same frequency, so the form of the solution will be Phi of X into cos omega T minus alpha, so now I substitute that and compute the kinetic energy, so we can get the kinetic energy as a term inside the brace and cos square omega T minus alpha and similarly for this is the potential energy this is kinetic energy which is given here. Now the maximum value of V of T is reached when this cos square omega T minus alpha takes its maximum value and that is 1 because cosine function is lie, cosine function lies between plus minus 1, so cos square omega T minus alpha would lie between 0 and 1, so similarly T of T the maximum value will be taken when this sin square omega T minus alpha is 1, so the maximum potential energy is given by this, maximum kinetic energy is given by this and since the system is conservative these two quantities must be equal to each other and this leads to the definition of the Rayleigh's quotient as shown here. Now to proceed further we have to make a choice on Phi of X, this is the definition of Rayleigh's quotient for the problem, now we have to make the choice for the Rayleigh's quotient. What we can do is as I pointed out already if this mass is not there and if this spring is not there we know that sine N pi X by L for N equal to 1, 2, 3, 4, etc. are the set of exact Eigen functions. Now for this case we can as well take the trial function to be sine pi X by L or I can take a polynomial like X into L minus X or I can apply a UDL on this beam and find out the deflected profile and use that as the trial function, so there are many options that we can follow. Now the choice of this Phi of X determines how good this Rayleigh's quotient is as an approximation to the first natural frequency, no matter which choice you make this inequality will be honored, that means R of Rayleigh's quotient will be always greater than or equal to the square of the first natural frequency, so now suppose if you take Phi of X to be sine pi X by L now you get the natural frequency to be given by this, you can carry out this integration they are very simple you will get this. Now clearly when K equal to 0 and M equal to 0 the natural you know the frequency should revert back to the exact solution because the mode shape is exact in that situation so that indeed happens we get for K equal to 0 and M equal to 0 the square of the natural frequency as this which is which leads to the estimate of first natural frequency to be this, this is in fact pi square 1 by L square EI by M which is the exact solution. Now on the other hand if you now take the trial function to be X into L minus X so it satisfy the geometric boundary condition at X equal to 0 and at X equal to L but if you compute the second derivative Phi double prime of X you will see that this trial function doesn't satisfy the natural boundary condition that is the bending moment at the supports is not 0 according to this trial function but this doesn't matter in this application of the method because to compute Rayleigh's quotient at the highest derivative of Phi that I need is Phi double prime, so there is no problem here as far as application is concerned, so if I do that this integration is straight forward and I get the Rayleigh's quotient in this case to be this, and if in this case if K equal to 0 and M equal to 0 if you compute you will get this as 10.95 into this quantity which you know is higher than 9.87 into the same quantity that we got when we use the exact solution. Now for K not equal to 0 and M not equal to 0 to gain an understanding how this method works you can assign some numerical values so I have taken span to be 3 meters, X modulus to be 210 Giga Pascal so density A and B and beam cross section being rectangular with this width and depth and point mass is half of the beam mass and stiffness of the spring is this, and if I now use two trial function Phi of X equal to sine Phi X by L and Phi of X is X into L minus X, for the case when M equal to 0 and K equal to 0 we get the exact solution is possible that also I have tabulated, and this trial function now as I already said coincides with the exact solution therefore we should expect the exact solution which is as it should be. For X into L minus X I get 364.63 as an estimate for 328.51, so this is still a bound nothing wrong but this is as you can see this may not be an acceptable approximation. Now for nonzero values of M and K this solution turns out to be 270.31 and this gives an answer 2.97.42. Now we don't know the exact solution, so now the question now arises if you are using Rayleigh's quotient how do you improve the result? Suppose if you had started with this you have got 297.42 and by some logic suppose you move to this shape function you get a much better answer, right, but this is fortuitous you happen to know the exact mode shape therefore you are able to get this a much lesser value for Rayleigh's quotient than a simple function like this, but how do we systematically lower the Rayleigh's quotient, right, so that takes us into the next question that is how to lower the value of Rayleigh's quotient, this leads to the method known as Rayleigh-Ritz method. Now the main objective of Rayleigh-Ritz method is to lower the value of Rayleigh's quotient, so let's start with Rayleigh's quotient for the beam I am using the beam as an example an inhomogeneous beam with flexural rigidity varying with in space and mass per unit length varying in space. What we do is we expand the trial function in a set of linearly independent functions psi n of X, okay, this psi n of X are taken to be known functions whereas these A n's are unknowns, so this psi n are a set of known linearly independent function which satisfy all the boundary conditions, so let us start with this and let us revise this revisit this question as we go along but for moment we will assume that it satisfy all the boundary conditions. Now A n's are a set of unknown constant which need to be determined, now the strategy is now as follows, now by introducing a series representation like this in my definition of Rayleigh's quotient I have introduced a set of undetermined parameters, so I have a handle now to lower the Rayleigh's quotient, so the strategy that we will take is we will minimize Rayleigh's quotient with respect to this unknown coefficient, so select A n such that this Rayleigh's quotient is minimized, so that leads to suppose after I represent phi of X in terms of these n parameters now this R now is a function of n parameters, it's no longer a functional it's a set of n, a function of set of n parameters let us call it as capital Omega square, and this is written like this, now the condition for optimality is dou R by dou AI must be equal to 0 for I equal to 1, 2, 3 and n, so this simple calculation can be done as follows we can call this denominator as A and numerator as B, so I am asked to find out dou R by dou AI, A is a function of A1, A2, A3, B is a function of A1, A2, A3 therefore I use the rule of differentiation, so I get 1 by B square where B is the denominator into B into dou A by dou AI minus A into dou B by dou AI that must be equal to 0, now since I am equating this to 0 this B square is inconsequential so I can forget that, so I get the condition for optimality to be this, now if I divide in this case both sides by B I get A by B, A by B is nothing but Omega square, so this will be therefore dou A by dou AI into Omega square dou B by dou AI equal to 0, so that is what we are doing here for, now if I want dou A by dou AI I have to differentiate this with respect to AI, AI is one of the terms here, so you differentiate this become 2 into this summation into psi I double prime of X, so that is what is written here, and this 2 is because of this 2, similarly I get the differentiation with respect to the denominator this is equal to 0, now I divide by B and I get Omega square, and this I have to do for I equal to 1, 2, 3, and N, so it leads to capital N number of equations, so we denote now KIN as EI psi N double prime psi I double prime DX, and MIN as a term here as M of X psi N psi I X DX, now consequently this equation can be written as N equal to 1, KIN AN minus Omega square N equal to 1, MIN AN equal to 0 for I running from 1 to N, so this is a set of capital N number of equations in the unknowns A1, A2, A3, A capital N, so this can be put in the matrix form I get KA is equal to Omega squared MA, so this K is a N cross N, capital N cross N square matrix, M is a square matrix of the same size, and if you see here the expression for KIN and MIN, KIN is same as KNI, and MIN is same as MNI that would mean K and M are symmetric N by N matrices, so this A is N cross 1 column, so this KA equal to Omega squared MA represents an algebraic eigenvalue problem, K we call it as generalized stiffness matrix, M as generalized mass matrix, K is symmetric, M is symmetric. Now actually to get a representation like this we should, since phi is a function of X we should really be writing this as N equal to 1 to infinity, now since we are terminating this expansion at lower case N equal to capital N this would mean AK equal to 0 for all K greater than N, so to achieve this we need to supply external forces to the system which increases the stiffness of the system, consequently the natural frequencies get overestimated. Now we can achieve some computational simplicity's if we select psi N of X to be orthogonal with respect to M of X as shown here this is chronicle delta function then what happens is this generalized mass matrix will be an identity matrix, so in which case the eigenvalue problem will be the standard eigenvalue problem KA equal to Omega squared A, now once you find the eigenvalues you will get associated eigenvectors and using those eigenvectors we can also construct the approximation to the eigenfunctions which is phi K of X, AKN summation over AKN psi N of X, so this method the basic aspiration was to lower the Rayleigh's quotient, and as a byproduct we actually what we have done is we have ended up with discretizing the system as a N degree of freedom system, and from that we are able to get approximations to not only the first natural frequency but up to the capital N number of natural frequency and associated more shapes, so this is a byproduct of the method, so this is a way of this can also be viewed as way of discretizing a continuous system as a multi degree freedom system, discrete multi degree freedom system. Now we will examine this solution bit more carefully, suppose if I retain only one term so I get a single degree freedom approximation, so from that model I can get only an estimate to the first natural frequency, on the other hand if I now retain two terms I will get a two degree freedom approximation therefore I will get approximation to the first natural frequency and a new information that I will get is an approximation to the second natural frequency, similarly a three term expansion will lead to estimate of first and second natural frequency and also the third natural frequency, so on and so forth, so a capital N degree of freedom system will give approximation to the first N natural frequency. Now we can show that if omega N are the true but unknown natural frequencies these omegas that is this omegas that I am showing here have this bounding property, that means this eigenvalues converge from the above, so I will show with a numerical example what it means is as a degree of freedom increases the natural frequency approach the true natural frequency from above, how do I show that maybe we can work through this example once again, suppose the same problem a beam with a point mass and a spring and the numerical data as we had used in the previous example, so what I will do is now I will take the trial function to be a linear combination of sine N pi X by L, okay, now just again let us quickly recall when K equal to 0 and M equal to 0 sine N pi X by L is the exact nth mode shape for the system, nth eigenfunction of the system, so with K and M not 0 this is no longer the exact solution but we can use that as a trial function, so that is what we are doing. Now we will get now the eigenvalue problem as shown here where K ij we can work through the steps and show that K ij is given by this, M ij is given by this, and delta ij again is a chronicle delta that we can get, now the fact that for the first two terms here we are getting chronicle delta is not surprising because sine of X is the exact mode shape for the beam problem, so sine of X is orthogonal to EI of X that is sine double prime into psi K double prime EI of X integral 0 to L DX is chronicle delta function, so that is known, so that is why we are getting chronicle deltas here, but of course K ij and M ij are not diagonal because of these terms involving the discrete spring and the point mass. Now if we now take only one term I get 270.309 as approximation, if I now take two terms I get now for the first natural frequency an estimate of 266.024 and an estimate for second natural frequency as well, so three terms I get a slightly different answer for first natural frequency and a new information on third natural, so on and so forth if you continue I will list here for first five modes, so for a 5 degree of freedom system these five numbers are respectively the estimates of the first five natural frequency, now if you look at this column you can see that it is converging from the top, similarly this is converging from the top, so this happens for all the I get value, so this is one of the properties of Rayleigh-Ritz method and provided you use linearly independent trial function, this is a very useful property. Now we have been talking now till now about finding approximations to the natural frequencies and more shapes, but we could as well ask the question now how to get an approximation to the force response itself, so how do we approach that, so there is a class of methods known as method of weighted residuals, so for sake purpose of illustration I will consider again a beam problem, a beam with inhomogeneous flexural rigidity, mass per unit length and damping properties driven by a distributed loading F of X, T, now let's assume that initial conditions are initial displacement is this, initial velocity is this, and appropriate boundary conditions involving geometric and natural boundary conditions have been specified, which could be simply supported cantilever, fix-fix so on and so forth, now in this class of methods the starting point is we assume the solution to be represented as a series N equal to 1 to capital N, A N of T phi N of X, this A N of T are the unknown functions of time to be determined, we don't know that, these are called generalized coordinates, this phi N of X are the known trial functions, now what condition these should satisfy we will see as we go along, so what I do now I substitute this assumed form of the solution into the governing field equation, and if I substitute here I will get this is the field equation and this is the assumed solution, now I have substituted here, what happens is this equation won't be exactly satisfied, we don't expect that this equation will be exactly satisfied and we will be left with an error, if indeed this is an exact solution this error will be 0 for all X and T, but that is not going to, I mean that's not the situation that we are discussing, so there will be a error or this is also known as residue, for the exact solution this error is 0 for all times and for all X in 0 to L, for an approximate solution as I said this is not so, now what is the strategy now, we have to still select A N of T which are not known, the strategy that we follow is we select these unknown functions such that a specified measure of the total error is minimized in some sense, we can't expect E of X, T to be 0 for all X, suppose you fix a time T and if we can find A N so that E of X, T is 0 for all X then that's the exact solution, but that is not what we are doing, we want to minimize a measure of this error, a measure of total error, but this measure itself has no unique interpretation, there is no unique sense of minimization because there is no unique sense of criterion of total error and how to minimize it, depending upon the sense of minimization we have different methods, method known as method of least squares what we do is this is the residue which is given as shown here and we define what is known as total mean square error which is epsilon of T which is integral 0 to L E square of X, T DX, we are squaring because we don't want positive and negative errors to cancel, so we could have taken the integral of error over 0 to L, but if for certain parts of the beam error is positive and for certain other parts of the beam it is negative we get a false impression that total error is 0, so we square that and find the total error, now we select A N of T so that this epsilon of T is minimized with respect to A N for N equal to 1 to N, so this leads to a set of capital N number of equations for the generalized coordinates A N of T, so if you do that so you will get 2E dou E by dou XT DX equal to 0 that 2 is irrelevant because on the right hand side is 0, so this is for 1 to N, so this can be written as in general WN X, T E of X, T DX equal to 0, where WN X, T is this function, okay for N equal to this, now if you write all these equations you can show that the resulting equation will be of the form MA double dot plus CA dot plus KA equal to P of T, this P of T will be of course integral of F of X, T multiplied by this quantity and integrated 0 to L and that need not have units of force, this function need not have units of force, so the interpretation of these matrices also needs to be carefully done, they need not have the traditional units and they may or may not be symmetric, okay, so that also we have to bear in mind. Now let us see if I can find the expression for that the weight, so what I am supposed to find is dou E by dou A N, right, so dou by dou N of the residue you can see that it is given by EI phi N double prime of X the entire thing differentiated twice, this is the weight function in this method that is WN, so the criterion for minimization of the error is given by this for N equal to 1 to N and this is independent of time, so W is although I wrote it as X, T it is actually W of X only, so this is MIJ is given by this, CIJ is given by this, PJ is given by this and KIJ is this, so you can see here MIJ is not same as MJI and KIJ is of course symmetric because of the way these terms appear here but it is not true for CIJ, and P of T is may not have units of force, okay, so these are some issues that you may have to bear in mind when you are solving and interpreting this solution. In another strategy what is known as method of collocation, what we do is we again look at this residue which is as before E of X, T which is this, we select this N of T such that this error is 0 for a given T this error is 0 at N chosen points along the beam, that means I select X1, X2, XN and demand that at all these capital N number of points error is exactly equal to 0, so that requirement leads to capital N number of equations for the generalized coordinates and that is the governing equation. Now we can see that this is equivalent to taking the weight function to be a Dirac-Delta function, so this is same as Dirac-Delta of X-XN E of X, T DX equal to 0, that means the statement that E of X of T 0 for XN equal to XN, N equal to 1 to N is equivalent to this. Now that is we can write this criteria as integral 0 to L WN X, T it's actually WN of X because there is no time here, this time can be deleted from this E of X, T DX equal to 0 where WN is this, now this leads to again a set of ordinary differential equations which are coupled and you can verify that these matrices will be unsymmetric, non-diagonal and this will be the equation, now the problem with this method is although the error is exactly 0 here we don't have an idea how small or large the error is at an X which doesn't coincide with any of these X, so this we don't get any idea on that so I have to be careful in using the collocation method. So the details of the matrices are obtained here, so this is the criterion for finding the generalized coordinates and MIJ turns out to be MXI Phi J XI and which is obviously not equal to MJI and similarly CIJ is obtained as this and KIJ is obtained as this, this is again not symmetric, PJ of T is F of XJ, T, so this is how we get this equation. A modified version of this is what is known as sub-domain collocation method, it is somewhat similar to collocation method but the criteria is that we again divide the beam into capital N number of segments but instead of demanding that at all these points the error is 0 we demand that over each of these segments the average error is 0, that is a total error over this segment is 0, total error over this segment is 0, so I write this as integral XN-1 to XN E of X, T DX equal to 0 for N from 1, 2 to capital N, where X naught is 0 XN is capital L, now this can be rewritten as integral is 0 to L, WN of X, T it is WN of X, E of X, T DX equal to 0 where WN is a box function written in terms of hue side step function, so for this method sub-domain collocation method the weight function is a box function, okay, so this is a modified version of collocation method, there is another method known as Galerkin's method this is going to be used often in our course, here the idea is the weight function is selected to be the trial function itself, so we are expanding V as in terms of phi N of X and the weight function we take it as phi N of X itself, so this leads to that means error is orthogonal to the trial function, so this leads to capital N number of equations and we will see now what is the form of this equation, but from this statement we can see that we are going to get a set of N ordinary differential equations, so how what exactly happens here, again let us start with the residue term, so the Galerkin's statement is this and you can see here if you implement this M I J will be M of X, Phi of X, Phi J of X DX and this is M J I, C I J is C of X, Phi I of X, Phi J X this is C J I and K I J is actually this we can see what happens to, this is a fourth derivative present here and this is there is no derivative on this we will see what happens to this term, but in any case this is the K I J, now P J of T is the forcing function, now what happens to K I J, so now this if you integrate by parts you do it twice, so upon each integration one of these derivative passes on to the weight function, so after you integrate by parts twice you get the integral 0 to L we will have Phi I of double prime of X, Phi J double prime of X, now if we take all trial functions to satisfy all the boundary conditions then these quantities inside the brace these two terms will be 0, because these encapsulate the class of admissible geometric and natural boundary conditions and this will be 0 and we get K I J to be given by this, and K I J in this case is symmetric, so in this method I am getting again equation of the form M A double dot plus C A dot plus K equal to P of T, but the advantage is MCK are symmetric and there are further advantages I will come to that as we go along it is something to do with the demands on differentiability of the trial functions, now how do we get initial conditions, we have got now the governing equation for the generalized coordinates we need to find A of 0 and A dot of 0, so how do we get that, so this is a assumed solution V of X, T is A N of T Phi N of X therefore at T equal to 0 I get this displacement at T equal to 0 the velocity is given by this, so the unknowns are A N dot of 0 and A N of 0, so this in collocation method again we can demand that the error is 0 for at N points and I get these equations for I equal to 1 to N and you can solve this equation to get the required initial condition, similarly this is the equation for the initial displacement initial velocity, in the least squares and Galerkin methods A of 0 is given by the inverse of this matrix into the integral over the initial condition as shown here, so this again requires some effort to evaluate the initial condition, so let me summarize what we have talked now, there is a least squares method, collocation method, Galerkin's method, subdomain collocation and petro Galerkin, there is one more thing known as petro Galerkin I didn't talk about it, here it is similar to Galerkin technique except that the weight functions are taken to be another set of cyanophics, this cyanophics is different from phi N of X, that technique is known as petro Galerkin, so in all these approaches you can see here one of the common term that is present in all these expressions is E of X, T, what distinguishes one method from the other is what multiplies this residue term in the integrand, for least squares it is this, for collocation it is this, for Galerkin it is phi N of X, for subdomain collocation it is this box function, for petro Galerkin it is psi N of X, so all these methods nevertheless can be put as a weight function into a residue equal to 0 for N equal to 1 to N, so this class of methods is known as method of weighted residues, now to implement these methods there are certain choices that we have to make, first set of N trial functions, then a set of N weight functions we have to make, choose, now issues when you are selecting trial function or the weight function order to which the functions need to be differentiable, and how well they satisfy the boundary conditions, and also how many terms need to be included, so we will revisit these questions but bear in mind these questions need to be addressed, so the general format is we have the field equation and we assume this solution and after applying the method of weighted residual in any one of this form I get this set of equation, of course the meaning of MCK and P would be different depending on which method you have used, so what does method of weighted residuals achieves, a partial differential equation governing the behavior of a continuous system has been replaced by an equivalent set of ordinary differential equations which are initial value problems with a view to obtain an approximate solution, so that's what we are doing, so this is a way of discretizing an infinite dimensional system by a finite dimensional system, so in fact finite element method also achieves this but in a somewhat different way as we will see shortly. Now the trial functions there is a scheme of classifying trial functions depending on the extent to which they satisfy the boundary conditions, see this is a field equation for example free vibration problem this is a field equation, and we have boundary conditions one on geometric quantities, another one that is a force quantities they are bending moment and shear force call it as set 1 and set 2, so in normal mode oscillation this is a solution that we assume and this is a Eigen value problem that we get. Now we say that a trial function is admissible if geometric boundary conditions are satisfied, we need not have to satisfy the natural boundary condition nor the field equation, by field equation I mean this, by comparison functions we mean those trial functions which satisfy both geometric and natural boundary conditions, but the field equation is not satisfied, on the other hand an Eigen function is a exact solution to the problem therefore it satisfies the geometric boundary conditions, natural boundary conditions as well as a field equation, so as you see here the least that we should expect from a trial function that it should be admissible, okay, so we will see what is the implication of this as we go along it is also tied up with the requirements on differentiability of the trial functions, so if a function is admissible the requirements is much less whereas if it is comparison function where you need to satisfy conditions on bending moment and shear force the requirements on differentiability increases and the class of functions that we can select shrinks, so ideally we would like to deal with admissible functions so that we have a huge class of functions to select from, as we go in this direction that choice of trial functions narrows, now let me just go back to the methods one by one and see what type of requirements each method imposes, in Rayleigh's quotient method I need to evaluate only the second derivative of the trial function, so they must be admissible, okay, and possess nonzero derivatives at least up to the second order, okay, in Rayleigh-Ritz method again the trial functions need to be at least admissible and possess nonzero derivatives up to the second order because if you see the implementation of Rayleigh's method I need to compute this Kij and Mij they don't require functions beyond second derivative, so trial functions need not you know we are happy to use trial functions which have only second order derivatives exist, okay, so this is required whereas in least squares and collocation trial functions need to be comparison functions, okay, you can verify that you need derivatives up to the fourth order, in Galerkin the trial functions need to be comparison functions, we will see some of this. Now there is yet another approach which is slightly different from what we have been discussing that is known as assumed mode method that employ Lagrange's equation, this is slightly different, so again let us consider for purpose of illustration the free vibration problem of an inhomogeneous beam, so we have the field equation and appropriate boundary conditions and initial conditions, so we have seen that solution of this differential equation is equivalent to minimizing this functional, okay, this we have seen, this is Hamilton's principle, so the idea of this assumed mode method is to approach the problem of minimizing the action integral at the stage of developing the approximation, that means I won't develop this partial differential equation but I will develop an approximate method based on the minimization of action integral itself, so what does that mean? That means approximate solution can be developed either by addressing the field equation or by optimizing the action integral, so either we can adopt this or this, so for purpose of illustration suppose if I take a beam lie simply supported beam the approximate solution is let us take it to be N equal to 1 to N of T phi N of X, where N of T is set of unknown generalized coordinates and phi N of X are known trial function taken to satisfy all the boundary conditions say, for example phi N of X is sine N pi X phi, now this is an expression for kinetic energy this is an expression for potential energy, so now I will form the Lagrange, Lagrangian for this, this is this in terms of these generalized coordinates and their derivatives, and on this I will apply the Lagrange equation N number of times, if I do this, this is simple calculation we can get the mass and the stiffness elements of mass and stiffness coefficients and we can get the equation M A double dot plus K equal to 0 with certain initial condition, this K IJ that is K NK is EI phi N double prime phi K double prime which is symmetric in this case, similarly if you look at M and K it is symmetric, so this is similar to what we got in Galerkin's method, so we can see here in this case M and K are N by N matrices and these are the generalized coordinates with N degrees of freedom both M and K are symmetric actually we can show that M is positive definite and K is a positive semi-definite, one of the limitations of all the approaches including this one is that these generalized coordinates, this vector of generalized coordinate does not have direct physical meaning, they need to be substituted into a series expansion and evaluate that series then only it emerges as a quantity to which we can assign a physical meaning like displacement or velocity things like that. Now clearly this M and K depend on choice of this trial function, so again in this approach also a partial differential equation governing the behavior of a continuous system has been replaced by an equivalent set of ODEs with a view to obtain approximate solution. Now in discussion of variational approaches to solution of partial differential equations there is a, the set of terminology that we should become now familiar, I will introduce them now, we talk about what is our strong form or a weighted residual form and a weak form of the governing equations, so what they are? Now the discussion on these forms of equation is basically in the context of finding approximate solution in this form V of X is N equal to 1 A N phi N of X, it is in this context that we talk about these different forms, in the strong form we use a field equation itself that means we write the governing partial differential equation and want to develop the approximate solution starting from this, so this field equation plus appropriate geometric and natural boundary conditions is used to develop approximations based on strong form. In the weighted residual form what we do is operationally what we do? F of X is taken to this side and we multiply by a weight function and demand that this integral is 0 is 0, so this is like virtual displacement in structural mechanics problems, okay, so I will come to that in due course. So the idea is we can select a set of N weight functions and obtain governing equations for the unknowns A N of T, A N, in this case I am considering a static problem it is A N so only, so the trial function phi N of X here you know if you are using this method you see in this method you would need V differentiated 4 times, that means you need D4 phi N by DX4, okay, whereas W of X there is no such requirement, it can be, it should be integrable but there is no requirements on its derivatives. Now in weak form what we do is we carry out integration by parts of this weighted residual statement and see if I can pass on some demands on differentiability of these trial functions to the weight function, what does that mean? Suppose let us do that integration by parts, I get, you do integration by parts twice and I get here if you carefully see one of the term is W double prime X EI of X V double prime, this is P V prime W prime K VW minus F of X into W of X DX equal to 0 plus these boundary terms. Now we know that EI V double prime is nothing but the shear force I can call it as F2, so this is W of L, this is W is here and this is shear force W of 0 bending moment at the support into W1 prime of L and this, plus this, this equation is known as the weak form of the governing equation. Now in this case if you carefully see the demand on differentiability of trial function is only up to the second order, the original equation was fourth order, now trial function need to be differentiable only up to the second order, but the weight function now need to possess derivatives up to second order, that means the demands on differentiability on V of X, the trial functions has been weakened and that has been passed on to the weight function, in that sense it is a weak form, that means demand and continuity goes down on the trial functions. So for purpose of illustration let us consider the beam be simply supported at X equal to 0 and X equal to L and it be subjected to bending moments M1 and M2, so the weak form of this equation is this, now M1 and M2 are the applied bending moments, so they will substitute now for that is M1 and M2 will appear there, but now the beam is simply supported therefore what I will do is I will demand that this W of L and W of 0 which is a weight function must be 0, at X equal to 0 and X equal to L, so if I do that I will get the weak form in this, okay, so that means the inhomogeneous boundary conditions now explicitly appear in the weak form, so now we have to select Phi of X and weight functions so that V of 0 is 0, V of L is 0, because that is a geometric boundary conditions because it is simply supported beam at X equal to 0 and X equal to L and this weight function also need to, we take that this weight function also need to satisfy this, it need not, this need not be 0 therefore this condition will also be satisfied, okay. So again let me emphasize that the trial functions are here in this approximation need to be differentiable only up to the order 2, whereas in the weighted residual form they need to be differentiable up to order 4 whereas there was no demand on weight function, but here the demand on differentiability of trial function and weight function are balanced, okay, so this is the advantage of the method, so if trial functions need to be only differentiable up to the second order the pool of function from which we can select the trial function becomes very large, okay, as compared to those function which need to be differentiated up to fourth order, okay, that is the main advantage of this method. So now there are certain other terminologies that we should be aware of, so if we go back to the problem of beam on elastic foundation carrying axial load I think this is a beam carrying axial loads we got these two terms in the variational form, so the coefficients of weight function we call them as secondary variables, okay, the coefficient of W of X is this and this is declared as a secondary variable. The primary variable V is written in the same form as what is appearing here, that is W of L means V of L, this is a primary variable, okay, so in this case for example this will be the secondary variable this is a primary variable, similarly at the other end the bending moment becomes a secondary variable, the primary variable is a slope, right, because this is the function that multiplies the gradient of the weight function, so that coefficient of weight function or its derivative are called secondary variable, this is the primary variable is in the same form as what the weight function appears but W is replaced by V, so this is the classification, so that means the dependent variable expressed in the same form as the weight function appearing in the boundary term is known as primary variable, now coefficients of weight functions are called secondary variable, so we will make now a few remarks, the weighted residual statement in the weighted residual statement as I already said the trial functions phi N of X must be differentiable up to the fourth order, W of X on the other hand which is the weight function can be any integrable function, so equations for the undetermined coefficients can be obtained by choosing a set of N weight function, so weight functions come from a very large pool whereas trial functions come from a very small, relatively smaller pool, in the weak form on the other hand we distribute the requirements on differentiability evenly between weight function and the trial function, thus the continuity requirements on the trial functions is weakened and hence the name weak form, we have larger pool of functions to select from for the trial functions in this case, the trial function therefore need to satisfy only the geometric boundary conditions as the natural boundary conditions are included in the weak form, for example I got a bending moment term that appeared explicitly in the weak form, right, so and also we talked about primary and secondary variables we will talk about this again in a later part of the course, the number of primary and secondary variables are the same whenever we consider even ordered you know partial differential equations which often is what we do, they appear in pairs, for example translation and shear force, slope and bending moment and so on and so forth, only one item in each pair can be specified at the boundary, for example at the simply supported end we can specify displacement to be 0, we cannot say displacement and shear force to assume certain values that is not on, similarly if you talk about slope bending moment is 0 and you cannot say slope is also 0, okay, so only one item in each pair can be specified at the bond. Now in the methods that we have discussed so far one of the common features that we have seen is that the trial functions are valid over the entire domain of the structure being studied, that means these are globally valid shape functions or trial functions, this is okay if you are talking about simple geometries like line elements that we are considering but for more complicated problems it is not easy to construct global trial functions, one more difficulty that we have is the generalized coordinates N of T they do not have direct physical meaning, unless I substitute in this series and sum it up I would not know only then I would arrive at a physically meaningful quantity, if I say A3 of T is at some time instant point 02 it doesn't specify much, you know it is a, we cannot interpret directly in any useful way. So we would like to deal with not simply line elements we would like to deal with trusses like this, frames like this and you can imagine constructing global trial function for a truss like this would not be easy it will be impossible. Similarly constructing a trial function for a much less, you know much more modest structure like this also is not easy. So in the next class we will see how we can address these limitations and we will start talking about the development of finite element method which aims to overcome the difficulties that I just mentioned, so with this we will conclude this lecture.