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