数值积分方法：复合梯形公式、复合辛普森公式和龙贝格算法

本文将介绍三种常用的数值积分方法：复合梯形公式、复合辛普森公式和龙贝格算法，并通过 MATLAB 代码示例展示如何实现这些方法。

1. 定义函数

首先，我们需要定义一个函数来计算被积函数的值。这里以函数 sqrt(x) * log(x) 为例：

function y = fun41(x)
y = sqrt(x) .* log(x);
end

2. 用复合梯形公式和复合辛普森公式计算数值积分

接下来，我们将分别使用复合梯形公式和复合辛普森公式计算函数 fun41(x) 在区间 [0, 1] 上的积分值。

2.1 复合梯形公式

复合梯形公式通过将积分区间等分为若干个小区间，并在每个小区间上使用梯形公式近似计算积分，最后将所有小区间的积分值累加得到最终的积分值。

clear;
clc;
h = 0.001;     % 步长
n = 1/h; 
t = 0; 
for i = 1:n-1
    t = t + fun41(i*h);
end
T = h/2 * (0 + 2*t + fun41(1));
T = vpa(T, 10)

2.2 复合辛普森公式

复合辛普森公式与复合梯形公式类似，但它在每个小区间上使用辛普森公式近似计算积分。

s1 = 0;
s2 = 0;
for i = 0:n-1
    s1 = s1 + fun41(h/2 + i*h);
end
for i = 1:n-1
    s2 = s2 + fun41(i*h);
end
S = h/6 * (0 + 4*s1 + 2*s2 + fun41(1));
S = vpa(S, 10)

3. 使用龙贝格算法计算数值积分

龙贝格算法是一种基于递推公式的数值积分方法，它利用复合梯形公式的计算结果，通过迭代的方式逐步提高积分精度。

clear;
clc;
m = 16;
h = 1;
T(1) = (0 + fun41(1)) * h/2
for i = 2:m
    h = h/2;
    n = 1/h;
    t = 0;
    for j = 1:2:n-1
        t = t + fun41(j*h);
    end
    T(i) = T(i-1)/2 + h*t; % 梯形公式
end
for i = 1:m-1
    for j = m:i+1
        T(j) = 4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1);
        % 通过不断的迭代求得 T(j)，即 T 表的对角线上的元素。
    end
end
vpa(T(m), 10)

代码逐句解析

定义函数
```
function y = fun41(x)
y = sqrt(x) .* log(x);
end
```
代码定义了一个名为 fun41 的函数，它接收一个实数 x 作为参数，并返回 sqrt(x) * log(x) 的值。
复合梯形公式
```
h = 0.001;     % 步长
n = 1/h; 
t = 0; 
for i = 1:n-1
    t = t + fun41(i*h);
end
T = h/2 * (0 + 2*t + fun41(1));
T = vpa(T, 10)
```
- h 定义了步长，表示将积分区间 [0, 1] 等分为 n 个小区间，每个小区间的长度为 h。
- n 表示小区间的数量。
- t 用于累加每个小区间梯形面积的中间部分。
- for 循环计算每个小区间中间节点的函数值，并累加到 t 中。
- 最后的 T 计算了整个积分区间的积分值。
- vpa(T, 10) 将结果保留 10 位有效数字。
复合辛普森公式
```
s1 = 0;
s2 = 0;
for i = 0:n-1
    s1 = s1 + fun41(h/2 + i*h);
end
for i = 1:n-1
    s2 = s2 + fun41(i*h);
end
S = h/6 * (0 + 4*s1 + 2*s2 + fun41(1));
S = vpa(S, 10)
```
- s1 和 s2 用于累加每个小区间辛普森公式中不同权重的函数值。
- for 循环计算每个小区间不同节点的函数值，并累加到 s1 和 s2 中。
- 最后的 S 计算了整个积分区间的积分值。
- vpa(S, 10) 将结果保留 10 位有效数字。
龙贝格算法
```
m = 16;
h = 1;
T(1) = (0 + fun41(1)) * h/2
for i = 2:m
    h = h/2;
    n = 1/h;
    t = 0;
    for j = 1:2:n-1
        t = t + fun41(j*h);
    end
    T(i) = T(i-1)/2 + h*t; % 梯形公式
end
for i = 1:m-1
    for j = m:i+1
        T(j) = 4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1);
        % 通过不断的迭代求得 T(j)，即 T 表的对角线上的元素。
    end
end
vpa(T(m), 10)
```
- m 表示迭代次数，即 T 表的行数。
- h 定义了初始步长，并随着迭代次数的增加而减半。
- 第一个 for 循环计算 T 表的第一列元素，也就是用复合梯形公式计算不同步长下的积分值。
- 第二个 for 循环计算 T 表的对角线上的元素，通过迭代公式不断提高积分精度。
- 最后的 vpa(T(m), 10) 将最后一次迭代得到的积分值保留 10 位有效数字。

通过以上代码示例，我们可以清楚地看到三种数值积分方法的实现过程。在实际应用中，选择哪种方法取决于具体的应用场景和精度要求。