0

我要找到 $I_m=\int_0^{2\pi} \prod_{k=1}^m cos(kx){}dx$,其中 $m=1,2,3\ldots$

简单的 SAGE 代码:

x=var('x')
f = lambda m,x : prod([cos(k*x) for k in range(1,m+1)])
for m in range(1,15+1):
    print m, numerical_integral(f(m,x), 0, 2*pi)[0],integrate(f(m,x),x,0,2*pi).n()

输出:

1 -1.47676658757e-16 0.000000000000000
2 -5.27735962315e-16 0.000000000000000
3 1.57079632679 1.57079632679490
4 0.785398163397 0.785398163397448
5 -2.60536121164e-16 0.000000000000000
6 -1.81559273097e-16 0.000000000000000
7 0.392699081699 0.392699081698724
8 0.343611696486 0.147262155637022
9 -1.72448482421e-16 0.294524311274043
10 -1.8747663502e-16 0.196349540849362
11 0.214757310304 0.312932080728671
12 0.190213617698 0.177941771394734
13 -1.30355375996e-16 0.208621387152447
14 -1.25168280013e-16 0.0859029241215959
15 0.138441766107 0.134223318939994

如您所见,数字答案是正确的,但是对于 $m=1,2,\ldots,7$ 而言,integrate(...) 的结果是正确的,然后存在一些错误。我们可以打印不定积分:

for m in range(7,11+1):
    print 'm=',m
    print 'Indef_I_m=',integrate(f(m,x),x)

和输出:

m = 7
Indef_I_m = 1/16*x + 1/16*sin(2*x) + 1/32*sin(4*x) + 7/384*sin(6*x) +
7/512*sin(8*x) + 3/320*sin(10*x) + 5/768*sin(12*x) + 5/896*sin(14*x) +
1/256*sin(16*x) + 1/384*sin(18*x) + 1/640*sin(20*x) + 1/704*sin(22*x) +
1/1536*sin(24*x) + 1/1664*sin(26*x) + 1/1792*sin(28*x)
m = 8
Indef_I_m = 3/128*x + 5/256*sin(2*x) + 1/32*sin(3*x) + 5/512*sin(4*x) +
5/768*sin(6*x) + 1/256*sin(8*x) + 1/256*sin(10*x) + 1/256*sin(12*x) +
1/256*sin(14*x) + 1/256*sin(16*x) + 7/2304*sin(18*x) + 3/1280*sin(20*x)
+ 5/2816*sin(22*x) + 1/768*sin(24*x) + 3/3328*sin(26*x) +
1/1792*sin(28*x) + 1/1920*sin(30*x) + 1/4096*sin(32*x) +
1/4352*sin(34*x) + 1/4608*sin(36*x) + 3/32*sin(x)
m = 9
Indef_I_m = 3/64*x + 3/128*sin(2*x) + 23/768*sin(3*x) + 3/256*sin(4*x) +
3/640*sin(5*x) + 1/128*sin(6*x) + 5/1792*sin(7*x) + 5/2304*sin(9*x) +
3/2816*sin(11*x) + 1/832*sin(13*x) + 1/1280*sin(15*x) + 3/4352*sin(17*x)
+ 5/4864*sin(19*x) + 1/1344*sin(21*x) + 3/2944*sin(23*x) +
7/6400*sin(25*x) + 1/1152*sin(27*x) + 3/3712*sin(29*x) +
5/7936*sin(31*x) + 1/2112*sin(33*x) + 3/8960*sin(35*x) +
1/4736*sin(37*x) + 1/4992*sin(39*x) + 1/10496*sin(41*x) +
1/11008*sin(43*x) + 1/11520*sin(45*x) + 23/256*sin(x)
m = 10
Indef_I_m = 1/32*x + 1/64*sin(2*x) + 17/512*sin(3*x) + 1/128*sin(4*x) +
7/2560*sin(5*x) + 1/192*sin(6*x) + 3/1792*sin(7*x) + 1/1152*sin(9*x) +
5/5632*sin(11*x) + 3/6656*sin(13*x) + 1/2560*sin(15*x) +
5/8704*sin(17*x) + 3/9728*sin(19*x) + 1/2688*sin(21*x) +
1/2944*sin(23*x) + 1/6400*sin(25*x) + 1/4608*sin(27*x) +
3/14848*sin(29*x) + 3/15872*sin(31*x) + 5/16896*sin(33*x) +
3/8960*sin(35*x) + 3/9472*sin(37*x) + 1/3328*sin(39*x) +
5/20992*sin(41*x) + 1/5504*sin(43*x) + 1/7680*sin(45*x) +
1/12032*sin(47*x) + 1/12544*sin(49*x) + 1/26112*sin(51*x) +
1/27136*sin(53*x) + 1/28160*sin(55*x) + 13/128*sin(x)
m = 11
Indef_I_m = 51/1024*x + 53/2048*sin(2*x) + 13/768*sin(3*x) + 53/4096*sin(4*x) +
13/1536*sin(6*x) + 1/2048*sin(8*x) + 1/2560*sin(10*x) + 1/3072*sin(12*x)
+ 5/14336*sin(14*x) + 1/4096*sin(16*x) + 5/18432*sin(18*x) +
1/4096*sin(20*x) + 1/5632*sin(22*x) + 5/24576*sin(24*x) +
5/26624*sin(26*x) + 5/28672*sin(28*x) + 1/5120*sin(30*x) +
3/16384*sin(32*x) + 5/34816*sin(34*x) + 1/9216*sin(36*x) +
5/38912*sin(38*x) + 1/10240*sin(40*x) + 1/10752*sin(42*x) +
3/22528*sin(44*x) + 3/23552*sin(46*x) + 1/8192*sin(48*x) +
3/25600*sin(50*x) + 5/53248*sin(52*x) + 1/13824*sin(54*x) +
3/57344*sin(56*x) + 1/29696*sin(58*x) + 1/30720*sin(60*x) +
1/63488*sin(62*x) + 1/65536*sin(64*x) + 1/67584*sin(66*x) +
13/256*sin(x)

所以对于 $m=7$ 的答案是正确的,通过 WolframAlpha 与 Indef_I_7比较 ,对于 $m=8$ 的答案是不正确的Indef_I_8 via WolframAlpha 应该有 Indef_I_8=$\frac{7x}{128}+\ldots$ 并且没有 $\ sin(x)$, $\sin(3x)$ 的总和,只有 $\sin(2k)$ for $k=1,2,3,\ldots 18$

对不起体积计算!问题是 - 这是符号集成中的错误,我说得对吗?

4

2 回答 2

0

好吧,显然设置algorithm='mathematica_free'解决了这个问题;这可能是 SAGE ( 'maxima') 使用的默认算法中的一个错误。

于 2013-06-17T19:03:18.383 回答
0

所以答案是——是的。这是 algorithm='maxima' 中的一个错误,所以使用 algorithm='mathematica_free' (定义新函数来找到定积分)或 simple_full() 用于 cos(kx) 的这种乘积而不是积分。

f(8,x).simplify_full().integral(x,0,2pi) == 7/64pi and 
(7/64*pi).n() == 0.343611696486384 is correct
于 2013-06-18T04:07:13.250 回答