# How to solve coupled complex-valued equations faster and get simplified solution

I have six coupled equations. How to solve these equations symbolically faster in Mathematica. In these equations, variables p,q,x,y are Reals and variables a,b,c are complex terms. These equations should be solved for {R11, R22, R33, R21, R31, R32}, but I want to only get the solution for R21.

eqns={ R11 + R22 + R33 == 1,
p Im[R21] == R22 x + R33 y,
q Im[R32] == R33 (y + z),
(I p R11)/2 + a R21 - (I p R22)/2 + (I q R31)/2 == 0,
(I q R21)/2 + b R31 - (I p R32)/2 == 0,
(I q R22)/2 - (I p R31)/2 + c R32 - (I q R33)/2 == 0};

term = {R11, R22, R33, R21, R31, R32};


I can do the series expansion for q and p. I try to solve these equations as a series solution of p,q given here...

Series[R21 /. Solve[eqns, term], {p, 0, 3}, {q, 0, 3}] // Timing


Do you have better solution? I also need to define that variables p,q,x,y are Reals and variables a,b,c are complex.

Edited to simplify the derivation and reduce the length of the final result.

If R21 and R32 are pure imaginary, then the equation is solved easily.

R21 /. Solve[eqns /. {Im[R21] -> -I R21, Im[R32] -> -I R32}, term][[1]]
(* (I p (2 b q^2 x + 2 b q^2 y - 4 b c x y - p^2 x y - 4 b c x z -
p^2 x z))/(6 b p^2 q^2 - 4 a b q^2 x - p^2 q^2 x - q^4 x -
8 b c p^2 y - 2 p^4 y - 4 a b q^2 y + 3 p^2 q^2 y - q^4 y +
8 a b c x y + 2 a p^2 x y + 2 c q^2 x y - 8 b c p^2 z - 2 p^4 z +
p^2 q^2 z + 8 a b c x z + 2 a p^2 x z + 2 c q^2 x z) *)


which indeed are pure imaginary, if all the parameters are real.

The general case also can be solved, although the result is lengthy. Since the challenge here arises from the second and third of eqns and R21 and R32, we begin by eliminating the other equations and variables.

Flatten[Solve[Delete[eqns, {{2}, {3}}], Delete[term, {{4}, {6}}]]];
s = Simplify[Delete[eqns, {{1}, {4}, {5}, {6}}] /. %]
(* {p Im[R21] == -((I ((q R21 - p R32) (p^2 (x - 2 y) + q^2 (x + y)) +
b (4 a q R21 (x + y) + p (-4 c R32 (x - 2 y) + 2 I q (x + y)))))/(6 b p q)),
q Im[R32] == -((I ((2 p^2 - q^2) (-q R21 + p R32) +
b (2 I p q + 4 a q R21 + 8 c p R32)) (y + z))/(6 b p q))} *)


Next, we Solve to obtain R32 and then its imaginary part, and substitute both into the remaining equation.

Collect[Solve[s[[1]], R32][[1, 1, 2]], {R21, Im[R21]}, Simplify];
Im[%];
Simplify[s[[2]] /. {R32 -> %%, Im[R32] -> %}]
(* q (Im[(q R21 (p^2 (x - 2 y) + 4 a b (x + y) + q^2 (x + y)))/
(p (4 b c (x - 2 y) + p^2 (x - 2 y) + q^2 (x + y)))] +
2 Re[(b q (x + y))/(4 b c (x - 2 y) + p^2 (x - 2 y) + q^2 (x + y))] -
6 Re[(b p q Im[R21])/(4 b c (x - 2 y) + p^2 (x - 2 y) + q^2 (x + y))]) ==
-(((y + z) ((-p^3 + 2 I a p^2 R21 + 2 I c q^2 R21 - 4 b c (p - 2 I a R21)) x + p^2
(8 b c + 2 p^2 - q^2) Im[R21]))/(p (4 b c (x - 2 y) + p^2 (x - 2 y) + q^2 (x + y)))) *)


Solve and Reduce seem unable to solve this equation as written, so instead we treat the real and imaginary parts of R21 as independent variables.

Simplify[% /. {Re[z_  Im[R21]] -> Re[z] R21i, Im[z_  R21] -> R21i Im[z] + R21r Re[z],
Im[R21] -> R21i} /. R21 -> R21r + I R21i]
Simplify[CoefficientArrays[%, {R21r, R21i}] // Normal];
Map[{R21r, R21i}.# &, Transpose[ReIm[%[[2]]]]] + ReIm[%[[1]]];
Simplify[%, (R21r | R21i | p | q | x | y | z) ∈ Reals];
Map[# == 0 &, %];


Finally, we solve for the real and imaginary parts and then reconstitute R21.

Simplify[(R21r + I R21i) /. Flatten[Solve[%, {R21r, R21i}]]]


The result, with a LeafCount of 1229, is a bit long to be reproduced here. The computations described above take only seconds, however. Also, it seems likely that the final result can be simplified somewhat, given enough effort.

• Final simplify gives me an error...Solve::naqs: (ReIm[-(((4 b c+p^2) x (y+z))/(Times[<<4>>]+Times[<<2>>]+Times[<<2>>]))+2 q Re[b q Plus[<<2>>] Power[<<2>>]]]==0)+(Transpose[{R21r,R21i}.ReIm[{Times[<<6>>]+Times[<<2>>],Times[<<4>>]+Times[<<8>>]+Times[<<6>>]+Times[<<7>>]+Times[<<2>>]+Times[<<3>>]}]]==0) is not a quantified system of equations and inequalities. >> – santosh Aug 23 '15 at 19:02
• Map[{R21r, R21i}.# &, Transpose[ReIm[%[[2]]]]] + ReIm[%[[1]]]; this line is not clear to me. What is ReIm[] here?? – santosh Aug 23 '15 at 19:10
• ReIm[] is the same as {Re[],Im[]}. What version of Mathematica are you using? – bbgodfrey Aug 23 '15 at 19:24
• I am using Mathematica 9 – santosh Aug 23 '15 at 20:24
• ReIm is new to Mathematica 10. Replace ReIm[%[[2]]] by {Re[%[[2]]], Im[%[[2]]]}. I just tried it successfully. Be sure to count brackets carefully! Best wishes. – bbgodfrey Aug 23 '15 at 21:13