数值积分辛普森法的推导以及MATLAB实现

数值积分在船体计算中很常用。以前只会用, 花点时间搞明白它近似计算的原理。(我没有完全搞懂, 本文可能有错误)

我个人理解的数值积分, 就是把原函数 \( f(x) \) 在 \( f(0) \) 处用麦克劳林级数展开, 用多项式来拟合原函数。多项式的次数为\( 1 \) 时, 就是梯形法; 为 \( 2 \) 时就是辛普森一法; 为 \( 3 \) 时就是辛普森二法。 多项式的次数越高, 拟合精度越高。

数值积分的一般形式

对于函数 \( f(x) \) 在区间 \( [a,b] \) 上的积分有

$$A = \int_0^L f(x) dx = \int_0^L y dx \tag{1} \label{1}$$

假定用 \( n \) 次多项式曲线对 \( f(x) \) 进行拟合。将该多项式曲线在 \( f(0) \) 处用麦克劳林级数展开, 则得

$$f(x) \approx \sum_{i = 0}^{n} \frac{1}{i!} f^{(i)}(0) x^i \tag{2} \label{2}$$

将式 \( \eqref{2} \) 代入式 \( \eqref{1} \) 可得

$$\begin{aligned} A &= \int_0^L f(x) dx \\\\ &\approx \frac{1}{0!} f(0) \int_0^L dx + \frac{1}{1!} f^{'}(0) \int_0^L x dx + \cdots + \frac{1}{n!} f^{(n)}(0) \int_0^Lx^ndx \\\\ &= \sum_{i = 0}^{n} \frac{L^{i+1} }{(i+1)!}f^{(i)}(0) \end{aligned} \tag{3} \label{3}$$

令面积 \( A \) 为

$$A = \sum_{i = 0}^{n} k_iy_i \tag{4} \label{4}$$

上式中的 \( y_i \) 可由式 \( \eqref{2} \) 得到

$$\begin{aligned} y_0 &= f(0) \\\\ y_1 &= f(0) + f^{'}(0)x_1 + \cdots + \frac{1}{n!} f^{(n)}(0)x_1^n \\\\ & \cdots \cdots \\\\ y_n &= f(0) + f^{'}(0)x_n + \cdots + \frac{1}{n!} f^{(n)}(0)x_n^n \end{aligned} \tag{5} \label{5}$$

将式 \( \eqref{5} \) 代入式 \( \eqref{4} \), 得 (注意:如果看不出来, 可以先在式 \( \eqref{5} \) 左右同时乘以\( k_i \) 然后再把 \( f^{(i)}(0) \) 提出来 )

$$A = (k_0 + k_1 + \cdots + k_n)f(0) + \\\\ \frac{1}{1!} (k_1x_1 + k_2x_2 + \cdots + k_nx_n) f^{'}(0) + \\\\ \cdots + \\\\ \frac{1}{n!} (k_1x_1^n + k_2x_2^n + \cdots + k_nx_n^n) f^{(n)}(0) \tag{6} \label{6}$$

将式 \( \eqref{6} \) 与式 \( \eqref{3} \) 相比较, 可以得到

$$k_0 + k_1 + \cdots + k_n = L \\\\ k_1x_1 + k_2x_2 + \cdots + k_nx_n = \frac{L^2}{2} \\\\ \cdots \\\\ k_1x_1^n + k_2x_2^n + \cdots + k_nx_n^n = \frac{L^{n+1}}{n+1} \\\\ \tag{7} \label{7}$$

至此得到了数值积分的一般形式。通过式 \( \eqref{7} \) 解得 \( k \) 的值, 再通过式 \( \eqref{4} \) 就可以求得面积 \( A \) , 即数值积分。

梯形法

取 \( n = 1 \) , 以直线段近似地代替 \( y = f(x) \) 的曲线, 坐标值只有 \( y_0 \) 和 \( y_1 \) , 区间长度 \( L = x_1 - x_0 \) , 且 \( x_1 = L \)。

由式 \( \eqref{7} \) , 可得

$$\left\{ \begin{array}{l} k_0 + k_1 = L \\\\ k_1x_1 = \frac{L^2}{2} \end{array} \right.$$

解得

$$\left\{ \begin{array}{l} k_0 = \frac{L}{2} \\\\ k_1 = \frac{L}{2} \end{array} \right.$$

代入式 \( \eqref{4} \) , 可得

$$A = \frac{L}{2} (y_0 + y_1)$$

辛普森一法

步骤和梯形法一样。取 \( n = 2 \) , 以二次抛物线近似地代替 \( y = f(x) \) 的曲线, 坐标值有 \( y_0 \) 、 \( y_1 \) 和 \( y_2 \) , 区间长度 \( L = x_2 - x_0 \) , 且 \( x_1 = L/2 , x_2 = L \)。

由式 \( \eqref{7} \) , 可得

$$\left\{ \begin{array}{l} k_0 + k_1 + k_2= L \\\\ k_1x_1 + k_2x_2 = \frac{L^2}{2} \\\\ k_1x_1^2 + k_2x_2^2 = \frac{L^3}{3} \end{array} \right.$$

解得

$$\left\{ \begin{array}{l} k_0 = \frac{L}{6} \\\\ k_1 = \frac{4L}{6} \\\\ k_2 = \frac{L}{6} \end{array} \right.$$

代入式 \( \eqref{4} \) , 可得

$$A = \frac{L}{6} (y_0 + 4y_1 + y_2)$$

一般使用的时候, 会把积分区间等分为偶数份, 每三个点用一次辛普森一法。

MATLAB 程序实现

MATLAB 中的 trapz 函数是梯形法的实现, 但是对于辛普森法没有对应的函数。quad 函数采用自适应辛普森法, 但是好像只能对函数对象进行积分, 而不能对离散的函数值进行积分。
所以我自己实现了一个 simpson 函数, 用法和 trapz 函数类似。但是多了个使用限制, 即必须把函数等分成偶数份 (输入的 y_vec 必须有奇数个值), 这个限制是为了方便生成辛普森系数。

使用方法:

clear;clc;
% 辛普森一法
res1 = 10 * simpson([1 2 3 4 5].^2);
% 梯形法
res2 = 10 * trapz([1 2 3 4 5].^2);
% 自适应数值积分
res3 = 10 * integral(@(x)x.^2 ,1,5);
fprintf('辛普森一法结果: %.2f \n', res1);
fprintf('梯形法结果: %.2f \n', res2);
fprintf('自适应数值积分结果: %.2f \n', res3);
% 执行结果: 
% 辛普森一法结果: 413.33 
% 梯形法结果: 420.00 
% 自适应数值积分结果: 413.33 

simpson 函数的实现:

function [ result ] = simpson( y_vec , d_num)
%SIMPSON 适用于离散数据的辛普森一法
%   y_vec : 等间距 y 值向量
%   d_num : 间距大小 (默认为 1)

%% 适用条件判断
if (nargin == 1)
    d_num = 1;
end
if length(y_vec) < 3
    error('y_vec 的长度太小了');
end
if mod(length(y_vec),2) == 0
    error('y_vec 的长度必须是奇数');
end
if d_num <= 0
    error('d_num 必须大于0');
end

%% 常用量
% y_vec 的长度
Len = length(y_vec);

%% 生成系数向量 M
M = zeros(1,Len);
if d_num == 3
    M = [1 4 1];
else
    M(1) = 1;
    M(2) = 4;
    for i = 1: (Len - 4)
        if mod(i,2) == 1
            M(2 + i) = 2;
        else
            M(2 + i) = 4;
        end
    end
    M(Len - 1) = 4;
    M(Len) = 1;
end

%% 计算结果
t = y_vec .* M;
result = (1/3) * sum(t(:)) * d_num;

end