Skip to main content

calculus and analysis - How can I Integrate large trigonometric expressions more efficiently?


I have a requirement to integrate some large symbolic expressions containing products of trigonometric functions. A single integration is averaging 7 minutes, and I have 2628 expressions to integrate, so it will take 13 days to complete unless I can speed it up.


Here is an example:


expr = c14 (a0+a1 Cos[n x]+a3 Cos[2 n x]+a5 Cos[3 n x]+a7 Cos[4 n x]+a2 Sin[n x]+a4 Sin[2 n x]+a6 Sin[3 n x]+a8 Sin[4 n x]) (e1 Cos[n x]+e3 Cos[2 n x]+e5 Cos[3 n x]+e7 Cos[4 n x]+e2 Sin[n x]+e4 Sin[2 n x]+e6 Sin[3 n x]+e8 Sin[4 n x])+c44 (e1 Cos[n x]+e3 Cos[2 n x]+e5 Cos[3 n x]+e7 Cos[4 n x]+e2 Sin[n x]+e4 Sin[2 n x]+e6 Sin[3 n x]+e8 Sin[4 n x])^2+c66 (j0+j1 Cos[n x]+j3 Cos[2 n x]+j5 Cos[3 n x]+j7 Cos[4 n x]+j2 Sin[n x]+j4 Sin[2 n x]+j6 Sin[3 n x]+j8 Sin[4 n x])^2+c11 (1/2 (h0+h1 Cos[n x]+h3 Cos[2 n x]+h2 Sin[n x]+h4 Sin[2 n x])^2+3/2 (a0+a1 Cos[n x]+a3 Cos[2 n x]+a5 Cos[3 n x]+a7 Cos[4 n x]+a2 Sin[n x]+a4 Sin[2 n x]+a6 Sin[3 n x]+a8 Sin[4 n x])^2+1/2 (d0+d1 Cos[n x]+d3 Cos[2 n x]+d5 Cos[3 n x]+d7 Cos[4 n x]+d2 Sin[n x]+d4 Sin[2 n x]+d6 Sin[3 n x]+d8 Sin[4 n x])^2)+c14 ((b0+b1 Cos[n x]+b3 Cos[2 n x]+b5 Cos[3 n x]+b7 Cos[4 n x]+b2 Sin[n x]+b4 Sin[2 n x]+b6 Sin[3 n x]+b8 Sin[4 n x]) (d0+d1 Cos[n x]+d3 Cos[2 n x]+d5 Cos[3 n x]+d7 Cos[4 n x]+d2 Sin[n x]+d4 Sin[2 n x]+d6 Sin[3 n x]+d8 Sin[4 n x])+2 (a0+a1 Cos[n x]+a3 Cos[2 n x]+a5 Cos[3 n x]+a7 Cos[4 n x]+a2 Sin[n x]+a4 Sin[2 n x]+a6 Sin[3 n x]+a8 Sin[4 n x]) (e1 Cos[n x]+e3 Cos[2 n x]+e5 Cos[3 n x]+e7 Cos[4 n x]+e2 Sin[n x]+e4 Sin[2 n x]+e6 Sin[3 n x]+e8 Sin[4 n x])+(h0+h1 Cos[n x]+h3 Cos[2 n x]+h2 Sin[n x]+h4 Sin[2 n x]) (g0+g1 Cos[n x]+g3 Cos[2 n x]+g5 Cos[3 n x]+g7 Cos[4 n x]+g2 Sin[n x]+g4 Sin[2 n x]+g6 Sin[3 n x]+g8 Sin[4 n x]))+c12 (1/2 (b0+b1 Cos[n x]+b3 Cos[2 n x]+b5 Cos[3 n x]+b7 Cos[4 n x]+b2 Sin[n x]+b4 Sin[2 n x]+b6 Sin[3 n x]+b8 Sin[4 n x])^2+1/2 (e1 Cos[n x]+e3 Cos[2 n x]+e5 Cos[3 n x]+e7 Cos[4 n x]+e2 Sin[n x]+e4 Sin[2 n x]+e6 Sin[3 n x]+e8 Sin[4 n x])^2+1/2 (g0+g1 Cos[n x]+g3 Cos[2 n x]+g5 Cos[3 n x]+g7 Cos[4 n x]+g2 Sin[n x]+g4 Sin[2 n x]+g6 Sin[3 n x]+g8 Sin[4 n x])^2)+c13 (1/2 (c0+c1 Cos[n x]+c3 Cos[2 n x]+c2 Sin[n x]+c4 Sin[2 n x])^2+1/2 (f0+f1 Cos[n x]+f3 Cos[2 n x]+f5 Cos[3 n x]+f7 Cos[4 n x]+f2 Sin[n x]+f4 Sin[2 n x]+f6 Sin[3 n x]+f8 Sin[4 n x])^2+1/2 (j0+j1 Cos[n x]+j3 Cos[2 n x]+j5 Cos[3 n x]+j7 Cos[4 n x]+j2 Sin[n x]+j4 Sin[2 n x]+j6 Sin[3 n x]+j8 Sin[4 n x])^2);


