Help at speeding up simplification/nested symbolic integration

Question

I have the following problem: I want to compute an integral of the form $$\int\limits_0^{t_g}\mathrm{d}t_2\int\limits_0^{t_2}\mathrm{d}t_1\left[A(t_1),A(t_2)\right],$$where $\left[B,C\right]=BC-CB$ is the commutator of two matrices (in my case both 3x3). All entries of the resulting commutator have a similiar form, so let us just focus on the first diagonal entry, which reads $$\left(\Omega_1(t_2)+e^{i\gamma t_2}\Omega_2(t_2)\right)\left(\Omega_1^*(t_1)+e^{-i\gamma t_1}\Omega_2^*(t_1)\right)-\text{h.c.}\quad (1)$$where h.c. means the hermitian conjugate. The functions $\Omega_i(t_j)$ are complex functions that are basically composed by sums and products of Sine/Cosine with different frequencies (please find the explicit definitions in the code block below).

Obviously Eq. (1) is just $2i \;{\rm Im[\ldots]}$ where $\ldots$ denotes the term in front of "h.c.". So I wanted to use this observation as a starting point for following symbolic computation. Nevertheless, extracting the imaginary part by using ComplexExpand[Im[...]] yields a term consisting of approximately 5100 summands that could "easily" be simplified by factoring out some stuff. Using Simplify for that job seems to take forever on my machine.

Even worse is the nested integration which I would like to compute symbolically in order to see if some conditions can be obtained from the result (results should be stored as a function to be used for further calculation).

So my questions are:

1. Is there a way to speed up simplification of the ComplexExpand part?
2. Is there a more elegant/faster way of performing the nested integration, apart from Integrate[Integrate[...],{t1,0,t2}],{t2,0,tg}]

And here is the code defining my functions

ex1Temp[t_, tg_] := (1 - Cos[(2*π*t)/tg])*(1 -A1*Cos[ωx1*(t - tg/2)]);
ex1[t_, tg_] :=ex1Temp[t, tg] + b21*D[ex1Temp[x, tg], {x, 2}] /. {x -> t};
ey1[t_, tg_] := (b11*D[ex1Temp[x, tg], {x, 1}] + b31*D[ex1Temp[x, tg], {x, 3}]) /. {x -> t};
O1[t_, tg_] := ex1[t, tg] + I*ey1[t, tg];
O1c[t_, tg_] := ex1[t, tg] - I*ey1[t, tg];

ex2Temp[t_,tg_] := (1 - Cos[(2*π*t)/tg])*(1 - A2*Cos[ωx2*(t - tg/2)]);
ex2[t_, tg_] := ex2Temp[t, tg] + b22*D[ex2Temp[x, tg], {x, 2}] /. {x -> t};
ey2[t_, tg_] := (b12*D[ex2Temp[x, tg], {x, 1}] + b32*D[ex2Temp[x, tg], {x, 3}]) /. {x -> t};
O2[t_, tg_] := ex2[t, tg] + I*ey2[t, tg];
O2c[t_, tg_] := ex2[t, tg] - I*ey2[t, tg];

comm00[t1_, t2_, tg_] = ComplexExpand[Im[(O1[t2, tg] + Exp[I*γ*t2]*O2[t2, tg])*(O1c[t1, tg] + Exp[-I*γ*t1]*O2c[t1, tg])]];

fun[tg_]=Integrate[comm00[t1,t2,tg],{t1,0,t2}],{t2,0,tg}]

Notice that $\Omega$ from my latex equations at the top are represented by O in the Mathematica code. Any help and/or tipps are well appreciated.

