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

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.

• Use exact integer for m, i.e., m={{1,1,3},...}, or use m=Rationalize@m Or see Chop May 1, 2020 at 5:56
• @OkkesDulgerci: Are you serious? Think of matrices of bigger sizes. Thank you anyway. The Chop command does not help at all. BTW, this is a slightly modified example from Mathematica documentation and I follow that documentation. May 1, 2020 at 5:59
• Yes, I am serious. cosm = MatrixFunction[Cos, m] // Chop or m = Rationalize@{{1., 1., 3.}, {0., 1., 0.}, {0., 1.0*I, 2.}}; cosm = MatrixFunction[Cos, m] // N work May 1, 2020 at 6:02
• @OkkesDulgerci: Please, pay your attention to the second elements of the first rows. Don't hesitate to ask for further explanation in need. May 1, 2020 at 6:04
• For what it's worth, MatrixFunction[Cos, m]+I*MatrixFunction[Sin, m] agrees with MatrixExp[I*m] so there is some consistency. And MatrixExp is numerically a fairly reliable function, to the extent that it can be at least. May 1, 2020 at 16:54

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 *)

• Thank you. Does the approach suggested by you work for matrices of size say 20? May 1, 2020 at 6:37
• Having tried that approach for a random matrix over complexes of size 20, I obtained a satisfactory result so I accept your answer. May 1, 2020 at 6:58
• EDIT previous comment mistakenly used MatrixPower when I meant to say MatrixFunction. I'm not suggesting that you use the Sum and MatrixPower. It was merely a means of checking the Mma result from MatrixFunction. The fact that the sum gave the same result provides additional confidence in the Mma MatrixFunction result. Try a similar test with Maple to see if its result checks out. If it doesn't, then there would be even more reason to trust the Mma MatrixFunction results. – May 1, 2020 at 7:04

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, 3 I Cos - I (3 Cos - (3 + I) Sin), -3 Cos + 3 Cos},
{0, Cos, 0}, {0, -I Cos + I Cos, Cos}}

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, Cos', 0}, {0, Cos, 0}, {0, 0, Cos}}.Inverse[sm]
{{Cos, -3 I Cos + 3 I Cos - (1 - 3 I) Sin, -3 Cos + 3 Cos},
{0, Cos, 0}, {0, -I Cos + I Cos, Cos}}

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 - 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.

• +1. Thank you. It's useful. May 2, 2020 at 18:30