Why the difference in the results of MatrixFunction between Mathematica and Maple?

Question

Mastering my Mathematica skills, I consider an example

m = {{1., 1., 3.}, {0., 1., 0.}, {0., 1.0*I, 2.}};
cosm = MatrixFunction[Cos, m]

$$\left( \begin{array}{ccc} 0.540302\, +0. i & -0.841471-0.344934 i & -2.86935+0. i \\ 0.\, +0. i & 0.540302\, +0. i & 0.\, +0. i \\ \text{2.7755575615628914$\grave{ }$*${}^{\wedge}$-17}+0. i & 0.\, -0.956449 i & -0.416147+0. i \\ \end{array} \right)$$ Just to compare, I do the same in Maple 2019.1 by

m := Matrix([[1., 1., 3.], [0., 1., 0.], [0., 1.0*I, 2.]]):
LinearAlgebra:-MatrixFunction(m, cos(x), x);

$$ \left[ \begin {array}{ccc} 0.5403023059&- 0.6666666665- 0.8693474277 \,i&- 2.869347427\\ 0.0& 0.5403023059& 0.0 \\ 0.0&- 0.9564491424\,i&- 0.4161468365\end {array} \right] $$ I don't understand the difference between the second element of the first row of both results. Of course, I looked in the documentation and found it too poor in both CASes.

Answer 1

m = {{1., 1., 3.}, {0., 1., 0.}, {0., 1.0*I, 2.}} // Rationalize;

cosm = MatrixFunction[Cos, m] // N

(* {{0.540302, -0.841471 - 0.344934 I, -2.86935}, {0., 0.540302, 0.}, {0., 
  0. - 0.956449 I, -0.416147}} *)

Checking using MatrixPower

coef[n_] = SeriesCoefficient[Cos[x], {x, 0, n}]

(* Piecewise[{{(I^n*(1 + (-1)^n))/(2*n!), n >= 0}}, 0] *)

cosm == Sum[coef[n]*MatrixPower[m, n], {n, 0, Infinity}] // N

(* True *)

Answer 2

Since Bob's answer has already been accepted, I will just leave a more extensive note on how to check for consistency. This is mostly a rehash of my other answers, so please refer to them for further details.

The starting point, as always, is to recall the (Cauchy integral-like) definition

$$f(\mathbf A) = \frac{1}{2\pi i} \oint_\gamma f(z)\, (z \mathbf I- \mathbf A)^{-1}\,\mathrm dz$$

where $\gamma$ is a closed contour enclosing the eigenvalues of $\mathbf A$, and where $f(z)$ is analytic within the contour.

First, this definition can be used to prove Daniel's identity in the comments, by letting $f(z)=\exp(i z)=\cos z+i\sin z$ and then decomposing accordingly.

Thus,

mat = {{1., 1., 3.}, {0., 1., 0.}, {0., 1.0 I, 2.}};

MatrixFunction[Cos, mat] // Chop
   {{0.5403023058681398, -0.8414709848078964 - 0.34493447282215695 I, -2.869347427245846},
    {0, 0.5403023058681398, 0},
    {0, 0. - 0.9564491424152821 I, -0.41614683654714246}}

(MatrixExp[I mat] + MatrixExp[-I mat])/2 // Chop
   {{0.5403023058681397, -0.8414709848078965 - 0.34493447282215695 I, -2.869347427245847},
    {0, 0.5403023058681397, 0},
    {0, 0. - 0.9564491424152821 I, -0.41614683654714235}}

---

Of course, the contour integral formula itself can be used for computation. I will temporarily consider the exact version of the OP's matrix for this demonstration:

mex = {{1, 1, 3}, {0, 1, 0}, {0, I, 2}};
eig = Eigenvalues[mex]
   {2, 1, 1}

At this juncture, we note that all the eigenvalues are real.

We can then employ the residue theorem and the Cauchy integral theorem to convert the evaluation to a sum of residues of the integrand over the eigenvalues of the given matrix:

Sum[Map[Residue[#, {z, λ}] &, 
        Cos[z] Inverse[z IdentityMatrix[Length[mex]] - mex], {2}],
    {λ, Union[eig]}]
   {{Cos[1], 3 I Cos[2] - I (3 Cos[1] - (3 + I) Sin[1]), -3 Cos[1] + 3 Cos[2]},
    {0, Cos[1], 0}, {0, -I Cos[1] + I Cos[2], Cos[2]}}

N[%]
   {{0.5403023058681398, -0.8414709848078965 - 0.3449344728221573 I, -2.8693474272458466},
    {0., 0.5403023058681398, 0.},
    {0., 0. - 0.9564491424152821 I, -0.4161468365471424}}

Compare this with the more popular evaluation method that uses the Jordan decomposition:

{sm, jm} = JordanDecomposition[mex]
   {{{1, 0, 3}, {0, 1/10 + 3 I/10, 0}, {0, 3/10 - I/10, 1}},
    {{1, 1, 0}, {0, 1, 0}, {0, 0, 2}}}

sm.{{Cos[1], Cos'[1], 0}, {0, Cos[1], 0}, {0, 0, Cos[2]}}.Inverse[sm]
   {{Cos[1], -3 I Cos[1] + 3 I Cos[2] - (1 - 3 I) Sin[1], -3 Cos[1] + 3 Cos[2]},
    {0, Cos[1], 0}, {0, -I Cos[1] + I Cos[2], Cos[2]}}

N[%]
   {{0.5403023058681398, -0.8414709848078965 - 0.3449344728221573 I, -2.8693474272458466},
    {0., 0.5403023058681398, 0.},
    {0., 0. - 0.9564491424152821 I, -0.4161468365471424}}

---

The contour integral formula also readily lends itself to numerical evaluation. Earlier, we noted that the eigenvalues of mat are real, so a convenient choice for the contour $\gamma$ is an axis-aligned rectangle enclosing the eigenvalues:

With[{ε = 1/20},
     contour = (Tuples[{MinMax[Eigenvalues[mat]] + {-ε, ε},
                        {-ε, ε}}].{1, I})[[{1, 3, 4, 2, 1}]]];

NIntegrate[] can then be used for the evaluation:

NIntegrate[Cos[z] Inverse[z IdentityMatrix[3] - mat],
           {z, Sequence @@ contour} // Evaluate]/(2 π I) // Chop
   {{0.5403023058681499, -0.8414709848079495 - 0.3449344728223653 I, -2.8693474272458817},
    {0, 0.5403023058681499, 0},
    {0, 0. - 0.9564491424152939*I, -0.41614683654714557}}

This evaluation will throw a few NIntegrate::izero, because some of the matrix elements are zero. Nevertheless, the result is consistent with all the other methods previously presented.