%复化辛普森公式
function T_n = C_t_r(f,n,a,b)
%f为函数,n为精度,a,b为区间端点
h=(b-a)/n;
%f=sqrt(x);
for i=0:n
    x(i+1)=a+i*h;
end
T_1=h/6*(f(a)+f(b));
T_2=0;
T_3=0;
for j=2:n
    F(j)=h/3*f(x(j));
    T_2=T_2+F(j);
end
for k=n+2:2*n+1
    F(k)=h/3*(f(x(k))+h/2);
    T_3=T_3+F(k);
end
%T_2=sum(F);
T_n=T_1+T_2+T_3;

复化辛普森公式求二重积分matlab源码
 
直接拷贝到matlab编辑器，傻瓜式操作。具体算法自行探究，网上都有，小编只提供代码。用的好的请加个关注，篱落～～成殇～～再次先行谢过。
 
%%%%%%%%%%     2020.6.5        %%%%%%%%%
%%%%%%%%%%复化Simpson公式求二重积分%%%%%%%%%
%%%%%%%%%%     Liu Deping      %%%%%%%%%
clear all;
%%%被积函数及积分上下限导入；
s=input('请输入函数表达式：f = ','s');
f=inline(s);
a = input('请输入积分变量x左边界a的值:');
b = input('请输入积分变量x右边界b的值:');
c = input('请输入积分变量y左边界c的值:');
d = input('请输入积分变量y右边界d的值:');
h1 = input('请输入积分变量x步长h1的值:');
h2 = input('请输入积分变量y步长h2的值:');
m=round((b-a)/h1);
n=round((d-c)/h2);
%%%系数表T，t[i,j]为复化梯形公式的系数；
T=zeros(m+1,n+1);
R=[1,4,1;4,16,4;1,4,1];
for p=0:m/2-1
    for q=0:n/2-1
        for i=1:3
            for j=1:3
                T(i+2*p,j+2*q)=T(i+2*p,j+2*q)+R(i,j);
            end
        end
    end
end
fprintf('系数表：');
T
%%%计算各二维节点函数值，并存放于矩阵F中；
F=zeros(m+1,n+1);
for i=0:m
    for j=0:n
       F(i+1,j+1)=f(a+i*h1,c+j*h2);    
    end
end 
%%%结果输出；
format long;
fprintf('各节点函数值；')
F
fprintf('复化Simpson公式计算结果；')
Snm=(b-a)*(d-c)/(9*m*n)*sum(sum(T.*F))
 
下面附上一个例题吧，感受一下代码的强大。
 
 只看第(2)问，利用辛普森公式计算积分。(关于第一问的代码可以去关注上一篇文章〈复化梯形公式公式求二重积分matlab源码〉)
 步骤:
 1.复制代码到matlab编辑器，点击运行;
 2.按照提示输入:
 
 3.得到结果(包括系数、节点函数值和积分结果):
 
 
 
 如果代码对你有帮助，请加一下关注，如果有问题可以私聊！最后，再次感谢喜欢我的小伙伴。

 （1）Matlab代码:
 f1:
 
if x==0
    y=1;
else
    %y=exp(-x^2);
    y=(sin(x))/x;
end
 
梯形公式：
 
a=0;
b=1;
tx=(b-a)/2*(f1(a)+f1(b));
tx=vpa(tx,5);
tx
 
simpson公式：
 
a=0;
b=1;
xps=(b-a)/6*(f1(a)+4*f1((a+b)/2)+f1(b));
xps=vpa(xps,5);
xps
 
（2）Matlab代码：
 复化梯形公式：
 
a=0;
b=1;
e=1e-6;
m=1;
h=b-a;
T(1)=h/2*(f1(a)+f1(b));
while 1
    h=h/2;
    sum=0;
    for j=1:2^(m-1)
        sum=sum+f1(a+(2*j-1)*h);
    end
    S(1)=T(1)/2+h*sum;
    for j=1:m
        S(j+1)=S(j)+(S(j)-T(j))/(4^j-1);
    end
    if(abs(S(m+1)-T(m))<=e)
        fprintf('求得最终解为：%f\n',S(m+1)); 
        fprintf('迭代次数为:%d\n',m);
        return
    end
    T=S;
    m=m+1;
