# Help at speeding up simplification/nested symbolic integration

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)$.

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 *)

• Thanks alot for your contribution! As I am reading on a mobile right now, it is hard to look deeper into your rules. Your time result looks really promising and I will look into details tomorrow and maybe ask one or two questions about your suggestions. – Lukas Dec 8 '14 at 15:58
• The speed improvement is absolutely brilliant! It only takes 9.5sec on my machine. Nevertheless I have a small issue with a simple test. testfun[t1_,t2_]=Sin[t1]+Cos[t2] fails with your method. To me, the problem here seems to be with the Count part - e.g. Count[Sin[t1],Sin[_],Infinity] returns 0 whereas Count[Sin[t1]+Cos[t2],Sin[_],Infinity] returns 1 as I expected. Why does it give 0 in the first case? However the method works fine with comm00 which really confuses me. Could you or someone else maybe elaborate why I fail at the simple test? – Lukas Dec 9 '14 at 10:45
• Here the result of testfun[t1,t2]/.trigInt[t1,0,t2]: t2 Cos[t2] + Integrate[Sin[t1],{t1,0,t2}]. And here tracing trigCheck[testfun[t1,t2]]: tinyurl.com/trace-trigCheck To me, there seems to be an issue with the Count and Exponent parts. – Lukas Dec 9 '14 at 11:05
• @Lukas Cases/Count start at level 1. To match a simple Sin[t1], we need to tell it to start at level 0. Fixed (and made slightly more efficient); see update. – Michael E2 Dec 9 '14 at 13:47
• Thanks for hinting at the levelspec issue. But I found another issue with your code: During the integration there are terms like t2*Sin[t2] that are not properly covered by the rules as they only lead to replacing the Sin. This can be shown by FreeQ[int,t2] with definitions above (10/12/14,8:42AM (MEZ)) which evaluates to False. I fixed it by adding another expand and two more rules as can be seen in my second edit. However as you did most of the job I'd like to accept your answer in case you add the specific lines. Thanks alot again! – Lukas Dec 10 '14 at 7:48

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.

• Thanks alot for your suggestion(s). I just implemented improvements and started the script. I will give feedback about the result! Edit: Sorry, I cannot upvote your answer due to a lack of reputation... – Lukas Dec 3 '14 at 15:23
• Update: So after running the code for a while it seems that it still takes a lot of time. Tomorrow I will start the calculation on another computer and see how long it really takes to do the calculation. Is there anything else (apart from Daniel Lichtblau's answer) I could do to increase speed? – Lukas Dec 3 '14 at 16:25
• Sorry for adding comments over comments, but I was just thinking about using ParallelMap instead of Map. The only thing I am wondering about is how the "mapped integration" is performed, i.e. do I have problems with parallelisation or not? Personally, I do not see any problems with e.g. side effects, but unfortunately I am not a master of parallelisation... Any comments on that? – Lukas Dec 4 '14 at 6:44
• I cannot claim to be terribly familiar with ParallelMap but offhand I don't see why it wouldn't just work. Maybe experiment on a smallish example, say the first 10 summands, and print the kernel of evaluation before doing the integral e.g.  ParallelMap[Print[$KernedID];Integrate[...]]. That should give an indication of whether it is working to expectations. – Daniel Lichtblau Dec 4 '14 at 15:10 • Printing $KernelID only gives the output 0. As a test function I used x^2 - x + 4*x^5 + Sin[x^2]*x - Sin[x^2]. Unfortunately it seems that running the job in parallel does not work (at least not that easily). Moreover my computer crashed, so I had to start the computation again to see how long your approach actually takes. – Lukas Dec 5 '14 at 9:36