数值积分——复化(Simpson)辛普森求积公式 | 北太天元 or matlab

文章对应的视频讲解 ：复化(Simpson)辛普森求积公式

复化Simpson求积公式

Simpson公式
$${\int }_a^{b}f(x)\mathrm{d}x\approx \frac{b-a}{6}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right).$$

将积分区间$[a,b]$进行$2m$（偶）等分，记$n=2m$，其中 $n+1$ 是节点总数，$m$ 是积分子区间的总数。

步长$h=\frac{b-a}{n}$，节点$x_k=a+kh,(k=0,1,2,\cdots ,n)$

$${\int }_a^{b}f(x)\mathrm{d}x=\sum _{i=0}^{m-1}{\int }_{x_{2i}}^{x_{2i+2}}f(x)\mathrm{d}x$$
在$ [x_{2i},x_{2i+2}] $上用Simpson公式
$${\int }_{x_{2i}}^{x_{2i+2}}f(x)\mathrm{d}x\approx \frac{h}{3}[f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})]$$

累加得复化Simpson求积公式
$$S_n(f)\approx \frac{h}{3}\left[f(a)+4\sum _{i=0}^{m-1}f(x_{2i+1})+2\sum _{i=1}^{m-1}f(x_{2i})+f(b)\right]$$

算法

$♡$ 复化Simpson积分：S = comp_simpson_integral(a,b,n,f)

1. 输入

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

2. 实现过程

   • 判断$n$是否是偶数，若不是，$+1$变为偶数
   • $m=\frac{n}{2}$
   • 计算出$[a,b]$ $n$等分后得到的$n+1$个节点，构成向量$x_0$
   • $y_0=f(x_0)$
   • 计算
     $$sumy1=\sum _{i=0}^{m-1}f(x_{2i+1})$$$$sumy2=\sum _{i=1}^{m-1}f(x_{2i})$$$$S=\frac{h}{3}\left[f(a)+4\ast sumy1+2\ast sumy2+f(b)\right]$$

3. 输出 $S$

北太天元 or matlab 实现

function S = comp_simpson_integral(a,b,n,f)
% 复化Simpson求积
% [a,b]
% n ：小区间的个数, 要求是偶数
% f：定义好的函数
%
%   Version:            1.0
%   last modified:      07/14/2023
    if mod(n,2) != 0    % 判断n是否为偶数，如果不是，使其变为偶数
    	   n = n+1;
    end
    h = (b-a)/n;
    k = 0:1:n;
    xi = a + k * h;
    yi = f(xi);
    m = n/2;
    i1 = 0:1:m-1;
        sumy1 = sum(yi(2*i1+1 +1)); % f(x_{2i+1})求和
    i2 = 1:1:m-1;
        sumy2 = sum(yi(2*i2 +1)); % f(x_{2i})求和
    S = (yi(1) + 4*sumy1 + 2*sumy2 + yi(n+1)) * h/3;
end

保存为 comp_simpson_integral.m 文件

数值算例

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

clc;clear all;format long;
f = @(x) 4./(1+x.^2);

N = [10 20 40 80 160];
delta = zeros(1,5);
delta1 =delta;
delta2 = delta;
k = 1;
for n = N
    S = comp_simpson_integral(0,1,n,f);
    T = comp_tra_integral(0,1,n,f);
    delta1(k) = abs(pi - S);
    delta2(k) = abs(pi - T);
    k++;
end
    figure(1)
    	plot(N,delta1,'-or');
    title('复化Simpson误差随n的变化')
    figure(2)
    	plot(N,delta2,'-*b');
    title('复化梯形误差随n的变化')

运行后得到


对比这两个图像，发现，n= 20 时，两者的误差远远不在一个量级

复化Simpson下，下降速度明显，优势巨大