end
 
复化辛普森公式：
 
a=0;
b=1;
e=1e-6;
N=20;
n=0;
h0=(b-a)/2;
S0=h0/3*(f1(a)+f1(b)+4*f1((a+b)/2));
while n<N
    n=n+1;
    h1=h0/2;
    sum1=0;
    sum2=0;
    for k=1:2^(n-1)
        sum1=sum1+f1(a+(4*k-2)*h1);
    end
    for k=1:2^n
        sum2=sum2+f1(a+(2*k-1)*h1);
    end
    S1=S0/2+h1/3*(-2*sum1+4*sum2);
    if(abs(S1-S0)<e)
        fprintf('求得最终解为：%f\n',S1); 
        fprintf('迭代次数为:%d\n',n);
        return
    end
      S0=S1;
      h0=h1;
end
 
龙贝格算法：
 
a=0;
b=1;
e=1e-6;
h=b-a;
T(1)=h.*(f1(a)+f1(b))./2;
m=1;
while(1)
    h=h/2;
    sum=0;
    for j=1:2^(m-1)
        sum=sum+f1(a+h*(2*j-1));
    end
    S(1)=T(1)/2+h*sum;
    for j=1:m
        S(j+1)=S(j)+(S(j)-T(j))/(4^j-1);
    end
    if(abs(S(m+1)-T(m))<=e)
        break;
    end
    T=S;
    m=m+1;
end
fprintf('求得解为：%f\n',S(m+1));
fprintf('迭代次数：%d\n',m);
 
自适应求积方法：
 
format short
a=0;
b=1;
e=1e-6;
p=[a,b];
p0=p;
ep=[e];
m=0;
q=0;
I=0;
while(1)
    n1=length(ep);
    n=length(p0);
    if n==1
        break;
    end
    h=p0(2)-p0(1);
    s0=h/6*(f1(p0(1))+4*f1(p0(1)+h/2)+f1(p0(1)+h));
    s1=h/12*(f1(p0(1))+4*f1(p0(1)+h/4)+f1(p0(1)+h/2));
    s2=h/12*(f1(p0(1)+h/2)+4*f1(p0(1)+h*3/4)+f1(p0(1)+h));
    if abs(s0-s1-s2)<=15*ep(1)
        I=I+s1+s2;
        p0=p0(2:n);
        if n1>=2;
            ep=ep(2:n1);
        end
        q=q+1;
    else
        m=m+1;
        p0=[p0(1),p0(1)+h/2,p0(2:n)];
        if n1==1
            ep=[ep(1)/2,ep(1)/2];
        else
            ep=[ep(1)/2,ep(1)/2,ep(2:n1)];
        end
        if q==0
            p=p0;
        else
            p=[p(1:q),p0];
        end
    end
end
fprintf('求得的解为：%f\n',I);
fprintf('分层数为%d\n',m);
p=p';
disp('分点为');p

 

 

  
 
 
数值积分的基本思想
 
由积分中值定理可知，在积分区间 
 

   内存在一点 
 
 

   ，成立
 
 

 

  
 
 

  
 
 

 

   式的几何意义即为：底为 
 
 

   而高为 
 
 

   的矩形的面积恰等于所求曲边梯形的面积 
 
 

   。因此，要想求出 
 
 

   式左端积分，我们只需要知道三个值： 
 
 

   即可，这里 
 
 

   是显然的，问题在于 
 
 

   的具体位置一般是不清楚的，从而 
 
 

   未知。我们暂且将 
 
 

   称为区间 
 
 

   上的平均高度，我们的目标就是寻求一种求出平均高度 
 
 

   的算法，这样我们的问题也就解决了。
 
 

 
