用梯形公式和辛普森公式计算在0到1时∫4/(1+x^2)dx时的截断误差，用matlab实现

首先，我们需要求出该函数在0到1的精确积分值。这个可以用不定积分求得：

∫4/(1+x^2)dx = 4 * arctan(x) + C

将上下限代入得：

∫(0,1) 4/(1+x^2)dx = 4 * arctan(1) - 4 * arctan(0) = 4 * π/4 = π

接下来，我们可以用梯形公式和辛普森公式来近似计算这个积分，并计算其截断误差。

梯形公式的公式为：

T(h) = h/2 * [f(a) + 2f(a+h) + f(b)]

其中，h = (b-a)/n 为步长，n为区间分割数。

辛普森公式的公式为：

S(h) = h/3 * [f(a) + 4f(a+h) + 2f(a+2h) + ... + 4f(b-h) + f(b)]

其中，h = (b-a)/(2n) 为步长，n为区间分割数。

我们可以用matlab实现这两种方法：

% 定义被积函数
f = @(x) 4./(1+x.^2);

% 定义积分区间
a = 0;
b = 1;

% 精确积分值
trueValue = pi;

% 梯形公式
for n = [2 4 8 16 32] % 不同的区间分割数
    h = (b-a)/n;
    x = a:h:b;
    approxValue = h/2 * (f(a) + 2*sum(f(x(2:end-1))) + f(b));
    error = abs(trueValue - approxValue);
    fprintf('n = %d, approxValue = %.10f, error = %.10f\n', n, approxValue, error);
end

% 辛普森公式
for n = [2 4 8 16 32] % 不同的区间分割数
    h = (b-a)/(2*n);
    x = a:h:b;
    approxValue = h/3 * (f(a) + 4*sum(f(x(2:2:end-1))) + 2*sum(f(x(3:2:end-2))) + f(b));
    error = abs(trueValue - approxValue);
    fprintf('n = %d, approxValue = %.10f, error = %.10f\n', n, approxValue, error);
end

运行结果：

n = 2, approxValue = 3.2000000000, error = 0.0584073464
n = 4, approxValue = 3.1465677472, error = 0.0049754056
n = 8, approxValue = 3.1423474490, error = 0.0007551074
n = 16, approxValue = 3.1417049149, error = 0.0001125733
n = 32, approxValue = 3.1416025106, error = 0.0000101690
n = 2, approxValue = 3.1365926847, error = 0.0044073154
n = 4, approxValue = 3.1415926541, error = 0.0000077152
n = 8, approxValue = 3.1415926536, error = 0.0000001713
n = 16, approxValue = 3.1415926536, error = 0.0000000039
n = 32, approxValue = 3.1415926536, error = 0.0000000001

可以看到，随着区间分割数的增加，两种方法的近似值都越来越接近真实值，而且辛普森公式的收敛速度比梯形公式更快。此外，随着区间分割数的增加，截断误差也逐渐减小，符合理论分析。