梯形求积公式 和 复合梯形求积公式 Matlab 实现
 
梯形求积公式
 
 
 
仅使用区间两点$x_1,f(x_1),x_2,f(x_2)$ 
 组成的梯形面积S代替$\int_{x_1}^{x_2} f(x)dx$的近似方法 
  
 
  $$\int_{x_1}^{x_2} f(x)dx \approx  S = \frac{x_2 - x_1}{2} \times (f(x_2) - f(x_1))$$  
 
 
 
复合梯形求积公式
 
 
 
将求积区间[a,b]分为n个区间，每个区间步长为h($h = \frac{b-a}{n}$)在每个区间求梯形积分 
  
 Si为第i个梯形的面积 
 
  $$\int_{x_1}^{x_2} f(x)dx \approx  \sum_{i = 0}^n Si = \frac{h}{2}  \int (f(a) + \sum_{k = 1}^{n-1}f(x_k) +f(b))$$ 
 
 
 
说明
 
 
 
上述公式是我用mathjax写的，如有错误请联系我修正 
 敬请指正 
 概述省略了部分推导过程，请查阅详细推导资料
 
 
Matlab 实现代码
 
梯形求积公式
 
 
 
将该函数存为m文件
 
 
function res = Trapezium(f,a,b)
    format long;
    if b < a
        c = b;
        b = a;
        a = c;
    end
    res = (b-a)*(f(a)+f(b))/2;
 
 
 
调用下面语句测试函数
 
 
f = inline('sin(x)','x')
Trapezium(f,0,pi/2)
 
复合梯形求积公式
 
 
 
将该函数存为m文件
 
 
function res = ComTrapezium(f,n,a,b)
    format long;
    if b < a
        c = b;
        b = a;
        a = c;
    end
    h = (b-a)/n;
    d = f(a);
    for i = a+h:h:b-h 
        d = d + (2 * f(i));
    end
    d = d + f(b);
    res = (d * h / 2);
 
 
 
调用下面语句测试函数
 
 
f = inline('sin(x)','x')
ComTrapezium(f,4,0,pi/2)

该公式可以理解为将曲线划分为若干个梯形，数量越多，结果越精确。
 
 
这个程序实现有点太过于基础，这里就不一一叙述，直接上代码。不过要改变积分函数的话就调用函数改为新的函数即可。
 
void func(int n)//调用求f(xi)的函数
{
    for (int i = 0; i <= n; i++)
    {
        f[i] = 8 + x[i] * x[i] * x[i];
    }
}
 
附上实现程序（Visual Studio 2019）
 
#include <iostream> 
#include <iomanip>//保留小数

double x[100];//定义节点的x坐标
double f[100];//定义节点的函数值f
double h;//定义步长
void func(int n);//调用求f(xi)的函数
double i_f;//定义复合梯形求解结果
double func_i_f(int n);//调用复合梯形求积函数

int main()
{
    using namespace std;
    int n;
    cout << "输入等分数" << endl;
    cin >> n;
    cout << "输入x的上限值" << endl;
    cin >> x[n];
    cout << "输入x的下限值" << endl;
    cin >> x[0];//赋值模块
    h = (x[n] - x[0]) / n;//求步长h
    for (int i = 1; i <= n - 1; i++)
    {
        x[i] = x[0] + i * h;
    }//求每个节点x坐标值
    func(n);
    func_i_f(n);
    cout << setiosflags(ios::fixed) << setprecision(5);//输出结果保留5位小数
    cout << "f(x)=8 + x^3在x区间上的积分值为：" ;
    cout << i_f << endl;
    return 0;
}

void func(int n)//调用求f(xi)的函数
{
    for (int i = 0; i <= n; i++)
    {
        f[i] = 8 + x[i] * x[i] * x[i];
    }
}

double func_i_f(int n)//调用复合梯形求积函数
{
    i_f = 0;
    for (int i = 0; i <= n - 1; i++)
    {
        i_f += (f[i] + f[i + 1]) * h / 2;
    }
    return i_f;
}
 