首先我们想到的是用区间两端点的“高度” 
 

   与 
 
 

   的算术平均值作为平均高度 
 
 

   的近似值，从而导出下述求积公式
 
 

 

  
 
 

  
 
 

 

   式即为我们我们熟知的
 
 梯形公式。写出MATLAB函数为
 

  
function
 
 
而若用区间中点 
 

   的“高度” 
 
 

   近似取代平均高度 
 
 

   ，则可导出
 
 中矩形公式（简称
 矩形公式）
 

 

  
 
 

  
 
 
其MATLAB函数为
 

  
function
 
 
总结：更一般地，我们可以在区间 
 

   上适当选取某些节点 
 
 

   ，然后用 
 
 

   的加权平均得到平均高度 
 
 

   的近似值，这样构造出来的求积公式具有下列通式
 
 

 

  
 
 

  
 
 
其中 
 

   称为
 
 求积节点， 
 

   称为
 
 求积系数，也称作伴随节点 
 

   的
 
 权。
 
这种数值积分方法通常称为机械求积，其优势在于将积分求值问题归结为被积函数值的计算，很适合在计算机上使用。
 
代数精度
 
由于数值积分方法是近似方法，为保证精度，自然希望我们的求积公式能够对“尽可能多”的函数准确成立，这便提出了代数精度的概念。
 

  如果某个求积公式对于次数不超过 
  

    的多项式均能够准确成立，但对于 
  
  

    次多项式就不准确成立，则称该求积公式具有
  
  
  

    次代数精度
  
  （或
  代数精确度）。
 
 
前面提到的梯形公式 
 

   和矩形公式 
 
 

   均具有一次代数精度。
 
 
求积公式的收敛性与稳定性
 

  在求积公式 
  

    中，若 
  
  

   
  
  

    其中 
  
  

    ，则称求积公式 
  
  

    是
  
  收敛的。
  
对任给 
  

    ，若 
  
  

    ，只要 
  
  

    就有 
  
  

   
  
  

    成立，其中 
  
  

    为计算 
  
  

    时产生误差 
  
  

    后实际得到的值，即 
  
  

    ，则称求积公式 
  
  

    是
  
  稳定的。
 
 
定理 
 

  
 
 若求积公式 
 

  
 
 中系数 
 

  
 
 ，则此求积公式是稳定的。
 
牛顿-柯特斯(Newton-Cotes)公式
 
将积分区间 
 

   划分为 
 
 

   等份，此时步长为 
 
 

   ，选取等距节点 
 
 

   构造出的插值型求积公式
 
 

 

  
 
 

  
 
 
称为牛顿-柯特斯(Newton-Cotes)公式，其中 
 

   称为
 
 柯特斯系数。可以通过下式确定
 

 

  
 
 

  
 
 
当 
   

     时， 
   
   

     这时的求积公式便是之前的梯形公式 
   
   

     。
   
当 
   

     时， 
   
   

    
   
   

     此时对应的求积公式便是
   
辛普森(Simpson)公式： 
   

    
   
   

     写成MATLAB函数为
   
 

  
function
 
 
当 
   

     时的牛顿-柯特斯公式则特别称为
   
柯特斯公式，其形式为 
   

    
   
   

     这里 
   
   

     。
   
 
其MATLAB代码为
 

  
function
 
 
复合求积公式
 
由于牛顿-柯特斯公式在 
 

   时不具有稳定性，故不可能再通过提高阶的方法来提高求积精度。
 
 复合求积法便是通过把积分区间分成若干个子区间（通常是等分），再在每个子区间上使用低阶求积公式的方法来提高精度的。
 
其实细心的同学可以发现，前面我所展示的MATLAB代码使用的便是对每个求积公式使用复合法完成的。常用的有复合梯形公式和复合辛普森公式。
 
龙贝格(Romberg)求积公式
 
