数值分析实验

实验四 数值积分与数值微分

专业班级：信计131班 姓名：段雨博 学号：2013014907 一、实验目的

1、熟悉matlab 编程。

2、学习数值积分程序设计算法。

3、通过上机进一步领悟用复合梯形、复合辛普森公式，以及用龙贝格求积方法计算积分的原理。

二、实验题目 P137

1

、用不同数值方法计算积分

4

9

xdx =-⎰

。

（1）取不同的步长h .分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？ （2）用龙贝格求积计算完成问题（1）。 三、实验原理与理论基础

1.1复合梯形公式及其复合辛普森求解

[]()()()11

101()()222n n n k k k k k h h T f x f x f a f x f b --+==⎡⎤=+=++⎢⎥⎣⎦

∑∑

误差关于h 的函数：()()2

12

n b a R f

h f η-''=-

复合辛普森公式：()()()()11

1/201426n n n k k k k h S f a f x f x f b --+==⎡⎤

=+++⎢⎥⎣⎦

∑∑

误差关于h 的函数：(

)()4

41802n n b a h R f I S f η-⎛⎫=-=- ⎪⎝⎭

1.2龙贝格求积算法：

龙贝格求积公式是梯形法的递推化，也称为逐次分半加速法，它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种计算积分的方法，同时它有在不断增加计算量的前提下提高误差的精度的特点。 计算过程如下：

（1）取0,k h b a ==-，求：

()()()[]()00.,.2

h

T f a f b k a b =

+→⎡⎤⎣⎦令k 1记为区间的二分次数 （2）求梯形值02k b a T -⎛⎫ ⎪⎝⎭

即按递推公式12102122n n n k k h T T f x -+=⎛⎫=+ ⎪⎝⎭∑计算0k

T .

（3）求加速值，按公式()

()()

111444141

m m k k k m

m m m m T T T +--=---逐个求出T 表的地k 行其余各

元素()()1,2,,k j j

T j k -=L

（4）若()

()001k

k T T ε--<（预先给定的精度）

，则终止计算，并取()

()0;1k T I k k ≈+→否则令转（2）继续计算。

上表指出了计算过程，第二列2

k h =

给出了子区间长度，i 表示第i 步外推。 可以证明，如果()f x 充分光滑，那么T 表每一列的元素及对角线元素均收敛到所求的积分

I 。

四、实验内容 程序设计如下：

1、复合梯形法M 文件：

function t=natrapz(fname,a,b,n) h=(b-a)/n;

fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0.001*h); t=h*(0.5*(fa+fb)+sum(f));

2、复合辛普森法M 文件：

function t=natrapz(fname,a,b,n) h=(b-a)/n;

fa=feval(fname,a);fb=feval(fname,b);f1=feval(fname,a+h:h:b-h+0.001*h);

f2=feval(fname,a+h/2:h:b-h+0.001*h); t=h/6*(fa+fb+2*sum(f1)+4*sum(f2));

3、龙贝格算法M 文件：

function [I,step]=Roberg(f,a,b,eps) if (nargin==3) eps=1.0e-4; end ; M=1;

tol=10;

k=0;

T=zeros(1,1);

h=b-a;

T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym( f)),b));

while tol>eps

k=k+1;

h=h/2;

Q=0;

for i=1:M

x=a+h*(2*i-1);

Q=Q+subs(sym(f),findsym(sym(f)),x);

end

T(k+1,1)=T(k,1)/2+h*Q;

M=2*M;

for j=1:k

T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);

end

tol=abs(T(k+1,j+1)-T(k,j));

end

I=T(k+1,k+1);

step=k;

1、复合梯形法运行：

>> format long;natrapz(inline('sqrt(x).*log(x)'),eps,1,10),format short;

ans =

-0.417062831779470

>> format long;natrapz(inline('sqrt(x).*log(x)'),eps,1,100),format short;

ans =

-0.443117908008157

>> format long;natrapz(inline('sqrt(x).*log(x)'),eps,1,1000),format short;

ans =

-0.444387538997162

2、复合辛普森法运行：

>> format long;natrapzz(inline('sqrt(x).*log(x)'),eps,1,10),format short;

ans =

-0.435297890074689

>> format long;natrapzz(inline('sqrt(x).*log(x)'),eps,1,100),format short;

ans =

-0.444161178415673

>> format long;natrapzz(inline('sqrt(x).*log(x)'),eps,1,1000),format short;

ans =

-0.444434117614180

3、龙贝格算法运行：