- Email: [email protected]
Implicit integration with coordinate partitioning
Applied Mathematics and Computation 111 (2000) 7±31 www.elsevier.nl/locate/amc
Implicit integration with coordinate partitioning Joseph F. McGrath
*,1,
Rajiv Rampalli
Mechanical Dynamics, Inc., Ann Arbor, MI, USA
Abstract The dynamics of a mechanical system are governed by dierential and algebraic (constraint) equations. If these DAE's are reduced by coordinate partitioning to a system of ®rst order dierential equations for the independent unknowns, then numerical methods for equations in explicit form can be used to compute the solution. We have developed techniques for the evaluation of the sparse Jacobian matrices for solving the partitioned equations with integrators based on backward dierentiation formulas, an important method for the simulation of the inherently sti mechanical systems. Ó 1999 Elsevier Science Inc. All rights reserved. AMS classi®cation: 34A50; 34A65; 65C20; 65L05 Keywords: Coordinate partitioning; Dynamics; Dierential and algebraic equations; Mechanical system simulation
1. Introduction and summary Coordinate partitioning is a numerical procedure for transforming the system of dierential and algebraic equations that governs mechanical motion into an initial value problem consisting of ordinary dierential equations. The reformulation de®ned by the method explicitly evaluates the time-derivatives of the independent coordinates in the problem. Thus, the unknowns can be integrated to obtain the solution by computer programs for systems of explicit dierential equations. The integrators may employ either explicit or implicit * Corresponding author. Silicon Graphics, Inc., 102 Central Park Square, Los Alamos, NM 87544, USA. Tel.: 505 662 4000; fax: 505 662 9457; e-mail: [email protected]. 1 Adjunct member of the faculty, Department of Mathematics, University of Michigan.
0096-3003/00/$ - see front matter Ó 2000 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 9 8 ) 1 0 0 2 3 - 1
8
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
methods. However, more information than just the time-derivatives of the unknowns must be generated for the techniques implemented in the implicit integrators. The algorithm in the implicit integrators based on backward dierentiation formulas needs the Jacobian matrix of the time-derivatives of the unknowns with respect to the unknowns themselves to advance the solution from one time step to the next. As described in this paper, it turns out that, in addition to evaluating the time-derivatives, the coordinate partitioning algorithm provides the machinery for the formation of the required Jacobian matrix. Speci®cally, the matrices that are used in the coordinate partitioning function to calculate the derivatives are the same as the coecient matrices needed to compute the various components of the Jacobian matrix for the implicit integrator. Background material and the context for the discussion of the coordinate partitioning algorithm are included in Sections 2 and 3. A review of the coordinate partitioning algorithm is contained in Section 4 where a basic example is developed to illustrate the derivation of the Jacobian matrix for an implicit integrator. The general concepts involved in implicit integration with coordinate partitioning are covered by the end of Section 4. For the reader interested in the full details, Section 5 contains the derivation of the Jacobian matrix for the general case. Future developments for acceleration-dependent forces and non holonomic constraints are mentioned in Section 6. 2. Explicit dierential equations A system of ordinary dierential equations is classi®ed as explicit if it has the form w_ f w; t 1 as opposed to a system in the form _ w; t 0; f w; where the derivatives are implicitly de®ned. In either case, the array 2 3 f1 6 f2 7 6 7 f 6 . 7 4 .. 5
2
3
fN
consists of N functions of the N unknowns in the array 2 3 w1 6 w2 7 6 7 w 6 .. 7: 4 . 5 wN
4
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
For the implicit equations (2), f is a function of both w and 2 3 w_ 1 6 w_ 2 7 6 7 w_ 6 .. 7: 4 . 5
9
5
w_ N
Finally, each component of f may also be a function of time t. The following schematic diagram of a function evaluation illustrates the explicit form (1) of a system of dierential equations. The essential
characteristic is the functional nature of the system. No information about the calculations themselves is contained in the output from an evaluation of the dierential equations. The coordinate partitioning algorithm satis®es the above criteria for a system of explicit dierential equations: The implicit system (2) is converted into an explicit system by concealing the implicit nature of the equations within the algorithmic box in the diagram and by computing the derivatives w_ internally as a function of the input values of w and t.
3. Methods for explicit initial value problems An important category of methods for solving initial value problems (IVP's) needs the system of equations to be in the explicit form (1). The various methods in this category of solvers include both explicit and implicit integrators. An explicit integrator uses current and previous information about the solution to compute its values at the next time step in a ®xed number of function evaluations. In contrast, an implicit integrator employs an iterative
10
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
algorithm (implying an unknown, variable number of function evaluations to achieve convergence) to solve the equations in either form (1) or (2) at the future time. For example, the DEPAC software package [1,2] contains the following three solvers: DDERKF DDEABM DDEBDF
The single step, explicit, fourth order Runge±Kutta±Fehlberg algorithm with ®fth order corrections; A multi-step, variable order, explicit Adams±Bashforth±Moulton integrator (predict, evaluate, correct, evaluate); A multi-step, variable order, implicit method based on backward dierentiation formulas.
Each one of the integrators on the above list expects the problem to be de®ned in the explicit form (1). The DDERKF and DDEABM integrators employ explicit algorithms while an implicit method is implemented in the DDEBDF integrator. Other software products, such as the DASSL code, use implicit methods to solve IVP's in the implicit form (2) [3,4]. To compute the solution, an explicit integrator requires only the evaluation of the derivatives in the form of Eq. (1) as illustrated by the schematic diagram on page 3. The algorithm does not need any details about the calculations themselves. However, an implicit integrator makes use of additional information about the behavior of the function f . In particular, the implicit integrators DDEBDF and DASSL need the Jacobian matrix 2 of1 of1 of1 3 . . . ow ow1 ow2 N 7 6 of of2 of2 7 6 2 6 ow1 ow2 . . . owN 7 7 6 6 7 6 . .. 7 6 .. . 5 4 ofN ofN ofN . . . ow ow1 ow2 N for the vector function f in Eq. (1). The matrix (6) provides the implicit algorithm with quanti®ed information about the sensitivity of each component of the output to each component of the input for an explicit system of equations. 4. Coordinate partitioning The system of n second order dierential equations and m algebraic equations _ q; t M q ÿ UTq k F q; U q; t 0
7
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
11
governs the dynamics of a mechanical system. Eq. (7) do not form an explicit system of ordinary dierential equations in their present form. However, as explained in the following paragraphs, the coordinate partitioning algorithm performs a numerical transformation of the system (7) into an explicit initial value problem [5]. In the above system of dierential and algebraic equations (DAE's), the ncomponent array q consists of the translational and rotational coordinates that describe the current con®guration of the mechanism. The m Lagrange multipliers in the array k correspond to the m constraint equations represented by U. The coordinate partitioning software has three components: 1. The algorithm selects a p-element vector of independent coordinates from the vector q where p 6 n is equal to the number of degrees of freedom (DOF's) in the mechanical system governed by the dierential and algebraic equations (7); 2. It evaluates f the vector-function of the N 2p unknowns consisting of the p independent coordinates plus their ®rst derivatives; 3. It computes all of the other variables in the problem including the dependent variables, their derivatives, and the Lagrange multipliers. For item 2, the derivatives w_ f w; t are calculated by successively solving three systems of equations, one or two of which are non linear. In addition to the derivatives, implicit integrators usually require the Jacobian matrix Eq. (6) of the derivatives of the unknowns with respect to the unknowns themselves. For the coordinate partitioning algorithm, the entries in the Jacobian matrix are computed by solving two pairs of linear equations. The coecient matrices for the four systems of equations are the same as two of the matrices used for the evaluation of the function f . The expressions for the constant terms (the right-hand sides) involve additional derivations of a non trivial level of complexity. Another unfortunate complication is the reason for the rather technical nature of the subject matter in Section 5 of this report. There is no transformation of coordinates involved in the coordinate partitioning algorithm. The p independent coordinates in the system (7) are simply a selection of the entries in q. If m < n and the m equations that make up U are independent, then the number of independent variables is p n ÿ m and the dynamic motion of the system has p DOF's. If m n and the constraint equations are independent, then p 0 and the mechanism can only undergo kinematic motion. If m > n and more than n of the constraint equations are independent, the system is overconstrained. The latter two cases are outside of the scope of the current discussion. The coordinate partitioning algorithm is used only for the simulation of dynamic motion, that is, p > 0.
12
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
4.1. The point-mass pendulum A complete derivation of the coordinate partitioning equations for a basic mechanism with one moving part and the corresponding Jacobian matrix for a BDF integrator is contained in the following paragraphs. The eort spent developing the details for this non general case serves two purposes: It lays the groundwork needed by the authors to construct the general formulation to be developed in Section 5 and it provides an explanation of the method which is not obscured by the more abstract notation. Thus, by reading Section 4, it is possible for a mathematician or computational engineer to understand how the method works without mastering all of the details contained in Section 5. A mechanism made up of a single moving part attached by a revolute joint to a ®xed part (to ground) provides a good illustration of the coordinate partitioning algorithm. This two dimensional, single-DOF problem is governed by the system of dierential and algebraic equations (DAE's) _ x; y; t mx ÿ xk F1 _x; y; _ x; y; t m y ÿ yk F2 _x; y; 1 2 x y 2 ÿ 1 0 2
8
where · The point-mass m is at location x; y a ®xed, unit distance from the origin; · F1 and F2 are the horizontal and vertical components of the applied forces, respectively; · The unknown k is a Lagrange multiplier. For example, if F1 0 and F2 ÿmg where g is the acceleration due to gravity, then the DAE's (8) govern the motion of a simple pendulum with point-mass m. The constraint equation in the system (8) is scaled by 1=2 for the sake of conformity with the standard representation Eq. (7) for DAE's governing the dynamics of a mechanism. 4.2. The coordinate partitioning function Suppose the partitioning algorithm has selected x as the independent coordinate leaving y as the dependent coordinate 2 for the problem (8). Then, the input to the coordinate partitioning function is x_ w 9 x
2
For the current example, we will simply assume the existence of a satisfactory partitioning algorithm. A description of the process is included in Section 5.2.
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
13
and the expected output consists of the computed components of x w_ f w; t 10 x_ which is in the explicit form (1). The algorithm for computing x and x_ in terms of x_ and x is de®ned in the following schematic diagram where, for the sake of clarity, the ®nal calculations of the two output values are enclosed in boxes:
The D, V, and A/k equations in this basic example form a nonlinear system of seven equations in the seven unknowns 2 3 x 6y 7 6 7 6 x_ 7 6 7 6 y_ 7: 11 6 7 6 x 7 6 7 4 y 5 k
14
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
The equations for the models of mechanical systems can be uncoupled in a manner which improves the eciency of the calculations. That is, the solution can be computed sequentially in the coordinate partitioning code. The output from the solution to the ®rst block of equations D is part of the input data for the subsystem V. The cumulative results from solving the subsystems D and V along with the values of x_ x_ in and x xin make up the input for solving the subsystem A/k. In general, as illustrated by this example, if values are speci®ed for the independent displacements and their ®rst derivatives, the solutions to each of the subsystems, D, V, and A/k, yield the dependent displacements, the dependent velocities, and the accelerations and constraint forces (Lagrange multipliers), in that order. Hence, there are two important and very valuable advantages to the coordinate partitioning method: 1. While the integrator advances the independent displacements and velocities, the coordinate partitioning calculations simultaneously provide consistent values for all of the other entries in the state vector for the mechanism; and 2. The solution satis®es not only the dierential equations and the constraint equations but also the ®rst and second derivatives of the constraint equations. In summary, the coordinate partitioning method yields a complete state vector for the mechanism consisting of the displacements, the ®rst and second derivatives, and the constraint forces at each integration time step. In the current single-DOF example, this state vector is the column matrix (11). The Jacobian and coecient matrices used to solve the systems of equations in the function evaluation de®ned by the schematic diagram shown above are x y 12 1 0 for the dependent displacements D and the velocities V and 2 3 m 0 ÿx 6 7 40 m ÿy5 x y 0
13
for the accelerations and constraint forces A/k. Recall that the ®rst component of the algorithm, as described on page 5, calls for the partitioning of the coordinates. Although we have not said anything yet about the selection process, the roles of x and y could be reversed with y as the independent variable, then the Jacobian matrix (12) for the dependent displacements would be x y 14 0 1 instead of the matrix (12). The coecient matrix (13) would not change.
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
15
4.3. The Jacobian matrix for the BDF integrator Numerically sti initial value problems, that is, initial value problems with widely separated time constants in the various components of the solution, can be solved with implicit integration techniques such as the methods based on backward dierentiation formulas (BDF's). The implicit BDF integrator for explicit dierential equations in the form (1) requires the Jacobian matrix (6). Thus, for the coordinate partitioning algorithm, the additional Jacobian matrix 3 " # 2 of1 o_x of2 o_x
of1 ox of2 ox
oxout o_x
oxout oxin
out
o_xout oxin
4 o_x in o_xin
5
15
is needed to solve the problem (8). From the de®nition of the coordinate partitioning function in the schematic diagram on page 7, we can see that the output value of x_ is completely independent of the input value of x. Indeed, x_ is simply set equal to its own input value. Hence, ox_ out 1 ox_ in
16
ox_ out 0: oxin
17
and
which takes care of the second row of the 2 ´ 2 Jacobian matrix (15). Although the derivation of the entries in the ®rst row is not quite as easy as the second row of the matrix (15), the calculations can still be done quite ef®ciently as will be shown in the following discussion. In fact, the evaluation of the Jacobian matrix (15) can be achieved with just the coecient matrices (12) and (13) needed for the evaluation of the function itself. The signi®cance of this theorem (which will be proven for the general case in Section 5) lies, of course, with the fact that the (LU, QR, or perhaps another) decomposition of the matrices (12) and (13) need not be repeated for the evaluation of the requisite Jacobian matrix (15). In the following work, partial dierentiation of the subsystems V and A/k with respect to x_ will yield a system of equations for the ®rst entry in the ®rst row of the matrix (15). To get a system of equations for the second entry, we will need the partial derivatives of the subsystems V and A/k with respect to x.
16
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
4.3.1. First row, ®rst entry Since y is dependent only upon x and x is an independent variable, we have ox oy 0: ox_ ox_
18
Hence, partial dierentiation of the ®rst equation in the system V from the diagram is subsection 4.2 with respect to x_ yields xy
oy_ 0 ox_
_ x easily enough. which can be solved for oy=o_ If Eqs. (16) and (19) are collected into the linear system " #2 o_x 3 " # 0 x y 4 o_x 5 ; oy_ 1 1 0 o_x
19
20
it might appear that we have complicated matters. However, in order to illustrate the more general case involving multiple dependent and multiple independent coordinates, the dierentiation of the system V must be written in this more elaborate fashion, rather than as the individual equations that apply to the single-DOF example Eq. (8) currently under consideration. In particular, it is signi®cant that the coecient matrix for the system (20) is the Jacobian matrix (12) for the dependent displacements. Partial dierentiation of the system A/k with respect to x_ yields the system of linear equations 2 32 ox 3 2 oF1 oF1 oy_ 3 oy_ o_x m 0 ÿx o_x o_x 7 7 6 6 76 6 oF2 oF2 oy_ 7 o y7 6 0 m ÿ y 76 21 6 6 7 4 54 o_x 5 4 o_x oy_ o_x 7 5 ok x y 0 ÿ 2_x 2y_ oo_yx_ o_x where the constant terms on the right-hand side and the coecients have already been determined. One component of the solution to Eq. (21) is of1 ox : 22 ox_ ox_ We have now evaluated three quarters of the Jacobian matrix (15). Note that, the coecient matrix for the system (21) is the coecient matrix (13) for the accelerations and Lagrange multipliers. 4.3.2. First row, second entry Partial dierentiation of the ®rst equation in each of the systems D and V from the diagram on page 7 with respect to x yields
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
17
y
oy ÿx ox
23
y_
oy oy_ y ÿ_x; ox ox
24
and
respectively. Eqs. (23) and (24) can be solved directly ®rst for oy=ox and then _ for oy=ox. However, once again for the sake of the more general formulation, the values of these quantities have to be obtained by sequentially solving the system " #" ox # " # 0 x y ox 25 oy 1 1 0 ox and, using Eq. (17), the system # " #" o_x # " x y ÿ _x y_ oy ox ox : 1 0 ooxy_ 0
26
The computed value of oy=ox is needed to solve the Eq. (26). However, it turns out that the system (25) is super¯uous. A comparison of the two systems of linear equations (20) and (25) proves that oy_ oy 27 ox_ ox which, as will be shown in subsection 5.3.1, is true in general. Consequently, the work load can be lightened by eliminating the system (25) and using the results from Eqs. (20) and (27) to solve the system (26). Finally, partial dierentiation of the system A/k with respect to x yields 2 32 ox 3 2 k oF1 oF1 oy 3 m 0 ÿx ox oy ox ox 7 7 6 oy oF 6 76 o y7 2 6 7 6 28 4 0 m ÿ y 54 ox 5 4 k ox ox oFoy2 oy ox 5 ok x y 0 ÿ x y oy 2y_ oy_ ox ox
ox
for which the solution requires the results from Eqs. (26) and (27). One of the components of the solution to the system of equations (28) is of1 ox : 29 ox ox Once again note that the coecient matrices in the systems (25), (26), and (28) are the matrices (12) and (13) from the coordinate partitioning function for the time-derivatives. Substituting Eqs. (22), (16), (29) and (17) into the de®nition Eq. (15) yields the Jacobian matrix
18
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
" of
1
o_x of2 o_x
of1 ox of2 ox
"
#
oxout o_xin
1
oxout oxin
0
# 30
where the ®rst row has been computed by solving Eqs. (21) and (28). The evaluation of the Jacobian matrix (15) is now complete. An implicit integrator, such as the DDEBDF program, can use the coordinate partitioning function de®ned in the diagram in subsection 4.2 and the Jacobian matrix (30) to compute the solution to the problem (8). 5. The general case The constraint equations in the system (7), repeated here for convenience, _ q; t M q ÿ UTq k F q;
31 U q; t 0 reduce the number of degrees of freedom in the motion of the mechanism by imposing interdependencies between some of the variables. Assume that there are n dierential equations and m < n constraint equations (the algebraic equations) in the system (31). Also, assume that the m constraints are mutually independent in which case there are p n ÿ m mechanical degrees of freedom in the problem. Let B be the p n matrix of 0's and 1's that selects the p independent variables from the array q of the n translational and rotational coordinates in the problem. Thus, for the explicit initial value problem (10), also repeated here for convenience, w_ f w; t; 32 the array of unknowns is r_ Bq_ w 33 r Bq where r is the p-component column matrix of independent coordinates. The remaining m n ÿ p dependent coordinates will constitute the s array. For the moment, we take this partitioning process for granted. A few words about the formation of B are included toward the end of Section 5.2. In the present context, a couple of examples are sucient to explain the relationship Eq. (33). If y and w are the independent variables in the system of coordinates 2 3 x 6y7 6 7 6z7 7 34 q6 6 w 7; 6 7 4h5 /
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
then
1 0 0 0 0 1 y r Bq ; w B
and
0 0
0 0
0 ; 0
2 3 x 6z7 7 s6 4 h 5: /
19
35 36
37
For the second illustration, the single-DOF example discussed in Section 4 has B 1 0; r x, and s y. 5.1. The general coordinate partitioning function For the general case Eq. (31), the following schematic diagram de®nes the coordinate partitioning function where, for the sake of clarity, the ®nal expressions for the output values are enclosed in boxes:
20
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
Note that the same two signi®cant advantages listed in subsection 4.2 for the basic example pertain to the algorithm de®ned in this ®gure. That is, the coordinate partitioning method yields the values of the complete state vector for the mechanism including the displacements, their ®rst and second derivatives, and the constraint forces, all of which satisfy the DAE's (31) plus the ®rst and second derivatives of the constraint equations. 5.1.1. Notation In the notation used in the above schematic diagram and throughout Section 5 of this report, · U represents the m constraint equations in the form U q; t 0; · Uq is the m n (generally sparse) Jacobian matrix of U with respect to the n unknowns in the array q; · Ut is the m-component column matrix of the partial derivatives of U with respect to time t; · U_ q means the m by n matrix of total time-derivatives of each of the components of Uq ; and · U_ t is the m-component column matrix of total time-derivatives of each of the components of Ut . Also, for any n-component column matrix q and p-component column matrix r, the expression oq=or means the n by p rectangular array 2 oq
1
or1
6 6 oq2 6 or1 6 oq 6 oq3 6 or1 or 6 6 6 .. 6 . 4 oqn or1
oq1 or2
...
oq1 orp
oq2 or2
...
oq2 orp
oq3 or2
... ..
oqn or2
.
...
3
7 7 7 7 oq3 7 7 orp 7: 7 7 7 5
38
oqn orp
That is, oq=or is the Jacobian matrix of q with respect to r. 5.1.2. Matrix derivatives The de®nition (38) gives meaning to the partial dierentiation of one column matrix with respect to another. Thus, the quantity o=oq Uq q, for example, is well de®ned because Uq q is an m 1 matrix. By an extension to the product rule for dierentiation, the equation o oUq oq Uq q q Uq oq oq oq
39
forces a consistent de®nition for the ®rst term on the right-hand side where the other two terms in the equation are already well de®ned. According to
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
21
Eq. (38), the i; j component of the m n matrix on the left-hand side of equation Eq. (39) is
o Uq q oq
i;j
o
P
n oUi k1 oqk qk
oqj
n n X X o2 U i oUi oqk qk oq oq oq j k k oqj k1 k1
40
since Uq i;k
oUi oqk
41
and Ui is the ith component of U. By comparing Eqs. (39) and (40), we get
n X oUq o2 U i q q oq oqj oqk k i;j k1
42
for i 1; 2; . . . ; m and j 1; 2; . . . ; n. Note that we do not attempt to de®ne oUq =oq as a stand-alone quantity. Similarly, the i; j element of the m n matrix oUq =oqq_ is
oUq q_ oq
i;j
n X o2 U i q_ : oqj oqk k k1
43
These results will be useful in Section 5.3.3 where expressions are derived for the constant terms in the equations for the entries in the Jacobian matrix (46). In technical terms, the ith row of the matrix (42) or (43) is the transpose of _ rethe product of the Hessian matrix of the scalar function Ui times q or q, spectively. 5.2. The coordinate partitioning matrices The Jacobian matrix used in the Newton±Raphson algorithm to solve the system D for the displacements q is the n by n array Uq : 44 JD B The unknown dependent displacements are extracted from the solution q with the matrix B as shown in the diagram at the beginning of subsection 5.1. The expression (44) is also the coecient matrix for solving the system of linear equations V to get the ®rst derivatives of the dependent displacements.
22
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
Although the derivation of the matrix B is not within the scope of this report, for the sake of completeness, we mention the following possibility [5]. The coordinates can be partitioned by assigning the unknowns corresponding to the pivot elements in the LU decomposition with partial pivoting of the m by n rectangular Jacobian matrix Uq to the dependent set. The remaining entries in q are the independent coordinates. The n m by n m array JA
M Uq
ÿ UTq 0
45
is the coecient matrix for solving the system A/k for the accelerations and Lagrange multipliers. The matrices JD and JA are the generalizations of the matrices (12) and (13), respectively. 5.3. The Jacobian matrix for the implicit integrator The evaluation of the Jacobian matrix (6) for the BDF integrator is divided into four square submatrices just as before for the single-DOF example discussed in Section 4.3. That is, the Jacobian matrix used by the implicit algorithm to integrate the explicit coordinate partitioning function can be written as 2 3 o rout
o rout orin
out
or_out orin
6 or_in 4 or_ or_in
7 5
46
where each entry is a p p matrix. The input and output values of the p 1 _ and r contain the independent coordinates and their derivatives as arrays r; r, selected with Eq. (33). In the diagram of the coordinate partitioning function on page 13, the evaluation of r_ is independent of r. In fact, r_ is simply set equal to its own input value. Hence, or_out I or_in
47
or_out 0; orin
48
and
where I and 0 are the p p identity and zero matrices, respectively. Eqs. (47) and (48) supply the entries for half of the Jacobian matrix (46).
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
23
5.3.1. Simplifying relationships The derivation of the remaining entries in the Jacobian matrix is not as easy as the analysis needed above for Eqs. (47) and (48). Consequently, it is helpful to pause at this point to develop some simplifying relationships. The following results will be used to derive linear equations which can be solved for the coecients in the top half of the Jacobian matrix (46): (1) Since the values of the dependent coordinates s are determined by the system D in terms of r according to the diagram at the beginning of subsection 5.1 and because r and r_ are mutually independent in the initial value problem, we have or 0 or_
49
os 0: or_
50
and
Therefore, oq 0: or_
51
(2) Consider any dierentiable, scalar function of time and of the displacement variables q such as w w q; t:
52
The time-derivative of w is ow ow ow ow q_ 1 q_ 2 q_ n w_ oq1 oq2 oqn ot
53
from which we can see that ow_ ow oq_j oqj
54
for every j 1; 2; . . . ; n. Thus, we have oq_ oq or_ or
55
as well as other applications of the same principle (54) in the following derivations. The result Eq. (54) was applied previously to get Eq. (27). (3) From the de®nitions of r and q at the beginning of Section 5 and with Eq. (55), we see that
24
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
B
oq_ oq B I; or_ or
56
where I is the p p identity and B is de®ned by the relationship (33) which _ maps q into r as well as q_ into r. Eq. (56) is merely a restatement of some fundamental principles: It says that the derivative of each independent coordinate in q with respect to the corresponding element of r is equal to 1. Also, the derivatives of the independent coordinates in q with respect to all of the other elements of r are equal to 0. Since the rectangular array B is not an invertible matrix, Eq. (33) cannot imply anything about the dependent coordinates in q which are simply thrown away in the multiplication B times q. (4) The entries in the Jacobian matrix (46) and some of the terms in the following analyses involve partial derivatives of time-derivatives of the generalized coordinates q with respect to lower or the same order time-derivatives of q. If the coordinates were entirely independent, such matrices would all be diagonal. Of course, it is precisely the existence of the constraint equations which gives meaning to quantities such as o q=oq. Some algebraic manipulations of the map (33) of q into r yield the equalities or oq oq B B BT or or oq
57
and similar relationships between the various time-derivatives of q and r. It is helpful to work out the details of Eq. (57) for the example de®ned by Eqs. (34) and (35). The conclusions reached in this section will be incorporated as needed to obtain the derivatives with respect to r and r_ in the following developments. 5.3.2. Equations for the entries in the Jacobian matrix The results established in the previous paragraphs of Section 5 are sucient for the proof of the following statement: Theorem 1. The evaluation of the Jacobian matrix (46) can be achieved with just the coecient matrices (44) and (45) already decomposed for the coordinate partitioning function itself. Proof. Partial dierentiation of the system A/k, shown in the schematic diagram at the beginning of subsection 5.1., with respect to q_ and with respect to q yields the linear equations
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
"
2
M
ÿ UTq
Uq
0
and
M Uq
25
3
oF #2 oq 3 oq_ oq_ 6 7 7 4 56 4 _ _ 5 oU ok ÿ oq_ q q_ U_ q ooqq__ ooUq_ t oq_
58
2 3 oUT " # q oF o q k oq oq ÿ UTq 7 6 oq 4 _q _t 5 ok oU o U oq_ q oU 0 _ _ q q U ÿ oq q oq
oq
oq
59
oq
where, according to the de®nition (38), the unknowns o q=oq_ and o q=oq are n by n matrices and the unknowns ok=oq_ and ok=oq are m n rectangular arrays. The remaining rows of the Jacobian matrix (46), that is, the top p rows, form a proper subset of, possibly much smaller than, the solutions to the systems of linear equations (58) and (59). The super¯uous columns can be removed before computing the solutions by multiplying the equations from the right by BT to get 2 3 oF T " #" # B o q oq_ T M ÿ Uq 7 6 or_ 4 5 _q _t oU oq_ ok oU _ Uq 0 ÿ oq_ q_ Uq oq_ oq_ BT or_ 2 3 oF T B oq_ 5 60 4 oUq T t _ U_ q oq ÿ oq qB oU BT oq or and "
M Uq
ÿ
UTq
0
#" # o q or ok or
2
oUT q oq
kBT oF BT oq
3
6 7 7 6 _ 4 _ oU oU U t BT 5 ÿ oqq q oq q q_ U_ q ooqq_ ooq 2
oUT q oq
kBT oF BT oq
3
6 7 7 6 _ 4 _ oU oU U t BT 5 _ T U_ q oorq_ ooq ÿ oqq qBT oq q qB
61
where the relationships (54) and (57) have come into play. The solutions to these equations contain the remaining two p p components orout oq B or_in or_
62
and orout oq B orin or of the Jacobian matrix (46).
63
26
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
The coecient matrix for the two systems of equations (60) and (61) is inherited from the evaluation of the coordinate partitioning function. Thus, once the constant terms have been evaluated, the equations can be solved to form the Jacobian matrix to be used for the implicit integration of the function in the form (32). The integration of the coordinate partitioning function produces the solution to the dynamics problem (31). It remains to compute the missing ingredients of the constant terms for the systems (60) and (61). _ partial dierentiation of the system V with Since U does not depend on q, respect to r_ yields the linear equations " # " # 0 Uq oq_ ; 64 B or_ I where we have made use of the extended product rule discussed in Sec_ r_ oq=or in the system of equations (64) form tion 5.1.2. The unknowns oq=o an n p rectangular array. Eq. (64) can be solved in a straightforward manner for the p columns of corresponding to the p columns of constant unknowns 0 terms in the array . The matrix oq=or provides the missing data for the I constant terms in Eq. (60). Partial dierentiation of the system V with respect to q yields the equations "
Uq B
#
2 3 oUq oUt _ q ÿ oq_ 4 oq oq 5: oq 0
65
As done in a previous procedure, multiplying from the right by BT transforms Eq. (65) into 2 3 " # oU t Uq oq_ _ T oU ÿ oqq qB BT oq 5; 4 66 B or 0 _ where the unknowns oq=or form an n p rectangular array which provides the missing data for the constant terms in system (61). The proof of Theorem 1 is complete. _ r_ oq=or. Eq. (64) have p columns of unknowns in the n p matrix oq=o 0 The solutions correspond to the p columns in the constant term, , the rightI hand side of the equation, which is also an n p matrix. Similarly, there are p columns of unknowns in the n ´ p system (66) and in each of the n m p systems (60) and (61). In all cases, each column of the unknowns is the solution to the system with the corresponding column of constant terms from the right-
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
27
hand side. However, solving the latter three systems (66), (60), and (61) will require a little more work than needed for Eq. (64) because the constant terms are more complicated. The solutions to the two pairs of systems of linear equations (64) with (60) and (66) with (61) provide the values of the entries in the top half of the Jacobian matrix (46). There are two steps, in the terminology used in the ®eld of numerical linear algebra, to the solution procedure for such systems of equations: 1. Decomposition (or factorization) of the matrix of coecients which is done just once for a given array of values and 2. Solving (the forward reduction and backward substitution calculations) for the unknowns that correspond to each column array of the constant terms. Although the ®rst step is generally quite computationally intensive, much more so than the second, the ecacy of the implicit method is quite favorable for integrating the coordinate partitioning function. The four systems of equations listed above share the two coecient matrices (44) and (45) which were already required for the evaluation of the function as shown in the schematic diagram at the beginning of subsection 5.1. Hence, the ®rst step of the solution process, the decompositions of the matrices, need not be repeated for the current application. Even though each of the four systems has to be solved a number of times, the additional work is not large compared to the function evaluation. 5.3.3. Expressions for the constant terms The basic elements of the constant terms for Eqs. (66), (60) and (61) must already exist for the evaluation of the coordinate partitioning function which consists of the sequential solutions to the non linear and linear systems D, V, and A/k in the diagram at the beginning of subsection 5.1. It remains to construct the formulas for the constant terms from the basic components. While the complexity of the constant terms appears daunting at ®rst, a number of simpli®cations apply to the constraint equations associated with the typical mechanical system (31). In particular, the standard mechanical joints have no time dependency in which case the corresponding derivatives of U with respect to t are equal to zero. For the sake of completeness, however, we conclude Section 5 in the following paragraphs with the fully general expressions for all of the components of the constant terms. For Eq. (66), we need two formulas: (1) The matrix oUq =oqq_ is given by Eq. (43) derived in Section 5.1.2. (2) The partial time-derivative of U is
28
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
2
oU1 ot
3
7 6 6 oU2 7 6 ot 7 7 6 Ut 6 7; 6 .. 7 6 . 7 5 4
67
oUm ot
where Ui is the ith element in the column matrix U of expressions for the constraint functions. Hence, the i; j component of oUt =oq is
oUt oq
i;j
o2 Ui oqj ot
68
for i 1; 2; . . . ; m and j 1; 2; . . . ; n. Of course, if the ith component of U does not depend on t, then oUi =ot 0. Indeed, we have Ut 0 in many applications. _ T and oUt =oqBT should not be computed as The products oUq =oqqB matrix-multiplies. Rather, multiplication from the right by BT amounts to the selection of the p columns corresponding to the p independent variables from the m ´ n matrices oUq =oqq_ and oUt =oq. There are three expressions involving the constraints U in the constant terms for Eq. (60): _ as already noted above in item 1 for Eq. (66), was 1. The expression oUq =oqq, derived in Section 5.1.2. It is is given by Eq. (43). 2. The m n Jacobian matrix Uq consists of the partial derivatives of U with respect to the displacement variables q. Thus, as stated in Eq. (41), the i; j component of Uq is oUi =oqj where Ui is the ith element in the column matrix U of the expressions for the constraints. Since the matrix U_ q is the term by term time-derivative of Uq , the i; j component of this m n array is oU_ i ; U_ q i;j oqj
69
where it is assumed that the order of the time-derivatives o=ot and d=dt is interchangeable. The matrix-product U_ q oq=or is computed by multi_ r_ to plying the matrix de®ned by Eq. (69) times the solution oq=or oq=o Eq. (64). If the ith component of U does not depend on t, then the entire ith row of U_ q is equal to zero. That is, oU_ i 0 oqj
70
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
29
for j 1; 2; . . . ; n. However, if the ith constraint is time-dependent, then oU_ i o2 Ui o2 Ui o2 U i o2 Ui q_ 1 q_ 2 . . . q_ n oqj oq1 oqj oq2 oqj oqn oqj otoqj
71
for each j 1; 2; . . . ; n. That is, the dot indicates the total time derivative of the expression.3. The m by n matrix oUt =oq is de®ned above by Eq. (68) in item 2 for Eq. (66).As described above at the end of the paragraph on the constant terms for Eq. (66), multiplication from the right by BT constitutes the selection of p columns from the various m n matrices. To get all of the constant terms involving U in Eq. (61), we have to develop ®ve more formulas: 1. The construction of an expression analogous to Eqs. (42) and (43) yields oUTq k oq
! i;j
n X o2 Uk kk oqj oqi k1
72
for the i; j component of the n n square matrix oUTq =oqk. q is de®ned by Eq. (42) which was derived in Sec2. The m n matrix oUq =oq tion 5.1.2. 3. As in item 1, we need to construct an expression analogous to Eqs. (42) and _ Thus, using the de®nition (69) of the matrix U_ q , (43) for the term oU_ q =oqq. we obtain oU_ q q_ oq
! i;j
n X o2 U_ i q_ oqj oqk k k1
73
for i 1; 2; . . . ; m and j 1; 2; . . . ; n. Once again, if the ith component of U _q oU 0 for j 1; 2; . . . ; n. Otherwise, we is independent of time, then oq q_ i;j
have to substitute the rather messy expression n oU_ i X o2 Ui o2 Ui q_ l oqk oql oqk otoqk l1
74
into Eq. (73). _ is computed by multiplying the matrix de4. The matrix-product U_ q oq=or ®ned by Eq. (69) times the solution to Eq. (66). 5. The column matrix Ut is de®ned by Eq. (67). The time-derivative of Ut , still a column matrix with m components, is given by
30
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
2
_1 oU ot
3
7 6 6 oU_ 2 7 7 6 ot 7 U_ t 6 6 . 7 6 . 7 4 . 5 _m oU ot
2 o2 U
1
oq1 ot
q_ 1
6 o2 U 6 2 q_ 1 6 oq ot 6 1 6 4 o2 Um q_ oq1 ot 1
o2 U1 q_ oq2 ot 2
...
... .. . ...
o2 U2 q_ oqn ot n
o2 U1 ot2 o2 U2 ot2
o2 Um q_ oqn ot n
o2 Um ot2
o2 U2 q_ oq2 ot 2
o2 Um q_ oq2 ot 2
o2 U1 q_ oqn ot n
3 7 7 7 7; 7 5
75
where it is assumed that the order of the time-derivatives o=ot and d=dt is interchangeable. Of course, the components of U_ t corresponding to the timeindependent entries in U are equal to zero. As for the previous two cases, multiplication from the right by BT constitutes the selection of p columns from the n n matrix in item 1 and p columns from the m ´ n matrices in the other items. Finally, the partial derivatives of the applied forces oF =oq_ and oF =oq are de®ned in a manner which is analogous to Eq. (38) in Section 5.1.1. Multiplication from the right by BT constitutes the selection of p columns from these two n n matrices. 6. Future developments If the forces in the system of DAE's (31) are acceleration-dependent, which occurs, for example, when a mechanism is modeled with resistance to non zero accelerations, then the equations take the form _ q; t; q; q; M q ÿ UTq k F U q; t 0;
76
where some components of the forces F are functions of the second derivatives q of the coordinates. Most of the analysis in Section 5 still applies to the problem (76). The systems D and V in the coordinate partitioning function on page 13 remain the same. However, the system A/k becomes quite dierent. It may even be non linear. Another kind of complication develops when the system of DAE's (31) takes the form _ q; t; M q ÿ UTq k F q; _ q; t 0; U q;
77
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
31
where the algebraic equation depend on the velocities, a condition which characterizes non holonomic constraints. The system of equations D in the coordinate partitioning function is no longer independent of the system V. Hence, the equation D and V have to be solved simultaneously, for which the analysis presented in Section 5 will have to be extensively reworked. Finally, both acceleration-dependent forces and non holonomic constraints can be incorporated in the same mechanical system, which also means the previous results no longer apply. The development of the details for the formation of the Jacobian matrices corresponding to problems (76) and (77) and to problems with both complications awaits a future report. References [1] L.F. Shampine, H.A. Watts, DEPAC-Design of a user- oriented package of ODE solvers, Sand79-2374, Sandia National Laboratories, 1979. [2] L.F. Shampine, Numerical Solutions of Ordinary Dierential Equations, Chapman & Hall, New York, 1994. [3] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Dierential-Algebraic Equations, North Holland, New York, 1989. [4] C.W. Gear, Numerical Initial Value Problems in Ordinary Dierential Equations, PrenticeHall, Englewood Clis, NJ, 1971. [5] R.A. Wehage, E.J. Haug, Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic system, J. Mechanical Design 104 (1982) 247±255.
Implicit integration with coordinate partitioning Joseph F. McGrath
*,1,
Rajiv Rampalli
Mechanical Dynamics, Inc., Ann Arbor, MI, USA
Abstract The dynamics of a mechanical system are governed by dierential and algebraic (constraint) equations. If these DAE's are reduced by coordinate partitioning to a system of ®rst order dierential equations for the independent unknowns, then numerical methods for equations in explicit form can be used to compute the solution. We have developed techniques for the evaluation of the sparse Jacobian matrices for solving the partitioned equations with integrators based on backward dierentiation formulas, an important method for the simulation of the inherently sti mechanical systems. Ó 1999 Elsevier Science Inc. All rights reserved. AMS classi®cation: 34A50; 34A65; 65C20; 65L05 Keywords: Coordinate partitioning; Dynamics; Dierential and algebraic equations; Mechanical system simulation
1. Introduction and summary Coordinate partitioning is a numerical procedure for transforming the system of dierential and algebraic equations that governs mechanical motion into an initial value problem consisting of ordinary dierential equations. The reformulation de®ned by the method explicitly evaluates the time-derivatives of the independent coordinates in the problem. Thus, the unknowns can be integrated to obtain the solution by computer programs for systems of explicit dierential equations. The integrators may employ either explicit or implicit * Corresponding author. Silicon Graphics, Inc., 102 Central Park Square, Los Alamos, NM 87544, USA. Tel.: 505 662 4000; fax: 505 662 9457; e-mail: [email protected]. 1 Adjunct member of the faculty, Department of Mathematics, University of Michigan.
0096-3003/00/$ - see front matter Ó 2000 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 9 8 ) 1 0 0 2 3 - 1
8
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
methods. However, more information than just the time-derivatives of the unknowns must be generated for the techniques implemented in the implicit integrators. The algorithm in the implicit integrators based on backward dierentiation formulas needs the Jacobian matrix of the time-derivatives of the unknowns with respect to the unknowns themselves to advance the solution from one time step to the next. As described in this paper, it turns out that, in addition to evaluating the time-derivatives, the coordinate partitioning algorithm provides the machinery for the formation of the required Jacobian matrix. Speci®cally, the matrices that are used in the coordinate partitioning function to calculate the derivatives are the same as the coecient matrices needed to compute the various components of the Jacobian matrix for the implicit integrator. Background material and the context for the discussion of the coordinate partitioning algorithm are included in Sections 2 and 3. A review of the coordinate partitioning algorithm is contained in Section 4 where a basic example is developed to illustrate the derivation of the Jacobian matrix for an implicit integrator. The general concepts involved in implicit integration with coordinate partitioning are covered by the end of Section 4. For the reader interested in the full details, Section 5 contains the derivation of the Jacobian matrix for the general case. Future developments for acceleration-dependent forces and non holonomic constraints are mentioned in Section 6. 2. Explicit dierential equations A system of ordinary dierential equations is classi®ed as explicit if it has the form w_ f w; t 1 as opposed to a system in the form _ w; t 0; f w; where the derivatives are implicitly de®ned. In either case, the array 2 3 f1 6 f2 7 6 7 f 6 . 7 4 .. 5
2
3
fN
consists of N functions of the N unknowns in the array 2 3 w1 6 w2 7 6 7 w 6 .. 7: 4 . 5 wN
4
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
For the implicit equations (2), f is a function of both w and 2 3 w_ 1 6 w_ 2 7 6 7 w_ 6 .. 7: 4 . 5
9
5
w_ N
Finally, each component of f may also be a function of time t. The following schematic diagram of a function evaluation illustrates the explicit form (1) of a system of dierential equations. The essential
characteristic is the functional nature of the system. No information about the calculations themselves is contained in the output from an evaluation of the dierential equations. The coordinate partitioning algorithm satis®es the above criteria for a system of explicit dierential equations: The implicit system (2) is converted into an explicit system by concealing the implicit nature of the equations within the algorithmic box in the diagram and by computing the derivatives w_ internally as a function of the input values of w and t.
3. Methods for explicit initial value problems An important category of methods for solving initial value problems (IVP's) needs the system of equations to be in the explicit form (1). The various methods in this category of solvers include both explicit and implicit integrators. An explicit integrator uses current and previous information about the solution to compute its values at the next time step in a ®xed number of function evaluations. In contrast, an implicit integrator employs an iterative
10
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
algorithm (implying an unknown, variable number of function evaluations to achieve convergence) to solve the equations in either form (1) or (2) at the future time. For example, the DEPAC software package [1,2] contains the following three solvers: DDERKF DDEABM DDEBDF
The single step, explicit, fourth order Runge±Kutta±Fehlberg algorithm with ®fth order corrections; A multi-step, variable order, explicit Adams±Bashforth±Moulton integrator (predict, evaluate, correct, evaluate); A multi-step, variable order, implicit method based on backward dierentiation formulas.
Each one of the integrators on the above list expects the problem to be de®ned in the explicit form (1). The DDERKF and DDEABM integrators employ explicit algorithms while an implicit method is implemented in the DDEBDF integrator. Other software products, such as the DASSL code, use implicit methods to solve IVP's in the implicit form (2) [3,4]. To compute the solution, an explicit integrator requires only the evaluation of the derivatives in the form of Eq. (1) as illustrated by the schematic diagram on page 3. The algorithm does not need any details about the calculations themselves. However, an implicit integrator makes use of additional information about the behavior of the function f . In particular, the implicit integrators DDEBDF and DASSL need the Jacobian matrix 2 of1 of1 of1 3 . . . ow ow1 ow2 N 7 6 of of2 of2 7 6 2 6 ow1 ow2 . . . owN 7 7 6 6 7 6 . .. 7 6 .. . 5 4 ofN ofN ofN . . . ow ow1 ow2 N for the vector function f in Eq. (1). The matrix (6) provides the implicit algorithm with quanti®ed information about the sensitivity of each component of the output to each component of the input for an explicit system of equations. 4. Coordinate partitioning The system of n second order dierential equations and m algebraic equations _ q; t M q ÿ UTq k F q; U q; t 0
7
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
11
governs the dynamics of a mechanical system. Eq. (7) do not form an explicit system of ordinary dierential equations in their present form. However, as explained in the following paragraphs, the coordinate partitioning algorithm performs a numerical transformation of the system (7) into an explicit initial value problem [5]. In the above system of dierential and algebraic equations (DAE's), the ncomponent array q consists of the translational and rotational coordinates that describe the current con®guration of the mechanism. The m Lagrange multipliers in the array k correspond to the m constraint equations represented by U. The coordinate partitioning software has three components: 1. The algorithm selects a p-element vector of independent coordinates from the vector q where p 6 n is equal to the number of degrees of freedom (DOF's) in the mechanical system governed by the dierential and algebraic equations (7); 2. It evaluates f the vector-function of the N 2p unknowns consisting of the p independent coordinates plus their ®rst derivatives; 3. It computes all of the other variables in the problem including the dependent variables, their derivatives, and the Lagrange multipliers. For item 2, the derivatives w_ f w; t are calculated by successively solving three systems of equations, one or two of which are non linear. In addition to the derivatives, implicit integrators usually require the Jacobian matrix Eq. (6) of the derivatives of the unknowns with respect to the unknowns themselves. For the coordinate partitioning algorithm, the entries in the Jacobian matrix are computed by solving two pairs of linear equations. The coecient matrices for the four systems of equations are the same as two of the matrices used for the evaluation of the function f . The expressions for the constant terms (the right-hand sides) involve additional derivations of a non trivial level of complexity. Another unfortunate complication is the reason for the rather technical nature of the subject matter in Section 5 of this report. There is no transformation of coordinates involved in the coordinate partitioning algorithm. The p independent coordinates in the system (7) are simply a selection of the entries in q. If m < n and the m equations that make up U are independent, then the number of independent variables is p n ÿ m and the dynamic motion of the system has p DOF's. If m n and the constraint equations are independent, then p 0 and the mechanism can only undergo kinematic motion. If m > n and more than n of the constraint equations are independent, the system is overconstrained. The latter two cases are outside of the scope of the current discussion. The coordinate partitioning algorithm is used only for the simulation of dynamic motion, that is, p > 0.
12
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
4.1. The point-mass pendulum A complete derivation of the coordinate partitioning equations for a basic mechanism with one moving part and the corresponding Jacobian matrix for a BDF integrator is contained in the following paragraphs. The eort spent developing the details for this non general case serves two purposes: It lays the groundwork needed by the authors to construct the general formulation to be developed in Section 5 and it provides an explanation of the method which is not obscured by the more abstract notation. Thus, by reading Section 4, it is possible for a mathematician or computational engineer to understand how the method works without mastering all of the details contained in Section 5. A mechanism made up of a single moving part attached by a revolute joint to a ®xed part (to ground) provides a good illustration of the coordinate partitioning algorithm. This two dimensional, single-DOF problem is governed by the system of dierential and algebraic equations (DAE's) _ x; y; t mx ÿ xk F1 _x; y; _ x; y; t m y ÿ yk F2 _x; y; 1 2 x y 2 ÿ 1 0 2
8
where · The point-mass m is at location x; y a ®xed, unit distance from the origin; · F1 and F2 are the horizontal and vertical components of the applied forces, respectively; · The unknown k is a Lagrange multiplier. For example, if F1 0 and F2 ÿmg where g is the acceleration due to gravity, then the DAE's (8) govern the motion of a simple pendulum with point-mass m. The constraint equation in the system (8) is scaled by 1=2 for the sake of conformity with the standard representation Eq. (7) for DAE's governing the dynamics of a mechanism. 4.2. The coordinate partitioning function Suppose the partitioning algorithm has selected x as the independent coordinate leaving y as the dependent coordinate 2 for the problem (8). Then, the input to the coordinate partitioning function is x_ w 9 x
2
For the current example, we will simply assume the existence of a satisfactory partitioning algorithm. A description of the process is included in Section 5.2.
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
13
and the expected output consists of the computed components of x w_ f w; t 10 x_ which is in the explicit form (1). The algorithm for computing x and x_ in terms of x_ and x is de®ned in the following schematic diagram where, for the sake of clarity, the ®nal calculations of the two output values are enclosed in boxes:
The D, V, and A/k equations in this basic example form a nonlinear system of seven equations in the seven unknowns 2 3 x 6y 7 6 7 6 x_ 7 6 7 6 y_ 7: 11 6 7 6 x 7 6 7 4 y 5 k
14
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
The equations for the models of mechanical systems can be uncoupled in a manner which improves the eciency of the calculations. That is, the solution can be computed sequentially in the coordinate partitioning code. The output from the solution to the ®rst block of equations D is part of the input data for the subsystem V. The cumulative results from solving the subsystems D and V along with the values of x_ x_ in and x xin make up the input for solving the subsystem A/k. In general, as illustrated by this example, if values are speci®ed for the independent displacements and their ®rst derivatives, the solutions to each of the subsystems, D, V, and A/k, yield the dependent displacements, the dependent velocities, and the accelerations and constraint forces (Lagrange multipliers), in that order. Hence, there are two important and very valuable advantages to the coordinate partitioning method: 1. While the integrator advances the independent displacements and velocities, the coordinate partitioning calculations simultaneously provide consistent values for all of the other entries in the state vector for the mechanism; and 2. The solution satis®es not only the dierential equations and the constraint equations but also the ®rst and second derivatives of the constraint equations. In summary, the coordinate partitioning method yields a complete state vector for the mechanism consisting of the displacements, the ®rst and second derivatives, and the constraint forces at each integration time step. In the current single-DOF example, this state vector is the column matrix (11). The Jacobian and coecient matrices used to solve the systems of equations in the function evaluation de®ned by the schematic diagram shown above are x y 12 1 0 for the dependent displacements D and the velocities V and 2 3 m 0 ÿx 6 7 40 m ÿy5 x y 0
13
for the accelerations and constraint forces A/k. Recall that the ®rst component of the algorithm, as described on page 5, calls for the partitioning of the coordinates. Although we have not said anything yet about the selection process, the roles of x and y could be reversed with y as the independent variable, then the Jacobian matrix (12) for the dependent displacements would be x y 14 0 1 instead of the matrix (12). The coecient matrix (13) would not change.
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
15
4.3. The Jacobian matrix for the BDF integrator Numerically sti initial value problems, that is, initial value problems with widely separated time constants in the various components of the solution, can be solved with implicit integration techniques such as the methods based on backward dierentiation formulas (BDF's). The implicit BDF integrator for explicit dierential equations in the form (1) requires the Jacobian matrix (6). Thus, for the coordinate partitioning algorithm, the additional Jacobian matrix 3 " # 2 of1 o_x of2 o_x
of1 ox of2 ox
oxout o_x
oxout oxin
out
o_xout oxin
4 o_x in o_xin
5
15
is needed to solve the problem (8). From the de®nition of the coordinate partitioning function in the schematic diagram on page 7, we can see that the output value of x_ is completely independent of the input value of x. Indeed, x_ is simply set equal to its own input value. Hence, ox_ out 1 ox_ in
16
ox_ out 0: oxin
17
and
which takes care of the second row of the 2 ´ 2 Jacobian matrix (15). Although the derivation of the entries in the ®rst row is not quite as easy as the second row of the matrix (15), the calculations can still be done quite ef®ciently as will be shown in the following discussion. In fact, the evaluation of the Jacobian matrix (15) can be achieved with just the coecient matrices (12) and (13) needed for the evaluation of the function itself. The signi®cance of this theorem (which will be proven for the general case in Section 5) lies, of course, with the fact that the (LU, QR, or perhaps another) decomposition of the matrices (12) and (13) need not be repeated for the evaluation of the requisite Jacobian matrix (15). In the following work, partial dierentiation of the subsystems V and A/k with respect to x_ will yield a system of equations for the ®rst entry in the ®rst row of the matrix (15). To get a system of equations for the second entry, we will need the partial derivatives of the subsystems V and A/k with respect to x.
16
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
4.3.1. First row, ®rst entry Since y is dependent only upon x and x is an independent variable, we have ox oy 0: ox_ ox_
18
Hence, partial dierentiation of the ®rst equation in the system V from the diagram is subsection 4.2 with respect to x_ yields xy
oy_ 0 ox_
_ x easily enough. which can be solved for oy=o_ If Eqs. (16) and (19) are collected into the linear system " #2 o_x 3 " # 0 x y 4 o_x 5 ; oy_ 1 1 0 o_x
19
20
it might appear that we have complicated matters. However, in order to illustrate the more general case involving multiple dependent and multiple independent coordinates, the dierentiation of the system V must be written in this more elaborate fashion, rather than as the individual equations that apply to the single-DOF example Eq. (8) currently under consideration. In particular, it is signi®cant that the coecient matrix for the system (20) is the Jacobian matrix (12) for the dependent displacements. Partial dierentiation of the system A/k with respect to x_ yields the system of linear equations 2 32 ox 3 2 oF1 oF1 oy_ 3 oy_ o_x m 0 ÿx o_x o_x 7 7 6 6 76 6 oF2 oF2 oy_ 7 o y7 6 0 m ÿ y 76 21 6 6 7 4 54 o_x 5 4 o_x oy_ o_x 7 5 ok x y 0 ÿ 2_x 2y_ oo_yx_ o_x where the constant terms on the right-hand side and the coecients have already been determined. One component of the solution to Eq. (21) is of1 ox : 22 ox_ ox_ We have now evaluated three quarters of the Jacobian matrix (15). Note that, the coecient matrix for the system (21) is the coecient matrix (13) for the accelerations and Lagrange multipliers. 4.3.2. First row, second entry Partial dierentiation of the ®rst equation in each of the systems D and V from the diagram on page 7 with respect to x yields
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
17
y
oy ÿx ox
23
y_
oy oy_ y ÿ_x; ox ox
24
and
respectively. Eqs. (23) and (24) can be solved directly ®rst for oy=ox and then _ for oy=ox. However, once again for the sake of the more general formulation, the values of these quantities have to be obtained by sequentially solving the system " #" ox # " # 0 x y ox 25 oy 1 1 0 ox and, using Eq. (17), the system # " #" o_x # " x y ÿ _x y_ oy ox ox : 1 0 ooxy_ 0
26
The computed value of oy=ox is needed to solve the Eq. (26). However, it turns out that the system (25) is super¯uous. A comparison of the two systems of linear equations (20) and (25) proves that oy_ oy 27 ox_ ox which, as will be shown in subsection 5.3.1, is true in general. Consequently, the work load can be lightened by eliminating the system (25) and using the results from Eqs. (20) and (27) to solve the system (26). Finally, partial dierentiation of the system A/k with respect to x yields 2 32 ox 3 2 k oF1 oF1 oy 3 m 0 ÿx ox oy ox ox 7 7 6 oy oF 6 76 o y7 2 6 7 6 28 4 0 m ÿ y 54 ox 5 4 k ox ox oFoy2 oy ox 5 ok x y 0 ÿ x y oy 2y_ oy_ ox ox
ox
for which the solution requires the results from Eqs. (26) and (27). One of the components of the solution to the system of equations (28) is of1 ox : 29 ox ox Once again note that the coecient matrices in the systems (25), (26), and (28) are the matrices (12) and (13) from the coordinate partitioning function for the time-derivatives. Substituting Eqs. (22), (16), (29) and (17) into the de®nition Eq. (15) yields the Jacobian matrix
18
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
" of
1
o_x of2 o_x
of1 ox of2 ox
"
#
oxout o_xin
1
oxout oxin
0
# 30
where the ®rst row has been computed by solving Eqs. (21) and (28). The evaluation of the Jacobian matrix (15) is now complete. An implicit integrator, such as the DDEBDF program, can use the coordinate partitioning function de®ned in the diagram in subsection 4.2 and the Jacobian matrix (30) to compute the solution to the problem (8). 5. The general case The constraint equations in the system (7), repeated here for convenience, _ q; t M q ÿ UTq k F q;
31 U q; t 0 reduce the number of degrees of freedom in the motion of the mechanism by imposing interdependencies between some of the variables. Assume that there are n dierential equations and m < n constraint equations (the algebraic equations) in the system (31). Also, assume that the m constraints are mutually independent in which case there are p n ÿ m mechanical degrees of freedom in the problem. Let B be the p n matrix of 0's and 1's that selects the p independent variables from the array q of the n translational and rotational coordinates in the problem. Thus, for the explicit initial value problem (10), also repeated here for convenience, w_ f w; t; 32 the array of unknowns is r_ Bq_ w 33 r Bq where r is the p-component column matrix of independent coordinates. The remaining m n ÿ p dependent coordinates will constitute the s array. For the moment, we take this partitioning process for granted. A few words about the formation of B are included toward the end of Section 5.2. In the present context, a couple of examples are sucient to explain the relationship Eq. (33). If y and w are the independent variables in the system of coordinates 2 3 x 6y7 6 7 6z7 7 34 q6 6 w 7; 6 7 4h5 /
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
then
1 0 0 0 0 1 y r Bq ; w B
and
0 0
0 0
0 ; 0
2 3 x 6z7 7 s6 4 h 5: /
19
35 36
37
For the second illustration, the single-DOF example discussed in Section 4 has B 1 0; r x, and s y. 5.1. The general coordinate partitioning function For the general case Eq. (31), the following schematic diagram de®nes the coordinate partitioning function where, for the sake of clarity, the ®nal expressions for the output values are enclosed in boxes:
20
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
Note that the same two signi®cant advantages listed in subsection 4.2 for the basic example pertain to the algorithm de®ned in this ®gure. That is, the coordinate partitioning method yields the values of the complete state vector for the mechanism including the displacements, their ®rst and second derivatives, and the constraint forces, all of which satisfy the DAE's (31) plus the ®rst and second derivatives of the constraint equations. 5.1.1. Notation In the notation used in the above schematic diagram and throughout Section 5 of this report, · U represents the m constraint equations in the form U q; t 0; · Uq is the m n (generally sparse) Jacobian matrix of U with respect to the n unknowns in the array q; · Ut is the m-component column matrix of the partial derivatives of U with respect to time t; · U_ q means the m by n matrix of total time-derivatives of each of the components of Uq ; and · U_ t is the m-component column matrix of total time-derivatives of each of the components of Ut . Also, for any n-component column matrix q and p-component column matrix r, the expression oq=or means the n by p rectangular array 2 oq
1
or1
6 6 oq2 6 or1 6 oq 6 oq3 6 or1 or 6 6 6 .. 6 . 4 oqn or1
oq1 or2
...
oq1 orp
oq2 or2
...
oq2 orp
oq3 or2
... ..
oqn or2
.
...
3
7 7 7 7 oq3 7 7 orp 7: 7 7 7 5
38
oqn orp
That is, oq=or is the Jacobian matrix of q with respect to r. 5.1.2. Matrix derivatives The de®nition (38) gives meaning to the partial dierentiation of one column matrix with respect to another. Thus, the quantity o=oq Uq q, for example, is well de®ned because Uq q is an m 1 matrix. By an extension to the product rule for dierentiation, the equation o oUq oq Uq q q Uq oq oq oq
39
forces a consistent de®nition for the ®rst term on the right-hand side where the other two terms in the equation are already well de®ned. According to
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
21
Eq. (38), the i; j component of the m n matrix on the left-hand side of equation Eq. (39) is
o Uq q oq
i;j
o
P
n oUi k1 oqk qk
oqj
n n X X o2 U i oUi oqk qk oq oq oq j k k oqj k1 k1
40
since Uq i;k
oUi oqk
41
and Ui is the ith component of U. By comparing Eqs. (39) and (40), we get
n X oUq o2 U i q q oq oqj oqk k i;j k1
42
for i 1; 2; . . . ; m and j 1; 2; . . . ; n. Note that we do not attempt to de®ne oUq =oq as a stand-alone quantity. Similarly, the i; j element of the m n matrix oUq =oqq_ is
oUq q_ oq
i;j
n X o2 U i q_ : oqj oqk k k1
43
These results will be useful in Section 5.3.3 where expressions are derived for the constant terms in the equations for the entries in the Jacobian matrix (46). In technical terms, the ith row of the matrix (42) or (43) is the transpose of _ rethe product of the Hessian matrix of the scalar function Ui times q or q, spectively. 5.2. The coordinate partitioning matrices The Jacobian matrix used in the Newton±Raphson algorithm to solve the system D for the displacements q is the n by n array Uq : 44 JD B The unknown dependent displacements are extracted from the solution q with the matrix B as shown in the diagram at the beginning of subsection 5.1. The expression (44) is also the coecient matrix for solving the system of linear equations V to get the ®rst derivatives of the dependent displacements.
22
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
Although the derivation of the matrix B is not within the scope of this report, for the sake of completeness, we mention the following possibility [5]. The coordinates can be partitioned by assigning the unknowns corresponding to the pivot elements in the LU decomposition with partial pivoting of the m by n rectangular Jacobian matrix Uq to the dependent set. The remaining entries in q are the independent coordinates. The n m by n m array JA
M Uq
ÿ UTq 0
45
is the coecient matrix for solving the system A/k for the accelerations and Lagrange multipliers. The matrices JD and JA are the generalizations of the matrices (12) and (13), respectively. 5.3. The Jacobian matrix for the implicit integrator The evaluation of the Jacobian matrix (6) for the BDF integrator is divided into four square submatrices just as before for the single-DOF example discussed in Section 4.3. That is, the Jacobian matrix used by the implicit algorithm to integrate the explicit coordinate partitioning function can be written as 2 3 o rout
o rout orin
out
or_out orin
6 or_in 4 or_ or_in
7 5
46
where each entry is a p p matrix. The input and output values of the p 1 _ and r contain the independent coordinates and their derivatives as arrays r; r, selected with Eq. (33). In the diagram of the coordinate partitioning function on page 13, the evaluation of r_ is independent of r. In fact, r_ is simply set equal to its own input value. Hence, or_out I or_in
47
or_out 0; orin
48
and
where I and 0 are the p p identity and zero matrices, respectively. Eqs. (47) and (48) supply the entries for half of the Jacobian matrix (46).
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
23
5.3.1. Simplifying relationships The derivation of the remaining entries in the Jacobian matrix is not as easy as the analysis needed above for Eqs. (47) and (48). Consequently, it is helpful to pause at this point to develop some simplifying relationships. The following results will be used to derive linear equations which can be solved for the coecients in the top half of the Jacobian matrix (46): (1) Since the values of the dependent coordinates s are determined by the system D in terms of r according to the diagram at the beginning of subsection 5.1 and because r and r_ are mutually independent in the initial value problem, we have or 0 or_
49
os 0: or_
50
and
Therefore, oq 0: or_
51
(2) Consider any dierentiable, scalar function of time and of the displacement variables q such as w w q; t:
52
The time-derivative of w is ow ow ow ow q_ 1 q_ 2 q_ n w_ oq1 oq2 oqn ot
53
from which we can see that ow_ ow oq_j oqj
54
for every j 1; 2; . . . ; n. Thus, we have oq_ oq or_ or
55
as well as other applications of the same principle (54) in the following derivations. The result Eq. (54) was applied previously to get Eq. (27). (3) From the de®nitions of r and q at the beginning of Section 5 and with Eq. (55), we see that
24
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
B
oq_ oq B I; or_ or
56
where I is the p p identity and B is de®ned by the relationship (33) which _ maps q into r as well as q_ into r. Eq. (56) is merely a restatement of some fundamental principles: It says that the derivative of each independent coordinate in q with respect to the corresponding element of r is equal to 1. Also, the derivatives of the independent coordinates in q with respect to all of the other elements of r are equal to 0. Since the rectangular array B is not an invertible matrix, Eq. (33) cannot imply anything about the dependent coordinates in q which are simply thrown away in the multiplication B times q. (4) The entries in the Jacobian matrix (46) and some of the terms in the following analyses involve partial derivatives of time-derivatives of the generalized coordinates q with respect to lower or the same order time-derivatives of q. If the coordinates were entirely independent, such matrices would all be diagonal. Of course, it is precisely the existence of the constraint equations which gives meaning to quantities such as o q=oq. Some algebraic manipulations of the map (33) of q into r yield the equalities or oq oq B B BT or or oq
57
and similar relationships between the various time-derivatives of q and r. It is helpful to work out the details of Eq. (57) for the example de®ned by Eqs. (34) and (35). The conclusions reached in this section will be incorporated as needed to obtain the derivatives with respect to r and r_ in the following developments. 5.3.2. Equations for the entries in the Jacobian matrix The results established in the previous paragraphs of Section 5 are sucient for the proof of the following statement: Theorem 1. The evaluation of the Jacobian matrix (46) can be achieved with just the coecient matrices (44) and (45) already decomposed for the coordinate partitioning function itself. Proof. Partial dierentiation of the system A/k, shown in the schematic diagram at the beginning of subsection 5.1., with respect to q_ and with respect to q yields the linear equations
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
"
2
M
ÿ UTq
Uq
0
and
M Uq
25
3
oF #2 oq 3 oq_ oq_ 6 7 7 4 56 4 _ _ 5 oU ok ÿ oq_ q q_ U_ q ooqq__ ooUq_ t oq_
58
2 3 oUT " # q oF o q k oq oq ÿ UTq 7 6 oq 4 _q _t 5 ok oU o U oq_ q oU 0 _ _ q q U ÿ oq q oq
oq
oq
59
oq
where, according to the de®nition (38), the unknowns o q=oq_ and o q=oq are n by n matrices and the unknowns ok=oq_ and ok=oq are m n rectangular arrays. The remaining rows of the Jacobian matrix (46), that is, the top p rows, form a proper subset of, possibly much smaller than, the solutions to the systems of linear equations (58) and (59). The super¯uous columns can be removed before computing the solutions by multiplying the equations from the right by BT to get 2 3 oF T " #" # B o q oq_ T M ÿ Uq 7 6 or_ 4 5 _q _t oU oq_ ok oU _ Uq 0 ÿ oq_ q_ Uq oq_ oq_ BT or_ 2 3 oF T B oq_ 5 60 4 oUq T t _ U_ q oq ÿ oq qB oU BT oq or and "
M Uq
ÿ
UTq
0
#" # o q or ok or
2
oUT q oq
kBT oF BT oq
3
6 7 7 6 _ 4 _ oU oU U t BT 5 ÿ oqq q oq q q_ U_ q ooqq_ ooq 2
oUT q oq
kBT oF BT oq
3
6 7 7 6 _ 4 _ oU oU U t BT 5 _ T U_ q oorq_ ooq ÿ oqq qBT oq q qB
61
where the relationships (54) and (57) have come into play. The solutions to these equations contain the remaining two p p components orout oq B or_in or_
62
and orout oq B orin or of the Jacobian matrix (46).
63
26
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
The coecient matrix for the two systems of equations (60) and (61) is inherited from the evaluation of the coordinate partitioning function. Thus, once the constant terms have been evaluated, the equations can be solved to form the Jacobian matrix to be used for the implicit integration of the function in the form (32). The integration of the coordinate partitioning function produces the solution to the dynamics problem (31). It remains to compute the missing ingredients of the constant terms for the systems (60) and (61). _ partial dierentiation of the system V with Since U does not depend on q, respect to r_ yields the linear equations " # " # 0 Uq oq_ ; 64 B or_ I where we have made use of the extended product rule discussed in Sec_ r_ oq=or in the system of equations (64) form tion 5.1.2. The unknowns oq=o an n p rectangular array. Eq. (64) can be solved in a straightforward manner for the p columns of corresponding to the p columns of constant unknowns 0 terms in the array . The matrix oq=or provides the missing data for the I constant terms in Eq. (60). Partial dierentiation of the system V with respect to q yields the equations "
Uq B
#
2 3 oUq oUt _ q ÿ oq_ 4 oq oq 5: oq 0
65
As done in a previous procedure, multiplying from the right by BT transforms Eq. (65) into 2 3 " # oU t Uq oq_ _ T oU ÿ oqq qB BT oq 5; 4 66 B or 0 _ where the unknowns oq=or form an n p rectangular array which provides the missing data for the constant terms in system (61). The proof of Theorem 1 is complete. _ r_ oq=or. Eq. (64) have p columns of unknowns in the n p matrix oq=o 0 The solutions correspond to the p columns in the constant term, , the rightI hand side of the equation, which is also an n p matrix. Similarly, there are p columns of unknowns in the n ´ p system (66) and in each of the n m p systems (60) and (61). In all cases, each column of the unknowns is the solution to the system with the corresponding column of constant terms from the right-
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
27
hand side. However, solving the latter three systems (66), (60), and (61) will require a little more work than needed for Eq. (64) because the constant terms are more complicated. The solutions to the two pairs of systems of linear equations (64) with (60) and (66) with (61) provide the values of the entries in the top half of the Jacobian matrix (46). There are two steps, in the terminology used in the ®eld of numerical linear algebra, to the solution procedure for such systems of equations: 1. Decomposition (or factorization) of the matrix of coecients which is done just once for a given array of values and 2. Solving (the forward reduction and backward substitution calculations) for the unknowns that correspond to each column array of the constant terms. Although the ®rst step is generally quite computationally intensive, much more so than the second, the ecacy of the implicit method is quite favorable for integrating the coordinate partitioning function. The four systems of equations listed above share the two coecient matrices (44) and (45) which were already required for the evaluation of the function as shown in the schematic diagram at the beginning of subsection 5.1. Hence, the ®rst step of the solution process, the decompositions of the matrices, need not be repeated for the current application. Even though each of the four systems has to be solved a number of times, the additional work is not large compared to the function evaluation. 5.3.3. Expressions for the constant terms The basic elements of the constant terms for Eqs. (66), (60) and (61) must already exist for the evaluation of the coordinate partitioning function which consists of the sequential solutions to the non linear and linear systems D, V, and A/k in the diagram at the beginning of subsection 5.1. It remains to construct the formulas for the constant terms from the basic components. While the complexity of the constant terms appears daunting at ®rst, a number of simpli®cations apply to the constraint equations associated with the typical mechanical system (31). In particular, the standard mechanical joints have no time dependency in which case the corresponding derivatives of U with respect to t are equal to zero. For the sake of completeness, however, we conclude Section 5 in the following paragraphs with the fully general expressions for all of the components of the constant terms. For Eq. (66), we need two formulas: (1) The matrix oUq =oqq_ is given by Eq. (43) derived in Section 5.1.2. (2) The partial time-derivative of U is
28
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
2
oU1 ot
3
7 6 6 oU2 7 6 ot 7 7 6 Ut 6 7; 6 .. 7 6 . 7 5 4
67
oUm ot
where Ui is the ith element in the column matrix U of expressions for the constraint functions. Hence, the i; j component of oUt =oq is
oUt oq
i;j
o2 Ui oqj ot
68
for i 1; 2; . . . ; m and j 1; 2; . . . ; n. Of course, if the ith component of U does not depend on t, then oUi =ot 0. Indeed, we have Ut 0 in many applications. _ T and oUt =oqBT should not be computed as The products oUq =oqqB matrix-multiplies. Rather, multiplication from the right by BT amounts to the selection of the p columns corresponding to the p independent variables from the m ´ n matrices oUq =oqq_ and oUt =oq. There are three expressions involving the constraints U in the constant terms for Eq. (60): _ as already noted above in item 1 for Eq. (66), was 1. The expression oUq =oqq, derived in Section 5.1.2. It is is given by Eq. (43). 2. The m n Jacobian matrix Uq consists of the partial derivatives of U with respect to the displacement variables q. Thus, as stated in Eq. (41), the i; j component of Uq is oUi =oqj where Ui is the ith element in the column matrix U of the expressions for the constraints. Since the matrix U_ q is the term by term time-derivative of Uq , the i; j component of this m n array is oU_ i ; U_ q i;j oqj
69
where it is assumed that the order of the time-derivatives o=ot and d=dt is interchangeable. The matrix-product U_ q oq=or is computed by multi_ r_ to plying the matrix de®ned by Eq. (69) times the solution oq=or oq=o Eq. (64). If the ith component of U does not depend on t, then the entire ith row of U_ q is equal to zero. That is, oU_ i 0 oqj
70
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
29
for j 1; 2; . . . ; n. However, if the ith constraint is time-dependent, then oU_ i o2 Ui o2 Ui o2 U i o2 Ui q_ 1 q_ 2 . . . q_ n oqj oq1 oqj oq2 oqj oqn oqj otoqj
71
for each j 1; 2; . . . ; n. That is, the dot indicates the total time derivative of the expression.3. The m by n matrix oUt =oq is de®ned above by Eq. (68) in item 2 for Eq. (66).As described above at the end of the paragraph on the constant terms for Eq. (66), multiplication from the right by BT constitutes the selection of p columns from the various m n matrices. To get all of the constant terms involving U in Eq. (61), we have to develop ®ve more formulas: 1. The construction of an expression analogous to Eqs. (42) and (43) yields oUTq k oq
! i;j
n X o2 Uk kk oqj oqi k1
72
for the i; j component of the n n square matrix oUTq =oqk. q is de®ned by Eq. (42) which was derived in Sec2. The m n matrix oUq =oq tion 5.1.2. 3. As in item 1, we need to construct an expression analogous to Eqs. (42) and _ Thus, using the de®nition (69) of the matrix U_ q , (43) for the term oU_ q =oqq. we obtain oU_ q q_ oq
! i;j
n X o2 U_ i q_ oqj oqk k k1
73
for i 1; 2; . . . ; m and j 1; 2; . . . ; n. Once again, if the ith component of U _q oU 0 for j 1; 2; . . . ; n. Otherwise, we is independent of time, then oq q_ i;j
have to substitute the rather messy expression n oU_ i X o2 Ui o2 Ui q_ l oqk oql oqk otoqk l1
74
into Eq. (73). _ is computed by multiplying the matrix de4. The matrix-product U_ q oq=or ®ned by Eq. (69) times the solution to Eq. (66). 5. The column matrix Ut is de®ned by Eq. (67). The time-derivative of Ut , still a column matrix with m components, is given by
30
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
2
_1 oU ot
3
7 6 6 oU_ 2 7 7 6 ot 7 U_ t 6 6 . 7 6 . 7 4 . 5 _m oU ot
2 o2 U
1
oq1 ot
q_ 1
6 o2 U 6 2 q_ 1 6 oq ot 6 1 6 4 o2 Um q_ oq1 ot 1
o2 U1 q_ oq2 ot 2
...
... .. . ...
o2 U2 q_ oqn ot n
o2 U1 ot2 o2 U2 ot2
o2 Um q_ oqn ot n
o2 Um ot2
o2 U2 q_ oq2 ot 2
o2 Um q_ oq2 ot 2
o2 U1 q_ oqn ot n
3 7 7 7 7; 7 5
75
where it is assumed that the order of the time-derivatives o=ot and d=dt is interchangeable. Of course, the components of U_ t corresponding to the timeindependent entries in U are equal to zero. As for the previous two cases, multiplication from the right by BT constitutes the selection of p columns from the n n matrix in item 1 and p columns from the m ´ n matrices in the other items. Finally, the partial derivatives of the applied forces oF =oq_ and oF =oq are de®ned in a manner which is analogous to Eq. (38) in Section 5.1.1. Multiplication from the right by BT constitutes the selection of p columns from these two n n matrices. 6. Future developments If the forces in the system of DAE's (31) are acceleration-dependent, which occurs, for example, when a mechanism is modeled with resistance to non zero accelerations, then the equations take the form _ q; t; q; q; M q ÿ UTq k F U q; t 0;
76
where some components of the forces F are functions of the second derivatives q of the coordinates. Most of the analysis in Section 5 still applies to the problem (76). The systems D and V in the coordinate partitioning function on page 13 remain the same. However, the system A/k becomes quite dierent. It may even be non linear. Another kind of complication develops when the system of DAE's (31) takes the form _ q; t; M q ÿ UTq k F q; _ q; t 0; U q;
77
J.F. McGrath, R. Rampalli / Appl. Math. Comput. 111 (2000) 7±31
31
where the algebraic equation depend on the velocities, a condition which characterizes non holonomic constraints. The system of equations D in the coordinate partitioning function is no longer independent of the system V. Hence, the equation D and V have to be solved simultaneously, for which the analysis presented in Section 5 will have to be extensively reworked. Finally, both acceleration-dependent forces and non holonomic constraints can be incorporated in the same mechanical system, which also means the previous results no longer apply. The development of the details for the formation of the Jacobian matrices corresponding to problems (76) and (77) and to problems with both complications awaits a future report. References [1] L.F. Shampine, H.A. Watts, DEPAC-Design of a user- oriented package of ODE solvers, Sand79-2374, Sandia National Laboratories, 1979. [2] L.F. Shampine, Numerical Solutions of Ordinary Dierential Equations, Chapman & Hall, New York, 1994. [3] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Dierential-Algebraic Equations, North Holland, New York, 1989. [4] C.W. Gear, Numerical Initial Value Problems in Ordinary Dierential Equations, PrenticeHall, Englewood Clis, NJ, 1971. [5] R.A. Wehage, E.J. Haug, Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic system, J. Mechanical Design 104 (1982) 247±255.