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

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

复化Simpson求积公式

Simpson公式

∫ a b f ( x ) d x ≈ b − a 6 ( f ( a ) + 4 f ( a + b 2 ) + f ( b ) ) . \int_a^bf(x)\mathrm{d}x\approx\frac{b-a}{6}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right). ∫ab​f(x)dx≈6b−a​(f(a)+4f(2a+b​)+f(b)).

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

步长 h = b − a n h=\frac{b-a}{n} h=nb−a​，节点 x k = a + k h , ( k = 0 , 1 , 2 , ⋯   , n ) x_{k}=a+kh,(k=0,1,2,\cdots,n) xk​=a+kh,(k=0,1,2,⋯,n)

∫ a b f ( x ) d x = ∑ i = 0 m − 1 ∫ x 2 i x 2 i + 2 f ( x ) d x \int_{a}^{b}f(x)\mathrm{d}x=\sum_{i=0}^{m-1}\int_{x_{2i}}^{x_{2i+2}}f(x)\mathrm{d}x ∫ab​f(x)dx=i=0∑m−1​∫x2i​x2i+2​​f(x)dx

在$ [x_{2i},x_{2i+2}] $上用Simpson公式

∫ x 2 i x 2 i + 2 f ( x ) d x ≈ h 3 [ f ( x 2 i ) + 4 f ( x 2 i + 1 ) + f ( x 2 i + 2 ) ] \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})] ∫x2i​x2i+2​​f(x)dx≈3h​[f(x2i​)+4f(x2i+1​)+f(x2i+2​)]

累加得复化Simpson求积公式

S n ( f ) ≈ h 3 [ f ( a ) + 4 ∑ i = 0 m − 1 f ( x 2 i + 1 ) + 2 ∑ i = 1 m − 1 f ( x 2 i ) + f ( b ) ] 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] Sn​(f)≈3h​[f(a)+4i=0∑m−1​f(x2i+1​)+2i=1∑m−1​f(x2i​)+f(b)]

算法

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

1. 输入

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

   • 实现过程

     • 判断 n n n是否是偶数，若不是， + 1 +1 +1变为偶数
     • m = n 2 m = \frac{n}{2} m=2n​
     • 计算出 [ a , b ] [a, b] [a,b] n n n等分后得到的 n + 1 n + 1 n+1个节点，构成向量 x 0 x_0 x0​
     • y 0 = f ( x 0 ) y_0 = f(x_0) y0​=f(x0​)
     • 计算

       s u m y 1 = ∑ i = 0 m − 1 f ( x 2 i + 1 ) sumy1 = \sum_{i=0}^{m-1}f(x_{2i+1}) sumy1=i=0∑m−1​f(x2i+1​) s u m y 2 = ∑ i = 1 m − 1 f ( x 2 i ) sumy2 = \sum_{i=1}^{m-1}f(x_{2i}) sumy2=i=1∑m−1​f(x2i​) S = h 3 [ f ( a ) + 4 ∗ s u m y 1 + 2 ∗ s u m y 2 + f ( b ) ] S =\frac{h}{3}\left[f(a)+4*sumy1+2*sumy2+f(b)\right] S=3h​[f(a)+4∗sumy1+2∗sumy2+f(b)]

     • 输出 S S 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 文件

数值算例

用数值积分法近似计算

π = 4 ∫ 0 1 1 1 + x 2 d x \pi = 4\int_0^1 \frac{1}{1+x^2}\mathrm{d}x π=4∫01​1+x21​dx编写复化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下，下降速度明显，优势巨大