数值积分——复化梯形求积公式 | 北太天元

复化求积法的思想：
将区间$[a,b]$进行$n$等分，步长$h=\frac{b-a}{n}$，等分点$x_k=a+kh,k=0,1,2,\cdots ,n$,
先在每个子区间$[x_k,x_{k+1}]$上采用低阶的数值求积公式求得近似积分值$I_k$,
再将它们累加并以和${\sum }_{k=0}^{n-1}{I}_k$作为积分$\mathrm{I}={\int }_{\mathrm{a}}^{\mathrm{b}}\mathrm{f}(\mathrm{x})\mathrm{dx}$的近似值.

复化梯形求积公式

将积分区间$[a,b]$进行$n$等分，步长为$h=\frac{b-a}{n}$，节点$x_k=a+kh,(k=0,1,2,\cdots ,n)$
$$I(f)={\int }_a^{b}f(x)\mathrm{d}x=\sum _{i=0}^{n-1}{\int }_{x_i}^{x_{i+1}}f(x)\mathrm{d}x$$

在$[x_i,x_{i+1}]$上使用梯形公式
$${\int }_{x_i}^{x_{i+1}}f(x)\mathrm{d}x=\frac{h}{2}[f(x_i)+f(x_{i+1})]-f^{\prime \prime }({\xi }_i)\frac{h^3}{12}\approx \frac{h}{2}[f(x_i)+f(x_{i+1})]$$

有
$$I(f)=\sum _{i=0}^{n-1}\left\{\frac{h}{2}[f(x_i)+f(x_{i+1})]-f^{\prime \prime }({\xi }_i)\frac{h^3}{12}\right\}\approx h\left[\frac{1}{2}f(a)+\sum _{i=1}^{n-1}f(x_i)+\frac{1}{2}f(b)\right]$$

得复化梯形公式
$$T(h)=T_n(f)=h\left[\frac{1}{2}f(a)+\sum _{i=1}^{n-1}f(x_i)+\frac{1}{2}f(b)\right]\approx {\int }_a^{b}f(x)\mathrm{d}x$$

---

算法

$♡$ 复化梯形求积公式：T = comp_tra_integral(a,b,n,f)

1. 输入

   • $[a,b]$
   • $n$：将$[a,b]$ $n$等分
   • $f$：已经定义好的函数，支持向量运算

2. 实现步骤

   • 计算出$[a,b]$ $n$等分后得到的$n+1$个节点，构成向量$x_0$
   • $y_0=f(x_0)$
   • 代入
     $${\int }_a^{b}f(x)\mathrm{d}x\approx h\left[\frac{1}{2}f(a)+\sum _{i=1}^{n-1}f(x_i)+\frac{1}{2}f(b)\right]=T$$

3. 输出

   • $T$：通过复化梯形求积公式得到的积分近似值

北太天元源程序

function T = comp_tra_integral(a,b,n,f)
% [a,b]
% n ：小区间的个数
% f：定义好的函数
%
%   Version:            1.0
%   last modified:      07/10/2023
    h = (b-a)/n;
    k = 0:1:n;
    xi = a + k * h;
    yi = f(xi);
    sumy = sum(yi(2:1:n));
    T = (yi(1)/2 + sumy + yi(n+1)/2)*h;
end

保存为 comp_tra_integral.m文件

数值算例

用数值积分法近似计算
$$\pi =4{\int }_0^1\frac{1}{1+x^2}\mathrm{d}x$$
编写复化梯形公式的实现程序，分别取剖分段数 $ n = 10, 20, 40, 80, 160, $ 计算积分值与 $\pi$ 的误差并作图;

% 复化梯形求积例子
%   last modified:      07/11/2023
%   file need: comp_tra_integral.m
%%
clc;clear all;format long;
f = @(x) 4./(1+x.^2);

N = [10 20 40 80 160];
delta = zeros(1,5);
k = 1;
for n = N
    T = comp_tra_integral(0,1,n,f);
    delta(k) = abs(pi-T);
    k++;
end
    	plot(N,delta,'b');
disp(delta);

运行后得到

可以发现，复化求积法下对积分的计算是非常稳定的，避开了Newton-Cotes 公式下随阶数升高带来的不稳定性。