A Direct Eigenanalysis of Multibody System in Equilibrium

Abstract

This paper presents a direct eigenanalysis procedure for multibody system in equilibrium. The first kind Lagrange’s equation of the dynamics of multibody system is a set of differential algebraic equations, and the equations can be used to solve the equilibrium of the system. The vibration of the system about the equilibrium can be described by the linearization of the governing equation with the generalized coordinates and the multipliers as the perturbed variables. But the multiplier variables and the generalize coordinates are not in the same dimension. As a result, the system matrices in the perturbed vibration equations are badly conditioned, and a direct application of the mature eigensolvers does not guarantee a correct solution to the corresponding eigenvalue problem. This paper discusses the condition number of the problem and proposes a method for preconditioning the system matrices, then the corresponding eigenvalue problem of the multibody system about equilibrium can be smoothly solved with standard eigensolver such as ARPACK. In addition, a necessary frequency shift technology is also presented in the paper. The importance of matrix conditioning and the effectiveness of the presented method for preconditioning are demonstrated with numerical examples.

1. Introduction

Modal analysis of multibody system [1, 2] has important application in structure analysis and modal reduction method [3–6]. The eigenanalysis methods for constrained multibody system can be mainly divided into two categories, namely, the eigenanalysis method by transforming the equations of motion from DAEs to ODEs by selecting independent coordinates [7–9] and the direct eigenanalysis method [10–12]. Compared with structural systems, the governing equation of constrained multibody system is differential algebraic equations, instead of the ordinary differential systems. In the reduction methods, the differential algebraic equations of multibody system are reduced into ordinary differential equations for independent general coordinates, through looking for a set of orthogonal basis of constraint Jacobian matrix to eliminate the constraints from the system governing equations. The reduction method is accurate, efficient, and stable provided that the degree of freedom of the system and the number of constraints are small. However, the reduction method is considered not suitable for large scale problems due to the spoilage of sparse matrices [13].

In the direct methods, the nonlinear differential and algebraic equations are simultaneously linearized about a given state of the system, that is, the algebraic constraint equations are kept during the eigenanalysis and the sparsity of system matrices is preserved. Depending on how the problem is actually formulated, the enlarged mass matrix may become structurally singular (eigenvalues associated to constraints become infinite); other approaches may result in a structurally singular enlarged stiffness matrix (eigenvalues associated to constraints become 0). Usually a frequency shift is applied, and structurally nonsingular enlarged matrices are resulted. But in most cases, matrices’ condition number of the frequency shifted system is very large and numerical singular due to the constraints and a corresponding precondition is needed. The root cause for numerical singularity can be understood by the fact that the dynamic equations and the Lagrange multipliers carry a dimension of force while the constraint equations and the general coordinate a dimension of length or angular. The literature [14, 15] suggested multiplying a parameter to the constraints’ corresponding rows and columns of system matrices as a preconditioning. The time step integration takes a parameter associate with time step and integration format to do the precondition. But the eigenanalysis about a static equilibrium has nothing to do with time integration, and we need to find an appropriate parameter.

For the direct eigenanalysis of constrained multibody system in static equilibrium, this paper proposes a procedure and a necessary preconditioning method of system matrices to guarantee correct eigensolution. This paper is organized as follows. Section 2 schematically describes the formulation of constrained multibody system based on Lagrange equation with multipliers and the direct eigenvalue problem about equilibrium. Section 3 presented a detailed condition number analysis of the system matrices and a preconditioning scaling technique. Section 4 introduces the selection of eigensolver and numerical examples.

2. Governing Equations of Multibody System and the Direct Eigenvalue Problem

