圆形区域函数数值积分未收敛至中心值的优化问询

我有一个函数u(x,y)，无显式表达式，由若干算术和三角函数定义，具备连续性等良好性质。我尝试用数值方法计算该函数在圆形区域（给定圆心、半径）和正方形区域（给定中心、边长）内的平均值（积分除以区域面积），所用Python代码如下：

import numpy as np

# 假设u函数已定义，示例为常数1时可替换为u=lambda x,y:1
def u(x, y):
    # 此处替换为实际函数实现
    return 1.0

# Integral in a column, x constant, y in ]y_c - l/2, y_c + l/2[, using n points:
def F(x, y_c, l, n):
    def f(b):
        return u(x,b)

    h = l / n

    y0 = y_c - l / 2
    yn = y_c + l / 2
    s = 0

    for i in range(1, n):
        yi = y0 + h * i
        s += f(yi)

    return h * (f(y0) + f(yn)) / 2 + h * s


# Integral in the circle of center (x_c, y_c) and radius r, using aprox n points:
def int_circ(x_c, y_c, r, n):
    def f(x, l):
        return F(x, y_c, l, int(np.sqrt(n)))

    h = 2 * r / np.sqrt(n)
    x0 = x_c - r
    xn = x_c + r
    s = 0

    for i in range(1, int(np.sqrt(n))):
        xi = x0 + h * i
        li = 2 * np.sqrt(r ** 2 - (xi - x_c) ** 2)
        s += f(xi, li)

    l = 2 * np.sqrt(r ** 2 - (r - h / 4) ** 2)

    return h / 2 * (f(x0 + h / 2, l) + f(xn - h / 2, l)) + h * s


# Integral in the square of center (x_c, y_c) and side 2r, using aprox n points:
def int_sqr(x_c, y_c, r, n):
    def f(x):
        return F(x, y_c, 2 * r, int(np.sqrt(n)))

    h = 2 * r / np.sqrt(n)
    x0 = x_c - r
    xn = x_c + r
    s = 0

    for i in range(1, int(np.sqrt(n))):
        xi = x0 + h * i
        s += f(xi)

    return h * (f(x0) + f(xn)) / 2 + h * s

分析结果发现，缩小正方形边长时，结果能如预期收敛到中心处的函数值；但圆形区域在n=10000时，结果与预期值仍有约10^-4的相对误差，即使将u替换为常数1，该问题依然存在。增大n能提升精度，但我认为应有更优方案。我尝试过多种圆形积分算法，均存在区域面积计算偏差问题，当前算法是最接近预期的。

优化方案 #

1. 修正当前切片积分的边界逻辑错误 #

你的int_circ函数在处理左右边界切片时，计算y方向长度的公式存在错误：

l = 2 * np.sqrt(r ** 2 - (r - h / 4) ** 2)

圆的最左/最右端（x0 = x_c - r、xn = x_c + r）对应的y方向长度应为0，而非自定义偏移值，这是常数函数积分仍有误差的核心原因。

修正后的int_circ函数（对齐梯形法则标准形式）：

def int_circ(x_c, y_c, r, n):
    def f(x, l):
        return F(x, y_c, l, int(np.sqrt(n))) if l > 1e-12 else 0.0  # 避免极小值计算误差

    n_x = int(np.sqrt(n))
    h = 2 * r / n_x
    x0 = x_c - r
    xn = x_c + r
    s = 0

    for i in range(1, n_x):
        xi = x0 + h * i
        dx = xi - x_c
        li = 2 * np.sqrt(r ** 2 - dx ** 2)
        s += f(xi, li)

    # 标准梯形法则：首尾项取半，中间项全加
    return h * (f(x0, 0) + f(xn, 0)) / 2 + h * s

2. 改用极坐标变换（高效且精度更高） #

圆形区域用极坐标变换能避免切片长度的复杂计算，积分元转化为ρ dρ dθ，收敛速度更快：

def int_circ_polar(x_c, y_c, r, n):
    n_theta = int(np.sqrt(n))
    n_r = int(np.sqrt(n))
    dtheta = 2 * np.pi / n_theta
    dr = r / n_r
    total = 0.0

    for i in range(n_theta):
        theta = dtheta * (i + 0.5)  # 中点采样提升精度
        cos_theta = np.cos(theta)
        sin_theta = np.sin(theta)
        r_sum = 0.0
        for j in range(n_r):
            rho = dr * (j + 0.5)
            x = x_c + rho * cos_theta
            y = y_c + rho * sin_theta
            r_sum += u(x, y) * rho
        total += r_sum * dr * dtheta

    return total

当u为常数1时，该方法计算的积分值会精确趋近于πr²，无额外偏差。

3. 使用自适应积分库函数 #

利用scipy的自适应二重积分函数，无需手动实现数值逻辑，精度和稳定性更优：

from scipy.integrate import dblquad

def int_circ_scipy(x_c, y_c, r):
    # 定义圆形区域的y边界
    def y_lower(x):
        return y_c - np.sqrt(r**2 - (x - x_c)**2)
    def y_upper(x):
        return y_c + np.sqrt(r**2 - (x - x_c)**2)
    # 计算二重积分，返回积分值和误差估计
    integral, error = dblquad(u, x_c - r, x_c + r, y_lower, y_upper)
    return integral

内容的提问来源于stack exchange，提问作者Martim Pinto Paiva