Python实现高斯-勒让德求积法求解含函数限二重积分

在不依赖任何第三方数值计算库的前提下，使用Python通过高斯-勒让德求积（Gauss–Legendre Quadrature）数值方法求解二重积分时，若积分上下限为函数形式，初始版本算法无法正常运行，初始实现代码如下：

def integrate(a: float, b: float, n: int, f_xy: callable, upper_func: callable, lower_func: callable) -> float:

    if n < 1:
        raise("n < 1 is invalid.")

    w, t = get_wt(n)

    e1_x = (b-a)/2
    e2_x = (a+b)/2
    sum = 0

    for i in range(n):

        x_i = e1_x*t[i]+e2_x

        if type(upper_func) == int:
            d = upper_func
        else:
            d = upper_func(x_i)

        if type(lower_func) == int:
            c = lower_func
        else:
            c = lower_func(x_i)

        e1_y = (d - c)/2
        e2_y = (c + d)/2

        som = 0
        for j in range(n):
            y_i = e1_y*t[j]+e2_y
            som += w[j]*f_xy(x_i, y_i)

        sum += w[i]*som

    result = (1/4)*(b-a)*(d-c)*sum

    return result

代码中get_wt()函数的作用是返回高斯-勒让德求积对应的权重数组与节点数组。

对原有逻辑做少量调整即可实现支持函数型上下限的二重积分计算，修正后的代码如下：

from typing import Union, Callable
from inspect import isfunction

def double_integrate(a: float, b: float, n: int, f_xy: Union[Callable, float],
                     up_func: Union[Callable, float], low_func: Union[Callable, float]) -> float:
     
    w, t = get_wt(n)

    e1_x = (b-a)/2
    e2_x = (a+b)/2

    sum = 0

    for i in range(n):

        som = 0

        x_i = e1_x*t[i]+e2_x

        if not isfunction(up_func):
            d = up_func
        else:
            d = up_func(x_i)

        if not isfunction(low_func):
            c = low_func
        else:
            c = low_func(x_i)

        e1_y = (d - c)/2
        e2_y = (c + d)/2

        for j in range(n):
            y_i = e1_y*t[j]+e2_y

            if isfunction(f_xy):
                som += w[j]*f_xy(x_i, y_i)
            else:
                som += w[j]*f_xy

        sum += w[i]*e1_y*som

    result = e1_x*sum

    return result

原代码的两处核心错误 #

• 缩放系数计算位置错误：原代码将y方向的区间缩放系数(d-c)放在所有循环外统一计算，但函数型上下限的y区间长度随每个x采样节点动态变化，放在外层会导致所有x节点共用循环最后一次得到的d、c值，计算结果完全偏离正确值。修正后将y方向的缩放系数e1_y放在每个x节点的循环内，乘入对应节点的内层求和结果中。
• 类型判断逻辑有缺陷：原代码仅判断上下限是否为int类型，若传入浮点数格式的常数上下限，会被误判为函数触发调用报错；修正后改为判断入参是否为函数类型，同时兼容整数、浮点数格式的常数上下限，也额外支持了常数形式的被积函数传入，无需额外封装返回常数的匿名函数。

内容的提问来源于stack exchange，提问作者user118799