web123456

Composite trapezoidal formula, composite simpson's formula matlab

 

 

1. Calculation of integrals using Newton-Cotes formulas of orders 1 to 4

                         

Program:

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. is known and can therefore be approximated by numerical integration.

(1) Approximations calculated using the composite trapezoidal formula and the composite Simpson formula by taking and, respectively;

Program:

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

%Calculating integrals using the compound trapezoidal formula

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

Compound Simpson:

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

%Calculating integrals using the composite Simpson's formula

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) Divide the interval [0,1] into equal parts and use the approximation calculated by the composite trapezoidal formula and the composite Simpson's formula to ask how many equal parts of the interval [0,1] need to be divided if the error is required not to be exceeded;

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

Composite Simpson's formula:

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) Choosing different, for both composite product formulas, try to characterize the error as of thefunction (math.), and compare the accuracy of the two methods.

Composite product formula:

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-')

 

Composite Simpson's formula:

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);

 

Composite Simpson's formula is more accurate as can be seen from the image comparison.

(4) Does there exist a certain value that, when less than this value and then continuing to decrease, the calculation no longer improves? Why?

be

Composite product formula:

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-')

 

 

Composite Simpson's formula:

 

It can be noticed through the images that the accuracy no longer changes significantly when h<0.025 is reached.

3. Calculation of integrals using the three-point and five-point Gauss-Legendre formulas, respectively

                         

Program:

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);

 

Three points:

syms x

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

a=0;b=1;

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

a =

0.359187

Five points:

a =

0.359141