Halo

A magic place for coding

0%

数值积分解法

介绍

   这篇博客讨论的是一个数值积分问题,对于一些难于求积的函数,使用 ** 牛顿 - 莱布尼茨 ** 公式显然是不科学的。因此对于这类问题,我们可以从 ** 积分中值定理 ** 出发,使用矩形或梯形的面积去近似积分值。
   这里主要讨论的是 ** 复合梯形公式 ** 和 ** 复合辛普森公式 **,统称为 ** 复合求积法 **。这种通过把积分区间细分成若干个子区间(通常是等分),再在每个子区间上使用低阶求积公式,从而提高了计算精度。

分析

1. 复合梯形公式

   该公式核心思想是对于细分后的每一个子区间,使用梯形公式求积。假设将区间 $[a,b]$ 等分为 $n$ 个子区间,分点 $x_k = a + kh,h = \frac {b-a}{n},k=0,1,…,n$,公式的原理就是对于每一块小梯形而言,上底为 $f (x_k)$,下底为 $f (x_{k+1})$,高为 $h$。
   复合梯形公式如下:
$$
I = \int^b_af (x) dx = \sum^{n-1}{k=0}\int^{x{k+1}}{x_k} f (x) dx = \frac {h}{2}\sum^{n-1}_{k=0}[f (x_k) + f (x_{k+1})],
$$
而在计算机编程时,通常使用的是下面这个形式:
$$
T_n = \frac {h}{2}\sum^{n-1}_{k=0}[f (x_k) + f (x_{k+1})] = \frac {h}{2}[f (a)+2\sum^{n-1}
{k=1} f (x_k) + f (b)]
$$

2. 复合辛普森公式

   复合辛普森公式的思路和复合梯形公式的思路一致,差别在于它的低阶求积公式使用了辛普森公式。辛普森公式是经过加权改造的梯形公式,集中在中间的点的权值更高,因此相比起简单的梯形公式,辛普森公式拥有更高的数值精度。
   将区间 &[a,b]& 分为 & n & 等分,在每个字区间 $[x_k, x_{k+1}]$ 上采用辛普森公式,记 $x_{k+1/2} = x_k + \frac {1}{2} h$,得:
$$
I = \int^b_af (x) dx = \sum^{n-1}{k=0}\int^{x{k+1}}{x_k} f (x) dx \ = \frac {h}{6}\sum^{n-1}_{k=0}[f (x_k) + 4f (x_{k+1/2}) + f (x{k+1})]
$$
在计算机编程时,通常写成下列形式:
$$
S_n = \frac {h}{6}[f (a) + 4\sum^{n-1}{k=0} f (x{k+1/2}) + 2\sum^{n-1}_{k=1} f (x_k) + f (b)]
$$

代码实现

1. 复合梯形公式

   根据上述推到的复合梯形公式,编写 matlab 代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
% 复合梯形公式 
function [result] = CompositeTrapezoid(a, b, n)
if (b < a)
c = b;
b = a;
a =c;
end
h = (b - a) /n; % 计算步长
result = myFun (a) + myFun (b);
for k = 1:n-1
x = a + k * h;
result = result + 2 * myFun (x);
end
result = (h / 2) * result;
end

function [y] = myFun(x)
if (x == 0)
y = 1;
else
y = sin(x) /x;
end
end

2. 复合辛普森公式

   根据上述推导的复合辛普森公式,编写 matlab 代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
% 复合辛普森公式 
function [result] = CompositeSimpson(a, b, n)
if (b < a)
c = b;
b = a;
a = c;
end
h = (b - a) /n; % 计算步长
result = myFun (a) + myFun (b);
for k = 1:n-1
x = a + k * h;
result = result + 2 * myFun (x);
end
for k = 1:n
x = a + (k-1)*h + 1/2 * h;
result = result + 4 * myFun (x);
end
result = result * (h / 6);
end

function [y] = myFun(x)
if (x == 0)
y = 1;
else
y = sin(x) /x;
end
end

小结

   采样点的数目分别为 5,9,17,33,分别进行简单的数值实验,得出结果如下图:

  1. 复合梯形公式
    ![数值积分结果分析 1](https://raw.githubusercontent.com/leungyukshing/Numerical-Computation-Methods/master/% E6%95% B0% E5%80% BC% E8% AE% A1% E7% AE%97% E7% AC% AC% E4% BA%8C% E6% AC% A1% E5% AE%9E% E9% AA%8C/Images/% E6%95% B0% E5%80% BC% E7% A7% AF% E5%88%86% EF% BC%88% E6% A2% AF% E5% BD% A2% EF% BC%89.png)
  2. 复合辛普森公式
    ![数值积分结果分析 2](https://raw.githubusercontent.com/leungyukshing/Numerical-Computation-Methods/master/% E6%95% B0% E5%80% BC% E8% AE% A1% E7% AE%97% E7% AC% AC% E4% BA%8C% E6% AC% A1% E5% AE%9E% E9% AA%8C/Images/% E6%95% B0% E5%80% BC% E7% A7% AF% E5%88%86% EF% BC%88% E8% BE%9B% E6%99% AE% E6% A3% AE% EF% BC%89.png)

** 分析 **:对于同一种方法来讲,采样点越多,积分值越准确;对比两种方法来说,复合辛普森公式更加准确,原因是它对中间的值分配更高的权重,使得数值积分更加准确。
   当然,除此之外还有著名的 ** 龙贝格求积公式 ** 和 ** 高斯求积公式 ** 等,这篇博客只涉及到很小的一部分,以后有机会继续补充,谢谢!


参考资料:

  1. 数值分析(第 5 版) 李庆扬,王能超,易大义 * 编 *

Welcome to my other publishing channels