Matlab数值积分（2）

实验目的：

　　掌握理查森外推法

实验要求：

1. 给出理查森外推算法

2. 用Matlab实现理查森外推算法

3. 用Matlab实现自适应积分算法

实验内容：

　　1. 理查森外推算法，数学知识：利用Richardson外推对逐次分半，若记则有由Richardson外推方法，可得到左式的误差为考虑舍入误差，m不能取得太大。

　　算法描述：

　　（1）命名函数；（2）如果输入的未知数少于4个，默认精度0.001；（3）描述T表示矩阵坐标；（4）依次赋值计算T表第一列；（5）根据Richardson公式求出T表矩阵的值；（6）若达到精度则运算结束，若未达到则循环计算；（7）输出T表，得出的值就是导数值。

　　

　　2.用Matlab实现理查森外推算法（见实验步骤）。

　　3. 用Matlab实现自适应积分算法，被积函数在整个积分区间上的变化是不均衡到，在[a,b]分成的若干子区间中，有些变化缓慢，有些则变化大。为了使计算结果达到预期精度，可以使用Simpson求积公式。（此处只展示MATLAB程序）。

实验步骤：

　　1.理查森外推算法代码：

function t=romberg(fname,a,b,ep)
if nargin<4
    ep=1e-5;
end
i=1;j=1;h=b-a;
T(i,1)=h/2*(feval(fname,a)+feval(fname,b));
T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h/2+0.001*h))*h/2;
T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);
while abs(T(i+1,i+1)-T(i,i))>ep
    i=i+1;h=h/2;
    T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+0.001*h))*h/2;
    for j=1:i
        T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);
    end
end
T
t=T(i+1,j+1);

　　运行示例：

　　

　　2.自适应积分算法代码：

function I=squad1(fun,a,b,ep)
if nargin<4
    ep=1e-5;
end
N=1;
h=b-a;
T1=h/2*(feval(fun,a)+feval(fun,b));
S0=T1;
while 1
    h=h/2;
    T2=T1/2;
    for k=1:N
        T2=T2+h*feval(fun,a+(2*k-1)*h);
    end
    I=(4*T2-T1)/3;
    if abs(I-S0)<ep
        break;
    end
    N=2*N;
    T1=T2;
    S0=I;
end

　　运行：

　　

小结：

　　在编写数学类的程序时，必须要熟读相关的数学知识。