作为研究生的入门课，数值计算的大作业算是所有研究生开学的重要编程作业。

   针对数值积分的编程，我把梯形、辛普森与龙贝格方法在MATLAB中编写的程序放在文章最后了，需要的同学自取。

PS：附录中的程序并不是一整个M文件，同学参考的时候一定要注意将程序拆解为如下图所示的几个m文件，并在同一工作目录下运行

下文为作业详解

题目：计算积分，要求绝对误差限ε≤10-6次方 ，完成以下工作：

    1.使用复合梯形公式计算，输出等分次数及积分的近似值

梯形公式数学原理：

     复合梯形法要想控制积分精度，可以采用如下的方法：设积分区间已经划分为n个子区间，这时再把区间划分更细，给出新的积分结果，如果通过复化梯形法计算的积分与积分准确值的绝对误差比给定的误差限小的话，则停止细化，否则继续增加积分区间。

 具体td_fun程序在附录，运行如下命令

[s_1,num_1]=td_fun(0,pi,1e-6);

得到运行结果

2. 使用复合 Simpson 公式计算，输出等分次数及积分的近似值

辛普森数学公式：

复合辛普森公式

具体S_fun程序也在附录，运行如下命令

[s_1,num_1]=S_fun(0,pi,1e-6);

得到运行结果

 结论：由上图 两个运行结果所示，在达到相同运算精度的条件下，复合辛普森公式达到的速度更快，运行效果更好。

3.使用龙贝格算法计算，列表显示积分的近似值

龙贝格数学原理：

     龙贝格求积公式也称为逐次分半加速法。是数值计算方法之一，用以求解数值积分。是在梯形公式、辛普森公式和柯特斯公式之间关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法，在不增加计算量的前提下提高了误差的精度

    在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。

执行如下命令

K=rb_fun(0,pi,1e-6);

 本题运算结果如下

 结论：由上图所示，实际进行的外推步骤很少就已经达到了所求精度，因此龙贝格算法比之前两个数值求积效果更好，收敛更快

附录：Matlab程序(复合梯度法，辛普森公式法，龙贝格法)

复合梯度法

function  f=chang_t_f(x)
 f=exp(x)*cos(x); %定义积分函数



function [jifen,num] = td_fun(a,b,tol)%a,b为积分区间
%tol为积分精度，a,b为积分区间
n=1;%初始等分区间数
h=b-a;%初始的步长
fun_1=- exp(pi)/2 - 1/2;%matlab计算出的准确的积分值
fun_2=(chang_t_f(a)+chang_t_f(b))*(h/2); %调用方程函数
while abs(fun_2-fun_1)>tol 
 n=n+1;   h=(b-a)/n;   fun_2=0;
for i=0:n-1
      x=a+i*h;    x1=x+h; 
fun_2=fun_2+(h/2)*(chang_t_f(x)+chang_t_f(x1));%复化梯形公式
    end
end
%如果利用复化梯形公式得到的积分值与准确值的绝对误差小于给定的精度则停止细化区间
jifen=fun_2 %积分结果      num=n; %区分划分细度

辛普森公式法

function f=chang_xps_f(x)
f=exp(x)*cos(x);定义积分函数



function [jifen,num] = s_fun(a,b,tol)
%a,b为积分区间，tol为积分精度
n=1;h=b-a;
int_1=- exp(pi)/2-1/2;
int_2=(chang_xps_f(a)+chang_xps_f(b))*(h/2); %调用方程函数
%如果利用复化辛普森公式得到的积分值与准确值的绝对误差小于给定的精度则停止细化区间
while abs(int_2-int_1)>tol
    n=n+1;    h=(b-a)/n;   int_2=0;
   for i=0:n-1
        x=a+i*h;   x1=x+h;
  int_2=int_2+(h/6)*(chang_xps_f(x)+4*chang_xps_f((x+x1)/2)+chang_xps_f(x1));%复化辛普森公式
    end
end
jifen=int_2; %积分结果
num=n; %区间划分细度

龙贝格法​​​​​​​

function f=romberg_f(x)
f=exp(x).*cos(x);定义积分函数

 function s=romberg(a,b,tol) 
%romberg积分方法函数
%a,b为积分区间
%tol为误差限
%s为最后积分值
n=1;   h=b-a;
delt=1; %设置设计误差初值
p=- exp(pi)/2 - 1/2;%matlab计算的积分值
x=a;   k=0;
R=zeros(4,4);%定义R为4×4的矩阵
R(1,1)=h*(romberg_f(a)+romberg_f(b))/2;%矩阵第一行第一列的数值
while delt>tol%%如果绝对误差小于给定的精度则进入循环
    k=k+1;   h=h/2;
    s=0;
    for j=1:n
        x=a+h*(2*j-1);
        s=s+romberg_f(x);
    end
    R(k+1,1)=R(k,1)/2+h*s;
    n=2*n;
    for i=1:k
        R(k+1,i+1)=((4^i)*R(k+1,i)-R(k,i))/(4^i-1);%龙贝格算法公式
    end
    delt=abs(R(k+1,k+1)-p);
end
s=R(k+1,k+1);