The dynamic equations of a general multibody system are usually established with Lagrange’s equations with multipliers [16]. The system equations are 𝑑𝑑𝑡𝜕𝜕𝐿𝐱̇−𝜕𝜕𝐿𝐱=𝜕𝜕𝚽𝐱^𝑇𝚽(𝐱)=𝟎,𝝀+𝐟,(2.1) where 𝐱∈ℝ^𝑛 is the vector of general coordinates, 𝝀∈ℝ^𝑚 is the vector of Lagrange multipliers, 𝐿 is the Lagrange function associated with (𝐱̇,𝐱), 𝐟 is the vector of general nonconservative forces which are functions of (𝐱̇,𝐱,𝑡), Φ is the vector of constraint functions associated with (𝐱), and (𝜕Φ/𝜕𝐱)^𝑇𝝀 is the vector of generalized constraint forces. Denote𝐅(𝐱̈,𝐱̇,𝐱,𝑡)=−𝑑𝐆(𝐱𝑑)𝑡=𝜕𝜕𝐿𝐱𝜕̇𝚽+(𝐱𝜕𝜕)𝐿𝐱+𝐟(𝐱̈,𝐱̇,𝐱,𝑡),𝜕𝐱^𝑇.(2.2) Then (2.1) becomes 𝐅(𝐱̈,𝐱̇,𝐱𝚽,𝑡()𝐱+)=𝐆𝟎(𝐱.)𝝀=𝟎,(2.3) The equilibrium state (𝐱_𝟎,𝝀_𝟎,𝑡₀) satisfies𝐅𝟎,𝟎,𝐱_𝟎𝚽,𝑡𝐱+𝐆𝐱_𝟎𝝀_𝟎=𝟎,_𝟎=𝟎.(2.4) Assume small vibration about the equilibrium (𝛿𝐱̈,𝛿𝐱̇,𝐱₀+𝛿𝐱,𝝀₀+𝛿𝝀,𝑡), then 𝐅𝛿𝐱̈,𝛿𝐱̇,𝐱_𝟎+𝛿𝐱𝚽,𝑡𝐱+𝐆𝐱₀+𝛿𝐱𝝀₀+𝛿𝝀=𝟎,_𝟎+𝛿𝐱=𝟎.(2.5) Carry out the Taylor expansion and reserve linear terms, we have 𝐌𝐱₀𝛿𝐲̈+𝐂𝐱₀𝛿𝐲̇+𝐊𝐱₀,𝝀₀𝛿𝐲=𝟎,(2.6) where 𝐌(𝐱₀)=^{𝐌(𝐱0)𝟎𝟎𝟎}, 𝐂(𝐱₀)^{𝐂(𝐱0)𝟎𝟎𝟎}, 𝐊(𝐱₀,𝝀₀)=^{𝐊(𝐱0,𝝀0)𝐆(𝐱0)𝐆}_{𝑇(𝐱0)𝟎}, 𝐌(𝐱₀)=(𝜕𝐅/𝜕𝐱̈)(𝐱₀), 𝐂(𝐱₀)=(𝜕𝐅/𝜕𝐱̇)(𝐱₀), 𝐊(𝐱₀,𝝀₀)=(𝜕(𝐅+𝐆𝝀)/𝜕𝐱)(𝐱₀,𝝀₀), 𝛿𝐲=^{𝛿𝐱𝛿𝝀}.

Write the solution as 𝛿𝐲(𝑡)=𝑒^𝑟𝑡𝐕, where 𝑟 is one of the eigenvalues, 𝐕 is the corresponding eigenvector, and the eigenvalue problem is 𝑟²𝐌+𝑟𝐂+𝐊𝐕=𝟎.(2.7) A frequency shift is applied to (2.7) by let 𝑟=𝑟−𝛼, 𝛼∈ℝ⁺, and then 𝑟²𝐌+𝑟𝐂+𝐊𝐕=𝟎,(2.8) where 𝐌=𝐌, 𝐊=𝛼²𝐌+𝛼𝐂+𝐊, 𝐂=2𝛼𝐌+𝐂. The frequency shift makes us calculate the eigenvalues about 𝛼 and makes sure 𝐊’s reversibility (the probability of singular is zero when 𝛼≠0) as the system’s energy must be in the form of kinetic energy or deformation energy. The eigenvalue problem can be transformed into the state space form as the following: 𝑟𝐌𝐘=𝐊𝐘,(2.9) where 𝐘=^𝐲_𝑟𝐲, 𝐌=_{−𝐂−𝐌𝐈𝟎}, 𝐊=_{𝐊𝟎𝟎𝐈}.

3. Ill-Condition in the Matrices of the Perturbed Vibration Equations and the Preconditioning

Consider the shifted eigenvalue problem (2.8) in the undamped case 𝑟²𝐌+𝐊𝐕=𝟎.(3.1) It is more convenient to consider the standard eigenvalue problem than the original generalized one 𝐊⁻¹𝐌+𝑟^−𝟐𝐈𝐕=𝟎.(3.2) As the inverse of matrix 𝐊 is used, cases in large condition number of 𝐊 could make the direct solution of 𝐊⁻¹𝐌 result in erroneous eigensolution. As shown in Section 4, the eigenvalue problems (2.8) obtained by the linearization of the Lagrange equations with multipliers for practical multibody systems usually encounter serious condition number problem. The reason mainly lies in the different dimensions of the multiplier variables and constraint equations with that of generalized coordinates and governing equations.

