I have equations that are products of sines and cosines that I want to integrate. As an example, the integrand may look something like this:

$-\cos \left(\frac{\sqrt{3} \pi x (p+q)}{A}\right) \sin \left(\frac{\pi y (p-q)}{A}\right) \sin \left(\frac{\sqrt{3} \pi x (r+s)}{A}\right) \sin \left(\frac{\pi y (r-s)}{A}\right)+\cos \left(\frac{\sqrt{3} \pi q x}{A}\right) \sin \left(\frac{\pi y (2 p+q)}{A}\right) \sin \left(\frac{\sqrt{3} \pi x (r+s)}{A}\right) \sin \left(\frac{\pi y (r-s)}{A}\right)-\cos \left(\frac{\sqrt{3} \pi p x}{A}\right) \sin \left(\frac{\pi y (p+2 q)}{A}\right) \sin \left(\frac{\sqrt{3} \pi x (r+s)}{A}\right) \sin \left(\frac{\pi y (r-s)}{A}\right)-\sin \left(\frac{\sqrt{3} \pi s x}{A}\right) \cos \left(\frac{\sqrt{3} \pi x (p+q)}{A}\right) \sin \left(\frac{\pi y (p-q)}{A}\right) \sin \left(\frac{\pi y (2 r+s)}{A}\right)+\cos \left(\frac{\sqrt{3} \pi q x}{A}\right) \sin \left(\frac{\sqrt{3} \pi s x}{A}\right) \sin \left(\frac{\pi y (2 p+q)}{A}\right) \sin \left(\frac{\pi y (2 r+s)}{A}\right)-\cos \left(\frac{\sqrt{3} \pi p x}{A}\right) \sin \left(\frac{\sqrt{3} \pi s x}{A}\right) \sin \left(\frac{\pi y (p+2 q)}{A}\right) \sin \left(\frac{\pi y (2 r+s)}{A}\right)+\sin \left(\frac{\sqrt{3} \pi r x}{A}\right) \cos \left(\frac{\sqrt{3} \pi x (p+q)}{A}\right) \sin \left(\frac{\pi y (p-q)}{A}\right) \sin \left(\frac{\pi y (r+2 s)}{A}\right)-\cos \left(\frac{\sqrt{3} \pi q x}{A}\right) \sin \left(\frac{\sqrt{3} \pi r x}{A}\right) \sin \left(\frac{\pi y (2 p+q)}{A}\right) \sin \left(\frac{\pi y (r+2 s)}{A}\right)+\cos \left(\frac{\sqrt{3} \pi p x}{A}\right) \sin \left(\frac{\sqrt{3} \pi r x}{A}\right) \sin \left(\frac{\pi y (p+2 q)}{A}\right) \sin \left(\frac{\pi y (r+2 s)}{A}\right)$

$x$ and $y$ are the only variables - everything else is a constant. Fortunately, as pointed out here: https://math.stackexchange.com/a/2425064/441529 It can be simplified by repeated application of the product-to-sum formula. Clearly this will be a very arduous task by hand, but a breeze for an appropriately programmed computer.

I have been trying to get Mathematica to help me, with some luck and also some issues. I have been applying the following rules:

CS = Cos[a_]*Sin[b_] :> (1/2)*(Sin[a + b] - Sin[a - b]);
SC = Sin[a_]*Cos[b_] :> (1/2)*(Sin[a + b] + Sin[a - b]);
CC = Cos[a_]*Cos[b_] :> (1/2)*(Cos[a + b] + Cos[a - b]);
SS = Sin[a_]*Sin[b_] :> (1/2)*(Cos[a - b] - Cos[a + b]);


Which work well to expand the trig functions to get someting like this:

$-\frac{1}{8} \sin \left(\frac{\sqrt{3} \pi x (p+q)}{A}-\frac{\pi y (p-q)}{A}-\frac{\sqrt{3} \pi x (r+s)}{A}-\frac{\pi y (r-s)}{A}\right)+\frac{1}{8} \sin \left(\frac{\sqrt{3} \pi x (p+q)}{A}-\frac{\pi y (p-q)}{A}+\frac{\sqrt{3} \pi x (r+s)}{A}-\frac{\pi y (r-s)}{A}\right)+\frac{1}{8} \sin \left(\frac{\sqrt{3} \pi x (p+q)}{A}+\frac{\pi y (p-q)}{A}-\frac{\sqrt{3} \pi x (r+s)}{A}-\frac{\pi y (r-s)}{A}\right)-\frac{1}{8} \sin \left(\frac{\sqrt{3} \pi x (p+q)}{A}+\frac{\pi y (p-q)}{A}+\frac{\sqrt{3} \pi x (r+s)}{A}-\frac{\pi y (r-s)}{A}\right)+\frac{1}{8} \sin \left(-\frac{\pi y (2 p+q)}{A}+\frac{\sqrt{3} \pi q x}{A}-\frac{\sqrt{3} \pi x (r+s)}{A}-\frac{\pi y (r-s)}{A}\right)-\frac{1}{8} \sin \left(-\frac{\pi y (2 p+q)}{A}+\frac{\sqrt{3} \pi q x}{A}+\frac{\sqrt{3} \pi x (r+s)}{A}-\frac{\pi y (r-s)}{A}\right)-\frac{1}{8} \sin \left(\frac{\pi y (2 p+q)}{A}+\frac{\sqrt{3} \pi q x}{A}-\frac{\sqrt{3} \pi x (r+s)}{A}-\frac{\pi y (r-s)}{A}\right)+\frac{1}{8} \sin \left(\frac{\pi y (2 p+q)}{A}+\frac{\sqrt{3} \pi q x}{A}+\frac{\sqrt{3} \pi x (r+s)}{A}-\frac{\pi y (r-s)}{A}\right)\text{...}$

