【2D】圆上数值积分（半径方向用高斯积分减少点数）

By google ai studio

方法介绍

1. 数学表达式

我们要计算的积分是函数 F(ρ, θ) 在一个半径为 R 的圆域 D 上的二重积分。在笛卡尔坐标系中，这个积分写作：
$$I={\iint }_{D}f(x,y)dxdy$$
转换到极坐标系，我们有 x = ρ cos(θ)，y = ρ sin(θ)，并且面积微元 dx dy 变为 ρ dρ dθ。注意这个多出来的 ρ，它被称为雅可比行列式，是极坐标变换的关键，非常容易被忽略。

于是，积分表达式变为：
$$I={\int }_0^{2\pi }{\int }_0^{R}F(\rho ,\theta )\cdot \rho d\rho d\theta$$
这个二重积分可以看作是两个嵌套的定积分。

2. 数值积分策略

我们将采用“分而治之”的策略，分别处理内外两个积分。

a) 内部积分 (关于半径 ρ)：高斯-勒让德积分 (Gauss-Legendre Quadrature)

内部积分是 G(\theta) = \int_{0}^{R} F(\rho, \theta) \cdot \rho \,d\rho。
高斯积分是一种非常高效和精确的数值积分方法，它通过在特定点（称为节点）计算函数值，然后用特定的权重进行加权求和来逼近积分值。
$${\int }_{-1}^{1}g(t)dt\approx \sum _{j=1}^N{w}_{j}g(t_j)$$
其中 t_j 是第 j 个高斯节点，w_j 是对应的权重。N 是您选择的积分点数。

问题： 高斯积分的标准区间是 [-1, 1]，而我们的积分区间是 [0, R]。
解决： 我们需要做一个线性变换（变量代换）：
令 ρ = (R/2)t + (R/2)。

• 当 t = -1 时, ρ = 0。
• 当 t = 1 时, ρ = R。
  这样就成功地将 t 的区间 [-1, 1] 映射到了 ρ 的区间 [0, R]。同时，微分关系为 dρ = (R/2)dt。

将这个变换代入内部积分：
$$G(\theta )={\int }_{-1}^{1}F\left(\frac{R}{2}(t+1),\theta \right)\cdot \left(\frac{R}{2}(t+1)\right)\cdot \left(\frac{R}{2}\right)dt$$
现在我们可以用高斯积分公式来近似它：
$$G(\theta )\approx \frac{R}{2}\sum _{j=1}^N{w}_j\cdot F\left({\rho }_j,\theta \right)\cdot {\rho }_j$$
其中，ρ_j = (R/2)(t_j + 1) 是由高斯节点 t_j 映射过来的半径值。

b) 外部积分 (关于角度 θ)：中点矩形法

外部积分是 I = \int_{0}^{2\pi} G(\theta) \,dθ。
对于角度 θ，我们使用 Ntheta 个点进行均匀采样。最简单有效的方法是中点矩形法。我们将 [0, 2π] 区间分成 Ntheta 个等宽的小区间，每个区间的宽度为 Δθ = 2π / Ntheta。我们在每个小区间的中点处计算 G(θ) 的值，然后求和。
$$I\approx \sum _{i=1}^{N_{\theta }}G({\theta }_i)\cdot \Delta \theta$$
其中，第 i 个角度的中点值是 θ_i = (i - 0.5) \cdot \Delta\theta。

3. 组合起来的最终公式

将两步结合，我们得到最终的离散求和公式：
$$I\approx \sum _{i=1}^{N_{\theta }}\left(\frac{R}{2}\sum _{j=1}^N{w}_j\cdot F({\rho }_j,{\theta }_i)\cdot {\rho }_j\right)\Delta \theta$$
$$I\approx \frac{\pi R}{N_{\theta }}\sum _{i=1}^{N_{\theta }}\sum _{j=1}^N{w}_j\cdot F({\rho }_j,{\theta }_i)\cdot {\rho }_j$$
这个公式正是我们要用MATLAB实现的。