Inspired by the eigenanalysis method through transforming the equations of motion from DAEs to ODEs by selecting independent coordinates, we divide the matrix into blocks and calculate its inverse by Gaussian elimination. Let 𝛼²𝐌+𝐊=𝑐𝐃 where 𝑐=‖𝛼²𝐌+𝐊‖_∞. Without loss of generality, let the first 𝑚 (the number of constraints) rows of 𝐆 be linearly independent and be denoted as 𝐆₁∈ℝ^𝑚×𝑚. Divide the matrix as 𝐊=𝑐𝐃¹¹𝑐𝐃₁₂𝐆₁𝑐𝐃₂₁𝑐𝐃₂₂𝐆₂𝐆^𝑇1𝐆^𝑇2𝟎,(3.3) where 𝐃₁₁∈ℝ^𝑚×𝑚, 𝐃₂₂∈ℝ^{(𝑛−𝑚)×(𝑛−𝑚)}, 𝐃₁₂∈ℝ^{𝑚×(𝑛−𝑚)}, 𝐃₂₁∈ℝ^{(𝑛−𝑚)×𝑚}, 𝐆₂∈ℝ^{(𝑛−𝑚)×𝑚}, 𝑛 is the number of freedom of general coordinates.

Then 𝐊⁻¹=𝑐⁻¹⎡⎢⎢⎣𝐆−^1−𝑇𝐆𝐏^𝑇2₋₁𝐏⁻¹𝐆₂𝐆¹⁻¹−𝐆^1−𝑇𝐆^𝑇2𝐏⁻¹𝑐𝐆^1−𝑇+𝐆^1−𝑇𝐆^𝑇2𝐏⁻¹𝐑𝑐𝐆₁₋₁𝐆₂𝐆¹⁻¹𝐏⁻¹−𝑐𝐏⁻¹𝐑+𝐐𝐏⁻¹𝐆₂𝐆¹⁻¹−𝑐𝐐𝐏⁻¹𝑐²𝐐𝐏⁻¹𝐑−𝐍⎤⎥⎥⎦,(3.4) where 𝐍=𝐆¹⁻¹𝐃₁₁𝐆^1−𝑇, 𝐏=𝐃₂₂+𝐆₂𝐍𝐆^𝑇2−𝐃₂₁𝐆^1−𝑇𝐆^𝑇2−𝐆₂𝐆¹⁻¹𝐃₁₂, 𝐐=𝐆¹⁻¹𝐃₁₂−𝐍𝐆^𝑇2, 𝐑=𝐃₂₁𝐆^1−𝑇−𝐆₂𝐍. See appendix for detail of Gaussian elimination.

Assume ‖𝐍‖_∞≈‖‖𝐏⁻¹‖‖_∞≈‖𝐐‖_∞≈‖𝐑‖_∞=𝑂(1).(3.5) Then it shows ‖𝐊⁻¹‖_∞=𝑂(𝑐) and by taking condition number cond(⋅)=‖⋅‖_∞‖⋅⁻¹‖_∞cond𝐊≈𝑂𝑐².(3.6)

Equation (3.5) can be expected when the elements of eigenvector of the underlying ODE system have about the same magnitude, that is, the elastic parameters of the physical system are about the same magnitude in physical description, and this is often the case. The elastic parameters mean the Young’ modulus of flexible bodies, stiffness of spring forces, and so forth.

