Matrix Trigonometry

Jelly, 23 22 bytes

³æ*÷!
®Ḥ‘©r0Ç€s2_@/µÐL

Try it online!

Background

This approach directly computes the Taylor series for sine and cosine, i.e.,

It increases the number of initial terms of both series until the result no longer changes, so its accuracy is only limited by the precision of the floating point type.

How it works

®Ḥ‘©r0Ç€s2_@/µÐL  Main link, Argument: A (matrix)

             µÐL  Loop; apply the chain until the results are no longer unique.
                  Return the last unique result.
®                   Yield the value of the register (initially zero).
 Ḥ                  Unhalve/double it.
  ‘©                Increment and copy the result (n) to the register.
    r0              Range; yield [n, ..., 0].
      ǀ            Apply the helper link to each k in the range.
        s2          Split the results into chunks of length 2. Since n is always
                    odd, this yields [[Ç(n), Ç(n-1)], ..., [Ç(1), Ç(0)]].
          _@/       Reduce the columns of the result by swapped subtraction,
                    yielding [Ç(1) - Ç(3) + ... Ç(n), Ç(0) - Ç(2) + ... Ç(n - 1)].


³æ*÷!             Helper link. Argument: k (integer)

³                 Yield the first command-line argument (A).
 æ*               Elevate A to the k-th power.
    !             Yield the factorial of k.
   ÷              Divide the left result by the right one.

Julia, 33 19 bytes

A->reim(expm(A*im))

This is a function that accepts a 2-dimensional array of floats and returns a tuple of such arrays correponding to the cosine and sine, respectively. Note that this is the reverse of the order given in the test cases, in which sine is listed first.

For a real-valued matrix A, we have

and

That is, the sine and cosine of A correspond to the imaginary and real parts of the matrix exponential e^iA. See Functions of Matrices (Higham, 2008).

Try it online! (includes all test cases)

Saved 14 bytes thanks to Dennis!

Mathematica, 27 bytes

{Im@#,Re@#}&@MatrixExp[I#]&

Based on @Rainer P.'s solution.

Takes the square matrix A as an argument and outputs a list containing {sin(A), cos(A)}.

The input is formatted with N to get a numerical value instead of a long exact formula and Column to display the results of sin(A) and cos(A) as separate matrices instead of a nested list.

Calculating the values separately requires 38 bytes

{(#2-#)I,+##}/2&@@MatrixExp/@{I#,-I#}&