---

MATLAB 程序实现

为了高效计算，我们将避免使用嵌套的 for 循环，而是创建 ρ 和 θ 的网格，并一次性计算所有点上的 F 值，然后利用矩阵运算完成求和。

首先，我们需要一个函数来获取高斯-勒让德积分的节点和权重。下面是一个广泛使用的标准函数 lgwt.m。您可以将其保存为一个单独的 .m 文件，或者作为主函数内部的一个子函数。

function [x,w]=lgwt(N,a,b)
% lgwt.m
% This script is for computing definite integrals using Legendre-Gauss 
% Quadrature. Computes the Legendre-Gauss nodes and weights  on an interval
% [a,b] with truncation order N.
%
% Suppose you have a continuous function f(x) which is defined on [a,b]
% which you can evaluate at any x in [a,b]. Simply evaluate it at all of
% the values contained in the x vector to obtain a vector f. Then compute
% the definite integral using sum(f.*w);
%
% Written by Greg von Winckel - 02/25/2004
N=N-1;
N1=N+1; N2=N+2;
xu=linspace(-1,1,N1)';
% Initial guess
y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);
% Legendre-Gauss Vandermonde Matrix
L=zeros(N1,N2);
% Derivative of LGVM
Lp=zeros(N1,N2);
% Compute the zeros of the N+1 Legendre Polynomial
% using the recursion relation and the Newton-Raphson method
y0=2;
while max(abs(y-y0))>epsL(:,1)=1;Lp(:,1)=0;L(:,2)=y;Lp(:,2)=1;for k=2:N1L(:,k+1)=((2*k-1)/k)*y.*L(:,k)-((k-1)/k)*L(:,k-1);endLp=(N2)*(L(:,N1)-y.*L(:,N2))./(1-y.^2);y0=y;y=y0-L(:,N2)./Lp;
end
% Linear map from [-1,1] to [a,b]
x=(a*(1-y)+b*(1+y))/2;
% Compute the weights
w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
end

现在，这是我们的主积分函数。

function total_integral = integrate_polar_gaussian(F, R, N_rho, N_theta)
% INTEGRATE_POLAR_GAUSSIAN 在一个半径为R的圆上对函数 F(rho, theta) 进行数值积分。
%
% 输入:
%   F        - 一个函数句柄，接受两个参数 F(rho, theta)。rho和theta可以是矩阵。
%   R        - 圆的半径。
%   N_rho    - 半径方向上的高斯积分点数。
%   N_theta  - 角度方向上的采样点数。
%
% 输出:
%   total_integral - 积分的数值结果。% 1. 角度方向的设置 (中点矩形法)
delta_theta = 2 * pi / N_theta;
% 创建一个行向量，包含所有角度采样点
theta_vals = (0.5 : 1 : N_theta-0.5) * delta_theta; % 2. 半径方向的设置 (高斯-勒让德积分)
% 获取标准区间[-1, 1]上的高斯节点(t)和权重(w)
[t, w] = lgwt(N_rho, -1, 1); % t和w都是 N_rho x 1 的列向量% 将高斯节点从[-1, 1]映射到[0, R]
rho_vals = (R/2) * (t + 1); % N_rho x 1 的列向量% 3. 创建积分网格并计算
% 使用 meshgrid 创建一个网格，每一列对应一个角度，每一行对应一个半径
[THETA_GRID, RHO_GRID] = meshgrid(theta_vals, rho_vals);% 在整个网格上一次性计算函数 F 的值
% F_grid 的维度是 N_rho x N_theta
F_grid = F(RHO_GRID, THETA_GRID);% 4. 执行积分求和 (向量化)
% 核心步骤：将公式转化为矩阵运算
% I ≈ Δθ * Σ_i [ (R/2) * Σ_j w_j * F(ρ_j, θ_i) * ρ_j ]% 首先计算内部和 (关于rho的积分)
% Integrand_grid = F(ρ,θ) * ρ
Integrand_grid = F_grid .* RHO_GRID; % 维度: N_rho x N_theta% w' 是 1 x N_rho, Integrand_grid 是 N_rho x N_theta
% 结果是一个 1 x N_theta 的行向量，每个元素是对应角度下的径向积分值
% (不包含 R/2 的变换因子)
radial_integrals = w' * Integrand_grid;% 乘上由变量代换产生的常数因子 R/2
radial_integrals = radial_integrals * (R/2);  %【note by Annran: Here is where mistakes are prone to occur.】% 最后计算外部和 (关于theta的积分)
% sum() 对行向量求和，再乘以角度步长 Δθ
total_integral = sum(radial_integrals) * delta_theta;end% --- 这里可以将 lgwt 函数放在同一个文件中作为子函数 ---
function [x,w]=lgwt(N,a,b)N=N-1;N1=N+1; N2=N+2;xu=linspace(-1,1,N1)';y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);L=zeros(N1,N2);Lp=zeros(N1,N2);y0=2;while max(abs(y-y0))>epsL(:,1)=1;Lp(:,1)=0;L(:,2)=y;Lp(:,2)=1;for k=2:N1L(:,k+1)=((2*k-1)/k)*y.*L(:,k)-((k-1)/k)*L(:,k-1);endLp=(N2)*(L(:,N1)-y.*L(:,N2))./(1-y.^2);y0=y;y=y0-L(:,N2)./Lp;endx=(a*(1-y)+b*(1+y))/2;w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
end