The condition number of the original system can be changed by scaling the multiplier variables and constraint equations. Let 𝝀=𝑘⁻¹𝝀 and multiply the constraint equations by 𝑘 in (2.1) we have 𝑑𝑑𝑡𝜕𝜕𝐿𝐱̇−𝜕𝜕𝐿𝐱=𝑘𝜕𝜕𝚽𝐱^𝑇𝑘𝚽(𝐱)=𝟎.𝝀+𝐟,(3.7) The corresponding block matrices are𝐊=𝑐𝐃¹¹𝑐𝐃₁₂𝑘𝐆₁𝑐𝐃₂₁𝑐𝐃₂₂𝑘𝐆₂𝑘𝐆^𝑇1𝑘𝐆^𝑇2𝟎,𝐊⁻¹=𝑐⁻¹⎡⎢⎢⎣𝐆−^1−𝑇𝐆𝐏^𝑇2₋₁𝐏⁻¹𝐆₂𝐆¹⁻¹−𝐆^1−𝑇𝐆^𝑇2𝐏⁻¹𝑘⁻¹𝑐𝐆^1−𝑇+𝐆^1−𝑇𝐆^𝑇2𝐏⁻¹𝐑𝑘₋₁𝐆₂𝐆¹⁻¹𝐏⁻¹−𝑘⁻¹𝑐𝐏⁻¹𝐑𝑐𝐆¹⁻¹+𝐐𝐏⁻¹𝐆₂𝐆¹⁻¹−𝑘⁻¹𝑐𝐐𝐏⁻¹𝑘⁻²𝑐²𝐐𝐏⁻¹𝐑−𝐍⎤⎥⎥⎦.(3.8) Then the condition number analysis shows cond𝐊≈𝑂max(𝑘,𝑐)×max𝑘⁻¹𝑐,𝑘⁻²𝑐².(3.9) Choosing 𝑘≈𝑐, the estimation of condition number decreases from about 𝑂(𝑐²) to about 𝑂(1).

The elements of matrix 𝐊 stand for the elastic coefficients of the system and the elements of 𝐆 and 𝐆^𝑇work similarly as the stiffness of constraints. Taking 𝐊 as an enlarged stiffness matrix, the scaling is equivalent to changing the stiffness of constraints from 𝑂(1) to the same magnitude of physical stiffness.

Note that the condition cannot be improved by only scaling multipliers or by only scaling constraint equations. It is suggested that any different dimension in the procedure of modeling could lead to bad condition of numerical calculation.

The complete process for direct eigenanalysis of constrained deformable multibody system is thus as follows.

Step 1. Calculate the static equilibrium (𝐱₀,𝝀₀) from (2.4).

Step 2. Calculate 𝐌, 𝐂, 𝐊 at (𝐱₀,𝝀₀).

Step 3. Precondition 𝐊 with 𝑐 which is the maximum absolute element of 𝐊.

Step 4. Frequency shift with parameter 𝛼(1∼10) with (2.8).

Step 5. Enlarge the matrix into state space as (2.9).

Step 6. Solve the eigenvalues and eigenvectors form (2.9) with eigensolver such as ARPACK.

Step 7. Add 𝛼to the eigenvalues and multiply the eigenvector components of Lagrange multipliers with 𝑐.

4. Selection of Eigensolver and Numerical Examples

The famous iteration method Implicitly Restarted Arnoldi Method (IRAM) [17, 18] is used in the numerical solution of eigenvalue problem (2.9) in the paper. Arnoldi method is one of the Krylov subspace methods [19] which calculate the largest few eigenvalues of a matrix. The presented numerical examples have verified the efficiency of IRAM in practical calculation. The famous parallel sparse linear system solver PARDISO [20, 21] is used to provide 𝐊⁻¹ of (2.9) in the paper. A cantilever and a four-bar linkage are calculated to illustrate the effect of the condition number of the system matrix and the effectiveness of the presented preconditioning method. It is shown that the numerical results are physically meaningless without scaling and the precondition presented in the paper is very effective.

4.1. A Cantilever

In this example, beam element and 3-dimensional solid element are separately used to model the cantilever in the multibody frame, as shown in Figure 1. The parameters of the cantilever are shown in Table 1. The models are computed with the presented direct eigenanalysis method in two cases, that is, with and without axial tension, and the results are compared with analytical results based on the mathematical physics, see Tables 2 and 3.

Table 1 

Parameters for the cantilever.

Table 2 

Eigenvalues comparison for cantilever (𝐹=0 Newton).

Table 3 

Eigenvalues comparison for cantilever (𝐹=1𝑒10 Newton).

Figure 1 

Two models of cantilever.