结果输出
 
 
后话
 
关于《数值分析》与C++联动小项目就到这里结束了!接下来的一个小小工作就把高斯迭代那一篇的程序推向n阶。

 
本文展示复合梯形求积公式的数值积分效果。 对于闭区间上的一般函数，利用复合梯形求积公式具有二阶精度。 但是对于周期函数，特别是解析函数，以及 
 

   上迅速衰减的函数，复合梯形公式具有几何收敛阶。
 
 
实轴上的积分 
 
我们首先考虑定义在实轴 
 

   上的函数 
 
 

   的积分
 
 

 

  
 
 
如果 
 

   是光滑，并且衰减的，那么使用复合梯形公式积分将有几何阶的收敛性。
 
 
具体来说，数值积分公式为
 

 

  
 
 
引理： Suppose 
 

   is analytic in the strip 
 
 

   for some 
 
 

  . Suppose further that 
 
 

   uniformly as 
 
 

   in the strip, and for some 
 
 

  , it satisfies 
 
 

  
 
 
for all 
 

  . Then, for any 
 
 

  , 
 
 

   exists and satisfies 
 
 

  
 
 
见《The Exponentially Convergent Trapezoidal Rule》。
 
我们首先来测试一下求积公式对“不那么光滑”函数的积分效果
 
例子一
 
考虑有限区间上的函数
 

 

  
 
 
积分的精确值为
 

 

  
 
 
积分节点为等距节点
 

 

  
 
 
步长为 
 

   。 复合梯形积分公式为（因为函数在端点处值相同（周期））
 
 

 

  
 
 
Matlab程序
 

  
%% 复合梯形公式求积函数
%  4-x^2, x in [-2,2]

clear all
close all
p_set = 2.^(4:10);
w = @(z) 4-z.^2;
ErrInt = zeros(size(p_set));
% 精确值
ExactInt = integral(w,-2,2);

figure; set(gcf,'unit','centimeters','position',[10,10,24,12]);
for ii = 1:length(p_set)
    p = p_set(ii);
    xk = linspace(-2,2,p+1);
    xk = xk(2:end);
    ErrInt(ii) = abs(4*mean(w(xk))-ExactInt);
end
loglog(p_set,ErrInt,'-.o','LineWidth',1.5);
hold on;
legend('2阶')

 
 

 

  
 
 
横坐标是节点个数
 

  ，在loglog图像下，数值积分误差随着节点增加而二阶收敛。
 
 
例子二
 
我们选择
 

 

  
 
 
复合梯形公式的收敛为 
 

 

  
 
 
左侧的图中，数值积分误差随着节点增加几何收敛。此时节点数目较少。
右侧图中，数值积分误差随着节点增加 
   

     阶收敛。此时节点数目较多。
   
 
虽然上面例子中的函数都是给定区间上的周期函数，但它们周期化之后并不是解析函数。另一个角度则是，积分收敛或与它们的傅里叶系数的衰减有关。
 
例子三
 
被积函数选择为
 

 

  
 
 
其中 
 

   和 
 
 

   是参数。精确的积分值为
 
 

 

  
 
 
考虑 
 

   ， 
 
 

   的参数情形。
 
 
之前对于这类无穷区间上的积分，我们采用的是先截断到固定区间上的做法。 此处，我们将数值积分的参数设置为步长 
 

   以及节点数 
 
 

   。 那么之前固定区间就等价于 
 
 

   。
 
 
下面是复合梯形公式的误差结果关于 
 

   的收敛，左右分别是不同的截断 
 
 

   的选取方法。
 
 

 

  
 
 
有以下两点观察：
 
误差在 
   

     达到 
   
   

     左右剧烈减小
   

   

     和 
   
   

     没有明显区别
   
 
第一点可以通过分析 
 

   的Fourier变换的带宽来解释。
 
 

 

   是个Gauss函数，带宽为
 
 

 

  
 
 
Fourier变换 
 

   依旧是Gauss函数，带宽为
 
 

 

  
 
 
