复合梯形公式、复合辛普森公式 matlab

 

 

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

posted @ 2018-11-10 22:52  夜游星  阅读(16743)  评论(0编辑  收藏  举报