Python 5种经典的求积分方法

Dilkople

21720人浏览 · 2022-04-22 22:56:40

Dilkople · 2022-04-22 22:56:40 发布

不同的方法最大区别在于公式的不同/矩形区域的不同，依据公式而定。本文代码为不严格代码，最严格的代码为先判断连续性。当你需要对大量函数积分（懒得写函数连续性检测的时候）可以使用本文代码用try…except…不严格替代。

代码使用方式：修改func(x)中的函数，然后可以看代码中的示例代码进行修改。只需要输入函数，开始值，结束值，迭代次数/精度，方法函数名即可使用。

1. Simpson积分法:
$\int _ { a } ^ { b } f ( x ) d x \approx \frac { ( b - a ) ( f ( a ) + f ( b ) + 4 f ( \frac { a + b } { 2 } ) ) } { 6 }$
2. 右矩形积分公式
$\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b - a } { n } f ( a + \frac { b - a } { n } i )$
3. 左矩形积分公式
$\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b - a } { n } f ( a + \frac { b - a } { n } (i - 1) )$
4. 中矩形积分公式
$\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b - a } { n } f ( a + \frac { (b - a)} { 2 n }(2i - 1 ) )$
5. 梯形积分公式
$\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b - a } { 2n } (f ( a + \frac { b - a } { n }( i - 1) ) + f ( a + \frac { b - a } { n } i ))$

from math import *


def sgn(x):
    if x > 0:
        return 1
    elif x == 0:
        return 0
    else:
        return -1


def func(x):
    # return sgn(sin(pi / sin(x)))
    return x ** 2 + 2


def simpson(begin, between, i):
    a = begin + (i - 1) * between
    b = begin + i * between
    return between * (func(a) + func(b) + 4 * func((a + b) / 2)) / 6


def rightRectangle(begin, between, i):
    return between * func(begin + between * i)


def leftRectangle(begin, between, i):
    return between * func(begin + between * (i - 1))


def midRectangle(begin, between, i):
    return between * func(begin + between / 2 * (2 * i - 1))


def trapezoidRectangle(begin, between, i):
    return between / 2 * (func(begin + between * (i - 1)) + func(begin + between * i))


def cal_IntegralByEpsilon(begin, end, epsilon, method):
    n = 1
    result = 0
    preResult = float("inf")
    while abs(preResult - result) >= epsilon:
        preResult = result
        result = 0
        n *= 2
        between = (end - begin) / n
        for i in range(n):
            try:
                result += method(begin, between, i + 1)
            except:
                return "Integrated function has discontinuity or does not defined in current interval"
    return result


def cal_IntegralByN(begin, end, n, method):
    result = 0
    between = (end - begin) / n
    for i in range(n):
        try:
            result += method(begin, between, i + 1)
        except:
            return "Integrated function has discontinuity or does not defined in current interval"


if __name__ == "__main__":
    begin = 0.2
    end = 0.5
    epsilon = 0.0001
    print(cal_IntegralByEpsilon(begin, end, epsilon, simpson))
    print(cal_IntegralByEpsilon(begin, end, epsilon, rightRectangle))
    print(cal_IntegralByEpsilon(begin, end, epsilon, leftRectangle))
    print(cal_IntegralByEpsilon(begin, end, epsilon, midRectangle))
    print(cal_IntegralByEpsilon(begin, end, epsilon, trapezoidRectangle))

    begin = 0.2
    end = 0.5
    n = 10
    print(cal_IntegralByN(begin, end, n, simpson))
    print(cal_IntegralByN(begin, end, n, rightRectangle))
    print(cal_IntegralByN(begin, end, n, leftRectangle))
    print(cal_IntegralByN(begin, end, n, midRectangle))
    print(cal_IntegralByN(begin, end, n, trapezoidRectangle))