前面介绍的复合求积方法可提高求积精度，如若精度仍不够，则可通过将步长逐次减半的方式来提高精度。如对复合梯形公式可导出其递推公式
 

 

  
 
 

  
 
 
其中 
 

   表示在 
 
 

   基础上步长 
 
 

   减半后的积分值。
 
 
定理 
 

  
 
  设 
 

   ，则有
 
 
 

 

  
 
 
 

   其中
 
  
 

   ，
 
 系数 
 

   与 
 
 
 

   无关。
 
 
 
从 
 

   式可以看出， 
 
 

   是 
 
 

   阶，若用 
 
 

   代替 
 
 

   ，有
 
 

  
 
 

  
 
 
再用 
 

   乘 
 
 

   式再减去 
 
 

   式后再除以 
 
 

   得
 
 

 

  
 
 

  
 
 
这里 
 

   是与 
 
 

   无关的系数。
 
 
从 
 

   式可看出，用近似积分值 
 
 

   ，其误差阶为 
 
 

   ，显然误差阶是提高了。类似这种将计算 
 
 

   的近似值的误差阶由 
 
 

   提高到 
 
 

   的方法称为
 
 外推算法，也称为
 理查德森(Richardson)外推算法。
 

  这是“数值分析”中一个重要的技巧，只要真值与近似值的误差能够表示成 
  

    的幂级数，如 
  
  

    式所示，都可以使用外推算法，提高精度。
  
 
 
龙贝格(Romberg)算法
 

 

   次外推加速为
 
 

 

  
 
 

  
 
 
余项为
 

 

  
 
 

  
 
 
此方法常称为理查德森外推加速方法。
 
设用 
 

   表示二分 
 
 

   次后求得的梯形值，且以 
 
 

   表示序列 
 
 

   的 
 
 

   次加速值，则由递推公式 
 
 

   可得
 
 

 

  
 
 

  
 
 
上次则称为龙贝格求积算法，计算过程如下：
 
取 
   

     ，求 
   
   

     。令 
   
   

     ( 
   
   

     记区间 
   
   

     的二分次数)。
   
求梯形值 
   

     ，即按递推公式 
   
   

     计算 
   
   

     。
   
求加速值。
若 
   

     (预先给定的精度)，则终止计算，并且 
   
   

     ；否则令 
   
   

     转 
   
   

     继续计算。
   
 
MATLAB代码为
 

  
function
 
 
自适应积分方法
 
复合求积方法通常适用于被积函数变化不太大的积分，如果在求积区间中被积函数变化很大，有的部分函数值变化剧烈，另一部分却变化平缓。这时如果统一将区间等分再用复合求积公式计算积分将会导致计算量很大，我们想实现的是在满足误差要求的前提下，对变化剧烈部分将区间细分，而平缓部分则可使用大步长，也即针对被积函数在区间上不同情形采用不同的步长，使得再满足精度前提下积分计算工作量尽可能小。其算法技巧是在不同区间上预测被积函数变化的剧烈程度确定相应步长。这就是自适应积分方法的基本思想。
 
下面是一个基于复合辛普森法的自适应积分算法的MATLAB代码：
 

  
function
 
 
下面我们计算积分 
 

   ，可输入下列语句调用上述MATLAB函数进行计算
 
 

  
clear

 
 
结果为
 

  
I 
 
 
高斯求积公式
 
下面研究带权积分 
 

   （ 
 
 

   为权函数）的求积公式
 
 

 

  
 
 

  
 
 
下面看定义
 

  如果求积公式 
  

    具有 
  
  

    次代数精度，则称其节点 
  
  

    为
  
  高斯点，相应公式 
  

    称为
  
  高斯型求积公式。
 
 
这里主要包括四种高斯型求积公式，即
 
高斯-勒让德求积公式
高斯-切比雪夫求积公式
高斯-拉盖尔求积公式
高斯-埃尔米特求积公式
 
多重积分
 
这里主要说一下多重积分计算的思路，归结为一句话，就是：多重积分化累次积分，再由里到外使用数值积分公式进行求积计算。