MATLAB数值分析之数值积分（一）

【实验原理】

一、复合辛普森方法

1）复合辛普森方法基本思路

把积分区间分成若干子区间，再在每个子区间上用低阶辛普森求积公式。

2）复合辛普森算法描述

二、龙贝格方法

1）基本思路

实质上是梯形公式的递推化，是在梯形公式、辛普森公式和柯特斯公式关系的基础上构造出的加速计算的方法，是一种外推算法。

2）算法描述

【实验内容】

一、回答下面的问题

1.什么是复合求积公式?写出事前误差估计和事后误差估计的公式。

2.当积分区间是奇数等份时，如何搭配使用n=2和n=3的辛普森公式，举例说明。

3.试说明高斯-勒让德积分、牛顿-柯特斯的区别和联系。

4.时间测试下，请对自己编写的龙贝格算法的执行效率进行分析，说明对分次数和计算执行时间的关系？

二、编程计算课后习题中规定题目，交回实验报告与计算

1.用不同数值方法计算积分

（1）取不同的步长h，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？

实验代码

%程序代码：复合梯形
function [t,dy]=tixing2(a,b,n)
format long
h=(b-a)/n;%等分为n份，计算步长
x=a:h:b;
x(1)=1;
y=sqrt(x).*log(x);%定义函数y
x(1)=0;%由于lnx会导致计算的瑕疵，故人工为瑕疵点赋值。计算得y1=0
t=0;
for k=1:(b-a)/h
     t=t+y(k)+y(k+1);%复合梯形公式
end
t=t*h/2;
dy=t+4/9; %误差
end

实验结果

运行结果

实验代码

%程序代码：复合辛普森
function [t,dy]=sampusen(a,b,n)
h=(b-a)/n; %等分为n份，计算步长
x=a:h:b;
x(1)=1;
y=sqrt(x).*log(x); %定义函数y
z=sqrt(x+h/2).*log(x+h/2);
x(1)=0; %由于lnx会导致计算的瑕疵，故人工为瑕疵点赋值。计算得y1=0
t=0;
for k=1:(b-a)/h
     t=t+y(k)+y(k+1)+4*z(k);%复合辛普森公式
end
t=t*h/6;
dy=t+4/9; %误差
end

实验结果

运行结果

（2）用龙贝格求积计算问题（1）

实验代码

clear;clc;
a=0;b=1;
e=1e-4;
n=10;  %算到R2才能用事后误差估计
R=zeros(1,4);%矩阵维数可以动态增加，初始设为1行4列，R的第1~4列分别存储T,S,C,R
R(1,1)=0;%人工为瑕疵点赋值
%求复合梯形公式，最高对分到2^n次，够算|R2n-Rn|<e，用事后误差估计
    for k=1:n-1%梯形值的递推过程
        h=(b-a)/2^k;
        sum=0;
        for i=1:2^(k-1)
            sum=sum+f(a+(2*i-1)*h);
        end
        R(k+1,1)=R(k,1)/2+h*sum;
    end
    for j=1:3%加速值
        fac=1/(4^j-1);
        for k=1:n-j
            R(k+j,j+1)=(1+fac)*R(k+j,j)-fac*R(k+j-1,j);
        end
    end
    format long
    format compact
    %--------------------判断精度---------------------
    e1=abs(R(n,4)-R(n-1,4)),%输出截断误差
    if e1<e
        I=R(n,4), %输出积分近似值
        break;
    end
    %--------------------判断精度---------------------

实验结果

运行结果

（3）用自适应辛普森积分，使得精度达到10-4

实验代码

function y=zishiyingsam(a,b,e)

h=(b-a)/2;

y0=sampusen(a,b,50); %自适应辛普森仍与n有较大关系

y1=sampusen(a,a+h,50); %(a+b)/2

y2=sampusen(a+h,b,50);

if(abs(y2+y1-y0)>e)

 y1=zishiyingsam(a,a+h,e/2);%递归过程

 y2=zishiyingsam(a+h,b,e/2);

 y=y1+y2;

else

 y=y1+y2+(1/15)*(y1+y2-y0);

end

实验结果

运行结果

---

相关文章：
数值分析之带初值的常微分方程数值解法（二）