文章对应的视频讲解 ：复化(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积分：<code>S = comp_simpson_integral(a,b,n,f)

输入

[

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下，下降速度明显，优势巨大