1. 用1阶至4阶Newton-Cotes公式计算积分

程序：

function I = NewtonCotes(f,a,b,type)

%

syms t;

t=findsym(sym(f));

I=0;

switch type

case 1,

I=((b-a)/2)*(subs(sym(f),t,a)+subs(sym(f),t,b));

case 2,

I=((b-a)/6)*(subs(sym(f),t,a)+4*subs(sym(f),t,(a+b)/2)+...

subs(sym(f),t,b));

case 3,

I=((b-a)/8)*(subs(sym(f),t,a)+3*subs(sym(f),t,(2*a+b)/3)+...

3*subs(sym(f),t,(a+2*b)/3)+subs(sym(f),t,b));

case 4,

I=((b-a)/90)*(7*subs(sym(f),t,a)+...

32*subs(sym(f),t,(3*a+b)/4)+...

12*subs(sym(f),t,(a+b)/2)+...

32*subs(sym(f),t,(a+3*b)/4)+7*subs(sym(f),t,b));

case 5,

I=((b-a)/288)*(19*subs(sym(f),t,a)+...

75*subs(sym(f),t,(4*a+b)/5)+...

50*subs(sym(f),t,(3*a+2*b)/5)+...

50*subs(sym(f),t,(2*a+3*b)/5)+...

75*subs(sym(f),t,(a+4*b)/5)+19*subs(sym(f),t,b));

case 6,

I=((b-a)/840)*(41*subs(sym(f),t,a)+...

216*subs(sym(f),t,(5*a+b)/6)+...

27*subs(sym(f),t,(2*a+b)/3)+...

272*subs(sym(f),t,(a+b)/2)+...

27*subs(sym(f),t,(a+2*b)/3)+...

216*subs(sym(f),t,(a+5*b)/6)+...

41*subs(sym(f),t,b));

case 7,

I=((b-a)/17280)*(751*subs(sym(f),t,a)+...

3577*subs(sym(f),t,(6*a+b)/7)+...

1323*subs(sym(f),t,(5*a+2*b)/7)+...

2989*subs(sym(f),t,(4*a+3*b)/7)+...

2989*subs(sym(f),t,(3*a+4*b)/7)+...

1323*subs(sym(f),t,(2*a+5*b)/7)+...

3577*subs(sym(f),t,(a+6*b)/7)+751*subs(sym(f),t,b));

end

syms x

f=exp(-x).*sin(x);

a=0;b=2*pi;

I = NewtonCotes(f,a,b,1)

N=1:

I =

0

N=2:

I =

0

N=3:

I =

(pi*((3*3^(1/2)*exp(-(2*pi)/3))/2 - (3*3^(1/2)*exp(-(4*pi)/3))/2))/4

N=4:

I =

(pi*(32*exp(-pi/2) - 32*exp(-(3*pi)/2)))/45

2. 已知，因此可以通过数值积分计算的近似值。

(1)分别取和，利用复合梯形公式和复合Simpson公式计算的近似值；

程序：

function Y= CombineTraprl(f,a,b,h)

%用复合梯形公式计算积分

syms t;

t= findsym(sym(f));

n=(b-a)/h;

I1= subs(sym(f),t,a);

l=0;

for k=1:n-1

xk=a+h*k;

l=l+2*subs(sym(f),t,xk);

end

Y=(h/2)*(I1+l+subs(sym(f),t,b));

syms x

f=4/(1+x^2);

a=0;b=1;

y= CombineTraprl(f,a,b,0.1);

vpa(y,6)

h=0.1:

ans =

3.13993

H=0.2:

ans =

1.04498

复合辛普森：

function Y= CombineSimpson(f,a,b,h)

%用复合辛普森公式计算积分

syms t;

t= findsym(sym(f));

n=(b-a)/h;

I1= subs(sym(f),t,a);

l=0;

for k=1:n-1

xk=a+h*k;

l=l+2*subs(sym(f),t,xk);

end

l2=0;

for k=1:n-1

xk2=a+h*(k+1)/2;

l2=l2+4*subs(sym(f),t,xk2);

end

Y=(h/6)*(I1+l+l2+subs(sym(f),t,b));

H=0.1:

ans =

3.22605

H=0.2:

ans =

2.93353

(2)把区间[0,1] 等分，利用复合梯形公式和复合Simpson公式计算的近似值，若要求误差不超过，问需要把区间[0,1]划分成多少等份；

function n=trap(f,a,b)

syms t;

t= findsym(sym(f));

I=zeros(1,500);

I(1)=((b-a)/2)*(subs(sym(f),t,a)+subs(sym(f),t,b));

I(2)=((b-a)/4)*(subs(sym(f),t,a)+2*subs(sym(f),t,(b-a)/2)+subs(sym(f),t,b));

k=3;

while((I(k-1)-I(k-2))>1/2*10^(-6))

l=0;

for i=1:k-1

xi=a+(b-a)/k*i;

l=l+2*subs(sym(f),t,xi);

end

I(k)=((b-a)/(2*k))*(subs(sym(f),t,a)+l+subs(sym(f),t,b));

k=k+1;

end

n=k-1;

syms x;

f=4./(1+x.^2);

a=0;b=1;

n=trap(f,a,b)

n =

88

复合辛普森公式：

function n=Simpson(f,a,b)

syms t;

t= findsym(sym(f));

I=zeros(1,500);

I(1)=((b-a)/6)*(subs(sym(f),t,a)+4*subs(sym(f),t,(b-a)/2)+subs(sym(f),t,b));

