（一） 实验目的
熟悉并掌握数值积分的方法，重要训练复化梯形公式，复化simpson 公式以及romberg 积分。
（二） 问题描述
问题三数值积分椭圆周长的计算。
考虑椭圆22221x y a b
+=,为计算其周长,只要计算其第一象限的长度即可.
用参数方程可以表示为cos (0/2)sin x a t t y b t π=⎧≤≤⎨=⎩
,
计算公式为/0π⎰
为计算方便,我们可以令1a =,即计算下面的积分
/
0π⎰/0π=⎰
(/0π⎰/0a π=⎰可以归结为上面的形式)
采用复化梯形公式，复化Simpson 公式以及Romberg 积分的方法计算积分
/
0()I b π=⎰
给出通用程序，该通用程序可以计算任何一个函数在任意一个区间在给定的精度下的数值积分。
程序输出为计算出的数值积分值以及计算函数值的次数。
（三） 算法介绍
首先利用给出的各迭代公式，设计程序。
在matlab 对话框中输入要计算的函数，给出区间和精度。
复化梯形的迭代公式为：
∫f(x)dx b a =h/2[f (a )+2∑f(x j )n−1j=1+f(b)];
复化simpson 迭代公式为：
∫f(x)dx b a =h/3[f (a )+2∑f(x 2j )(n 2)−1j=1+4∑f(x 2j−1)(n 2
)j=1+f(b)]; Romberg 迭代公式为：
R k,j =R k,j−1+R k,j−1—R k−1,j−1
4j−1−1。
（四） 程序
对于复化梯形公式和复化simpson 公式，我们放在jifenn.m 中。
（%标记后的程序可用来把b 看为变量时的算法实现）
%复化梯形公式
function y=jifenn(f,n,a,b) （说明：f 表示任一函数，n 精度，a ，b 为区间） fi=f(a)+f(b);
h=(b-a)/n;
d=1;
%function f=jifen(n,a,b,c)
%syms t
%y=sqrt(1+(c^2-1)*cos(t)^2);
%ya=subs(y,t,a);
%yb=subs(y,t,b);
%fi=ya+yb;
for i=1:n-1
x=a+i*h;
fi=fi+2*f(x);
d=d+1;
%yx=subs(y,t,x);
%fi=fi+2*yx;
end
f4=h/2*fi,d
%复化simposon 公式
f1=0;
f2=0;
dd=1;
for i=1:n-1
dd=dd+1;
if rem(i,2)~=0;
x1=a+i*h;
f1=f1+f(x1);
else rem(i,2)==0;
x2=a+i*h;
f2=f2+f(x2) ;
end
end
f3=(h/3)*(f(a)+4*f1+2*f2+f(b)),dd
对于romberg积分，建立romberg.m文件。
function y=romberg(f,n,a,b) （说明：f表示任一函数，n精度，a，b为区间）z=zeros(n,n);
h=b-a;
z(1,1)=(h/2)*(f(a)+f(b));
f1=0;
for i=2:n
for k=1:2^(i-2)
f1=f1+f(a+(k-0.5)*h);
end
z(i,1)=0.5*z(i-1,1)+0.5*h*f1;
h=h/2;
f1=0;
for j=2:i
z(i,j)=z(i,j-1)+(z(i,j-1)-z(i-1,j-1))/(4^(j-1)-1);
end
end
z,n
（五）运行结果
对于复化梯形公式和复化simpson公式，我们运行下列语句并得到结果：
>> fun=inline('sqrt(1+(0.5^2-1).*cos(t).^2)');
>> jifenn(fun,8,0,pi/2)
f4 =
1.2111
d =
8
f3 =
1.2111
dd =
8
>> 1.2111*4
ans =
4.8444
>> 1.2111*4
ans =
4.8444
（说明：在本题中将椭圆中的未知量a取为1，b取为0.5。
f4为复化梯形公式得到的椭圆周长，f3为复化simpson公式得到的椭圆周长）。
对于romberg，运行下列语句并最终得到结果为：
>> fun=inline('sqrt(1+(0.5^2-1).*cos(t).^2)');
>> romberg(fun,8,0,pi/0.5)
z =
3.1416 0 0 0 0 0 0 0
3.1416 3.1416 0 0 0 0 0 0
4.7124
5.2360 5.3756 0 0 0 0 0 4.8398 4.8823 4.8587 4.8505 0 0 0 0 4.8442 4.8457 4.8432 4.8430 4.8429 0 0 0 4.8442 4.8442 4.8441 4.8441 4.8442 4.8442 0 0 4.8442 4.8442 4.8442 4.8442 4.8442 4.8442 4.8442 0 4.8442 4.8442 4.8442 4.8442 4.8442 4.8442 4.8442 4.8442
ans =
4.8442
n =
8
（说明：其中最终结果为4.8442）。
（六）结果分析
我们计算了当椭圆长轴为1，短轴为0.5时的周长。
通过上述三种方法的计算可以看到，结果相差不大。
根据椭圆周长的一个计算公式L=2πa(长轴)+4（b−a）我们可以得到L=4.283。
因此三种方法都较好的接近真值。
（七）心得体会
应该熟练掌握这三种方法，才能在编程时正确快速的写出迭代公式。
同时在一种思想的前提下可以寻找多种方法实现算法，互相验证。