(The full thing is >10 times as long as this snippet! All terms have the same general form: - $\frac{1}{8} \sin (C x+D y)$)

This should be easily integrated now, but Mathematica's Integrate[] function won't solve the above for me, despite having no trouble with this equivalent integral:

Integrate[Sin[k*a + j*b], {a, a1, a2}, {b, b1, b2}]
= (-Sin[b1 j + a1 k] + Sin[b2 j + a1 k] + Sin[b1 j + a2 k] - Sin[b2 j + a2 k])/(j k)


So my questions is either:

How can I coax Integrate[] into providing a solution?

or

How can I write a rule to apply to the above to produce the result of the integral?

I have tried both the above with no luck so far. Thank you!

Edit: Here's the Mathematica code if you would like to run it: There are three types of wavefunction:

A1[x_, y_, p_, q_] := Cos[q*Sqrt[3]*Pi*x/A]*Sin[(2*p + q)*Pi*y/A] -    Cos[p*Sqrt[3]*Pi*x/A]*Sin[(2*q + p)*Pi*y/A] - Cos[(p + q)*Sqrt[3]*Pi*x/A]*Sin[(p - q)*Pi*y/A];
A2[x_, y_, p1_, q1_] :=   Sin[q1*Sqrt[3]*Pi*x/A]*Sin[(2*p1 + q1)*Pi*y/A] -   Sin[p1*Sqrt[3]*Pi*x/A]*Sin[(2*q1 + p1)*Pi*y/A] + Sin[(p1 + q1)*Sqrt[3]*Pi*x/A]*Sin[(p1 - q1)*Pi*y/A];
EE[x_, y_, p_, q_] := A2[x, y, p, q] + A1[x, y, p, q]*Sqrt[-1];


and the integrands are products of 2 such wavefunctions - e.g.

A1[x, y, p, q]*A2[x, y, r, s]


I then proceed to Expand[] and apply the above rules to the functions. They rapidly get far too long to copy/paste here.

(I start with A1[x, y, p, q]*A2[x, y, r, s]. I first use Expand[] once, then (the solution) /.CS, then /.SS, then Expand[] again and /.CS one more time, to arrive at the snippet above.)

EDIT 2:

I think I may have cracked it.

COL = Sin[a_] :> Sin[Collect[Collect[a, x], y]]
IN = Sin[k_ x + j_ y] :> (-Sin[b1 j + a1 k] + Sin[b2 j + a1 k] + Sin[b1 j + a2 k] - Sin[b2 j + a2 k])/(j k)


However the resulting formula is enormous. So large it caused my computer to crash when I attempted to simplify it!

• Could you edit your question to include code for the integrand that we can copy/paste and run? Moreover, can you post the exact code you used to expand the trig functions and do the integration? You don't need to include the full output, but a sample (like you included, but it would be better in code form) would be nice. Also, what happens with Integrate? Does it run without stopping, or does it give up and return unevaluated? Sep 11, 2017 at 15:42
• In addition to @jjc385's comment, you could use TrigReduce instead of manually coding the rules Sep 11, 2017 at 16:10
• @Mathe172 I tried this - but it's actually easier and faster to do it semi-manually (I.e. choose what rule to apply and when) to guide mathematica towards the form I want
– Ben
Sep 12, 2017 at 3:21
• @jjc385 I've just added the input code. I start with A1[x, y, p, q]*A2[x, y, r, s]. I first use Expand[] once, then (the solution) /.CS, then /.SS, then Expand[] again and /.CS one more time, to arrive at the snippet above.
– Ben
Sep 12, 2017 at 3:34
• For the integral I have tried everything from Integrate[blah[x, y], {x, 0, 1}, {y, 0, 1}] to Assuming[{x > 0, y > 0, x \[Element] Reals, y \[Element] Reals, a > 0, a \[Element] Reals, A \[Element] Reals, A > 0, p > 0, q > 0, p \[Element] Integers, q \[Element] Integers, p1 > 0, q1 > 0, p1 \[Element] Integers, q1 \[Element] Integers}, Integrate[ A1[x, y, p, q]*A2[x, y, p1, q1], {x, y} \[Element] SSSTriangle[a, a, a]]]
– Ben
Sep 12, 2017 at 3:37

