一、实验目的及题目
1. 实验目的：
(1) 学习用复化辛普森公式及龙贝格方法求解积分并掌握这种方法。
(2)了解这些辛普森公式及龙贝格方法的概念，参考课本写出用复化辛普森算法以及龙贝格方法计算目标题目的程序，在matlab 中实现，并用matlab 内置的函数计算出结果，并提出存在的问题。
2. 题目：
利用复化辛普森公式和龙贝格方法计算下列积分： (1)dx e x ?-5
.002
(2)dx x x ?20
2sin )2sin(cos π
二、实验用仪器设备、器材或软件环境 计算机、matlab 软件。
三、实验原理、程序框图、程序代码 1.实验原理：
根据微积分学基本定理，若被积函数f(x)在区间[a,b]上连续，只要能找到f(x)的一个原函数F(x)，便可利用牛顿-莱布尼茨公式求得积分值。但会经常遇到如下问题：找不到用初等函数，找到了原函数，但因表达式过于复杂而不便计算等等。此时则不能用牛顿-莱布尼茨公式，因此有必要研究如下公式。 1)复化求积公式及原理
由于高阶插值的不稳定性，为了提高计算积分的精度，可把积分区间分为若干个
小区间，将()I f 写成这些小区间上的积分之和，然后对每一个小区间上的积分应用到辛普森公式，或柯特斯公式，并把每个小区间上的结果累加，所得到的求积公式就称为复化求积公式。
辛普森公式的数值积分公式为：
?
+++-≈
b
a
b f b a f a f a b dx x f )]()2
(
4)([6
)(

%下面除了龙贝格外，其他均以此fun函数作为被积函数
%梯形、辛普森、科特斯：都是已知积分上下限和分段数，求解积分近似值。
function a = fun(x)
    a = cos(x*x);
end

一、复化梯形：

function Tn = ComplexTrap(xL, xR, n)
% xL,xR：积分区间左右端点，n：分段数
% Tn ：复化梯形积分结果
h = (xR-xL)/n;%求步长
Tn = 0;
x = xL:h:xR;
for k = 1:n
    Tn = Tn + fun(x(k)) + fun(x(k+1));
end 
Tn = Tn*h/2;
Tn = vpa(Tn,6);
end 

 

二、复化辛普森：

function Sn = ComplexSimpson(xL, xR, n)
% xL,xR：积分区间左右端点，n：分段数
% Sn ：复化辛普森积分结果
h = (xR-xL)/n;%求步长
Sn = 0;
x = xL:h:xR;
for k = 2:n+1
    Sn = Sn + 4*fun((x(k)-0.5)) + 2*fun(x(k)); 
end 
Sn = Sn + fun(xL) - fun(xR);
Sn = vpa(Sn,7);
Sn = Sn*h/6;
Sn = vpa(Sn,6);
end 

三、复化科斯特：

function Cn = ComplexCotes(xL, xR,n)
% 判断积分区间个数是否是2的倍数，满足则进行计算，否则打印提示
if mod(n,4) ~= 0
        print('等分区间数错误！');return
else
    % 获取步长h
    h = (xR-xL)/n;
    %累积计算
    Cn = 0;
    for i = 1:n-3
        Cn = Cn + 7*fun(xL+h*(i-1))...
                 +32*fun(xL+h*i)...
                 +12*fun(xL+h*(i+1))...
                 +32*fun(xL+h*(i+2))...
                 +7*fun(xL+h*(i+3));
    end % 循环结束
    Cn = Cn*h/90;
    Cn = vpa(Cn,6);
end % if判断结束
end % 函数结束

四、龙贝格：

function I = romberg(fun,a,b,e)
% 使用龙贝格(Romberg数值求解公式)
% 输入例如：I=romberg(@(x)x^(3/2),0,1,0.000001)
% 判断输入参数是否足够
if nargin~=4
    error('请输入需要求积分的f、上界和下界以及误差e')
end
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a; %上下界间距
T(1,1)=h/2*(fun(a)+fun(b));
d=b-a; %误差间距
while e<=d
    k=k+1;
    h=h/2;
    sum=0;
    % 计算第一列T
    for i=1:n
        sum=sum+fun(a+(2*i-1)*h);
    end
    T(k+1,1)=T(k)/2+h*sum;
    % 迭代
    for j=1:k
        T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
    end
    n=n*2;
    d=abs(T(k+1,k+1)-T(k,k));
end
T 
I=T(k+1,k+1);

 

%复化辛普森公式
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.得到结果(包括系数、节点函数值和积分结果):







如果代码对你有帮助，请加一下关注，如果有问题可以私聊！最后，再次感谢喜欢我的小伙伴。