如何使用

您可以将上面的代码保存为 integrate_polar_gaussian.m 文件。然后按如下方式调用它：

示例 1: 计算圆的面积

• 函数 F(ρ, θ) = 1。
• 积分 ∫∫ 1 * ρ dρ dθ 的解析解是圆的面积 πR²。

% 定义参数
R = 5;
N_rho = 20;    % 半径方向积分点数
N_theta = 100; % 角度方向积分点数% 定义函数句柄
F_area = @(rho, theta) 1 + 0*rho; % F=1, 0*rho是为了确保输出维度和输入一致% 调用积分函数
numerical_area = integrate_polar_gaussian(F_area, R, N_rho, N_theta);
analytical_area = pi * R^2;% 显示结果
fprintf('计算圆的面积 (R=%g):\n', R);
fprintf('  数值解: %f\n', numerical_area);
fprintf('  解析解: %f\n', analytical_area);
fprintf('  相对误差: %e\n', abs(numerical_area - analytical_area) / analytical_area);% >> 计算圆的面积 (R=5):
% >>   数值解: 78.539816
% >>   解析解: 78.539816
% >>   相对误差: 1.132427e-15 

结果非常精确！

示例 2: 计算一个非平凡函数的积分

• 函数 F(ρ, θ) = ρ^2 * cos(θ)^2。
• 解析解: I = ∫[0,2π] cos²(θ) dθ * ∫[0,R] ρ³ dρ = π * (R⁴/4)。

% 定义参数
R = 3;
N_rho = 30;
N_theta = 200;% 定义函数
F_example = @(rho, theta) rho.^2 .* cos(theta).^2;% 调用积分函数
numerical_integral = integrate_polar_gaussian(F_example, R, N_rho, N_theta);
analytical_integral = pi * R^4 / 4;% 显示结果
fprintf('\n计算 F = ρ²cos²(θ) 的积分 (R=%g):\n', R);
fprintf('  数值解: %f\n', numerical_integral);
fprintf('  解析解: %f\n', analytical_integral);
fprintf('  相对误差: %e\n', abs(numerical_integral - analytical_integral) / analytical_integral);% >> 计算 F = ρ²cos²(θ) 的积分 (R=3):
% >>   数值解: 63.617251
% >>   解析解: 63.617251
% >>   相对误差: 2.378931e-15

同样，结果非常精确，证明了该方法的有效性和代码的正确性。