Microsoft Word - numero_29_art_12 A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 128 Focussed on: Computational Mechanics and Mechanics of Materials in Italy Geometric numerical integrators based on the Magnus expansion in bifurcation problems for non-linear elastic solids A. Castellano, P. Foti, A. Fraddosio, S. Marzano, M. D. Piccioni Politecnico di Bari – Dipartimento di Scienze dell’Ingegneria Civile e dell’Architettura anna.castellano@poliba.it, pilade.foti@poliba.it, a.fraddosio@poliba.it, s.marzano@poliba.it, m.d.piccioni@poliba.it ABSTRACT. We illustrate a procedure based on the Magnus expansion for studying mechanical problems which lead to non-autonomous systems of linear ODE’s. The effectiveness of the Magnus method is enlighten by the analysis of a bifurcation problem in the framework of three-dimensional non-linear elasticity. In particular, for an isotropic compressible elastic tube subject to an azimuthal shear primary deformation we study the possibility of axially periodic twist-like bifurcations. The approximate matricant of the resulting differential problem and the first singular value of the bifurcating load corresponding to a non-trivial bifurcation are determined by employing a simplified version of the Magnus method, characterized by a truncation of the Magnus series after the second term. KEYWORDS. Nonlinear elasticity; Bifurcation; Geometric numerical integrators; Magnus expansion. INTRODUCTION n the framework of three-dimensional non-linear elasticity, a number of “small on large” bifurcation problems lead to the analysis of a system of linear second order ODE’s with varying coefficients. It is a common practice to reduce the system of linear second order ODE’s to a simpler non-autonomous first order linear ODE system. Nevertheless, the resulting differential system may be somewhat complex and only numerically tractable; thus, it is crucial to adopt an adequate strategy for obtaining an accurate numerical solution for determining the critical load and the corresponding bifurcation deformation field. For example, approximate expressions of the matricant of non-autonomous linear ODE systems, based either on the multiplicative integral of Volterra or on the truncated Peano expansion, have been recently proposed for bifurcation problems analyzed by means of the Stroh approach within the so-called sextic formalism [1-3]). In Section 2, we discuss an alternative numerical method – based on the Magnus expansion [4] – which furnishes an approximate exponential representations of the matricant of first order linear ODE systems. The Magnus expansion belongs to the so-called geometric numerical integrators, based on Lie-Group methods [5 -7] for appropriate references). To our knowledge, the applications of this method are not widespread in continuum mechanics, but it has been shown in [8] that for problems involving singularities, bifurcations and wave propagation the geometric numerical integrators may be useful and accurate. In particular, the Magnus expansion may be very efficient for a number of bifurcation analyses, since it leads to the determination of approximate solutions that preserve at any order of approximation the same qualitative properties of the exact (but unknown) solution; moreover, this method exhibits an improved accuracy with respect to other frequently used numerical schemes. As an application of the Magnus method, in An application of the Magnus method in solid mechanics: periodic twist-like deformations bifurcating from the azimuthal shear of a circular tube we investigate the possibility for a compressible hollow cylinder subject to I A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 129 azimuthal shear to support axially periodic twist-like bifurcation modes similar to the Taylor-Couette axially periodic cellular patterns observed when the flow of a viscous fluid confined in the gap between two rotating concentric cylinders becomes unstable. Here, the hollow cylinder is made of a Levinson-Burgess isotropic compressible elastic material. By restricting our attention to a class of incremental displacements characterized by three unknown scalar functions of the radial coordinate and having an axial periodic structure, we then study the related incremental boundary-value problem. This leads to a somewhat complex non-autonomous homogeneous system of six first order linear ODE’s with homogeneous boundary conditions, whose analysis is performed by a numerical approach. We show that a procedure based on the Magnus method outlined in An overview of geometric numerical integrators with special reference to the Magnus expansion reveals very effective either for approximating the matricant of the resulting differential problem or for determining the first singular value of the bifurcating load corresponding to a non-trivial twist-like solution. AN OVERVIEW OF GEOMETRIC NUMERICAL INTEGRATORS WITH SPECIAL REFERENCE TO THE MAGNUS EXPANSION number of problems in scientific and engineering areas, ranging from Quantum Mechanics, Optics, Electromagnetism and Molecular Physics to Control Theory and Bifurcation, require the resolution of systems of differential equations with varying coefficients. Besides some very special cases, for which a fundamental matrix solution is explicitly available, the analysis of the problem usually requires the construction of an approximate representation of the solution of the non-autonomous differential system. In the last decade, this issue has attracted many scientists and has also represented a field of intense research activity; in particular, recent studies have shown that in many cases the so-called geometric numerical integration methods may offer better approximation schemes compared to perturbative procedures or to standard numerical integrators. The Magnus method belongs to the class of geometric numerical integrators, based on Lie-Group methods (see [6, 7] for appropriate references), which, although not yet diffused in continuum mechanics, have recently found many applications in areas of molecular physics [9], electromagnetism [10] and celestial mechanics [11]. In this Section, we first give a brief overview of geometric numerical integrators and then we deal with the Magnus method as a particular case. Geometric numerical integration concerns the development of numerical approximations to the solution of an ODE system with the main goal of preserving at any order of approximation the qualitative properties of the exact (but unknown) solution and thereby of exhibiting an improved accuracy with respect to other numerical methods. Although such methods apply also for nonlinear matrix differential equations, for the purposes of bifurcation problems here we consider geometric numerical integrators based on the Magnus expansion for the special case of homogeneous non- autonomous first order linear ODE systems of the form         6x6R R R , R 0, 0   Y A Y Y I (1) where the unknown 6x6 matrix  RY , usually referred as the matricant of the Cauchy initial-value problem (1), is the 6- dimensional identity matrix when the real independent variable R is equal to zero. Generally speaking, these new numerical approaches involve Lie-Group methods, which are based on Lie groups and their associated Lie algebras. We recall that a Lie group corresponds to a differential manifold which is endowed with a group structure; the corresponding Lie algebra, defined as the tangent space to the Lie group at the identity, is a linear space endowed with the Lie bracket commutation operation and a natural mapping. Since the unknown matricant of a first order ODE system like (1) must be an invertible 6x6 matrix with variable entries, the setting for finding solutions of (1) is a 6x6 matrix Lie group. In this case, the above general definition of a Lie algebra implies that the 6x6 matrix  RA governing (1) is an element of a matrix Lie algebra. Moreover, in this setting the associated matrix Lie algebra is endowed with the matrix commutation law  , := A B AB BA (2) where  , A B is commonly called the matrix commutator. Finally, for these special matrix Lie groups the matrix exponential function plays the role of the natural mapping. We emphasize that in the theory of Lie groups the role of the exponential mapping is crucial. Indeed, within the context of differential geometry, the matrix exponential maps an element belonging to a matrix Lie algebra into an invertible matrix belonging to its corresponding Lie group; in other words, the exponential mapping allows to recapture the local A A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 130 group structure from the corresponding Lie algebra. Thus, for a first order ODE system like (1), it is natural to attempt to construct the solution as a matrix exponential with exponent given by a matrix which belongs to the associated Lie algebra. This peculiarity of the exponential mapping probably inspired the idea developed by Magnus [4], who proposed to express the exponent by a matrix series expansion made of infinite terms, each of them belonging to the Lie algebra corresponding to the Lie group of the unknown matricant. In detail, with reference to the differential problem defined in (1) the Magnus method assumes the matricant to be of exponential form with the exponent given by the following series expansion, called the Magnus expansion:            6x6 k k=1 R exp R 0 , 0 , R R       Y Y O    (3) Here we do not report the analytical procedures for determining the terms of the series expansion in (3). We address the reader to the seminal paper [5] where, in view of (1) and (3), and according to the notions of the first derivative of a matrix exponential map and of its inverse (see formulae (33) and (38) in [5]), it is clearly showed how it is possible to construct all the terms of the Magnus series. For our purposes, we report here only the first two terms of the series:           1R R 1 1 1 2 1 1 2 2 0 0 0 1 R d , R d , d 2             A A A  (4) where [. , .] is the matrix commutator (2) and  A is the governing matrix of (1). Notice that each of the remaining terms (here not reported) contains nested matrix commutators of the form (2) involving the matrix operator  RA . It is worth to note that only if the commutative condition        1 2 2 1 1 2R R R R for each R , RA A A A (5) holds, then, in view of (3), (1) and (4), the matricant of (1) takes the exponential form     R 0 R exp d          Y A (6) which is analogous to the classical exponential solution for a scalar linear ODE. Clearly, (6) yields the solution of (1) if the entries of the matrix A governing (1) are independent of R. Obviously, (5) holds only in special cases because the matrices in general do not commute (this is closely related to the structure of the commutation law (2) which characterizes matrix Lie algebras). A crucial aspect of the Magnus method related to its pertinence when seeking approximate solutions of differential problems like (1) emerges from the analysis of (3) and (4). By the above discussion, we recall that the matrix  RA in (1) belongs to a matrix Lie algebra for all R; thus, in view of (4) we see that any term of the Magnus expansion belongs to the same Lie algebra. This implies that any truncation of the Magnus series will also belong to the same Lie algebra and therefore the exponential map of any truncation will necessarily stay in the corresponding Lie group. This is basically the main reason why an approximated solution obtained by truncating the Magnus expansion at any order preserves the same qualitative features of the exact solution. Another major issue concerns the conditions on the matrix  RA which guarantee the convergence of the Magnus series. This problem has been studied in [12], where it is shown that a sufficient condition for the local convergence of the Magnus series in a certain interval  20, R is given by :   2R 2 0 R dR   A (7) where the integrand function in (7) is the spectral norm of  RA , i.e., the square root of the maximum eigenvalue of the product    TR RA A . Moreover, the fulfilment of the inequality (7) determines the range of the independent variable R for which it is possible to write the solution in the form (3). However, in many applicative cases it happens that the convergence condition (7) does not hold for the whole integration interval  20, R . In these cases, as usual for numerical A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 131 integration methods, it is common dividing the interval  20, R into N subintervals  i+1 iR , R such that the Magnus series converges in each of them. Given such a subinterval  i+1 iR , R , in view of (3) we may write:       i+1 i+1 iR exp R RY Y (8) wherein  i+1R is expressed by the infinite Magnus series. Consequently, the matricant of (1) evaluated on the whole domain of interest  20, R takes the form      N 2 i+1 i=1 R exp R        Y  (9) As it usually occurs when constructing approximate solutions, a crucial issue is that of choosing an appropriate order of truncation for  i+1R which guarantees a suitable level of accuracy. To this aim, we recall the efficient procedure suggested in [7], which may be summarized in the following three steps: first, the series  i+1R is truncated at a suitable order in each subinterval wherein the convergence is guaranteed. Then, the integrals appearing in each truncated series are approximated by appropriate quadrature rules. Finally, the exponential of the matrix  i+1R is evaluated. The qualified literature on this issue contains a number of proposals which correlate the level of accuracy of the method to the order of truncation of the Magnus series. We recall the paper [6] for an overview on this problem. Finally, the accurate calculation of the exponential of a matrix is still an open problem, as one may infer from the paper [13]. In the spirit of geometric numerical integration methods, the consequence of an inaccurate computation of the exponential matrix due, for example, by the replacement of the exponential map with one of its analytic (for example rational) approximants may lead to a solution not belonging to the desired Lie group, thus rendering the main advantages of geometric numerical integration ineffective. Here, we report two references on this issue, namely the splitting method proposed in [14] and the scaling and squaring with a Padè approximation method summarized in [13]. Notice that the command MatrixExp of the widely used software Mathematica employs the Padè approximation method for the calculation of the exponential matrix. We conclude by observing that the non-trivial evaluation of the exponential matrix may be circumvented when the unknown matrix of a differential problem belongs to special Lie groups. One of these cases involves a group commonly occurring in mechanical problems, namely the orthogonal group, that is the group of 3x3 matrices Q such that T 3x3 .QQ I Unlike many other Lie-group solvers, they do not require the evaluation of matrix exponentials because, as outlined in [15], the back translation from the Lie algebra to the corresponding Lie group may be obtained by using the Cayley transform. AN APPLICATION OF THE MAGNUS METHOD IN SOLID MECHANICS: PERIODIC TWIST-LIKE DEFORMATIONS BIFURCATING FROM THE AZIMUTHAL SHEAR OF A CIRCULAR TUBE n the present section, we introduce a bifurcation problem set within the context of the nonlinear theory of elastic solids. In particular, the last part of the following analysis contains an explicit application of the Magnus method based on the procedure outlined in the previous section. The fundamental azimuthal shear equilibrium deformation We consider a homogeneous, isotropic, compressible, elastic tube which occupies the region   1 2R, , Z 0 < R R R , 0 2 , 0 Z H          (10) in its natural reference configuration.  R, , Z denote the cylindrical coordinate of a point X in a cylindrical coordinate system with orthonormal basis  R Z, , e e e . The boundary of  is divided into two disjoint complementary parts:      , 1 1 2 2R, , Z R = R or R = R R, , Z Z = 0 or Z = H           (11) I A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 132 A deformation        : = r, , z R, , Z    f X x f X f  of  is assumed to be a smooth function with gradient    : F X f X . Since the tube is elastic and isotropic, the general form of the strain energy function is given by    W W I , II , IIIF B B B , where         22 2 22 T1 1I tr , II tr tr tr , III det det 2 2               B F F B B F F FF B FB B B (12) are the orthogonal principal invariants of T : B F F ; thus, the Piola stress takes the form       T1 2 3 DW W W I W III 2 2 2     S F F F F BF FB B (13) with 1 2 3 W W W W : , W : , W : I II III         B B B (14) We assume that the inner cylinder is kept fixed, whereas the outer is subject to a uniform angular displacement 0  around its axis. On the bases of  only tangential displacements are admitted. This leads to the following mixed boundary-value problem:  Div = in S 0F  (15)            1 R 2 2 1 Z 1 2 2 1 2, 0 at R = R = R cos 1 at R = R 0 at R = R = on R sin at R = R = 0 at R = R R                f X X e f X X e f X X e  (16)       2 = 0, = on      f X X e S e e 0F  (17) For the equilibrium problem (15)-(17), we consider the possibility of an azimuthal shear deformation f defined by  r R, R , z Z       (18) where  , the angular displacement field, is assumed to be a smooth function satisfying    1 2R 0, R 0      1on   (19) Consequently, in view of (18) we have     2 2R r rR R , R R                   I e I F e B e e e e e e (20) where the prime denotes differentiation with respect to R and  r z, , e ee is the deformed cylindrical orthonormal basis at  r, , z . Notice that the azimuthal shear is an isochoric deformation whose principal invariants (12) are given by 2 2 , I II 3 R II 1I     B B B (21) Moreover, (18) trivially satisfies the displacement boundary conditions (16)-(17)1 and, in view of (13) and (20), it is easily seen that also the traction boundary condition (17)2 holds. It remains to study the equilibrium field Eq. (15). The explicitly evaluation of (15) by means of (13) and (20) yields an over-determined system of two ODE for the single unknown function  R . In particular, here we consider the class C2 Levinson-Burgess strain energy function A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 133            21 2 1 2= 1 4 1I 3 1 J 3 2 1 2 III 1 2 III 1 2 1 2 W                       F B B B B (22) where and are the referential shear modulus and referential Poisson’s ratio, respectively, and  0, 1 is a dimensionless material parameter. For the constitutive class (22), the Piola stress takes the form        1 2 1 2-1 -T -T4 1 1 1 2 2 III 1 III 1 2                            B BS F F B F F (23) By following a result contained in [16], it is possible to generate classes of strain energy functions such that the azimuthal shear satisfies both the two equilibrium ODE; in particular, for the present case of a Levinson-Burgess material, deformations of the form (18) are allowed only if the following condition on the material parameters holds: 3 4   (24) Finally, in view of (23), (19)-(21) and (24), after a non-trivial calculation we obtain from (15) a single ODE, whose solution is      2 2 2 2 2 22 1 1 2R R R R R R 1 R,           (25) Periodic twist-like bifurcations We now consider the possibility of deformations which bifurcate from the primary equilibrium pure circular shear (18) as the loading parameter (angle)  increases from 0 (natural state). A condition (only necessary) for the occurrence of a local branch of bifurcating solutions at a fixed  is the existence of a non-trivial solution of the homogeneous linearized equilibrium problem which corresponds to an adjacent equilibrium state. Such solutions superposed to the primary azimuthal shear must satisfy the adjacent equilibrium Eq. (26) and the incremental boundary conditions (27)-(28), obtained by linearizing (15)-(17) around the fundamental deformation  x f X , which is conveniently assumed as the independent variable:     gradr, div in  0 uF  (26) 1 on u 0  (27)    gradr, z 2z z= 0, = on  u e e e 0 uF  (28) In particular,   3: u x  represents an incremental displacement field, “div” and “grad” are the divergence and gradient operators with respect to x, respectively, and  is the fourth-order instantaneous elasticity tensor, given by:           r, r, r, r, 2 T: D , W Lin          HH HF F F F (29) For the constitutive class (22), it follows from (29), (20)-(21) and (24) that                   r, T T1 1 T grad 1 grad grad grad 4 3 2 1 grad grad grad 4 1 2 2                        u B B u u u B I u I u     uF (30) where B is given by (20)2. We now restrict our analysis to the class           3 r 1 2 z 3 1 2 : : u v r cos z, u v r cos z, u v r sin z, n , n 0, 1, 2, .., r r H                       u u u  0   (31) A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 134 of periodic bifurcating displacements which are reminiscent of the instability patterns observed in the classical Taylor- Couette shear flow of viscous fluids. A more detailed discussions on such a choice, together with an exhaustive analysis of problems closely related to the one here presented, are developed in [17], [18]; in particular, in [17] one may also find stability issues based on [19] and partly on [20]. This allows us to skip here the description of certain analytical developments, for which we refer to the works above mentioned. Notice that (31) models the occurrence of an axially periodic cellular pattern in the gap between the inner and outer cylinders (n represents the number of possibly forming cells in the axial direction); in particular, inside each of the n possible forming cells, (31) describes a twist-like displacement perpendicular to the e -direction of primary annular shear, so that the vector lines of the incremental displacement u are similar to the streamlines of the twisting Taylor-like effects for fluids. We easily check from (31) that the displacement incremental boundary condition (28)1 trivially holds, whereas (27) holds provided the smooth scalar functions  1v r ,  2v r ,  3v r for r in  1 2R , R are such that                    1 1 1 1 2 2 21 2 3 2 3 21r : r r r rv , v , v v , v , v: r r r 0 0, , 0    v v 0 (32) For what concerns the traction boundary condition (28)2 and the field Eq. (26), we first observe by (31) that     1 1 1 r r 2 r 2 r 1 3 z z 1 r z 2 z 3 z r grad v v r v v r v v cos z sin z v v                             u e e e e e e e e e e e e e e e e (33) and by (20)2 that z z B e e ; thus, in view of (30)-(31), condition (28)2 immediately holds. Finally, by evaluating the divergence of (30) with the aid of (20)2 and (33), after some non-trivial calculation and also by employing separation of variables allowed by the choice (31) of the displacement field, we reduce the partial differential problem (26)-(28) to the following system of three homogeneous second-order ODE (see also (32)):                   T 1 2 1 2 r, r, r, , r, , r, , r, , r r r r r                           P v P R R v R Q v v v 0 0 (34) where  1 2 3v , : v , vv (we have omitted the dependence on r) and       rrrr rr r rr rrr rrzz rrr r r r rr rzz zrzr zrrz zr z 2 2 rzrz r rz z z r A r A 0 A A r A r, r A r A 0 , r, , A A r A 0 0 r A r A r A 0 1 1 A r A A r A A r r r, ,                                                          P R Q z 2 2 r zrz r r z z r zz 2 zz zzr zzzz 1 1 A r A A r A A r r A A r A                                    (35) are non-constant matrices determined by the components    ijhk h k i j hkij:    e e e e   of the symmetric fourth- order tensor  in the coordinate system (r, , z). Here, we do not report the explicit expression of such components, which may be obtained by (30), (20)2 and (26). We only observe that  r, P and  r, , Q are symmetric, since  is symmetric; moreover, for the constitutive class (22)  r, P is always invertible. The latter property of  r, P is crucial. Indeed, a common practice is that of transforming the set of three linear second order ordinary differential equations like (34) into a system of six linear first order ordinary differential equations, for A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 135 which well-known existence theorems and procedures for constructing the solutions are available in the literature. With this aim, once introduced the 3x3 matrices                 1 1 T : r, r, , r, , : r, r, r, , r, ,                        T P R Q K P P R R (36) we rewrite (34) in the form     1 2 1 2 , r r r r r          v Tv Kv v v 0 (37) Then, by setting 6x1 6x6 6x6 6x6 : , : : , :                                   v O I y A v T K I O O O M N O O I O (38) we easily see that (37) is equivalent to the linear homogeneous first order ODE boundary-value problem       1 2 1 2 6x1 r for r r r r r                  y A y My Ny 0 0 (39) Observe that in the notation above we have dropped the dependence of the basic parameters (the inner and outer radii, the height of the tube, the material moduli and the angular displacement  ), usually assumed as prescribed. We only maintain explicit the dependence on r, in order to emphasize that the linear ODE system is non-autonomous. Now, by following elementary results on systems of linear ODE's, we may write a solution of (39) in the form (see (38)1 and (37)2)           1 1 6x1 r r r r r           y Y y Y v 0 (40) where  1ry is an unknown constant vector and               1 2 1 2 1 1 2 6x66x6 r r r , with det r 0, for r < r < r and r r r                  I OU U Y Y Y O IU U (41) is a particular fundamental matrix solution for (39)1, usually called matricant. Then, substitution of (40) into (39)2 yields      1 2 1r r r             0 MY NY y 0 (42) which, in view of (38)3-4 , (41) and (40), is equivalent to the homogeneous system of three algebraic equations    2 2 1r r U v 0 (43) for which non-trivial solution  1rv 0  are possible if and only if  2 2det r 0U (44) A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 136 Finally, a solution of the primary ODE problem (34) can be written in the following form      2 1r r rv U v (45) which clearly satisfies the boundary conditions (34)2 (see (41) and (43)). The above analysis shows that the possibility of determining non-trivial solutions for the bifurcation differential problem (34) is strictly related to the knowledge of the 3x3 matrix  2 rU , which represents the “north-east” block of the 6x6 matricant (41) of the Cauchy initial-value problem         1 2 1 6x6 r r r for r r r r              Y A Y I O Y O I (46) Since (46) has the same structure of the Cauchy initial-value problem (1), in the following we illustrate a strategy based on the Magnus method for determining the solution of (46). A strategy for determining the bifurcating load through the Magnus method The theoretical preliminaries in An overview of geometric numerical integrators with special reference to the Magnus expansion combined with the bifurcation analysis of Subsections The fundamental azimuthal shear equilibrium deformation - Periodic twist-like bifurcations allow us to describe a strategy for determining the critical value cr of the primary angular displacement, that is the lowest value of the angle  for which the tube may support (during a loading process) a non-trivial twist-like bifurcation of the form (31) from the fundamental azimuthal shear deformation. The Magnus method is now the key ingredient for approximately determining the matricant of (46), whereas (44) plays the role of bifurcation condition. Following the steps outlined in the Section 2, we first fix the radii R1 and R2, the height H, the referential shear modulus  and the referential Poisson’s ratio  (recall that 3 / 4  by (24)). Then, we: - fix the number n of twisting cells in the axial direction (by (31) this is equivalent to fixing  ) starting from n =1; - fix the value of the loading parameter  starting from zero; - compute the matrices (35), which now depend only on r, and then compute  rA in (38) by means of (36); - divide the interval  1 2R , R into N subintervals i+1 ir , r such that the Magnus series converges in each of them (see the convergence condition (7)); - employ for simplicity a fourth-order method by truncating the Magnus series  i+1r after the second term in each subinterval (see (4)) and by approximating the integrals using a Gauss quadrature rule. This allow us to write                2 1 2 1 2 i+1 i i i i h h 3 r r r r , r 2 12        A A A A (47) where    1 2 2 1 i ii i R R1 3 1 3 r r h, r r h , with h = 2 6 2 6 N                   (48) are two Gauss points in correspondence of which  rA and its employed linear approximation agree. - compute the matrix exponential of  i+1r using the function MatrixExp in the software Mathematica; - evaluate the matricant of (46) at the external radius r2 = R2 by (8) and (9); - compute the “north-east” block of the matricant at r2 = R2 , that is the 3x3 matrix  2 2rU ; - compute the determinant of  2 2rU and check if the bifurcation condition (44) holds. Following this procedure, we define ncr as the smallest value of  which satisfies the bifurcation condition (44) corresponding to the fixed value of n. Of course, when one repeats the above procedure by varying the number of possible axial cells, the corresponding critical load ncr varies. We have performed a large number of computations which show that as n increases, the corresponding ncr determined by (44) decreases until a critical value of n above which (44) A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 137 definitively does not hold for any  . In correspondence, a load cr corresponding to crn is determined, and this is the critical value of the circular angle of shear which allows for the occurrence of a non-trivial axial periodic deformation of the form (31). Fig. 1 reports the bifurcation curve for a referential shear modulus 1  MPa, referential Poisson’s ratio 0.25  , 3 / 4  and H/R1=1. It shows the critical load cr as a function of the ratio 2 1R /R between the outer and the inner radius of the tube. We may easily infer, as expected, that twist-like bifurcation may be more easily supported in thin tubes. Figure 1: cr vs. R2/R1, for  = 1,  = 3/4,  = 0.25, H/R1=1. CONCLUSIONS e have studied a bifurcation problem concerning the possibility of periodic twits-like deformations superposed to a primary azimuthal shear state for a Levinson-Burgess isotropic elastic tube. The procedure proposed for determining approximate solutions of the resulting non-autonomous systems of linear ODE’s is based on the Magnus method, which represents an innovative numerical scheme belonging to the more general class of geometric numerical integrators. We have the feeling that these new methods, already employed in many scientific areas, may find interesting applications also in the field of continuum mechanics. In the future we are going to further explore the applicability and the advantages of the Magnus method by studying other bifurcation problems yielding non-autonomous ODE systems like, for example, a Taylor-Couette type bifurcation for a sheared incompressible annular tube modeled by an extended version of the Gent constitutive equation. ACKNOWLEDGMENTS he authors gratefully acknowledge the research projects MIUR-PRIN 2010-2011: “Dinamica, stabilità e controllo di strutture flessibili” and MIUR PON-REC: “MASSIME” – Sistemi di sicurezza meccatronici innovativi (cablati e wireless) per applicazioni ferroviarie, aerospaziali e robotiche. REFERENCES [1] Goriely, A., Vandiver, R., Destrade, M., Nonlinear Euler buckling, Proc. R. Soc. A, 464 (2008) 3003 – 3019. [2] Shuvalov, A. L., A sextic formalism for the three-dimensional elastodynamics of cylindrically anisotropic radially inhomogeneous materials, Proc. R. Soc. Lond., 459 (2003) 1611–1639. [3] Shuvalov, A. L., Poncelet, O., Deschamps, M, General formalism for plane guided waves in transversely inhomogeneous anisotropic plates, Wave Motion, 40 (2004) 413 – 426. W T A. Castellano et alii, Frattura ed Integrità Strutturale, 29 (2014) 128-138; DOI: 10.3221/IGF-ESIS.29.12 138 [4] Magnus, W., On the exponential solution of differential equations for a linear operator, Comm. Pure Appl. Math., VII (1954) 649 – 673. [5] Blanes, S., Casas, F., Oteo, J.A., Ros, J., The Magnus expansion and some of its applications, Physics Reports, 470 (2009) 151 – 238. [6] Iserles, A., Munthe-Kaas, H.Z., Norsett, S.P., Zanna, A., Lie-group methods, Acta Numer., 9 (2000) 215 – 365. [7] Iserles, A., Norsett, S.P., On the solution of linear differential equations in Lie groups, Phil. Trans. R. Soc. A, 357 (1999) 983 – 1019. [8] Budd, C. J., Piggott, M.D., The geometric integration of scale-invariant ordinary and partial differential equations, J. of Comput. and Appl. Mathematics, 128 (2001) 399 – 422. [9] Schek, I., Jortner, I., Sage, M.L., Application of the Magnus expansion for higher-order multiphoton excitation, Chem. Phys., 59 (1981) 11 – 27. [10] Botchev, M.A., Verwer, J.G.., Numerical integration of damped Maxwell equations, Technical Report CWI Amsterdam (2008). [11] Viswanath, D., Symbolic dynamics and periodic orbits of the Lorenz attractor, Nonlinearity, 16 (2003) 1035 – 1056. [12] Moan, P.C., Niesen, J., Convegence of the Magnus series, Found. Comput. Math., 8 (20089 291 – 301. [13] [13] Moler, C. B., Van Loan, C. F., Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Review, 45 (2003), 3 – 49. [14] Celledoni, E., Iserles, A., Approximating the exponential from a Lie algebra to a Lie group, Math. Comput., 69 (2000) 1457 – 1480. [15] Iserles, A., On Cayley-transform methods for the discretization of Lie-group equations, Found. Comput. Math., 1 (2000) 129 – 160. [16] Haughton D.M., Circular shearing of compressible elastic cylinders, Q. Jl. Mech. Appl. Math., 46 (1993) 471 – 486. [17] Fosdick R., Foti P., Fraddosio A., Marzano S., Shear driven planar Couette and Taylor-like instabilities for a class of compressible isotropic elastic solids, ZAMP, 61 (2010) 537 – 554. [18] Fosdick R., Foti P., Fraddosio A., Marzano S., Piccioni, M. D., Taylor-like bifurcations for a compressible isotropic elastic tube, Mathematics and Mechanics of Solids (MMS), published on line (2013) ISSN: 1081-2865, DOI: 10.1177/1081286513496576. [19] Fosdick, R., Foti, P., Fraddosio, A., Piccioni, M.D., Lower bound estimate of critical loads for isotropic elastic solids, Cont. Mech. and Thermodynamics (CMT), 22 (2010) 77 – 97. [20] Fosdick R., Foti P., Fraddosio A., Marzano S., Piccioni, M. D., A Lower Bound Estimate of the critical load in Bifurcation Analysis for Incompressible Elastic Solids, Mathematics and Mechanics of Solids (MMS), invited contribution on the special number in honor of K.R. Rajagopal, (2014). [21] Pease III, M.C.. Methods of Matrix Algebra, Academic Press, New York and London, (1965).