-1

I am trying to quickly solve the following problem:

f[r_] := Sum[(((-1)^n (2 r - 2 n - 7)!!)/(2^n n! (r - 2 n - 1)!))
             * x^(r - 2*n - 1), 
         {n, 0, r/2}]; 

Nw := Transpose[Table[f[j], {i, 1}, {j, 5, 200, 1}]]; 

X1 = Integrate[Nw . Transpose[Nw], {x, -1, 1}]

I can get the answer quickly with this code:

$starttime = AbsoluteTime[]; Quiet[LaunchKernels[]]; 
DIM = 50; 
Print["$Version = ", $Version, "  |||  ", 
      "Number of Kernels : ", Length[Kernels[]]]; 

Nw = Transpose[Table[f[j], {i, 1}, {j, 5, DIM, 1}]]; 
nw2 = Nw.Transpose[Nw]; 
Round[First[AbsoluteTiming[nw3 = ParallelMap[Expand, nw2]; ]]] 

intrule = (pol_Plus)?(PolynomialQ[#1, x]&) :> 
      (Select[pol, !FreeQ[#1, x] & ] /. 
         x^(n_.) /; n > -1 :> ((-1)^n + 1)/(n + 1)) + 2*(pol /. x -> 0)]); 

Round[First[AbsoluteTiming[X1 = ParallelTable[row /. intrule, {row, nw3}]; ]]]

X1 
Print["overall time needed in seconds: ", Round[AbsoluteTime[] - $starttime]];

But how can I manage this code if I need to solve the following problem, where a and b are known constants?

       X1 = a Integrate[Nw.Transpose[Nw], {x, -1, 0.235}]
          + b Integrate[Nw.Transpose[Nw], {x,  0.235,1}];
rcollyer
  • 10,475
  • 4
  • 48
  • 75
George Mills
  • 159
  • 6
  • 1
    This question is not particularly clear and is very poorly formatted, which makes it hard to read and thus discourages people from helping you. Can [you](http://stackoverflow.com/users/1031298) please try to fix it (this includes the title and text around the code)? -- Note that good questions encourage good answers and then everyone benefits! – Simon Dec 02 '11 at 11:00
  • Related questions: [SO/8021501](http://stackoverflow.com/q/8021501) & [SU/315337](http://superuser.com/questions/315337) – Simon Dec 02 '11 at 11:02
  • What is not clear Simon? I need firstly X1 = Integrate[Nw . Transpose[Nw], {x, -1, 1}]; Now I need X1 = a*Integrate[Nw . Transpose[Nw], {x, -1, 0.235}]+b*Integrate[Nw .Transpose[Nw], {x, 0.235,1}]; So try to obtain this two integrals with presented code, this is all – George Mills Dec 02 '11 at 11:13
  • It's understandable, just not clear. – Simon Dec 02 '11 at 11:22
  • 4
    Welcome back. I have a couple of questions. Would you mind if I did some formatting, including removing the `Print` statements? I'd like to include just the core code as it helps to see what the code is doing. Secondly, in the definition of `Nw` you use `:=` which re-evaluates every time `Nw` is called. Is this intentional, and why? Third, the `Table` in `Nw` only generates a single row (`i` is always 1), so do you have cases where there is more than one value for `i`? – rcollyer Dec 02 '11 at 13:57
  • @rcollyer: I say have at the reformatting. – Brett Champion Dec 02 '11 at 16:13
  • @BrettChampion, I plan on it. Sometime this afternoon, though. – rcollyer Dec 02 '11 at 17:25

1 Answers1

6

Here's a simple function to do definite integrals of polynomials

polyIntegrate[expr_List, {x_, x0_, x1_}] := polyIntegrate[#, {x, x0, x1}]&/@expr
polyIntegrate[expr_, {x_, x0_, x1_}] := Check[Total[# 
  Table[(x1^(1 + n) - x0^(1 + n))/(1 + n), {n, 0, Length[#] - 1}]
  ]&[CoefficientList[expr, x]], $Failed, {General::poly}]

On its range of applicability, this is about 100 times faster than using Integrate. This should be fast enough for your problem. If not, then it could be parallelized.

f[r_] := Sum[(((-1)^n*(2*r - 2*n - 7)!!)/(2^n*n!*(r - 2*n - 1)!))*
   x^(r - 2*n - 1), {n, 0, r/2}];
Nw = Transpose[Table[f[j], {i, 1}, {j, 5, 50, 1}]];

a*polyIntegrate[Nw.Transpose[Nw], {x, -1, 0.235}] + 
   b*polyIntegrate[Nw.Transpose[Nw], {x, 0.235, 1}] // Timing // Short

(* Returns: {7.9405,{{0.0097638 a+0.00293462 b,<<44>>,
             -0.0000123978 a+0.0000123978 b},<<44>>,{<<1>>}}} *)
Simon
  • 14,631
  • 4
  • 41
  • 101
  • I was about to admonish you for not putting in sufficient protections into `polyIntegrate` to deal with `expr == constant` and `expr == 1/x`, but apparently `polyIntegrate` handles the first one just fine and the second one never shows up, despite appearing like it would. btw, +1. – rcollyer Dec 02 '11 at 13:49
  • Thank you dear Simon, for any case how it can be parallelized? – George Mills Dec 02 '11 at 15:48
  • @Simon I just used old parallelize code, but improve it with poly in front and it is working very fast. Thank you very much! – George Mills Dec 02 '11 at 16:04
  • @rcollyer: `CoefficientList` will spit out an error if `expr` is not a polynomial in `x`. I've added a `Check` to this now, so that `polyIntegrate` will return `$Failed` when given a nonpolynomial. (I was too lazy to add this yesterday!) – Simon Dec 02 '11 at 23:07
  • 1
    @Simon, I noticed ... :P Alternatively, you could use `PolynomialQ[#,x]&`. – rcollyer Dec 03 '11 at 01:42