OPTIMAL PRECONDITIONERS FOR THE
SOLUTION OF CONSTRAINED MECHANICAL
SYSTEMS IN INDEX THREE FORM
POLItecnico
di MIlano
Carlo L. Bottasso, Lorenzo Trainelli
Politecnico di Milano, Italy
Daniel Dopico
University of La Coruña, Spain
ACMD 2006
The University of Tokyo
Komaba, Tokyo, August 1-4, 2006
Optimal Scaling of Index Three DAEs
Outline
• Introduction and background, with specific reference to
Newmark’s method;
• Conditioning and sensitivity to perturbations of index 3 DAEs;
• Optimal preconditioning by scaling;
• Conclusions:
Index 3 DAEs can be made as easy to integrate as well
behaved ODEs!
POLITECNICO di MILANO
Optimal Scaling of Index Three DAEs
Index 3 Multibody DAEs
Constrained mechanical (multibody) systems in index 3 form:
M v 0 = f (u; v; t) + G(u; t) ¸;
u0 = v;
0 = ©(u; t);
where
• Displacements:
• Velocities:
• Lagrange multipliers:
• Constraint Jacobian:
POLITECNICO di MILANO
u
v
¸
G := ©;T
u
Optimal Scaling of Index Three DAEs
Newmark’s Method
Newmark’s method in three-field form (Newmark 1959):
µ
¶
1
1
M (vn+1 ¡ vn ) = fn+1 + Gn+1 ¸n+1 ¡ 1 ¡
M an ;
°h
°
µ
µ
¶¶
µ
¶
¯
¯
¡ h2 1 ¡ 2 ¯ a ;
un+1 ¡ un = h
vn+1 + 1 ¡ vn
n
°
°
2
°
0 = ©n+1 ;
and similarly for HHT, generalized-α, etc.
Remark: the analysis is easily extended to other integration methods,
for example BDF-type integrators (Bottasso et al. 2006).
POLITECNICO di MILANO
Solving with Newton’s method at each time step:
J q = ¡b;
{
{
Multipliers
6
6
J =6 I
4
Displacements
0
0
0 right hand side is
and the
7
7
7
5
Equilibrium eqs.
Kinematic eqs.
{
Velocities
X
¡G
{
GT
1
U
°h
¡¯ h I
°
0
{
where
is
2 the iteration matrix 3
{
Optimal Scaling of Index Three DAEs
Linearization of Newmark’s Method
Constraint eqs.
X := ¡ (f;u )
n+1
Y := ¡ (f;v )
n+1
³
´
¡ (G ¸)
;u
;
U := M + ° h Y :
³
´
¡ v ) ¡ (f
¡ 1 Ma
1 M (v
+
G
¸
)
+
1
n+1
n+1³
n+1 ´´
n+1
B °h
³n
³ ° ´ n
b=B
@ un+1 ¡ un ¡ h ¯ vn+1 + 1 ¡ ¯ vn + h2 1 ¡ 2 ¯ an
°
°
2
°
©n+1
POLITECNICO di MILANO
n+1
1
C
C:
A
;
Optimal Scaling of Index Three DAEs
Perturbation Analysis
Model effects
arithmetics by means of the small
" of finite precision
lim b(h; 0) = 0
parameter :
!
¯
@b ¯¯
b(h; ") = b(h; 0) + "
@" ¯
h
0
+ O("2 );
(h;0)
¯
@J ¯¯
J(h; ") = J(h; 0) + "
+ O("2 );
¯
@"
¯ (h;0)
@q ¯¯
q(h; ") = q(h; 0) + "
+ O("2 ):
¯
@"
(h;0)
lim q(h; 0) = 0
h!0
Insert into ¯
J q = ¡b
and
¯ °neglect higher °order
° terms to °get
¯
¯ °
°
¯ lim q (h; ")¯ · ° lim J ¡1 (h; 0)°
¯ ! i
¯ ° !
°
h
0
POLITECNICO di MILANO
h
0
1
°
°
° lim b(h; ")°
° !
°
h
0
1
:
Optimal Scaling of Index Three DAEs
Perturbation Analysis
Interpretation of the main result
¯
¯ °
°
¯
¯ °
°
¯ lim q (h; ")¯ · ° lim J ¡1 (h; 0)°
¯ ! i
¯ ° !
°
h
J ¡1
0
h
b
0
1
°
°
°
°
° lim b(h; ")°
° !
°
h
0
:
1
h
If
and/or
depend on negative powers of , accuracy will
be very poor for small time steps (this will spoil the order of
convergence of the integration scheme).
POLITECNICO di MILANO
Optimal Scaling of Index Three DAEs
Perturbation Analysis of Newmark
2
3 form:
Apply perturbation
analysis¡to three-field
O(h0 ) O(h 1 ) O(h0 )
5;
0
lim J = 4 O(h0 ) O(h1 )
h!0
O(h0 )
0
0
2
3
O(h2 ) O(h0 )
O(h0 )
lim J ¡1 = 4 O(h1 ) O(h¡1 ) O(h¡1 ) 5 :
h!0
O(h0 ) O(h¡2 ) O(h¡2 )
8
9
¡
< O(h 1 ) =
O(h0 )
lim b =
:
:
;
!
h 0
O(h0 )
¯ Accuracy of¯ the solution:
¯
¯
¯
¯
¯
¯
¯ lim ¢u
¯ · O(h0 ); ¯ lim ¢v
¯ · O(h¡1 );
¯ !
¯
¯
n+1
n+1 ¯
!
h
0
h
0
¯
¯
¯
¯
¯ lim ¢¸
¯ · O(h¡2 ):
¯ !
n+1 ¯
h
0
Remark: severe loss of accuracy (ill conditioning) for small
POLITECNICO di MILANO
h
!
Optimal Scaling of Index Three DAEs
Perturbation Analysis of Newmark
Perturbation analysis for other possible forms of Newmark’s
method:
¢u
¢v
¢a
¢¸
(u; v; a; ¸)
O(h0 )
O(h¡1 )
O(h¡2 )
O(h¡2 )
POLITECNICO di MILANO
(u; v; ¸)
O(h0 )
O(h¡1 )
O(h¡2 )¤
O(h¡2 )
(u; ¸)
O(h0 )
O(h¡1 )¤
O(h¡2 )¤
O(h¡2 )
(v; ¸)
O(h0 )¤
O(h¡1 )
O(h¡2 )¤
O(h¡2 )
(a; ¸)
O(h0 )¤
O(h¡1 )¤
O(h¡2 )
O(h¡2 )
Optimal Scaling of Index Three DAEs
Index 3 Multibody DAEs
Solutions proposed in the literature:
• Index two form (velocity level constraints), but drift and need for
stabilization;
• Position and velocity level constraints (GGL of Gear et al. 1985,
Embedded Projection of Borri et al. 2004), but increased cost and
complexity;
• Scaling, simple and straightforward to implement (related work:
Petzold Lötstedt 1986, Cardona Geradin 1994, Negrut 2005).
POLITECNICO di MILANO
Optimal Scaling of Index Three DAEs
Left and Right Preconditioning
Scale equations and unknowns of Newton’s problem:
¹
J¹ q¹ = ¡b;
with
where
J¹ := DL J DR ;
DL
DR
q¹ := D¡1 q;
R
¹ := D b;
b
L
= left preconditioner, scales the equations;
= right preconditioner, scales the unknowns.
Remark: scale back unknowns once at convergence of Newton’s
iterations for the current time step.
POLITECNICO di MILANO
Optimal Scaling of Index Three DAEs
Optimal Preconditioning
Left scaling:
2
¯ h2 I
DL = 4 0
0
2
I
6 0
DR = 6
4
0
Right scaling: 6
0
°
I
¯h
3
0 0
h2
I 0 5
Scale equilibrium eqs. by
0 I
0
0
1
I
¯ h2
0
3
7
7 Scale velocities by h
7
5
Scale Lagrange multipliers by
¯ Accuracy of¯ the solution:
¯
¯
¯
¯
¯
¯
¯ lim ¢u
¯ · O(h0 ); ¯ lim ¢v
¯ · O(h0 );
¯ !
¯ !
n+1 ¯
n+1 ¯
h
0
h
0
¯
¯
¯
¯
¯ lim ¢¸
¯ · O(h0 ):
¯ !
n+1 ¯
h
0
Remark: accuracy (conditioning) independent of time step size, as for
well behaved ODEs! Same results for all other forms of the scheme.
POLITECNICO di MILANO
h2
Optimal Scaling of Index Three DAEs
Numerical Example
Andrews’ squeezing mechanism:
At each Newton iteration, solve ¡
J j q j = bj
Stop at first non-decreasing
kq j+1 kNewton
¸ kq j kcorrection, i.e.:
A measure of the tightest achievable convergence of Netwon method.
POLITECNICO di MILANO
Optimal Scaling of Index Three DAEs
Numerical Example
POLITECNICO di MILANO
Optimal Scaling of Index Three DAEs
Numerical Example
POLITECNICO di MILANO
Optimal Scaling of Index Three DAEs
Conclusions
• DAEs can be made as easy to solve numerically as well behaved ODEs;
• Recipe for Newmark’s scheme:
- Use any form of the algorithm (two field forms are more
computationally convenient);
- Use left and right scaling of primary unknowns;
- Recover scaled quantities at convergence.
• Simple to implement in existing codes.
POLITECNICO di MILANO