Yes, you can do trigonometry on matrices using the Taylor series, and identities will hold.

taylor.py

import math
import numpy as np

import numpy as np
from numpy.linalg import matrix_power

def pow(a, n):
    if isinstance(a, np.ndarray):
        return matrix_power(a, n)
    else:
        return a ** n

class TaylorExpansion:
    def __init__(self, func):
        self.func = func        

    def __call__(self, x, N=10):
        if isinstance(N, int):
            return sum(self.func(x, i) for i in range(N))
        else:
            raise TypeError('N must be an integer.')

        return NotImplemented

@TaylorExpansion
def exp(x, i):
    return pow(x, i) / math.factorial(i)

@TaylorExpansion
def log1p(x, i):
    if i == 0:
        return 0
    return (-1) ** (i+1) * pow(x, i) / i

@TaylorExpansion
def sin(x, i):
    k = 1 + 2 * i
    return (-1)**i * pow(x, k) / math.factorial(k)

@TaylorExpansion
def cos(x, i):
    return (-1)**i * pow(x, 2*i) / math.factorial(2*i)

def exponent(a, b):
    return exp(b * log1p(b - 1))

if __name__ == '__main__':
    A = np.array([[1, 2], [3, 4]])
    print(f'array A:\n{A}')
    print(f'sin(A):\n{sin(A)}')
    I_approx = sin(A)**2 + cos(A)**2
    print(f'sin^2 + cos^2:\n{I_approx}')