Adaptive Quadrature Algorithm using MATLAB (m file)
%
Adaptive quadrature algoritm
%
Find the integral of y=sin(x) from 0 to pi.
f = @(x)
sin(x);
aa = input('Enter
lower limit, a: ');
bb = input('Enter
upper limit, b: ');
n = input('Enter
maximum no. of levels, n: ');
eps = input('Enter
tolerance, tol: ');
cnt = 0;
app = 0;
i = 1;
tol = zeros(1,n);
a = zeros(1,n);
h = zeros(1,n);
fa = zeros(1,n);
fc = zeros(1,n);
fb = zeros(1,n);
s = zeros(1,n);
l = zeros(1,n);
fd = zeros(1,n);
fe = zeros(1,n);
v = zeros(1,7);
tol(i) = 10*eps;
a(i) = aa;
h(i) = 0.5*(bb-aa);
fa(i) = f(aa);
cnt = cnt+1;
fc(i) = f((aa+h(i)));
cnt = cnt+1;
fb(i) = f(bb);
cnt = cnt+1;
s(i) = h(i)*(fa(i)+4*fc(i)+fb(i))/3;
l(i) = 1;
while i > 0
fd = f((a(i)+0.5*h(i)));
cnt = cnt+1;
fe = f((a(i)+1.5*h(i)));
cnt = cnt+1;
s1 = h(i)*(fa(i)+4*fd+fc(i))/6;
s2 = h(i)*(fc(i)+4*fe+fb(i))/6;
v(1) = a(i);
v(2) = fa(i);
v(3) = fc(i);
v(4) = fb(i);
v(5) = h(i);
v(6) = tol(i);
v(7) = s(i);
lev = l(i);
i = i-1;
if abs(s1+s2-v(7)) < v(6)
app = app+(s1+s2);
else
if lev >= n
fprintf('Procedure
fails');
else
i = i+1;
a(i) = v(1)+v(5);
fa(i) = v(3);
fc(i) = fe;
fb(i) = v(4);
h(i) = 0.5*v(5);
tol(i) = 0.5*v(6);
s(i) = s2;
l(i) = lev+1;
i = i+1;
a(i) = v(1);
fa(i) = v(2);
fc(i) = fd;
fb(i) = v(3);
h(i) = h(i-1);
tol(i) = tol(i-1);
s(i) = s1;
l(i) = l(i-1);
end
end
end
fprintf('The
approximate integral is: %11.8f \n',
app);
OUTPUT:
>> Adaptive_quadrature_c
Enter lower limit, a: 0
Enter upper limit, b: pi
Enter maximum no. of levels, n: 6
Enter tolerance, tol: 0.001
The approximate integral is: 2.00026917
>>
.jpg)
Posts
a=0;
ReplyDeleteb=0.8;
n=2;
y= @(x) 0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5;
choice = menu ( 'choose a methode','trapizoidal','simpson1/3','simpson3/8','Gaussian Quadratute Algorithm');
switch choice
case 1
% trapziodal ;
h=(b-a)/n;
s=0.5*(y(a)+y(b));
for i=1 :n-1
s=s+y(a+i*h);
end
I=h*s
case 4
% Gaussian Quadratute Algorithm;
% Find the integral of f=0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5 from
% 0 to 0.8.
a = input('Enter lower limit, a: ');
b = input('Enter upper limit, b: ');
n = input('Enter the order, n: ');
if (n < 2 | n > 8)
fprintf('This code works only for order from 2 to 8\n');
else
a0 = 0.5*(b+a);
a1 = 0.5*(b-a);
if (n == 2)
c(1) = 1.0;
c(2) = c(1);
x_table(1) = -0.577350269;
x_table(2) = -x_table(1);
elseif (n == 3)
c(1) = 0.555555556;
c(2) = 0.888888889;
c(3) = c(1);
x_table(1) = -0.774596669;
x_table(2) = 0.0;
x_table(3) = -x_table(1);
elseif (n == 4)
c(1) = 0.347854845;
c(2) = 0.652145155;
c(3) = c(2);
c(4) = c(1);
x_table(1) = -0.861136312;
x_table(2) = -0.339981044;
x_table(3) = -x_table(2);
x_table(4) = -x_table(1);
elseif (n == 5)
c(1) = 0.236926885;
c(2) = 0.478628670;
c(3) = 0.568888889;
c(4) = c(2);
c(5) = c(1);
x_table(1) = -0.906179846;
x_table(2) = -0.538469310;
x_table(3) = 0.0;
x_table(4) = -x_table(2);
x_table(5) = -x_table(1);
elseif (n == 6)
c(1) = 0.171324492;
c(2) = 0.360761573;
c(3) = 0.467913935;
c(4) = c(3);
c(5) = c(2);
c(6) = c(1);
x_table(1) = -0.932469514;
x_table(2) = -0.661209386;
x_table(3) = -0.238619186;
x_table(4) = -x_table(3);
x_table(5) = -x_table(2);
x_table(6) = -x_table(1);
elseif (n == 7)
c(1) = 0.129484966;
c(2) = 0.279705391;
c(3) = 0.381830050;
c(4) = 0.417959184;
c(5) = c(3);
c(6) = c(2);
c(7) = c(1);
x_table(1) = -0.949107912;
x_table(2) = -0.741531186;
x_table(3) = -0.405845151;
x_table(4) = 0.0;
x_table(5) = -x_table(3);
x_table(6) = -x_table(2);
x_table(7) = -x_table(1);
elseif (n == 8)
c(1) = 0.101228536;
c(2) = 0.222381034;
c(3) = 0.313706645;
c(4) = 0.362683783;
c(5) = c(4);
c(6) = c(3);
c(7) = c(2);
c(8) = c(1);
x_table(1) = -0.960289856;
x_table(2) = -0.796666477;
x_table(3) = -0.525532409;
x_table(4) = -0.183434642;
x_table(5) = -x_table(4);
x_table(6) = -x_table(3);
x_table(7) = -x_table(2);
x_table(8) = -x_table(1);
end
integral = 0.0;
for i = 1:n
x(i) = a0 + a1*x_table(i);
f(i) = funkeval(x(i));
integral = integral + c(i)*f(i);
end
integral = integral*a1;
fprintf('The integral is: %10.8f\n\n',integral);
end
FUNCTION f(i) = funkeval(x(i));
f=0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5
end
>> UntitledRAQ
Enter lower limit, a: 0
Enter upper limit, b: 0.8
Enter maximum no. of levels, n: 4
Enter tolerance, tol: 0.001
The approximate integral is: 0.30329600
>> UntitledRAQ
Enter lower limit, a: 0
Enter upper limit, b: PI
Error using input
Undefined function or variable 'PI'.
Error in UntitledRAQ (line 8)
bb = input('Enter upper limit, b: ');
Enter upper limit, b: CLC
Error using input
Undefined function or variable 'CLC'.
Error in UntitledRAQ (line 8)
bb = input('Enter upper limit, b: ');
Enter upper limit, b: PI
Error using input
Undefined function or variable 'PI'.
Error in UntitledRAQ (line 8)
bb = input('Enter upper limit, b: ');
Enter upper limit, b: pi
Enter maximum no. of levels, n: 6
Enter tolerance, tol: 0.001
The approximate integral is: 2.00026917
>>