# How can I make Mathematica solve this set of four equations?

I use Mathematica version 12.1 on Windows 10.

I have four equations in four unknowns.

Is there a trick to make Mathematica solve these?

I tried Solve and NSolve and it is not able to do it. Either it hangs (waited long time) or it says can't solve.

I copied the same equations over to Maple, converted them to Maple using Maple's Mathematica translator, and Maple solved them immediately.

Here is the code

ClearAll[a,b,c,d];
T    = 2 Pi;
bt   = MatrixExp[{{a, b}, {c, d}}*T];
cmat = {{Sin, Cos}, {-Cos, Sin}};
eqs  = {cmat[[1, 1]] == bt[[1, 1]],
cmat[[1, 2]] == bt[[1, 2]],
cmat[[2, 1]] == bt[[2, 1]],
cmat[[2, 2]] == bt[[2, 2]]};
Solve[eqs, {a, b, c, d}] (* Hangs, and after about 30 minutes I gave up *)


If I do Solve[eqs, {a, b, c, d}, Reals] it gives

Solve::nsmet This system cannot be solved with the methods available to Solve

I also tried NSolve[eqs, {a, b, c, d}], but it seems to hang also.

Here it is in Maple:

restart;
T    := 2*Pi;
BT   := Matrix([[a,b],[c,d]])*T:
BT   := LinearAlgebra[MatrixExponential](BT);
cmat := Matrix([[sin(1), cos(1)], [-cos(1), sin(1)]]);
eqs  := [cmat[1,1]=BT[1,1],cmat[1,2]=BT[1,2],cmat[2,1]=BT[2,1],cmat[2,2]=BT[2,2]]:
sol  := solve(eqs,[a,b,c,d]);
evalf(sol)

[[a = 0.,
b = 0.09084505695 - (6.571202944*10^(-12))*I,
c = -0.09084505695 + (6.571202944*10^(-12))*I,
d = 0.]]


The answer time is less than one second.

At first I thought the matrix exponentials are different. But I compared them and they are the same. I then copied Mathematica's bt variable over to Maple, and used that, and Maple gave the same answer:

restart;
with(MmaTranslator); # Used to translate mma code to Maple
T    := 2*Pi;
cmat := Matrix([[sin(1), cos(1)], [-cos(1), sin(1)]]);

# This below is bt from Mathematica. This is the result of MatrixExp[{{a, b}, {c, d}}*T];
# done in Mathematica and copied here to see if Maple can solve it

BT    := FromMma({{-(((a - d - Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*E^((a + d - Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*Pi))/(2*Sqrt[a^2 + 4*b*c - 2*a*d + d^2])) +
((a - d + Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*E^((a + d + Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*Pi))/(2*Sqrt[a^2 + 4*b*c - 2*a*d + d^2]),
-((b*E^((a + d - Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*Pi))/Sqrt[a^2 + 4*b*c - 2*a*d + d^2]) + (b*E^((a + d + Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*Pi))/
Sqrt[a^2 + 4*b*c - 2*a*d + d^2]}, {-((c*E^((a + d - Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*Pi))/Sqrt[a^2 + 4*b*c - 2*a*d + d^2]) +
(c*E^((a + d + Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*Pi))/Sqrt[a^2 + 4*b*c - 2*a*d + d^2],
-(((-a + d - Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*E^((a + d - Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*Pi))/(2*Sqrt[a^2 + 4*b*c - 2*a*d + d^2])) +
((-a + d + Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*E^((a + d + Sqrt[a^2 + 4*b*c - 2*a*d + d^2])*Pi))/(2*Sqrt[a^2 + 4*b*c - 2*a*d + d^2])}}):

eqs  := [cmat[1,1]=BT[1,1],cmat[1,2]=BT[1,2],cmat[2,1]=BT[2,1],cmat[2,2]=BT[2,2]]:
sol  := solve(eqs,[a,b,c,d]);
evalf(sol)


gives

[[a = 0.,
b = 0.09084505695 - (6.571202944*10^(-12))*I,
c = -0.09084505695 + (6.571202944*10^(-12))*I,
d = 0.]]


What else should I try to make Mathematica solve these four equations?

The matrix logarithm in Alex's answer will give one out of many possible (complex!) solutions, in complete analogy with the scalar case.

One way to go about this is to simultaneously reduce cmat and the matrix within the exponential to the Jordan form:

{sm, jm} = JordanDecomposition[2 π {{a, b}, {c, d}}]
{sr, jr} = JordanDecomposition[{{Sin, Cos}, {-Cos, Sin}}] // FullSimplify


Conveniently, 1. both jm and jr are diagonal matrices; and 2. both sm and sr are normalized such that their second row is one. We then recall that the Jordan vectors of $$\mathbf A$$ and $$\exp(\mathbf A)$$ should be the same, so:

GroebnerBasis[Thread[First[sm] == First[sr]], {a, b, c, d}]
{b + c, a - d, -2 I c + Sqrt[a^2 + 4 b c - 2 a d + d^2]}


Immediately, we find that $$c=-b$$ and $$a=d$$. We can use this to simplify the next set of equations:

eq = Simplify[TrigToExp[Thread[Diagonal[jr] == Exp[Diagonal[jm]]] /.
{c -> -b, d -> a}], b < 0]
{-I + E^(I + 2 a π + 2 I b π) == 0,
E^(2 (a - I b) π) == -I E^I}


Feeding this to Solve[] (Solve[%, {a, b}] // FullSimplify) and then plugging the results into the original matrix will generate a set of parametrized solutions:

{{I u, 1/4 - 1/(2 π) + v},
{-((-2 + π + 4 π v)/(4 π)), I u}}


and

{{1/2 I (1 + 2 u), -((2 + π)/(4 π)) + v},
{(2 + π - 4 π v)/(4 π), 1/2 I (1 + 2 u)}}


where I have replaced the C[k] with simpler parameters for clarity. Here, u and v are integers. In particular, Alex's solution corresponds to the first set, with u = 0 and v = 0.

• Yes, you are right (+1). Nevertheless my answer is the same as Maple gives. So they probably also using not all set of solutions. May 5, 2020 at 9:34

There is a solution with using MatrixLog

ClearAll[a, b, c, d];
T = 2 Pi; m = T {{a, b}, {c, d}}; q = 1;
cmat = {{Sin[q], Cos[q]}, {-Cos[q], Sin[q]}};

NSolve[m == MatrixLog[cmat], {a, b, c, d}]

Out[]= {{a -> -1.76697*10^-17 + 0. I, b -> 0.0908451 + 0. I,
c -> -0.0908451 + 0. I, d -> 1.76697*10^-17 + 0. I}}


I have no idea why, but this seems to work. Starting with your code

T = 2 Pi;
bt = MatrixExp[{{a, b}, {c, d}}*T];
cmat = {{Sin, Cos}, {-Cos, Sin}};
eqs = {cmat[[1, 1]] == bt[[1, 1]], cmat[[1, 2]] == bt[[1, 2]],
cmat[[2, 1]] == bt[[2, 1]], cmat[[2, 2]] == bt[[2, 2]]};


and "eliminating" one of the variables

elim = FullSimplify @ Eliminate[eqs, a]


gives a new set of equations. I put "eliminate" in quotes because it doesn't actually do that, elim still has the variable a in it. However, plugging this to Solve

Solve[elim, {a, b, c, d}]


gives a bunch of solutions, one of which is the numerical solution you posted. Using Resolve instead, like

FullSimplify @ Reduce[elim, {a, b, c, d}]


gives a parametrized solution with $$a=d$$ and $$b=-c$$, like in J. M.'s elegant answer. 