Integrate[expr, {x, 0, 2 Pi}]

n is an integer, if that helps.


How can I compute these integrals more efficiently?


Take case with this example too:


a12 = Cos[n x] (c14 (a0+a1 Cos[n x]+a3 Cos[2 n x]+a5 Cos[3 n x]+a7 Cos[4 n x]+a2 Sin[n x]+a4 Sin[2 n x]+a6 Sin[3 n x]+a8 Sin[4 n x]) (e1 Cos[n x]+e3 Cos[2 n x]+e5 Cos[3 n x]+e7 Cos[4 n x]+e2 Sin[n x]+e4 Sin[2 n x]+e6 Sin[3 n x]+e8 Sin[4 n x])+c44 (e1 Cos[n x]+e3 Cos[2 n x]+e5 Cos[3 n x]+e7 Cos[4 n x]+e2 Sin[n x]+e4 Sin[2 n x]+e6 Sin[3 n x]+e8 Sin[4 n x])^2+c66 (j0+j1 Cos[n x]+j3 Cos[2 n x]+j5 Cos[3 n x]+j7 Cos[4 n x]+j2 Sin[n x]+j4 Sin[2 n x]+j6 Sin[3 n x]+j8 Sin[4 n x])^2+c11 (1/2 (h0+h1 Cos[n x]+h3 Cos[2 n x]+h2 Sin[n x]+h4 Sin[2 n x])^2+3/2 (a0+a1 Cos[n x]+a3 Cos[2 n x]+a5 Cos[3 n x]+a7 Cos[4 n x]+a2 Sin[n x]+a4 Sin[2 n x]+a6 Sin[3 n x]+a8 Sin[4 n x])^2+1/2 (d0+d1 Cos[n x]+d3 Cos[2 n x]+d5 Cos[3 n x]+d7 Cos[4 n x]+d2 Sin[n x]+d4 Sin[2 n x]+d6 Sin[3 n x]+d8 Sin[4 n x])^2)+c14 ((b0+b1 Cos[n x]+b3 Cos[2 n x]+b5 Cos[3 n x]+b7 Cos[4 n x]+b2 Sin[n x]+b4 Sin[2 n x]+b6 Sin[3 n x]+b8 Sin[4 n x]) (d0+d1 Cos[n x]+d3 Cos[2 n x]+d5 Cos[3 n x]+d7 Cos[4 n x]+d2 Sin[n x]+d4 Sin[2 n x]+d6 Sin[3 n x]+d8 Sin[4 n x])+2 (a0+a1 Cos[n x]+a3 Cos[2 n x]+a5 Cos[3 n x]+a7 Cos[4 n x]+a2 Sin[n x]+a4 Sin[2 n x]+a6 Sin[3 n x]+a8 Sin[4 n x]) (e1 Cos[n x]+e3 Cos[2 n x]+e5 Cos[3 n x]+e7 Cos[4 n x]+e2 Sin[n x]+e4 Sin[2 n x]+e6 Sin[3 n x]+e8 Sin[4 n x])+(h0+h1 Cos[n x]+h3 Cos[2 n x]+h2 Sin[n x]+h4 Sin[2 n x]) (g0+g1 Cos[n x]+g3 Cos[2 n x]+g5 Cos[3 n x]+g7 Cos[4 n x]+g2 Sin[n x]+g4 Sin[2 n x]+g6 Sin[3 n x]+g8 Sin[4 n x]))+c12 (1/2 (b0+b1 Cos[n x]+b3 Cos[2 n x]+b5 Cos[3 n x]+b7 Cos[4 n x]+b2 Sin[n x]+b4 Sin[2 n x]+b6 Sin[3 n x]+b8 Sin[4 n x])^2+1/2 (e1 Cos[n x]+e3 Cos[2 n x]+e5 Cos[3 n x]+e7 Cos[4 n x]+e2 Sin[n x]+e4 Sin[2 n x]+e6 Sin[3 n x]+e8 Sin[4 n x])^2+1/2 (g0+g1 Cos[n x]+g3 Cos[2 n x]+g5 Cos[3 n x]+g7 Cos[4 n x]+g2 Sin[n x]+g4 Sin[2 n x]+g6 Sin[3 n x]+g8 Sin[4 n x])^2)+c13 (1/2 (c0+c1 Cos[n x]+c3 Cos[2 n x]+c2 Sin[n x]+c4 Sin[2 n x])^2+1/2 (f0+f1 Cos[n x]+f3 Cos[2 n x]+f5 Cos[3 n x]+f7 Cos[4 n x]+f2 Sin[n x]+f4 Sin[2 n x]+f6 Sin[3 n x]+f8 Sin[4 n x])^2+1/2 (j0+j1 Cos[n x]+j3 Cos[2 n x]+j5 Cos[3 n x]+j7 Cos[4 n x]+j2 Sin[n x]+j4 Sin[2 n x]+j6 Sin[3 n x]+j8 Sin[4 n x])^2))

