使用辛普森法则计算错误的积分值

Question

我为 Simpson 的规则集成编写了以下代码来近似 sin。方程式在此attachment 中。我为偶数和奇数项编写了单独的循环，因为它们在附件中分组显示。

import math 

x1=0
x2=math.pi
N=6
delta_x=(x2-x1)/N
f=math.sin
sum=f(0)

#odd summation
for i in range (1,N+1):
    sum=sum+f(x1+(2*i+1)*delta_x)

sum=4*sum
print(sum)

#even summation

for i in range(2,N+1):
    sumeven=0
    sumeven=sumeven+f(x1+(2*i)*delta_x)
sumeven=2*sumeven

sumeven=sumeven+f(N)
print(sumeven)

integral=(delta_x/3)*(sum+sumeven) 
print(integral)   

但是当我打印这些值时，它给了我非常小的负数。

谁能看到我的代码有什么问题？

Answer 1

首先，学习一些基本的调试。您只为这段代码做了两行简单的print 行，这两行都没有跟踪更详细的问题。请参阅这个可爱的 debug 博客寻求帮助。我添加了一些工具，清理了一些样式，然后再次运行代码，更好地跟踪逻辑和数据流。诊断如下。

sum_odd = f(0)

# odd summation
for i in range (1, N+1):
    x_val = x1 + (2*i+1)*delta_x
    sum_odd = sum_odd + f(x_val)
    print (x_val/math.pi, "pi", f(x_val), sum_odd)

sum_odd = 4*sum_odd

print(sum_odd)

# even summation
for i in range(2, N+1):
    sum_even = 0
    x_val = x1 + (2*i)*delta_x
    sum_even = sum_even + f(x_val)
    print (x_val/math.pi, "pi", f(x_val), sum_even)

输出：

0.5 pi 1.0 1.0
0.8333333333333333 pi 0.5000000000000003 1.5000000000000004
1.1666666666666665 pi -0.4999999999999997 1.0000000000000007
1.5 pi -1.0 6.661338147750939e-16
1.8333333333333333 pi -0.5000000000000004 -0.4999999999999998
2.1666666666666665 pi 0.4999999999999993 -4.996003610813204e-16
-1.9984014443252818e-15
0.6666666666666666 pi 0.8660254037844387 0.8660254037844387
1.0 pi 1.2246467991473532e-16 1.2246467991473532e-16
1.3333333333333333 pi -0.8660254037844384 -0.8660254037844384
1.6666666666666665 pi -0.866025403784439 -0.866025403784439
2.0 pi -2.4492935982947064e-16 -2.4492935982947064e-16
-0.27941549819892636
-0.04876720424671586

诊断

你有两个直接的问题：

1. 您在循环中初始化 sum_even inside，这会丢弃该系列一半的先前计算。将其移到循环之前。
2. 您在指定范围内积分超过两次；仅此一项，将为您提供一个很小的值（理想情况下，这将是 0，即函数完整循环的积分）。

修复初始化，修复限制，你应该会看到好的结果，更像这样：

0.16666666666666666 pi 0.49999999999999994 0.49999999999999994
0.5 pi 1.0 1.5
0.8333333333333333 pi 0.5000000000000003 2.0000000000000004
8.000000000000002
0.3333333333333333 pi 0.8660254037844386 0.8660254037844386
0.6666666666666666 pi 0.8660254037844387 1.7320508075688772
1.0 pi 1.2246467991473532e-16 1.7320508075688774
3.184686116938829
1.9520959854268212

另外，我强烈建议您阅读有关调试和编码风格的教程；这些将帮助您在以后的工作中。我从经验中知道。 :-)