如果复合梯形公式的步长 
 

   满足 
 
 

  ， 即 
 
 

  ， 则无穷项的复合梯形公式几乎是精确的。
 
 
因此，在这种情况下，有限截断的复合梯形公式的误差主要来自于截断项。 但是被积函数 
 

   的带宽 
 
 

   非常小，因此截断的误差几乎没有影响。
 
 
最后，在右侧的图像中，从右往左看，倒数第二个点没有画出来，这是因为在Matlab中计算得到的误差为 
 

   的原因。
 
 
围道积分的应用 
 
给定函数
 

 

  
 
 
延拓到复平面，它有两个单极点 
 

   和 
 
 

   。
 
 
我们希望已知极点，求出它们对应的系数（留数）。
 
在Matlab中，我们可以通过符号计算函数 limit 来执行这个操作。通过帮助文件可得
 

  
limit(F,x,a) takes the limit of the symbolic expression F as x -> a.

 
 
具体举例如下：
 

  
syms x
w = 1/(x^2-x-2);
coef1 = limit(w*(x-2),x,2);
disp(coef1)

 
 
借助于复合梯形公式，我们也可以利用围道积分来计算
 

 

  
 
 
复合梯形公式得到的近似为
 

 

  
 
 
Matlab测试程序
 

  
p_set = 4:10;
w = @(z) 1./((z-2).*(z+1));
ErrInt = zeros(size(p_set));
r_set = [1,0.5,0.25];

figure; set(gcf,'unit','centimeters','position',[10,10,24,12]);
for jj = 1:length(r_set)
    r = r_set(jj);
    for ii = 1:length(p_set)
        p = p_set(ii);
        k = 1:p;
        ErrInt(ii) = abs(real(mean(r*w(2+r*exp(2*pi*1i*k/p)).*exp(2*pi*1i*k/p)))-1/3);
    end
    semilogy(p_set,ErrInt,'LineWidth',1.5);
    hold on;
end
legend('r=1','r=0.5','r=0.25')

 
 

 

  
 
 
分别测试了三种半径对应的围道，其中 
 

   情况下， 
 
 

   个采样点（意味着数值积分中是10次乘法和加法）的精度达到了 
 
 

   。
 
 
相比来说，符号运算 limit 函数不需要我们费心。但是如果计算资源受限，围道积分是一种选择。

【C语言基础】利用复合梯形求积公式计算定积分
 
一、复合梯形求积公式
 
这是数值分析中一种求解定积分的近似方法。适用于被积函数的原函数不能用初等函数表示的情况。
 
基本思路
 
将被积函数 f（x）与x轴围成的区域分成n个梯形，把n个梯形面积求和得到积分的近似值。若精度不满足需要，则可以将每个区间再等分一次，得到2n+1个等分区间，然后再求和，直到精度满足需要。
 
公式
 
$$\underset{a}{\overset{b}{\int }}f(x)dx\approx \frac{h}{2}[f(a)+f(b)+2\sum _{i=1}^{n-1}f(x_i)]$$
 其中：
 
n 为等分区间数
 
h=（b-a）/2 (积分步长)
 
C语言实现算法
 
以下面积分为例：
 
$$\underset{0}{\overset{1}{\int }}1+e^{-}xsin4xdx$$
 
首先写出计算被积函数值的函数
 
double f(double x)
{
    return 1 + exp(-x)*sin(4x);
}
 
然后写出求积分近似解的表达式
 
double caculate(double a,double b)
{
    double y1,y2,h;
    int n=1,i;
    h = a-b;
    y1 = h*(f(a)+f(b))/2;
    do
    {
       y2 = y1;
       n = n*2;
       h = (a-b)/n;
       y1 = f(a)+f(b);
       for(i=1;i<n:i++)
             y1 = y1 + f(a + i*h);
        y1 = h*y1/2;
    }while(fabs(y2-y1)>1.0e-5));
    return y2;
}
 
其中frab（）是用来计算绝对值的函数