I(2)=((b-a)/12)*(subs(sym(f),t,a)+4*subs(sym(f),t,(b-a)/4)+4*subs(sym(f),t,3*(b-a)/4)+2*subs(sym(f),t,(b-a)/2)+subs(sym(f),t,b));

k=3;

while((I(k-1)-I(k-2))>1/2*10^(-6))

l=0;

m=4*subs(sum(f),t,(a+((a+b)/(2*k))));

for i=1:k-1

xi=a+(b-a)/k*i;

l=l+2*subs(sym(f),t,xi);

end

for j=1:k-1

xj=a+(b-a)/(k*2)+(b-a)/k*j;

m=m+4*subs(sym(f),t,xj);

end

I(k)=((b-a)/(2*k))*(subs(sym(f),t,a)+l+m+subs(sym(f),t,b));

k=k+1;

end

n=k-1;

n =

5

(3)选择不同的，对两种复合求积公式，试将误差描述为的函数，并比较两种方法的精度。

复合求积公式：

function y=traprls(f,a,b,h)

syms t;

t= findsym(sym(f));

n=(b-a)/h;

l=0;

for k=1:n-1

xk=a+h*k;

l=l+2*subs(sym(f),t,xk);

end

I1=(h/2)*(subs(sym(f),t,a)+l+subs(sym(f),t,b));

h=(b-a)/(n-1);

n=(b-a)/h;

l=0;

for k=1:n-1

xk=a+h*k;

l=l+2*subs(sym(f),t,xk);

end

I2=(h/2)*(subs(sym(f),t,a)+l+subs(sym(f),t,b));

y=I2-I1;

y=abs(y);

y=vpa(y,8);

syms x;

f=4./(1+x.^2);

a=0;b=1;

h=0.01:0.05:0.5;

v=zeros(1,10);

for i=1:10

v(i)=traprls(f,a,b,h(i))

end

v

plot(h,v,‘r-‘)

复合辛普森公式：

function y=Simpsons(f,a,b,h)

syms t;

t= findsym(sym(f));

n=(b-a)/h;

l=0;

m=4*subs(sum(f),t,(a+h/2));

for k=1:n-1

xk=a++h*k;

l=l+2*subs(sym(f),t,xk);

end

for i=1:n-1

xi=a+h/2+h*i;

m=m+4*subs(sym(f),t,xi);

end

I1=(h/6)*(subs(sym(f),t,a)+l+m+subs(sym(f),t,b));

h=(b-a)/(n-1);

n=(b-a)/h;

l=0;

m=4*subs(sum(f),t,(a+h/2));

for k=1:n-1

xk=a++h*k;

l=l+2*subs(sym(f),t,xk);

end

for i=1:n-1

xi=a+h/2+h*i;

m=m+4*subs(sym(f),t,xi);

end

I2=(h/6)*(subs(sym(f),t,a)+l+m+subs(sym(f),t,b));

y=abs(I2-I1);

y=vpa(y,10);

通过图像对比可知，复合辛普森公式精度更高。

(4)是否存在某个值，当小于这个值之后，再继续减小，计算结果不再有改进？为什么？

是

复合求积公式：

syms x;

f=4./(1+x.^2);

a=0;b=1;

h=0.001:0.004:0.2;

v=zeros(1,10);

for i=1:50

v(i)=traprls(f,a,b,h(i));

end

plot(h,v,‘r-‘)

复合辛普森公式：

通过图像可以发现，当h<0.025后，精度不再有显著改变。

3. 分别用三点和五点Gauss-Legendre公式计算积分

程序：

function I = IntGaussLegen(f,a,b,n)

syms t;

t= findsym(sym(f));

ta = (b-a)/2;

tb = (a+b)/2;

switch n

case 0,

I=2*ta*subs(sym(f),t,tb);

case 1,

I=ta*(subs(sym(f),t,ta*0.5773503+tb)+...

subs(sym(f),t,-ta*0.5773503+tb));

case 2,

I=ta*(0.55555556*subs(sym(f),t,ta*0.7745967+tb)+...

0.55555556*subs(sym(f),t,-ta*0.7745967+tb)+...

0.88888889*subs(sym(f),t,tb));

case 3,

I=ta*(0.3478548*subs(sym(f),t,ta*0.8611363+tb)+...

0.3478548*subs(sym(f),t,-ta*0.8611363+tb)+...

0.6521452*subs(sym(f),t,ta*0.3398810+tb) +...

0.6521452*subs(sym(f),t,-ta*0.3398810+tb));

case 4,

I=ta*(0.2369269*subs(sym(f),t,ta*0.9061793+tb)+...

0.2369269*subs(sym(f),t,-ta*0.9061793+tb)+...

0.4786287*subs(sym(f),t,ta*0.5384693+tb) +...

0.4786287*subs(sym(f),t,-ta*0.5384693+tb)+...

0.5688889*subs(sym(f),t,tb));

case 5,

I=ta*(0.1713245*subs(sym(f),t,ta*0.9324695+tb)+...

0.1713245*subs(sym(f),t,-ta*0.9324695+tb)+...

0.3607616*subs(sym(f),t,ta*0.6612094+tb)+...

0.3607616*subs(sym(f),t,-ta*0.6612094+tb)+...

0.4679139*subs(sym(f),t,ta*0.2386292+tb)+...

0.4679139*subs(sym(f),t,-ta*0.2386292+tb));

end

I=simplify(I);

I=vpa(I,6);

三点：

syms x

f=x.*exp(x)./((1+x)^2);

a=0;b=1;

a=IntGaussLegen(f,a,b,2)

a =

0.359187

五点：

a =

0.359141