The governing equation and boundary conditions for the single-mode-free vibration of the cantilever under axial tension is as the following: 𝐸𝐼𝑑⁴𝑦𝑑𝑥⁴−𝐹𝑑²𝑦𝑑𝑥²𝑦||−𝜌𝜔²𝑦=0,₀=0𝑑𝑦𝑑𝑥|||₀=0𝑑²𝑦𝑑𝑥²||||_𝑙=0𝑑³𝑦𝑑𝑥³||||_𝑙=𝑎𝑑𝑦𝑑𝑥|||_𝑙,(4.1) where 𝑥 is the material coordinate, 𝑦 is the lateral vibration displacement mode, 𝐸 is the Young’s modulus, 𝐼 is the sectional moment of inertia, 𝐹 is the axial tension, 𝜌 is the linear density, and 𝜔 is the circular frequency of vibration. 𝑙 is the length of the cantilever and 𝑎=𝐹/𝐸𝐼. With the conventional mathematical procedure, characteristic equation for frequency 𝜔 is 𝑟²¹cosh𝑟₁𝑙+𝑟²²cos𝑟₂𝑙𝑟³¹s=inh𝑟₁𝑙−𝑟³²sin𝑟₂𝑙𝑟+𝑎−𝑟₁sinh𝑟₁𝑙−𝑟₂sin𝑟₂𝑙₂₁sinh𝑟₁𝑙+𝑟₁𝑟₂sin𝑟₂𝑙𝑟³¹cosh𝑟₁𝑙+𝑟₁𝑟²²cos𝑟₂𝑙+𝑎−𝑟₁cosh𝑟₁𝑙−𝑟₁cos𝑟₂𝑙,(4.2) where 𝑏=𝜌𝜔²/𝐸𝐼, 𝑟₁=(√𝑎²+4𝑏+𝑎)/2, and 𝑟₂=(√𝑎²+4𝑏−𝑎)/2.

It is shown in Tables 2 and 3 that the results without the scaling precondition are totally erroneous in both cases. In the two models, through computation 𝑐=1.39167𝑒+06 for 𝐹=0 and 𝑐=1.39252𝑒+06 for 𝐹=1𝑒+10. Precondition parameter 𝑘 is chosen the same as 𝑐 and the condition numbers in the two computed cases are decreased about 𝑂(𝑐²) times after the precondition. The results show that the presented precondition procedure could guarantee correct results for the direct eigenanalysis, and the unscaled models give out meaningless results by ARPACK as the condition numbers are too large for the nonphysical contribution of constraints. Note that model 2 is used to illustrate the precondition’s effectiveness for models including other flexible elements such as 3D solid elements, and the error of scaled model 2 to Euler beam analytical results is partly because the solid model does not agree with the plan cross section assumption of Euler-Bernoulli beam.

4.2. Four-Bar Linkage

The four-bar linkage is shown in Figure 2, and the parameters of the bar are shown in Table 4. The model is modeled with Cartesian coordinates, and the constraints in the model are all revolute joints. The beam element model is used in this example, as shown in Figure 3. The equilibrium under gravity is first calculated with dynamic relaxation method, [22] and the eigenvalues are calculated about the deformed configuration. The equilibrium configuration of the four-bar system is also shown in Figure 3; the computed results with and without the preconditioning scaling are shown in Table 5 compared with results of FEA software ABAQUS with the same model. Again the unscaled results in the case are all meaningless as they have nonphysical real parts. In this problem 𝑐=9.11047𝑒+13 and precondition parameter 𝑘 is chosen, the same value as 𝑐. The condition numbers are reduced about 𝑂(𝑐) times after the preconditioning, and correct results are obtained with the direct eigenanalysis method.

Table 4 

Parameters for the linkage.

Table 5 

Eigenvalues with or without scaling.

Figure 2 

Four-bar linkage discrete model.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Modes of the four-bar linkage.

5. Summary and Conclusions

In this paper, the bad condition of the eigenvalue problem formulated from the direct linearization of the differential-algebraic equations of multibody system is reviewed and then discussed. The root cause may be attributed to the different dimensions of constraint multipliers and equations with physical system. Through analysis, a simple scale of the constraint multiplier variables and equations according to the stiffness of the elastic elements in the system can reduce the condition number of the system matrix to be inverted. After the preconditioning scale, eigensolvers such as ARPACK can be directly applied to the preconditioned system and yield correct result. The soundness of proposed approach and applicability are illustrated by comparing the results of a cantilever and a four-bar linkage with those obtained from either analysis or ABAQUS.

Appendix

A. Gaussian elimination to get (3.4) from (3.3)