Answer



Your integrand(s), when expanded, consist of terms that are constants times zero, one or two factors of sine or cosine. The terms can be integrated symbolically and stored as rules. It will be convenient to find the average value and multiply by 2 Pi at the end. (This will automatically handle constant terms without a special rule.) Thanks to Simon Woods for a much simplified set of rules.


Clear[a, b];

iRules = {(Cos | Sin)[_]^2 -> 1/2, (Cos | Sin)[_] -> 0};

ia11 = 2 Pi Expand[a11] /. iRules; // AbsoluteTiming
(*  {0.010274, Null}  *)

It's a lot faster:


evalRules = {   (* assumes  n  is an  Integer *)
  Cos[(a_Integer: 1) n Pi] :> (-1)^a,
  Sin[(b_Integer: 1) n Pi] :> 0
  };

ja11 = Integrate[a11, {x, 0, 2 Pi}] /. evalRules; // AbsoluteTiming
(*  {177.270803, Null}  *)

And it gives the same result, provided the assumptions inherent in evalRules are correct:


ia11 == Expand[ja11] // Simplify
(*  True  *)


Update: Here's another way, about 20 times slower than iRules, but 1000 times faster than Integrate. It's a similar idea, but it should work on general products of sines and cosines. It assumes that sines and cosines are the only functions in the integrands.


ka11 = 2 Pi Expand[TrigToExp@a11] /. E^_ -> 0; // AbsoluteTiming

(*  {0.205734, Null}  *)

ka11 == Expand[ja11 /. evalRules] // Simplify
(*  True  *)

Or, sparked by a comment by the OP, more complicated combinations of sines and cosines can be fairly quickly simplified with TrigReduce.


2 Pi TrigReduce[a11] /. iRules; // AbsoluteTiming
(*  {0.086483, Null}  *)



Original rules


iRules =(*Dispatch@*)
 {Cos[(a_Integer: 1) n x] Sin[(a_Integer: 1) n x] ->
   Integrate[Cos[a n x] Sin[a n x], {x, 0, 2 Pi}]/(2 Pi),
  Cos[(a_Integer: 1) n x] Sin[(b_Integer: 1) n x] /; a == -b -> 
   Integrate[Cos[a n x] Sin[-a n x], {x, 0, 2 Pi}]/(2 Pi),
  Cos[(a_Integer: 1) n x] Sin[(b_Integer: 1) n x] -> 
   Integrate[Cos[a n x] Sin[b n x], {x, 0, 2 Pi}]/(2 Pi),
  Sin[(a_Integer: 1) n x] Sin[(b_Integer: 1) n x] -> 
   Integrate[Cos[a n x] Sin[b n x], {x, 0, 2 Pi}]/(2 Pi),

  Cos[(a_Integer: 1) n x] Cos[(b_Integer: 1) n x] -> 
   Integrate[Cos[a n x] Sin[b n x], {x, 0, 2 Pi}]/(2 Pi),
  Cos[(a_Integer: 1) n x]^2 -> 
   Integrate[Cos[a n x]^2, {x, 0, 2 Pi}]/(2 Pi),
  Sin[(b_Integer: 1) n x]^2 -> 
   Integrate[Sin[b n x]^2, {x, 0, 2 Pi}]/(2 Pi),
  Cos[(a_Integer: 1) n x] -> 
   Integrate[Cos[a n x], {x, 0, 2 Pi}]/(2 Pi),
  Sin[(b_Integer: 1) n x] -> 
   Integrate[Sin[b n x], {x, 0, 2 Pi}]/(2 Pi)};


Some remarks on the code: Note that the rule in iRules are Rule (->) not RuleDelayed (:>). This means that the integrals will be performed when iRules are defined instead of delayed until the application of iRules. Therefore the pattern symbols a and b must be cleared of any values. If that is not possible, then use other symbols or put in a Module (Module[{a, b}, iRules = {...}]).


Dispatch did not speed things up for me. That may change if the list of integration rules were to grow much larger.


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more