Edit 1 (based on Daniel Lichtblau's recommendations)

A good point was to apply rules to the different summands of the appearing integrals. I am working on this but need some help with generalization of the results. So for this Edit, please let us just focus on the following code snippet:

intRule = Integrate[a_ + b_, c_] :> Integrate[a, c] + Integrate[b, c];
trafo = Integrate[Cos[w_*t_]*ex1[t_, tg_], {t_, 0,t2_}] :> (1 - w^2*b21)*Integrate[Cos[w*t]*ex1Temp[t, tg], {t,0, t2}];

intRule is used to transform the integral over a sum to a sum of integrals. One of these summands is $\cos(\omega t)\epsilon_{x1}(t,t_g)$ which appears in my rule trafo. I am now facing a problem/uncertainty regarding efficiently applying the transformation rule(s). Inside trafo there is another integral which again be easily computed analytically, say we name the result F1[t2]. What is now a better way to "replace" that integral in rule trafo?

1. Define the function F1[t2] with Set (not SetDelayed) and change trafo to trafo=Integrate[Cos[w_*t_]*ex1[t_, tg_], {t_, 0,t2_}] :> (1 - w^2*b21)*F1[t2], or
2. Apply another rule in a second step that acts in a similar way as $#1$, but maybe performs better

Edit 2 (based on Michael E2's answer)

As there is a small issue in MichaelE2's answer (as of 8:49AM (MEZ), 10/12/14) with terms that go like t2*Sin[t2] the rules need to be changed to something like

trigInt[v_, a_, b_] := {
   integrand_ /;FreeQ[integrand, v] :> (b - a) integrand,(*constant rule*)
   integrand_Plus :> Replace[integrand, trigInt[v, a, b], 1],(*addition rule*)   
   integrand_ /; trigCheck[integrand] :>(*sin/cos rules*)
    Expand[Expand[integrand] /. {
        Sin[arg_] :> 1/Coefficient[arg, v]*First@Differences[-Cos[arg] /. {{v -> a}, {v -> b}}],
        Cos[arg_] :> 1/Coefficient[arg, v]*First@Differences[Sin[arg] /. {{v -> a}, {v -> b}}],
       v*Sin[arg_] :> 1/Coefficient[arg, v]*(First@Differences[-v*Cos[arg] +1/Coefficient[arg, v]*Sin[arg] /. {{v -> a}, {v -> b}}]),
       v*Cos[arg_] :> 1/Coefficient[arg, v]*(First@Differences[v*Sin[arg] +1/Coefficient[arg, v]*Cos[arg] /. {{v -> a}, {v -> b}}])}],
    integrand_ :> Inactive[Integrate][integrand, {v, a, b}]};(*failure rule*)

To generalize this one can think of using Gamma function to build up rules for $\int\mathrm{d}x\,x^n\sin(x)$.

Answer 1

Edit 2 - Fixed omitted integrand

Lukas pointed out that the first integration produces terms of the form x Sin[x], which were mishandled by the original rules. (See edit history, if curious.) I changed some things around a little. We have to Expand the result of the first integration before doing the second integration. Overall the speed is actually improved by a few seconds (was ~13, now ~9).

Here is a set of rules that handles a linear combination of sine/cosine functions and functions of the form x Sin[a x + b] and x Cos[a x + b], where a, b are constants and x is a variable. TrigReduce will convert the integrand comm00[t1, t2, tg] into such a form. The integrand is so large that I wanted to do some checking ahead of time. The function validQ[coeff, arg, v] checks that a term has a coefficient is free of the variable of integration v and that the argument of the trigonometric function arg is a linear polynomial in v.

validQ[coeff___, arg_, v_] := FreeQ[{coeff}, v] && Exponent[arg, v] == 1;
trigInt[v_, a_, b_] := {
   integrand_ /; FreeQ[integrand, v] :> (b - a) integrand,       (*constant rule*)
   integrand_Plus :> Replace[integrand, trigInt[v, a, b], 1],    (*addition rule*)
   coeff___ *v*Sin[arg_] /; validQ[coeff, arg, v] :>             (*v sin v rules*)
      coeff/Coefficient[arg, v] * (First @
      Differences[-v*Cos[arg] + 1/Coefficient[arg, v]*Sin[arg] /. {{v -> a}, {v -> b}}]),
  coeff___ *v*Cos[arg_] /; validQ[coeff, arg, v] :> 
   coeff/Coefficient[arg, v] * (First @
      Differences[ v*Sin[arg] + 1/Coefficient[arg, v]*Cos[arg] /. {{v -> a}, {v -> b}}]),
  coeff___ *Sin[arg_] /; validQ[coeff, arg, v] :>                (*sin/cos rules*)
   coeff/Coefficient[arg, v] * (First @ Differences[-Cos[arg] /. {{v -> a}, {v -> b}}]),
  coeff___ *Cos[arg_] /; validQ[coeff, arg, v] :> 
   coeff/Coefficient[arg, v] * (First @ Differences[Sin[arg] /. {{v -> a}, {v -> b}}]),
  integrand_ :> Inactive[Integrate][integrand, {v, a, b}]};      (*failure rule*)

TrigReduce@comm00[t1, t2, tg] itself about 5.6 sec. and produces 8443 terms, but applying the rules is rather fast relative to Integrate.

int3 = Expand[
     TrigReduce@comm00[t1, t2, tg] /. trigInt[t1, 0, t2]
     ] /. trigInt[t2, 0, tg]; // AbsoluteTiming
(* {9.093579, Null} *)

Check that all terms were integrated:

FreeQ[int, t1 | t2]
(* True *)

Check that it agrees with Lukas' version of these rules:

(* re-execute Lukas's definition of trigInt *)
intL = TrigReduce@comm00[t1,t2,tg] /. trigInt[t1, 0, t2] /. trigInt[t2, 0, tg]; // AbsoluteTiming
Expand[intL - int]
(*  {15.393432, Null}  *)
(*  0  *)

The integral has a lot of terms. Simplifying it might take as long as Integrate. I do not have a good idea how verify the solution further.

Length[int]
Length @ Expand[int]
(* 12898, 197904 *)

Answer 2

You might get a speed improvement by doing as follows.

(1) Change the expansion to give an explicit sum of products of trigs.

comm00[t1_, t2_, tg_] = 
 ComplexExpand[
  Im[Expand[
    ExpToTrig[(O1[t2, tg] + 
        Exp[I*γ*t2]*O2[t2, tg])*(O1c[t1, tg] + 
        Exp[-I*γ*t1]*O2c[t1, tg])]]]]

(2) Do the double integral without iteration, and give (plausible) assumptions on the upper bounds.

(3) Maybe integrate termwise.

--- edit ---

Putting all this together gives something like the following.

Map[Integrate[#, {t2, 0, tg}, {t1, 0, t2}, 
   Assumptions -> {tg > 0, 
     Element[{ωx1, ωx2, A1, A2, b11, b12, b21, b22, b31,
        b32, γ}, Reals]}] &, comm00[t1, t2, tg]]

To see where bottlenecks lie you can break up the integrand, print intermediate information, maybe write special-case definite integral code (I think you started in that direction; you'd want integral formulas for products such as Cos[(2 π t2)/tg] Cos[t1 γ] Sin[(2 π t1)/tg] or Cos[t1 ωx1-(tg ωx1)/2] Sin[(2 π t1)/tg]). Here is something specific.

c1 = comm00[t1, t2, tg]; (* A sum of 1550 terms *)

j = 0;
Map[(Print[j]; Print[#]; j++; 
   res[j] = 
    Integrate[#, {t2, 0, tg}, {t1, 0, t2}, 
     Assumptions -> {tg > 0, 
       Element[{ωx1, ωx2, A1, A2, b11, b12, b21, b22, 
         b31, b32, γ}, Reals]}]; Print[res[j]]; res[j]) &, c1]

It's taken a few of minutes to process 30 summands but I suspect some later ones get much slower (if not, it should complete in around 6 hours or so). At least the print statements will indicate which types are slow. Possibly those could be special-cased.