2
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ Use exact integer for m, i.e., m={{1,1,3},...}, or use m=Rationalize@m Or see Chop $\endgroup$ May 1, 2020 at 5:56
  • $\begingroup$ @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. $\endgroup$
    – user64494
    May 1, 2020 at 5:59
  • $\begingroup$ 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 $\endgroup$ May 1, 2020 at 6:02
  • $\begingroup$ @OkkesDulgerci: Please, pay your attention to the second elements of the first rows. Don't hesitate to ask for further explanation in need. $\endgroup$
    – user64494
    May 1, 2020 at 6:04
  • 2
    $\begingroup$ 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. $\endgroup$ May 1, 2020 at 16:54

2 Answers 2

1
$\begingroup$
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 *)
$\endgroup$
3
  • $\begingroup$ Thank you. Does the approach suggested by you work for matrices of size say 20? $\endgroup$
    – user64494
    May 1, 2020 at 6:37
  • $\begingroup$ Having tried that approach for a random matrix over complexes of size 20, I obtained a satisfactory result so I accept your answer. $\endgroup$
    – user64494
    May 1, 2020 at 6:58
  • $\begingroup$ 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. – $\endgroup$
    – Bob Hanlon
    May 1, 2020 at 7:04
3
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ +1. Thank you. It's useful. $\endgroup$
    – user64494
    May 2, 2020 at 18:30

Your Answer

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

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