python 上的辛普森整合

Question

我正在尝试使用 simpson 积分规则对 f(x) = 2x 从 0 到 1 进行数值积分，但一直出现大错误。期望的输出是 1，但是 python 的输出是 1.334。有人可以帮我找到解决这个问题的方法吗？谢谢。

import numpy as np

def f(x):
    return 2*x

def simpson(f,a,b,n):
    x = np.linspace(a,b,n)
    dx = (b-a)/n
    for i in np.arange(1,n):
        if i % 2 != 0:
            y = 4*f(x)
        elif i % 2 == 0:
            y = 2*f(x)
    return (f(a)+sum(y)+f(x)[-1])*dx/3

a = 0
b = 1
n = 1000
ans = simpson(f,a,b,n)
print(ans)

Answer 1

一切都错了。 x 是一个数组，每次调用 f(x) 时，您都是在整个数组上计算函数。由于 n 是偶数而 n-1 是奇数，所以最后一个循环中的 y 是 4*f(x) 并且从它的总和中计算出一些东西

那么n就是段数。点数为n+1。正确的实现是

def simpson(f,a,b,n):
    x = np.linspace(a,b,n+1)
    y = f(x)
    dx = x[1]-x[0]
    return (y[0]+4*sum(y[1::2])+2*sum(y[2:-1:2])+y[-1])*dx/3

simpson(lambda x:2*x, 0, 1, 1000)

然后正确 returns 1.000。如果 n 是偶数，您可能想添加一个测试，如果不是这样，则将其加一。

如果你真的想保持循环，你需要在循环内实际累加和。

def simpson(f,a,b,n):
    dx = (b-a)/n;
    res = 0; 
    for i in range(1,n): res += f(a+i*dx)*(2 if i%2==0 else 4); 
    return (f(a)+f(b) + res)*dx/3;

simpson(lambda x:2*x, 0, 1, 1000)

但是循环通常比向量化操作慢，所以如果你使用 numpy，请使用向量化操作。或者直接使用 scipy.integrate.simps.

python 上的辛普森整合

simpson integration on python

python-3.x

numerical-methods