𝑐𝐃¹¹𝑐𝐃₁₂𝐆₁𝑐𝐃₂₁𝑐𝐃₂₂𝐆₂𝐆^𝑇1𝐆^𝑇2𝟎∼𝐈^𝑚𝟎𝐈𝟎_𝑛−𝑚𝟎𝟎𝟎𝐈𝟎_𝑚._{1.Unitize1stand3rdlinebyleftmultiplying𝐆-11to1stlineand𝐆1−𝑇3rdline.2.Exchange1stand3rdline.}−−−−−𝐈−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−→_𝑚𝑐𝐃₂₁𝐆^1−𝑇𝐆^𝑇2𝟎𝑐𝐃₂₂𝐆₂𝑐𝐆¹⁻¹𝐃₁₁𝑐𝐆¹⁻¹𝐃₁₂𝐈_𝑚∼𝟎𝟎𝐆^1−𝑇𝐆𝟎₁₋₁𝐈_𝑛−𝑚𝟎𝟎𝟎_{1.Leftmultiply−𝑐𝐃21to1stlineandaddto2ndline.2.Leftmultiply−𝑐𝐆1−1𝐃11to1stlineandaddto3rdline.3.Leftmultiply−𝐆2to3rdlineandaddto2ndline.}.⎡⎢⎢⎣−−𝐈−−−−−−−−−−−−−−−−−−−−−−−−−−−−−→_𝑚𝟎𝑐𝐃₂₂𝐆^1−𝑇𝐆^𝑇2𝟎𝟎+𝐆₂𝐆𝑐¹⁻¹𝐆𝐃₁₁₁₋₁𝐆^1−𝑇𝐆^𝑇2−𝐃₂₁𝐆^1−𝑇𝐆^𝑇2−𝐆₂𝐆¹⁻¹𝐃₁₂𝟎𝐃₁₂−𝐆¹⁻¹𝐃₁₁𝐆^1−𝑇𝐆^𝑇2𝐈_𝑚∼⎡⎢⎢⎣𝟎𝟎𝐆_1−𝑇⎤⎥⎥⎦−𝐆𝐆₂𝐆¹⁻¹₁₋₁𝐈_𝑛−𝑚−𝑐𝐃₂₁𝐆^1−𝑇−𝐆₂𝐆¹⁻¹𝐃₁₁𝐆^1−𝑇𝟎−𝑐𝐆¹⁻¹𝐃₁₁𝐆^1−𝑇⎤⎥⎥⎦._{Let𝐍=𝐆1−1𝐃11𝐆1−𝑇,𝐏=𝐃22+𝐆2𝐍𝐆𝑇2−𝐃21𝐆1−𝑇𝐆𝑇2−𝐆2𝐆1−1𝐃12,𝐐=𝐆1−1𝐃12−𝐍𝐆𝑇2,𝐑=𝐃21𝐆1−𝑇−𝐆2𝐍.}−−𝐈−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−→_𝑚𝟎𝟎𝐆𝑐𝑐^1−𝑇𝐐𝐏𝐆^𝑇2𝐈𝟎𝟎_𝑚∼⎡⎢⎢⎣𝟎𝟎𝐆^1−𝑇𝐆𝐆₂𝐆¹⁻¹𝐈_𝑛−𝑚−𝑐𝐑₁₋₁𝟎−𝑐𝐍_{1.Unitize2ndlinebyleftmultiplying𝑐−1𝐏−1reversibilityofthewholematrixensuretheexistenceof𝐏−1.2.Leftmultiply−𝐆1−𝑇𝐆𝑇2to2ndlineandaddto1stline.3.Leftmultiply−𝑐𝐐to2ndlineandaddto3rdline.}⎤⎥⎥⎦.−−𝐈−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−→_𝑚𝟎𝐈𝟎_𝑛−𝑚𝟎𝟎𝟎𝐈𝟎_𝑚∼𝑐⁻¹⎡⎢⎢⎣𝐆−^1−𝑇𝐆𝐏^𝑇2₋₁𝐏⁻¹𝐆₂𝐆¹⁻¹−𝐆^1−𝑇𝐆^𝑇2𝐏⁻¹𝑐𝐆^1−𝑇+𝐆^1−𝑇𝐆^𝑇2𝐏⁻¹𝐑𝑐𝐆₁₋₁𝐆₂𝐆¹⁻¹𝐏⁻¹−𝑐𝐏⁻¹𝐑+𝐐𝐏⁻¹𝐆₂𝐆¹⁻¹−𝑐𝐐𝐏⁻¹𝑐²𝐐𝐏⁻¹𝐑−𝐍⎤⎥⎥⎦.(A.1)

Acknowledgment

Cheng Yang and Zhengru Zhang were supported by National NSF of China under Grant 11071124.