Consider two square matrices $A_1$ and $A_2$. Consider the following matrix involving matrix trigonometric functions:

\begin{equation} M_1(t)=\begin{bmatrix} \cos(tA_1) & t\mathrm{sinc}(t A_1) \\ -A_1\sin(tA_1) & \cos(tA_1) \end{bmatrix} \end{equation} and similarly $M_2(t)$ defined by changing $A_1$ to $A_2$. Using the double-angle identities, it can be shown that

\begin{align} \delta &= M_1(2t_1)-M_2(2t_2) \\ &= 2\begin{bmatrix} t_1\mathrm{sinc}(t_1A_1) & -t_2\mathrm{sinc}(t_2A_2) \\ \cos (t_1 A_1) & -\cos(t_2A_2) \end{bmatrix} \begin{bmatrix}-A_1\sin(t_1A_1) & \cos(t_1 A_1) \\ -A_2\sin(t_2A_2) & \cos(t_2 A_2) \end{bmatrix} \end{align}

which provides a factorization of the difference $\delta$.

This equality can be checked in MMA using random values for $A_1,A_2,t_1,t_2$:

M1[t_] := 
 ArrayFlatten[{{MatrixFunction[Cos, t*A1], 
    t*MatrixFunction[Sinc, t*A1]}, {-A1.MatrixFunction[Sin, t*A1], 
    MatrixFunction[Cos, t*A1]}}]
M2[t_] := 
 ArrayFlatten[{{MatrixFunction[Cos, t*A2], 
    t*MatrixFunction[Sinc, t*A2]}, {-A2.MatrixFunction[Sin, t*A2], 
    MatrixFunction[Cos, t*A2]}}]
delta := M1[2 t1] - M2[2 t2]
zero := delta - 
  2 ArrayFlatten[{{t1*MatrixFunction[Sinc, t1*A1], -t2*
        MatrixFunction[Sinc, t2*A2]}, {MatrixFunction[Cos, 
        t1*A1], -MatrixFunction[Cos, 
         t2*A2]}}].ArrayFlatten[{{-A1.MatrixFunction[Sin, t1*A1], 
       MatrixFunction[Cos, t1*A1]}, {-A2.MatrixFunction[Sin, t2*A2], 
       MatrixFunction[Cos, t2*A2]}}]

Block[{A1 = RandomReal[{-1, 1}, {2, 2}], 
   A2 = RandomReal[{-1, 1}, {2, 2}], t1 = RandomReal[10], 
   t2 = RandomReal[10]}, zero] // Chop

(* {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}} *)

My question is, would it have been possible to find this factorization using MMA?

  • $\begingroup$ There is an infinity of matrix products that equal δ. What is the criterion for selecting the one in the question? $\endgroup$
    – bbgodfrey
    Nov 10, 2017 at 21:38
  • $\begingroup$ @bbgodfrey In the question I ask to find this precise factorization. The reason is that $\det \delta = 1$ iff at least one of the two matrices of the factorized form is singular. Having $\delta = I_n \delta$ would not be useful. Here, it provides a factorization of $\det\delta$. Hope that answers your question. $\endgroup$
    – anderstood
    Nov 10, 2017 at 22:13
  • $\begingroup$ @bbgodfrey I doubt your were asking for that level of detail, but this factorization is associated with a very nice physical interpretation (in short, symmetry of some motions of impact oscillators), see this article, near equation (3.15) --- unfortunately I can't think of a short explanation... It is used to find such motions. $\endgroup$
    – anderstood
    Nov 10, 2017 at 22:19
  • $\begingroup$ Are $ A_1 $ and $ A_2 $ simultaneouly diagonalizable? $\endgroup$ Oct 7, 2018 at 7:43
  • $\begingroup$ @ΑλέξανδροςZεγγ unfortunately not... $\endgroup$
    – anderstood
    Oct 8, 2018 at 12:49

2 Answers 2


I did the problem in reverse and showed the matrix product (which I called δ2) reduces to the difference of matrices with

δ2 // FunctionExpand // Simplify 

Only then I could reduce the Sinc[A] as Sin[A]/A. The result is that of a difference of two matrices.

  • 2
    $\begingroup$ OK but my question was on factorizing, not expanding :) $\endgroup$
    – anderstood
    Nov 10, 2017 at 20:20

this is a question with some answer that might interest You:

mathematica program for PLUR decomposition of a symbolic matrix using full pivot.

Nobody needs it, nobody wants it following @daniel-lichtblau.

For numerical decomposition, You have the choices in Mathematica: AdvancedMatrixOperations.

The major problem atop is that



is not a function pair like Sin/ArcSin are. And that is even worse because this is the MatrixFunction of Sinc.

Simple example: rotation

{Q, R} = QRDecomposition[{{1, 1}, {-1, 1}}]
(* {{{1/Sqrt[2], -(1/Sqrt[2])}, {1/Sqrt[2], 1/Sqrt[2]}}, {{Sqrt[2],0}, {0, Sqrt[2]}}} *)


Solve[RotationMatrix[\[CurlyPhi] ] == Q, \[CurlyPhi]][[1]] /.C[1] -> 0 
(*{\[CurlyPhi] -> \[Pi]/4} *)

A simple example: regular matrix

S = MatrixPower[Transpose[M].M, 1/2]; R = Inverse[Transpose[M]].S(* rotationmatrix R^t.R=Identity*)

A variation:

m = {{1, 1}, {-1, 1}};
FullSimplify @ Solve[ConjugateTranspose[RotationMatrix[θ]].ScalingMatrix[{s1, s2}] == m, 
  {s1, s2, θ}, Reals] /. C[1] -> 0

As far as I understand the given source:

periodic solutions of n-dof autonomous vibro-impact oscillators with one lasting contact phase

There are three matrices:

M={{m1,0},{0,m2}} (*mass matrix*)


K={{k1+k2,-k2},{-k2,k2}} (coupling matrix of harmonic oscillators)

both from formular (2.5)



from (3.9).

This A matrix is then put into the MatrixExp to give (3.11) and (3.12). In both of them L is a scalar out of Reals. And that remains in all the following sections.

This seems together with the commend more or less the solution path to the proof of theorem 3.1 Please reformulate.

  • $\begingroup$ PLUR seems to be unrelated to the question that was asked. $\endgroup$ Jan 25, 2021 at 0:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.