**Expected behavior**: Using `simplify`

in the Symbolics.jl package, return the identity for the operation `A^-1 * A`

for a square matrix involving multiplication of trigonometric functions in the matrix elements (assuming `det(A) != 0`

).

**Actual behavior**: Using `simplify`

does not reduce the final expression to the identity, but substituting numerical values at the end of symbolic calculations does evaluate to the identity matrix.

I think the issue has to do with how `lu()`

is computed, but I’m not sure what it’s doing or how to fix it.

Is there a preferred way to do these kinds of operations? I’ve gone through the Symbolics.jl documentation and can’t find much about matrices that include trigonometric functions as elements.

Here’s a simple case with the 2x2 rotation matrix (I’ve tried other variations, too):

```
julia> using LinearAlgebra, Symbolics, SymbolicUtils
julia> @variables x::Real;
julia> I2 = [1 0; 0 1]; #identity matrix
julia> A = [cos(x) -sin(x); sin(x) cos(x)]
2×2 Matrix{Num}:
cos(x) -sin(x)
sin(x) cos(x)
julia> simplify.(inv(A) * A)
2×2 Matrix{Num}:
1 … sin(x)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1) - (sin(x)*(true - ((cos(x)^-1)*(sin(x)^2)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1)))*(cos(x)^-1))
0 cos(x)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1) + (cos(x)^-1)*(sin(x)^2)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1)
```

By eye, the result should simplify to `I2`

(except for those `true`

values). Trying to simplify `inv(A)`

first does not work. Directly accessing the LU factorization method and using left-division instead of `inv`

also does not work (i.e. solving `A^-1 = A / I`

. See Base.inv):

```
julia> simplify.(lu(A) \ I2 * A)
2×2 Matrix{Num}:
1 … sin(x)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1) - (sin(x)*(1 - ((cos(x)^-1)*(sin(x)^2)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1)))*(cos(x)^-1))
0 cos(x)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1) + (cos(x)^-1)*(sin(x)^2)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1)
```

Again, it evaluates the expression correctly, but can’t simplify further (even though it should be able to).

When setting `x = π/4`

(or any other value) using `substitute`

, the expression correctly reduces to the identity matrix, but some of the elements are converted to floats:

```
julia> B = inv(A) * A;
julia> answer = simplify.(substitute.(B, (Dict(x => π/4),)))
2×2 Matrix{Num}:
1 0.0
0 1.0
julia> answer == I2
true
```

I found one workaround, bypassing `inv`

altogether, but it only works for 2x2 matrices.

Using the Caley-Hamilton method (Wikipedia):

```
julia> A_inv = det(A) * (tr(A) * I2 - A)
2×2 Matrix{Num}:
cos(x)*(cos(x)^2 + sin(x)^2) sin(x)*(cos(x)^2 + sin(x)^2)
-sin(x)*(cos(x)^2 + sin(x)^2) cos(x)*(cos(x)^2 + sin(x)^2)
julia> simplify.(A_inv * A)
2×2 Matrix{Num}:
1 0
0 1
julia> simplify.(A_inv * A) == I2
true
```