matlab实现复化梯形公式,复化simpson公式以及romberg积分

（一） 实验目的

熟悉并掌握数值积分的方法，重要训练复化梯形公式，复化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）。