数值积分

$$ \int_{a}^{b} f(x) \, dx \approx \sum_{k=0}^{n} A_{k}f(x_k) $$

$A_k$为点$x_k$的权

简单的数值积分直接取权为1，n越大精度越高

def i_f(f, a, b, n):
    res = 0

    h = (b - a) / n

    for i in range(n):
        left = a + i * h
        res += (f(left) * h)

    return res

import math

print(i_f(f=math.sin, a=0, b=math.pi, n=10))
print(i_f(f=math.sin, a=0, b=math.pi, n=13))
print(i_f(f=math.sin, a=0, b=math.pi, n=16))
print(i_f(f=math.sin, a=0, b=math.pi, n=19))
print(i_f(f=math.sin, a=0, b=math.pi, n=22))
print(i_f(f=math.sin, a=0, b=math.pi, n=50))

1.9835235375094544
1.990257175347774
1.993570343772339
1.9954413183201944
1.9966002202692708
1.9993419830762613

复合辛普森求积公式

$$ I = \int_{a}^{b} f(x) \, dx = \frac{h}{6}\sum_{k=0}^{n-1}[f(x_k)+4f(x_{k+1/2})+f(x_{k+1})] + R_n(f) $$ 其中 $$ x_{k+1/2}=x_k+\frac{1}{2}h $$

def simpson_sum(f, a, b, n):
    h = (b - a) / n

    res = 0

    for i in range(n):
        xk = a + i * h
        xk1 = xk + h

        res += (f(xk) + 4 * f(xk + h / 2) + f(xk1))

    res *= (h / 6)
    return res

import math

print(simpson_sum(f=math.sin, a=0, b=math.pi, n=4))
print(simpson_sum(f=math.sin, a=0, b=math.pi, n=5))
print(simpson_sum(f=math.sin, a=0, b=math.pi, n=6))
print(simpson_sum(f=math.sin, a=0, b=math.pi, n=7))
print(simpson_sum(f=math.sin, a=0, b=math.pi, n=8))
print(simpson_sum(f=math.sin, a=0, b=math.pi, n=9))
print(simpson_sum(f=math.sin, a=0, b=math.pi, n=20))

2.0002691699483877
2.0001095173150043
2.0000526243411856
2.0000283435514685
2.0000165910479355
2.000010347705775
2.000000423093183