Matrix Trigonometry

13

3

Introduction

The two most common trigonometric functions, sine and cosine (or sin and cos for short), can be extended to be matrix-valued functions. One way to compute the matrix-valued analogs is as follows:

Consider these two important trigonometric identities:

trig identities

Using these identities, we can derive the following equations for sin and cos:

trig functions

The matrix exponential exists for all square matrices and is given by:

matrix exponential

where A0 is the identity matrix I with the same dimensions as A. Using the matrix exponential, these two trigonometric functions (and thus all the other trigonometric functions) can be evaluated as functions of matrices.

The Challenge

Given a square matrix A, output the values of sin(A) and cos(A).

Rules

  • Input and output may be in any convenient, reasonable format (2D array, your language's matrix format, etc.).
  • You may write a single program, two independent programs, a single function, or two functions. If you choose to write two functions, code may be shared between them (such as imports and helper functions).
  • The input matrix's values will always be integers.
  • Your solution may have accuracy issues as the result of floating-point imprecision. If your language had magical infinite-precision values, then your solution should work perfectly (ignoring the fact that it would require infinite time and/or memory). However, since those magical infinite-precision values don't exist, inaccuracies caused by limited precision are acceptable. This rule is in place to avoid complications resulting from requiring a specific amount of precision in the output.
  • Builtins which compute trigonometric functions for matrix arguments (including hyperbolic trig functions) are not allowed. Other matrix builtins (such as multiplication, exponentiation, diagonalization, decomposition, and the matrix exponential) are allowed.

Test Cases

Format: A -> sin(A), cos(A)

[[0]] -> [[0]], [[1]]
[[0, 2], [3, 5]] -> [[-0.761177343863758, 0.160587281888277], [0.240880922832416, -0.359709139143065]], [[0.600283445979886, 0.119962280223493], [0.179943420335240, 0.900189146538619]]
[[1, 0, 1], [0, 0, 0], [0, 1, 0]] -> [[0.841470984807897, -0.158529015192103, 0.841470984807897], [0, 0, 0], [0, 1, 0]], [[0.540302305868140, -0.459697694131860, -0.459697694131860], [0, 1, 0], [0, 0, 1]]
[[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1]] -> [[0.841470984807897, 0, 0, 0, 0], [0, 0.841470984807897, 0, 0, 0], [0, 0, 0.841470984807897, 0, 0], [0, 0, 0, 0.841470984807897, 0], [0, 0, 0, 0, 0.841470984807897]], [[0.540302305868140, 0, 0, 0, 0], [0, 0.540302305868140, 0, 0, 0], [0, 0, 0.540302305868140, 0, 0], [0, 0, 0, 0.540302305868140, 0], [0, 0, 0, 0, 0.540302305868140]]
[[-3, 2, -6], [3, 0, 4], [4, -2, 7]] -> [[-0.374786510963954, 0.135652884035570, -1.35191037980742], [1.14843105375406, 0.773644542790111, 1.21625749577185], [1.21625749577185, -0.135652884035570, 2.19338136461532]], [[4.13614256031450, -1.91289828483056, 5.50873853927692], [-2.63939111203107, 1.49675144828342, -3.59584025444636], [-3.59584025444636, 1.91289828483056, -4.96843623340878]]

Further Reading

This excellent question over at Math.SE includes some alternate derivations of the matrix-valued analogs of trigonometric functions.

Mego

Posted 2016-05-05T16:33:22.410

Reputation: 32 998

I got sin([[1, 0, 1], [0, 0, 0], [0, 1, 0]]) = {{0.841, -0.158, 0.841}, {0, 0, 0}, {0, 1, 0}} with Mathematica, can you check? – kennytm – 2016-05-05T17:37:14.387

1@kennytm That is what the test case shows. – Mego – 2016-05-05T17:38:01.250

So it is disallowed to sum some terms of the series, because adding a finite number of terms does not give an infinitely accurate value? – feersum – 2016-05-05T18:51:05.270

@feersum Correct. Unless it converges in a finite number of terms (which it doesn't), that's not a valid approach. – Mego – 2016-05-05T19:25:34.243

1@Mego Apparently all of the existing answers should be deleted then. – feersum – 2016-05-05T19:30:26.803

@feersum There are other ways to compute the matrix exponential other than just evaluating the series, such as diagonalization/Jordan decomposition. It's fair to assume that a builtin would use one of the more robust methods, rather than summing a finite number of terms of the sequence. – Mego – 2016-05-05T19:32:54.960

3@Mego It's completely unreasonable to think that all of the floating-point based builtins use an exact algorithm (or one that would be exact if floating-point operations were replaced with "real number" operations). – feersum – 2016-05-05T19:35:03.710

I'm not sure I'm interpreting this discussion correctly, but irrational numbers cannot be represented exactly as decimals or floating point numbers, unless you have both infinite memory to store the representation and infinite time to calculate its value. My answer keeps adding terms of the Taylor series until the result no longer changes. It gives an inaccurate result, but its accuracy increases as the precision of the floating point type increases. Is that allowed? – Dennis – 2016-05-05T19:41:25.307

@feersum I tried to make this as simple as possible when writing the challenge. Let me try stating it a different way: Use matrix exponential builtins if you want. Don't worry about inaccuracies due to limited precision with floating point values. Assume that, if your language had magical infinite-precision numeric values, the builtins would use them and everything would work like you expect. – Mego – 2016-05-05T19:46:36.237

@Dennis I'm going to allow that approach. What I specifically meant when I said that summing a finite number of terms was not valid was, if you sum some N terms of the series, where N is constant, it's not valid. If you sum terms of the series until the delta hits 0, that's valid (because N would be dependent on the precision of the numeric values, rather than a constant). Summing 5 terms and calling it good enough is not ok. Summing terms until they go below the machine/language floating point epsilon is ok. – Mego – 2016-05-05T19:50:03.003

Summing until the delta becomes zero would fail if the language did have infinite precision values as the program would not halt. – feersum – 2016-05-05T19:51:41.747

1@feersum I've addressed that in my latest edit: (ignoring the fact that it would require infinite time and/or memory) – Mego – 2016-05-05T19:52:38.377

Hrm. The current conditions sound like they forbid the technique of using the recurrence (cos(x) ,sin(x)) = (2 cos(x/2)^2 - 1, 2 sin(x/2) cos(x/2)) until the argument is small enough that you can just use cos(x) = 1 - x^2/2 and sin(x) = x. – None – 2016-05-05T23:11:22.647

(1) I've fixed that. (2) I'm rather sure that they do. (3) The notion of well-ordering has no relation to this topic. There are plenty of norms one can put on matrices. (4) The same technique works for negative real numbers. – None – 2016-05-05T23:19:31.243

Answers

6

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

sine

and

cosine

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

Try it online! (includes all test cases)

Saved 14 bytes thanks to Dennis!

Alex A.

Posted 2016-05-05T16:33:22.410

Reputation: 23 761

6

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.

Example

Calculating the values separately requires 38 bytes

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

miles

Posted 2016-05-05T16:33:22.410

Reputation: 15 654

6

Jelly, 23 22 bytes

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

Try it online!

Background

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

formula

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.

Dennis

Posted 2016-05-05T16:33:22.410

Reputation: 196 637

3

C++, 305 bytes

#include<cmath>
#include<iostream>
#include<vector>
int x,i=0, j;void p(std::vector<double> v){int x=sqrt(v.size());for(i=0;i<x;i++){for(j=0;j<x;j++) std::cout << v[x] << " ";std::cout << "\n";}}int main(){std::vector<double> s, c;while(std::cin >> x){s.push_back(sin(x));c.push_back(cos(x));}p(s);p(c);}

Input is a list of numbers that are a perfect square on stdin. Output is a pretty printed 2d array on stdout

HSchmale

Posted 2016-05-05T16:33:22.410

Reputation: 181

2

Matlab, 138 121 52 50 bytes

Since matrix exponentiation is allowed, (what I didnt notice first, d'oh) I do not need to define my helper funciton anymore, and the whole thing can be trivially solved:

A=input('')*i;a=expm(A);b=expm(-A);[(b-a)*i,a+b]/2

The input should be a matrix e.g. [1,2;4,5] or alternatively [[1,2];[3,4]]

An unexpected thing (in hindsight not so unexpected) is that the cosine and sine matrix still satisfy

I = sin(A)^2+cos(A)^2

flawr

Posted 2016-05-05T16:33:22.410

Reputation: 40 560

Isn't A^0 the same as eye(size(A))? – FryAmTheEggman – 2016-05-05T17:15:41.020

Oh, you're right, thank you! – flawr – 2016-05-05T17:17:02.587

2Why not use expm? – Luis Mendo – 2016-05-05T17:18:35.177

Matrix exponential is not forbidden. Only matrix sine/cosine/secant/cosecant. – Mego – 2016-05-05T17:22:48.957

2As per the identity: I should hope they satisfy that identity, considering that the scalar form was used to extend the functions to matrices! – Mego – 2016-05-05T17:27:37.700

1Well then the whole thing becomes almost trivial. – flawr – 2016-05-05T18:27:49.497

@LuisMendo Since sinm and cosm were nto allowed, I assumed that matrix exponentiation itself wasn't allowed either. – flawr – 2016-05-05T18:35:58.527

2

Matlab, 37 bytes

@(A){imag(expm(i*A));real(expm(i*A))}

Rainer P.

Posted 2016-05-05T16:33:22.410

Reputation: 2 457

0

Julia 0.4, 28 bytes

A->imag({E=expm(im*A),im*E})

Input is a matrix of floats, output is an array of matrices. Try it online!

Dennis

Posted 2016-05-05T16:33:22.410

Reputation: 196 637

0

Sage, 44 bytes

lambda A:map(exp(I*A).apply_map,(imag,real))

Try it online.

This anonymous function returns a list of 2 matrices corresponding to sin(A) and cos(A), respectively. exp(I*A) computes the matrix exponential for I*A (A with all elements multiplied by the imaginary unit), and matrix.apply_map(f) returns a matrix where f has been applied to all of its elements. By applying imag and real (the functions for getting the imaginary and real parts of a scalar value) to the matrices, we get the values of sin(A) and cos(A), thanks to Euler's famous identity (referenced in the challenge text).

Mego

Posted 2016-05-05T16:33:22.410

Reputation: 32 998