数值积分与数值微分实验报告

实验三 数值积分程序设计算法
1）实验目的
通过本次实验熟悉并掌握各种数值积分算法及如何在matlab 中通过设计程序实现这些算法，从而更好地解决实际中的问题。

2）实验题目
给出积分 dx x I ⎰-=3
221
1 1.用Simpson 公式和N=8的复合Simpson 公式求积分的近似值.
2.用复合梯形公式、复合抛物线公式、龙贝格公式求定积分，要求绝对误差为 710*2
1
-=ε，将计算结果与精确解做比较，并对计算结果进行分析。

3）实验原理与理论基础
Simpson 公式
)]()2(4)([6b f b a f a f a b S +++-=
复化梯形公式
将定积分⎰=b
a
dx x f I )(的积分区间],[b a 分隔为n 等分，各节点为n j jh a x j ,,1,0, =+= n
a b h -=复合梯形(Trapz)公式为 ])()(2)([211
∑-=++-=n j j n b f x f a f n a b T 如果将],[b a 分隔为2n 等分，而n a b h /)(-=不变，则 )]()(2)(2)([41021112b f x f x f a f n a b T n j j n j j n +++-=∑∑-=+-= 其中
h j a h x x j j )21(2121
++=+=+，)]()(2)(2)([4102
1112b f x f x f a f n a b T n j j n j j n +++-=∑∑-=+-= ∑-=-++-+=10)2)12((221n j n n
a b j a f n a b T n=1时，a b h -=，则)]()([2
1b f a f a b T +-=
)0(0T = )2
1(22112h a f a b T T +-+=)1(0T = 若12-=k n ，记)1(0-=k T T n ， ,2,1=k 12--=k a b h jh a x j +=12--+=k a b j a h x x j j 2121
+=+k a b j a 2
)12(-++=，则可得如下递推公式
)0(0T )]()([2b f a f a b +-= ∑-=--++-+-=120
001)2)12((2)1(21)(k j k k a b j a f a b k T k T k=1,2, 即为梯形递推公式。

由复化递推公式的余项)(3122n n n T T T I -≈- n n T T I 3
1342-≈ n n j j n T x f n a b T I 31))(221(341021--+≈∑-=+∑-=+-+≈1
02
1)(6)(431n j j n x f n a b T ≈I ])(4)(2)()((6102
111∑∑-=+-=+++-=n j j n j j x f x f b f a f n a b n S = 即为复化Simpson 公式。

Romberg 公式
由复合Cotes 公式的余项)(63122n n n C C C I -≈
-得n n C C I 63
163642-≈ )1(63
1)(636422--=k T k T 令)1(3-k T )1(63
1)(636422--=k T k T 由此综合可得)0(0T )]()([2b f a f a b +-= ∑-=--++-+-=120
001
)2)12((2)1(21)(k j k k a b j a f a b k T k T )1(1-k T )1(3
1)(3400--=k T k T )1(2-k T )1(15
1)(151611--=k T k T )1(3-k T )1(63
1)(636422--=k T k T ,2,1=k 4）实验内容
本次实验需要通过在matlab 中编程实现复化梯形及Romberg 等各种数值积分算法，从而更加熟练的掌握这几种算法，也通过在matlab 中的实现来比较这几种算法之间的优劣性，从而在实际应用中更好地选择算法，以利于解决实际中的问题。

5）实验结果
Simpson 算法 function z=simpson(a,b)
c=(a+b)/2;
z1=1/(a^2-1);z2=1/(b^2-1);z3=1/(c^2-1);
z=(b-a)*(z1+4*z3+z2)/6;
a=2;b=3;simpson(a,b)
ans =
0.2034
复合Simpson算法
function y=comsimpson(a,b,n)
z1=1/(a^2-1);z2=1/(b^2-1);
h=(b-a)/n;s1=0;x1=a+h/2;
s2=0;x2=a+h;
for i=0:1:(n-1)
x1=x1+h;
s1=s1+1/(x1^2-1);
end
for i=1:1:(n-1)
x2=x2+h;
s2=s2+1/(x2^2-1);
end
y=h*(z1+4*s1+2*s2+z2)/6;
>> a=2;b=3;n=8;comsimpson(a,b,n)
ans =
0.1804
复合梯形公式
function y=comti(a,b,e)%复合梯形公式求解
z1=1/(a^2-1);z2=1/(b^2-1);
c=(a+b)/2;z3=1/(c^2-1);
t=abs(z3-z2-z1);z4=z3;n=2;
while (t>e)
n=n+1;
h=(b-a)/n;
x=a;
s1=z1+z2;
for i=1:n-1
x=a+i*h;
s1=s1+1/(x^2-1);
end
s1=s1*h;
t=abs(s1-z4);z4=s1;
end
y=s1;
复合抛物线公式求解
%复合抛物线公式求解
function y=compwx(a,b,e)
s=1/(a^2-1)+1/(b^2-1);c=(b-a)/2+a;z=1/(c^2-1);
s3=(s+4*z)*(b-a)/6;
r=abs(s3-s*(b-a)/6);
n=0;
while(r<e)
n=n+1;
h=(b-a)/(2*n);
x=a;
s1=s;
for i=1:2:(2*n-1)
x=x+h;
s1=s1+4*1/(x^2-1);
x=x+h;
s1=s1+2*1/(x^2-1);
end
s2=s1*h/3;
r=abs(s2-s3);
s3=s2;
end
y=s3;
Romberg求积算法
%龙贝格算法
function[s,n]=Romberg(a,b,eps)
if nargin<3,eps=1e-6;
end
z1=1/(a^2-1);z2=1/(b^2-1);
s=10;s0=0;k=2;t(1,1)=(b-a)*(z1+z2)/2;
while (abs(s-s0)>eps)
h=(b-a)/2^(k-1);
w=0;
if(h~=0)
for i=1:(2^(k-1)-1)
x=a+i*h;
w=w+1/(x^2-1);
end
t(k,1)=h*(z1/2+w+z2/2);
for l=2:k
for i=1:(k-l+1)
t(i,l)=(4^(l-1)*t(i+1,l-1)-t(i,l-1))/(4^(l-1)-1);
end
end
s=t(1,k);
s0=(t(1,k-1));
k=k+1;
n=k;
else s=s0;
n=-k;
end
end
6）实验结果分析与小结
通过本次实验，在matlab中编写程序实现复合Simpson,复合梯形，复合抛物线及Romberg算法等，进一步熟悉了这几种数值积分算法，通过比较发现simpson算法不稳定，而梯形法简单但收敛慢，由梯形法递推而得到的Romberg算法具有更高的收敛速度。