I don't see any problems here. Your integrand is a sum of terms.

Every term can be integrated for its own easily. Here done for the first term.

Integrate[-Cos[(Sqrt[3] \[Pi] (p + q) x)/A] Sin[(
Sqrt[3] \[Pi] (r + s) x)/A] Sin[(\[Pi] (p - q) y)/
A] Sin[(\[Pi] (r - s) y)/A], {x, x1, x2}, {y, y1, y2}]

(*     -(A^2 (-(Cos[(Sqrt[3] \[Pi] (p + q - r - s) x1)/A]/(p + q - r - s)) +
Cos[(Sqrt[3] \[Pi] (p + q + r + s) x1)/A]/(p + q + r + s) +
Cos[(Sqrt[3] \[Pi] (p + q - r - s) x2)/A]/(p + q - r - s) -
Cos[(Sqrt[3] \[Pi] (p + q + r + s) x2)/A]/(
p + q + r +
s)) ((p - q - r + s) Sin[(\[Pi] (p - q + r - s) y1)/
A] + (-p + q - r + s) Sin[(\[Pi] (p - q - r + s) y1)/
A] + (-p + q + r - s) Sin[(\[Pi] (p - q + r - s) y2)/
A] + (p - q + r - s) Sin[(\[Pi] (p - q - r + s) y2)/
A]))/(4 Sqrt[3] \[Pi]^2 (p - q + r - s) (p - q - r + s))     *)


Now make a list of terms of the integrand (here done for the first three terms) and integrate each listelement alone and then sum up.

integrand = -Cos[(Sqrt[3] \[Pi] (p + q) x)/A] Sin[(
Sqrt[3] \[Pi] (r + s) x)/A] Sin[(\[Pi] (p - q) y)/
A] Sin[(\[Pi] (r - s) y)/A] +
Cos[(Sqrt[3] \[Pi] q x)/A] Sin[( \[Pi] (2 p + q) y)/A] Sin[(
Sqrt[3] \[Pi] (r + s) x)/A] Sin[(\[Pi] (r - s) y)/A] -
Cos[(Sqrt[3] \[Pi] p x)/A] Sin[( \[Pi] (p + 2 q) y)/A] Sin[(
Sqrt[3] \[Pi] (r + s) x)/A] Sin[(\[Pi] (r - s) y)/A];

Plus @@ (Integrate[#, {x, x1, x2}, {y, y1, y2}] & /@
List @@ integrand)

(* result too long to show here  *)


Appendix:

Since you told, that calculation takes very long

and noticing, that the terms of the integrand are all of the same form

Cos[a x]*Sin[b x]*Sin[c y]*Sin[d y]

there is a much faster way, applying a rule instead of integrating each time again.

intrule =
Cos[a_ x]*Sin[b_ x]*Sin[c_ y]*Sin[d_ y] ->
Integrate[
Cos[a x]*Sin[b x]*Sin[c y]*Sin[d y], {x, x1, x2}, {y, y1, y2}]

(*     Cos[x a_] Sin[x b_] Sin[y c_] Sin[
y d_] -> ((-b Cos[a x1] Cos[b x1] + b Cos[a x2] Cos[b x2] -
a Sin[a x1] Sin[b x1] +
a Sin[a x2] Sin[b x2]) (-d Cos[d y1] Sin[c y1] +
c Cos[c y1] Sin[d y1] + d Cos[d y2] Sin[c y2] -
c Cos[c y2] Sin[d y2]))/((a^2 - b^2) (c^2 - d^2))     *)


Test it with the first term of the integrand

i1 = Integrate[term1, {x, x1, x2}, {y, y1, y2}];

t1 = term1 /. intrule;

t1 - i1 // TrigExpand // FullSimplify

(*     0     *)


The code, now by orders of magnitude faster, is now

Plus @@ ((# /. intrule) & /@ List @@ integrand)

• It works! (But it is incredibly slow and caused my computer to crash twice! probably the laptop's fault more than anything else). Actually, the most useful (and faster) result was when I applied it to the fully expanded function, and changed the integral to over a triangle using Plus @@ (Integrate[#, {x, y} \[Element] SSSTriangle[a, a, a]] & /@ List @@ integrand)
– Ben
Sep 12, 2017 at 15:01
• I just appended an improvement to my former answer. Sep 13, 2017 at 12:53