EXACT ELASTIC STABILITY ANALYSIS BASED ON DYNAMIC STIFFNESS METHOD

The article is dedicated to the discussion on the exact dynamic stiffness matrix method applied to the problems of elastic stability of engineering structures. The detailed formulation of the member dynamic stiffness matrix for beams is presented along with the general guidelines on automatisation of the assembly of member dynamic stiffness matrices into the global matrix that corresponds to the whole structure. The advantage of the dynamic stiffness matrix in case of parametric studies is explained. The problem of computing the eigenvalues of transcendental matrix is addressed. The straightforward approach as well as a powerful WitrickWilliams algorithm are discussed in details. The general guidelines on programming the DS matrix method are given as well.


1.
Introduction Elastic stability problems faced by engineers these days feature extremely high levels of complexity. Usually, the structures of interest represent multilevel hierarchies of substructural assemblies. Finite element analysis has lately shown a considerable potential in solving such complex problems. However, the FE method has several disadvantages. The errors associated with the numerical approximations of the stresses and strains within a finite element using numerical shape functions are well known (Bathe 1996). Another weak point of the FE method is computational efficiency. Although, the effectiveness of programming of the FE procedures in commercial software has been led to perfection over the last thirty years the necessity to assemble large stiffness matrices and to perform costly operations like inversion is inevitable. These major drawbacks of the FE method give a strong background to look for any alternatives that can be used either in parallel or separately from the FE analysis.
This article is dedicated to the discussion of the exact Dynamic Stiffness method and its application to the problems of elastic stability. This method is not new and has been known for many decades. The exact date of the development of this method is difficult to trace. However, it is assumed to start from works by Professor John H. Argyris (Spalding 2014) who was the first to introduce the method of assembly structural components into a general stiffness matrix. The essence of the method is in writing exact stiffness relations for a single member of a structure in a general form and then assemble such equations into a matrix. Then the member matrices are assembled into a global structural stiffness matrix. The eigenvalues of the global matrix yield the critical buckling loads of a whole structure.
The member stiffness matrices derived from the exact stiffness relations are known as dynamic stiffness (DS) matrices (Williams et al. 2004). However, the name of such matrix was taken from the nomenclature widely accepted in the theory of vibrations, the approach is highly universal and can be readily applied to the problems of elastic stability. The member DS matrices as well as the assembled DS matrices are transcendental which means that their entries are transcendental functions of eigenvalues and define exact relations between nodal forces and nodal displacements of each member in a structure. In such a way the DS matrix method possesses considerable advantages compared to the finite element (FE) method in terms of computation efficiency and accuracy. However, the transcendental nature of a DS matrix significantly complicates the process of computing its eigenvalues. For this purpose an iterative algorithm capable of converging on any number of eigenvalues with any precision was developed by Wittrick and Williams and has their names (F. Williams and W. Wittrick 1970;Wittrick & Williams 1971).
This paper is dedicated to the detailed study of the DS matrix method in its application to the problems of elastic stability of structures. It contains numerous guidelines and tricks on automating and programming the method. We discuss here the formulation of the member DS matrix starting from member governing differential equation (Section 2). Then we proceed with the analysis of assembly of the global DS matrix and propose essential steps to program it. Finally, we discuss the methods of computing the eigenvalues of the DS matrix on a simple illustrative example. We study both the straightforward evaluation of the zeros of the determinant of the DS matrix and the Wittrick-Williams (W-W) algorithm. Although the W-W algorithm has been considerably studied in dynamics (Wittrick & Williams 1971; F. Williams and W. Wittrick 1970; Banerjee & Williams 1986) we believe that it requires a more in-depth analysis from the perspective of how it can be programmed for the case of elastic stability. We also believe that illustration of the application of the algorithm on a simple example is essential.

