2

I'm integrating a positive function f(t) times sin(t) from 0 to pi/5 and get -38.

Actually f is slightly negative for a short time (smallest value ~ -0.0005), but far from enough to explain this. Here is the code for the function: 

NormFact[k_, S_] := 
Sqrt[(4 S - 2 k - 1) (4 S - k) (4 S - k + 1) (k + 
  1) (k + 2)/(4*Pi*S)];

xi[k_, j_, S_] := (-1)^j*
 Binomial[k, j]*(4*S - k - 1)!/((j + 2)!*(4*S - k - j - 1)!);

Ofunc[t_, S_, coeffs_] := 
1 - Cos[t/2]^(4*S) + 
coeffs.Table[
 NormFact[k, S]*
  Sum[xi[k, j, S]*Sin[t/2]^(2*j + 2)*Cos[t/2]^(4*S - 2*j - 2), {j,
     0, k}], {k, 0, Length[coeffs] - 1}];

cfs = {-2.6155133, 0.9614036, 0.100279, -0.432464, 
0.39887624, -0.2508507, 0.10555472, -0.0007824, -0.0526054, 
0.0703274, -0.0688152, 0.05138949, -0.03833719, 
0.02062684, -0.0018939, -0.0027808};

And here's the result when I integrate:

In:= Integrate[Ofunc[t, 58.5, cfs]*Sin[t], {t, 0, Pi/5}]
Out:= -38.132 + 0. I

Can anyone see my mistake?

rcollyer
  • 33,976
  • 7
  • 92
  • 191
jorgen
  • 539
  • 2
  • 11
  • 2
    On v9.0.1, I get 37.8958. What version are you using? Run "ReleaseID" /. ("Kernel" /. SystemInformation["Small"]) – rcollyer Aug 08 '13 at 16:27
  • Hm! I get 8.0.1.0 (2063990, 2063802). Note: I tried restarting the program and get the same negative answer. – jorgen Aug 08 '13 at 16:30
  • 1
    I confirmed I get that answer on 8.0.1, but I don't get it on v9. So, it is likely a bug in v8. Although, I haven't tried v8.0.4. My only advice: upgrade. – rcollyer Aug 08 '13 at 16:36
  • Thanks for advice! Note: it seems to me the answer you get must also be wrong, since the function and sine are both too small on the interval to give 38 (try plotting it) – jorgen Aug 08 '13 at 16:37
  • 1
    Hmm, I hadn't paid attention to that. NIntegrate gives an answer I'm comfortable with (0.13956), and it works fine on v8.0.1, too. So, I'd use that, especially considering the coefficients in your function have vastly different scales, and may be confusing things. – rcollyer Aug 08 '13 at 16:48
  • Good idea, thanks! That works for me too. PS I can't upvote your comments, I guess I don't have enough rep on this account yet. I'll try to remember when I get there. – jorgen Aug 08 '13 at 17:41
  • I get around 0.13956 when I use exact exponents, and around 37.896 when using exponents with decimals. Offhand I do not know if this is from arithmetic or from Integrate going through slightly different somersaults. Given the coefficient sizes and oscillations, huge numerical error here is not entirely a surprise. – Daniel Lichtblau Aug 08 '13 at 18:06
  • When I moreover rationalize the coefficients I get a longish exact result. It numericizes to around 0.13956. – Daniel Lichtblau Aug 08 '13 at 18:27

0 Answers0