《计算方法》实验报告
姓名：
班级：
学号：
实验日期: 2011年10月26日
一、实验题目：
数值积分
二、实验目的：
1．熟悉matlab 编写及运行数值计算程序的方法。
2．进一步理解数值积分的基础理论。
3．进一步掌握应用不同的数值积分方法求解给定的积分并给出数据结果及误差分析。
三、实验内容：
1．分别用复合梯形求积公式及复合辛普森求积公式计算积分xdx x ln 10
⎰
,
要求计算精度达到410-,给出计算结果并比较两种方法的计算节点数. 2．用龙贝格求积方法计算积分dx x x ⎰+3
021,使误差不超过510-.
3．用3=n 的高斯-勒让德公式计算积分⎰3
1
sin x e x ,给出计算结果.
4．用辛普森公式(取2==M N ) 计算二重积分.5
.00
5
.00
dydx e x y ⎰
⎰
-
四、实验结果：
1.（1）复合梯形法:
将区间[a,b]划分为n 等份，分点n k n
a
b h kh a x k ,2,1,0,,=-=+=在每个区间[1,+k k x x ](k=0,1,2,···n-1)上采用梯形公式，则得
)()]()([2)()(1
11
1
f R x f x f h dx x f dx x f I n n k k k b
a
n k x x k k
++===∑⎰∑⎰
-=+-=+
故)]()(2)([21
1
b f x f a f h
T n k k n ++=∑-=称为复合梯形公式
计算步长和划分的区间
Eps=1E-4
h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1))) h1 =0.0600 N1=ceil(1/h1) N1 =17
用复合梯形需要计算17个结点。
复合梯形：
function T=trap(f,a,b,n) h=(b-a)/n;
T=0;
for k=1:(n-1); x0=a+h*k;
T=T+limit(f,x0); end;
T=h*(limit(f,a)+limit(f,b))/2+h*T; T=double(T);
（2）复合辛普森：
将区间[a,b]划分为n 等份，在每个区间[1,+k k x x ](k=0,1,2,···n-1)上采用辛普森公式，若记h x x k k 2
1
2/1+
=+，则得 )()]()(4)([6)()(1
12/11
1
f R x f x f x f h dx x f dx x f I n n k k k k b
a
n k x x k k
+++===∑⎰∑⎰
-=++-=+
记)]()(2)(4)([61
1
102/1b f x f x f a f h
S n k k n k k n +++=∑∑-=-=+为复合辛普森公式
h2=power(Eps/abs(-(1-0)/180*1/(1+4)),1/4) h2 =0.5477 N2=(1/h2) N2=ceil(1/h2) N2 =2
用复合辛普森需要计算2个结点。
复合辛普森：
function S=simpson(f,a,b,n) h=(b-a)/(2*n); s1=0; s2=0; for k=1:n
x0=a+h*(2*k-1); s1=s1+limit(f,x0); end
for k=1:(n-1)
x0=a+h*2*k; s2=s2+limit(f,x0); end
S=h*(limit(f,a)+limit(f,b)+4*s1+2*s2)/3; S=double(S);
2.龙贝格算法思想：
（1）将区间[a,b]划分为n 等分，分点为：（n x x x 10,）；根据梯形公式
)]()(2)([21
1
b f x f a f h
T n k k n ++=∑-=，求出n T ，再根据n T 和n T 2之间的递推公式
n T 2∑-=++=
1
2/1)(22n k k n x f h T 求出n T 2 （2）设m 为加速次数，k 为划分区间次数，则由加速公式：
)2,1(1
41144)(1)1(1)( =---=-+-k T T T
k m m
k m m m k m
求出第k 次划分，第m 次加速次数的梯形值)(k m T ，这样不断地循环，直到求出在满足精度条件下的某个)
(k m T 作为积分值
为止。
function [I,T]=romberg(f,a,b,n,Eps) if nargin<5 Eps=1E-5; end m=1; h=b-a; err=1; j=0; x=a;
T=zeros(4,4);
T(1,1)=h*(limit(f,a)+limit(f,b))/2; while((err>Eps)&&(j<n)||(j<4)) j=j+1; h=h/2; s=0;
for p=1:m;
x=a+h*(2*p-1); s=s+limit(f,x); end
T(j+1,1)=T(j,1)/2+h*s; m=2*m; for k=1:j
T(j+1,k+1)=T(j+1,k)+(T(j+1,k)-T(j,k))/(4^k-1); end
err=abs(T(j,j)-T(j+1,k+1)); end
I=T(j+1,j+1);
输入的命令：
syms x;
f=sym(x*(1+x^2)^(1/2)); [I,T]=romberg(f,0,3,4,1E-5)
以上为3.高斯-勒让德算法：
先做变换，将积分区间变换到[-1,1]，令2
2a
b t a b x ++-=则有 ⎰⎰
-++=1
1
)2(3
1
)2sin()sin(*dt t e dx x e t x
，再利用高斯-勒让德求积公式求解。
f=inline('exp(t+2).*sin(t+2)')
x=[-0.8611363,-0.33399810,0.33399810,0.8611363]; A=[0.3478548,0.6521452,0.6521452,0.3478548]; I=0;
for i=1:length(x)
I=I+feval(f,x(i)).*A(i); end
I=10.9676
即由高斯勒让德公式得到的结果为I=10.9676 4.辛普森公式计算二重积分：
将二次积分转化成两个一次积分做乘积的形式，分别对两个一次积分用辛普森公式，求积最终得到二次积分的结果。
function s=simprl(f,a,b,n) %f 是被积函数
%a,b 分别为积分的上下限 %n 是子区间的个数 %s 是梯形总面积 h=(b-a)/(2*n); s1=0; s2=0; for k=1:n
x=a+h*(2*k-1); s1=s1+feval('f',x); end
for k=1:(n-1) x=a+h*2*k;
s2=s2+feval('f',x); end
s=h*(feval('f',a)+feval('f',b)+4*s1+2*s2)/3; 第一两个 .m 文件 function y=f(x) y=exp(-x);
function y=g(x)
y=exp(-x);
输入simprl('exp(-x)',0,0.5,1)
ans =
0.393477815999854
输入simprl('exp(x)',0,0.5,1)
ans =
0.648735244787591
.5.00
5
.00
dydx e x y ⎰⎰
-=0.393477815999854*0.648735244787591=
0.255262927281152
五、实验体会或遇到问题：
在本次试验当中，遇到了很多问题，如在区间端点没有定义，这都需要考虑，将其中的原理整理清楚是最重要的。
有了这次课内实验，帮助自己把学的知识巩固了一遍。