Exact dynamic stiffness matrix for buckling of a beam member
Consider a standard Sturm-Liouville differential equation in a form (Arfken et  3) with appropriate boundary conditions defines the problem of elastic stability of a structural member represented by an elastic beam. It can be discretised according to the FE method and solved for eigenvalues of a generalised eigenvalue problem that define the critical buckling loads. Alternatively, this problem can be solved exactly employing the DS (Dynamic Stiffness) method. Here we proceed with the DS approach to the problem. We derive the dynamic stiffness matrix for the case of buckling of a beam member following the general guidelines published by Banerjee (Banerjee 1997 (2.5). Therefore, the system of four equations in terms of unknown coefficients is obtained. This system of equations relates displacements at the nodes of a member to the unknown constants. Written in a matrix form it reads: (2.14) This matrix is called the exact Dynamic Stiffness Matrix. The word "Dynamic" came from the origins of the DS matrix method from the theory of vibrations. Since then the method has been employed in various fields but the name has been kept unchanged. The entries of the matrix were computed using equation (2.14) and is presented in the simplified form using the notation from the paper by Williams and Wittrick (Williams & Wittrick 1983) who in turn formalised the elastic stability functions computed and tabulated by Livesley  (2.16) The member DS matrix shows some similarities with a well-known numeric member stiffness matrix employed in FE analysis. The FE stiffness matrix of a member usually contains two different types of degrees of freedom (DOF). The first type of DOF corresponds to the nodes of each element and the second type of DOF corresponds to the nodes of a whole member. Here by member we understand a single complete beam that can be discretized with any arbitrary number of elements. According to Wittrick and Williams (Wittrick & Williams 1971) in the limit case when the number of finite elements of a structure tends to infinity the so-called condensed dynamic stiffness matrix can be located within the FE numerical stiffness matrix. The condensed stiffness matrix corresponds only to member nodal displacements and approaches the DS matrix. Thanks to this fact, the entries of the DS matrix can be regarded as limits to which the series of the infinite numerical stiffness matrix converges.
The DS matrix is square, symmetric and positive definite similarly to the FE stiffness matrix because of the way how it is formulated from the self-adjoined operator equation. The most important property of this DS matrix is that it is exact because it simply represents the governing differential equation (2.3) without any discretisation or approximations. Thus, the eigenvalues of the DS matrix correspond to the exact values of the critical buckling loads. Another important property of the DS matrix is that it can be assembled into a global DS matrix following the procedure similar to the assembly of numerical stiffness matrix implemented in FE analysis.

3.
Assembly of the complete global dynamic stiffness matrix Consider a two-span beam compressed by an axial loading scathed in Fig. 3.1. Each span has different stiffness and is represented by a member shown in figure 2.1. We write stiffness equations as relations between nodal forces and nodal displacements for each member and assemble them into a system of equations: Here and are the vectors of nodal forces acting on each member and and are the vectors of nodal displacements respectively. The relations between nodal forces and nodal displacements for each member are defined with member DS matrices K and K . Any system of equations can be assembled and presented in a matrix form: .
(3.2) Here and are the assembled force and displacement vectors and is the global DS matrix of the structure. Equation (3.2) can be rewritten: . (3.3) Here and are the submatrices of the global DS matrix that correspond to the member DS matrices that are overlapped. Overlapping is explained by the fact that degrees of freedom of the coincident nodes (nodes 2 and 3 on Another important property of the global stiffness matrix is that it can be appended with other stiffness matrices that contain essential information about the structure. For instance, the two-span beam may be installed on the elastic foundation (see Fig. 3.2). In this case, the complete global DS matrix consists of two matrices: Here is a complete global DS matrix and K is a matrix containing stiffness of elastic foundation as well as stiffness of artificial coupling between members. In this case, the K and K matrices do not overlap because they do not have any coinciding degrees of freedom. On the other hand, the coupling matrix K provides the information on a coupling between members and ensures continuity if displacements. The matrix relation between nodal forces and nodal displacements thus reads as: (3.5) Here, submatrices K K and K K correspond to the member stiffness plus the stiffness of the elastic foundation, while K and K (off-diagonal terms of the global DS matrix) define stiffness of coupling between members. Matrix defined in equation (3.5) is a complete global DS matrix of a structure. In general, this matrix contains enough information to describe the majority of beam assemblies including plain frames with various clamping conditions.
After the complete DS matrix of a structure is assembled it can be studied for eigenvalues that yield the critical buckling loads of the structure. It is known (F. Williams and W. Wittrick 1970) that buckling occurs when compressing load reaches its critical values that correspond to the eigenvalues of the following eigenvalue problem: . (3.6) In matrix notation: .
Concluding the analysis presented in this section several key points should be underlined. The DS matrix is a matrix that contains exact information about each structural member stiffness and is derived from the governing differential equation. The DS matrices of each member can be assembled in a global structural DS matrix which contains exact information about the stiffness of a whole structure. The complete DS matrix contains exact information about the stiffness of an arbitrary complex structure provided it is assembled from members which member DS matrices are known. The eigenvalues of the DS matrices directly correspond to the critical buckling loads. The DS matrix is transcendental and contains all the eigenvalue spectrum of the continuous structure.
It is now apparent that DS matrix method is a powerful and universal method that can be used as an alternative to the FE method as any frame structures or, generally speaking, any structure that consists of assemblies of members for which member DS matrices are known can be addressed by this method. This method is exact; thus the common errors associated with the numerical discretization of a continuum are avoided. It is also computationally efficient as the global DS matrices are usually much smaller than those assembled in the FE analysis. It can be explained by the fact that DS matrix is defined for an entire member rather than for a finite element. For example, in the case of a beam member, the DS matrix has only four degrees of freedom that correspond to each node. In FE analysis the corresponding numerical stiffness matrix is considerably larger than the DS one because the beam should be meshed with multiple elements to reach acceptable accuracy. Dependence of the computation accuracy of the FE analysis on the number of elements is a drawback compared to the DS matrix method.
The DS method, however, has its limitations. By virtue of the fact that the DS matrix is exact and possesses intrinsic transcendental dependence on its eigenvalues, the computation of the eigenvalues of equation (3.7) has been an almost unresolvable problem for a long time. This fact along with the difficulties in the assembly of the DS matrices for complex structures have been the major obstacles in development and application of the DS matrix method for decades since it was first introduced. These two issues are discussed in details in the next sections.

4.
Eigenvalues of the exact dynamic stiffness matrix Computation of the eigenvalues of a DS matrix is a known problem (Wittrick & Williams 1971;Williams 1981;Wanxie et al. 1997). Before the computationally efficient Wittrick-Williams algorithm was developed the standard procedure for computing eigenvalues of transcendental matrices was quite inefficient and unreliable (Williams et al. 2002). With the increase in computational capacities of modern computers, the straightforward calculations of the eigenvalues have become more reliable and can be automated to some extent. However, the danger of missing some eigenvalues still exists. In the first subsection of this section, we discuss the issues with the assembly of the DS matrices and the straightforward approach to computing its eigenvalues. The second subsection is dedicated to the Wittrick-Williams algorithm, its strong and weak points.

Computation of eigenvalues of the parameter dependent DS matrix:
Consider the member DS matrix given by equations (2.15) and (2.16). It is apparent that practical assembly of a considerable number of DS matrices into a global DS matrix is associated with mathematical difficulties as the entries are not float point numbers but some transcendental functions. One solution to this issue would be to use a numerical form of the DS matrix and apply Wittrick-Williams algorithm (W-W algorithm). This option is discussed in details in the next subsection. The second option is reasonable when numeric form of the DS matrix is not convenient. Such situation is possible when an exact analytical dependence of an eigenvalue on some parameter is required.
If the DS matrix is relatively small, i.e. from one to five members, the modern symbolic mathematic software is capable of assembling it in almost no time. Each member of the DS matrix is encapsulated within a symbolic function that is then passed to a standard assembly process. This process is illustrated in Fig. 4.1.1 for the case of a two-span beam.
The functions and contain expressions for the entries of a member DS matrices given by equation (2.16). Thus, the global DS matrix contains the same functions or their combinations (overlap region) with altered labels defining their positions in the global matrix. It is often required to consider some other means of coupling between members such as springs. Then the coupling matrix should be assembled separately following the same procedure. Thereafter, the complete global DS matrix is assembled according to equations (3.4) and (3.5).
One of the most important strong points of the DS matrix method is that any structural parameter such as member length, stiffness etc. can be regarded as a parameter of the system in general. Because the way how DS matrix is formulated using the exact symbolic functions it can be easy adapted to perform parametric studies. Each of functions and can contain an arbitrary number of meaningful physical parameters that describe the structure. Then, it is possible to derive exact analytical dependencies of the eigenvalues of the structure on these parameters. In the case of parametric studies, the straightforward approach to computing the eigenvalues of the DS matrix is desirable.
Consider transcendental eigenvalue problem 3.7. Let us assume that the complete global DS matrix is parameter dependent. It is easy to come up with some practical reasoning for such matrix to exist recalling Fig. 3.1. For instance, assume, that the coupling spring stiffness is the governing parameter of the system. Thus, matrix would have some functional dependence on two parameters, e.g., its eigenvalue λ and a structural parameter . Then the eigenvalue problem 3.7 will be written as: , . This matrix numerically defines a surface of the determinant of the DS matrix with respect to its eigenvalue λ and the structural parameter . This surface should be intersected with the zero plane to obtain parametric eigenvalue curves . It can be done in a computationally efficient way according to the following steps. Firstly, some near-zero tolerance ε should be chosen 0. Secondly, the selection algorithm is run on the entries , , 1, with the condition | | . Therefore, only the entries that are close to zero are selected. The values , that correspond to are selected from and vectors. Therefore, a set of unique correspondences → for which | | 0 is defined and can be plotted. Such plots exhibit parameter dependent eigenvalue curves of a structure. In such a way, various important properties of the structure with respect to the change of governing parameters can be studied.

Wittrick-Williams algorithm:
While the straightforward algorithm for computing the eigenvalues of the DS matrices has some advantages it is not suitable to compute eigenvalues of large DS matrices. The reason is that it is concerned with large symbolical operations that become extremely costly for large DS matrices. Another drawback of the straightforward approach is that it cannot guarantee that exactly all eigenvalues are computed within a chosen interval and none of them is missed. Therefore, this approach requires much of a human's assistance and is limited in automatisation.
On the other hand, the Wittrick-Williams algorithm is known for its ability to locate all the eigenvalues of the DS matrix in the given interval with arbitrary precision without missing a single one (F. Williams & Williams 1971). Along with the ability to be automated the unprecedented computation efficiency makes W-W algorithm a good option in the case when numerous eigenvalues of a large DS matrix need to be precisely computed (Williams et al. 2004). However efficient and powerful, the W-W algorithm has not received enough attention during the decades since its creation. It is still an option for researchers rather than for large industrial applications due to the complexity of its formulation and relatively small amount of studies focused on its application. However, we believe that this tendency can be changed. When properly studied and programmed along with the straightforward method discussed in the previous subsection, this method can hugely benefit the DS matrix approach and make it a strong competitor to the FE analysis. Thus, this discussion is dedicated to the analysis of the W-W algorithm along with a detailed explanation of the steps that need to be undertaken to program it for the case of buckling. The authors of the algorithm have published a paper dedicated to its application for the case of buckling (Williams & Wittrick 1983) in 1983, however, we believe that this algorithm needs a thorough explanation in light of modern programming techniques.
Consider equation (3.7). In general, it is valid only if either 0 or (recall the previous section). The authors stated (Wittrick & Williams 1971) that both conditions (4.1.2) should be considered to ensure that all eigenvalues have been computed. Originally the algorithm was defined to compute the vibration frequencies of the frame structures; thus all the notation and terminology is oriented toward the nomenclature widely accepted in dynamics. However, the algorithm is focused on computing the eigenvalues of various DS matrices that came from various problems. Thus, in general, it can be applied to any eigenvalue problem that is described by equation (3.7).
The main idea of the algorithm is based on Rayleigh's theorem on the frequencies of constrained structures and on the Sturm's theorem on the number of real distinct roots of a polynomial. Firstly, we quote here the Rayleigh's theorem reformulated by the authors (Wittrick & Williams 1971).

Theorem 1. If one constraint is removed from a structure, the number of natural frequencies which lie below some fixed chosen frequency either remains unchanged or increases by one.
We underline here that the same statement is true for the case of buckling in the same way as for vibration. The critical buckling loads are obtained from the eigenvalues of the DS matrix for buckling in a similar way that frequencies are computed from the eigenvalues of the DS matrix for vibration.
Secondly, we define necessary background and definitions before stating Sturm's theorem according to (Miller Thomas 1941). Sturm sequence or Sturm chain of polynomials is the sequence , , , … , , 3. The product is negative close before the root of and positive close after.
The third property along with the first two provides an essential feature of the Sturm sequence. Just before the root of the product is negative, then at the root, the sign is lost since 0 and shortly after the root the product has the same sign as . The number of such lost signs in any specific semi-closed interval gives the exact number of roots of the polynomial . Finally, the Sturm's theorem can be stated. 3) The idea of the W-W algorithm is to define an arbitrary interval , and to count the number of all eigenvalues of the DS matrix within this range to multiplicity using two counters based on the Theorem 1 and Theorem 2 respectively. However, this procedure does not give the values for the eigenvalue, only the number of them in the interval. This issue is addressed by performing the iterative reduction of the initial range to a set of intervals each of which contains only one eigenvalue. Then, the eigenvalues can be easily computed up to any required precision within each interval using numerical root finding techniques such as Newton's method.
Consider a transcendental DS matrix and select any arbitrary trial magnitude of the eigenvalue . We will count the total number of eigenvalues that lie in the range 0, . It will be performed by counting how many times the determinant of DS matrix changes its sign in this interval. Such procedure is called sign count of a matrix (Wittrick & Williams 1971).
In order to perform the sign count the matrix should be reduced to its upper triangular form using Gauss elimination process without pivoting. It is an essential procedure as it transforms the matrix to a form that satisfies requirements of the Theorem 2. Wittrick and Williams have proven (Wittrick & Williams 1971) that if the sequence of minors of DS matrix is constructed in a form: This sequence possesses the property of the Sturm sequence from Theorem 2. It is essential to ensure that the reduction to the upper triangular form is performed without pivoting because otherwise the signs of the sequence will be altered leading to the erroneous results. Thus, the standard LU decomposition that is widely used instead of Gaussian elimination due to its efficiency is not recommended unless it is insured that the decomposition does not alter the signs of the sequence 4.2.5.
After the upper triangular form of the DS matrix has been obtained, the number of roots of its determinant can be computed in a straightforward way according to the Theorem 2. Substitution of into yields the numerical form of the matrix . Count the number of negative elements on the main diagonal to obtain the initial estimate of the number of the eigenvalues of the DS matrix in the chosen interval.
The sign count of a matrix does not give the number of all eigenvalues in the interval if the interval is sufficiently large. This can be easily demonstrated using member DS matrix presented in equation (2.15). It is a 4 4 matrix which becomes a 2 2 matrix if fixed-free boundary conditions are applied. Such matrix reads as: .
(4.2.6) Thus, the maximum number of negative elements on the main diagonal of such matrix is 2. However, the matrix contains an infinite number of eigenvalues. It is apparent that the sign count can only show the exact number of eigenvalues and only in a small interval that contains no more than two eigenvalues. The sign count of a matrix is capable of counting no more than eigenvalues as no more than negative elements exist on the main diagonal of its upper triangular form. We propose calling the interval ∈ , in which each element on the main diagonal of becomes negative only once a period of the DS matrix. Each period contains no more than eigenvalues and the periods are repeated up to infinity. To illustrate this we plotted the sign count of a The map of and vectors is plotted in Fig. 4.2.1 in red dots illustrating the process of converging on eigenvalues of DS matrix from above. Blue lines are drawn solely for the convenience of representation. It is apparent from the figure that the sign count "jumps" from 1 to 2 in a periodic manner. The only difference to this periodicity is exhibited in the first period 0.5 , 2 . In this period the sign count changes from 0 to 1 for 0.5 stating that a single eigenvalue exists in the interval 0.5 , 0.5 . Then, the sign count increases from 1 to 2 stating that another eigenvalue exists in the interval 1.5 , 1.5 . It is known from the textbook (Timoshenko & Gere 2012) that the exact magnitudes of the lowest two eigenvalues for the fixed-free beam are 0.5 and 1.5 . These points are shown in the figure with purple dots. Thus, we see that the sign count of the matrix easily defines small intervals that enclose the eigenvalues of the DS matrix in the first period. However, the sign count does not produce any reliable information in the second period and so forth.
In cases when higher eigenvalues of the buckling problem are required the second counter of the W-W algorithm should be applied along with the sign count. If the sign count gives the precise number of the eigenvalues within each period of the DS matrix the second counter should obviously be able to account for the number of periods and for the number of the eigenvalues located within each period. The formulation of such counter is based on the Theorem 1. According to the original paper by Wittrick and Williams (Wittrick & Williams 1971), the counter is equal to the number of eigenvalues that would lie below if the structure has been completely constrained (all possible degrees of freedom suppressed). In our case, equals to the number of eigenvalues of a fixed-fixed beam that are smaller in magnitude than . All the eigenvalues of a fixed-fixed beam are known (Timoshenko & Gere 2012) thus can be easily computed. In a case when there are many members assembled in one structure the counter should be simply multiplied by the number of members. If the members are of different kinds the counter for each kind of members should be computed and then all the counters summed up. The simplicity with which the multiple substructures are encountered defines the power of the W-W algorithm for the case of complex hierarchical structures with multiple substructures.
The W-W algorithm can now be stated (Wittrick & Williams 1971). When converging from above, i.e. performing iterative reduction of the trial eigenvalue the total number of the eigenvalues that lie below the chosen one is equal: where is the number of eigenvalues smaller than if the structure were fully constrained and is the sign count of the DS matrix .
The application of the algorithm is illustrated for the case of a single fixed-free member defined by equation (4.2.6) in Fig. 4.2.2. It is apparent from the second graph in Fig. 4.2.2 that the W-W algorithm has located the small intervals , that contains eigenvalues of the DS matrix. By changing initial interval and running the algorithm one obtains any desired number of eigenvalue estimations. If the precision of the magnitudes of eigenvalues required is higher than the chosen step, then either the step should be reduced or a root finding technique (Newton's method) applied within each small interval that contains eigenvalues.
Summarizing everything stated above the following set of steps to calculate eigenvalues of the assembled dynamic stiffness matrix is proposed:  Chose initial eigenvalue interval 0, .  Substitute the value into functions in equation (2.16) that define entries of matrix. This operation will transform the DS matrix into a numerical matrix .  Reduce matrix to the upper triangular form using the Gauss technique without pivoting.  Count the number of the negative elements on the main diagonal of . This number is the initial estimation of the number of eigenvalues in the chosen interval.  Compute correction to the initial count by evaluating the number of the fixed-fixed solutions of each member with magnitude smaller than . If the number of substructures is , then the correction for the whole structure is defined as .
 The sum of the initial estimation and the correction gives the exact number of eigenvalues with magnitudes smaller than . Here sub-subscript 1 denotes the first iteration.  Start the second iteration by reducing by a small increment . Thus, the new interval is written as 0, 0, .
 Repeat the procedure for the new eigenvalue interval. The new number of eigenvalues that lie within the new interval is computed. Then, the two numbers and are compared. If these two numbers coincide, then no eigenvalue of the dynamic stiffness matrix lie in the small interval . If, however, these two numbers are different, i.e.
, then eigenvalues lie in the small interval .  Start the third iteration from reducing the initial interval by the increment 2 and repeat the procedure described in the previous steps. Following iterations are performed in the same way until either 0, or 0. These conditions ensure that the whole initial interval was analysed.  In such a way, the whole initial interval is divided into small intervals . The number of eigenvalues within each interval is known. The intervals for which contain eigenvalues while those for which 0 do not contain any eigenvalues. Iterative root finding operation can be performed in each small interval that contains eigenvalues to locate eigenvalues up to any required precision.

Conclusions
The DS matrix method has been shown to be a powerful approach to study elastic stability of structures. The method is exact and provides analytical dependences between the eigenvalues or critical buckling loads of structures and the physical parameters of these structures. The approach was discussed in details at all stages. Firstly, we have shown how member DS matrices are derived from member governing differential equations. Secondly, we have explained the process of assembly of the complete DS matrix of a structure from the member DS matrices along with stiffness matrices that contain information on coupling between members or on elastic supports. Finally, we discussed the approaches to compute the arbitrary number of eigenvalues of DS matrices with any required precision.
We believe that the DS matrix approach can be a good substitution for the FE analysis in many cases. In buckling it can be successfully applied to exactly compute the critical buckling loads of complex hierarchical structures that consist of large number of members much faster than by FE method. The DS matrix method is also capable of giving the exact analytical dependences of the critical buckling loads on any number of physical parameters which is essential for parametric studies. There is no known engineering software that operates with the DS matrix method for buckling. However, we have shown that this method can be automated and made very universal. Member DS matrices are known and the assembly can be easily programmed following the guidelines given here. The computation of the eigenvalues of the complete global DS matrix according to the W-W algorithm can also be readily automated by the appropriate computer routine developed following the steps proposed in section 4.2.
The DS method can also be extended beyond the problems of vibrations or buckling. As long as the member DS matrix is known it can be used along with appropriately modified W-W algorithm to address various physical problems. This fact makes the DS matrix method as universal as the FE analysis.