Permalink
Browse files

initial commit for lab5

  • Loading branch information...
1 parent 56f4ffd commit caaded6d2df39f8ec4e17440d25f5860c947ae5c @uberj committed Nov 18, 2012
View
@@ -0,0 +1,4 @@
+all:
+ MATLABPATH=$MATLABPATH:./
+ export MATLABPATH
+ matlab -nosplash -nodesktop
View
Binary file not shown.
View
@@ -0,0 +1,25 @@
+function [val,bp,wf]=gaussint(a,b,n,fun)
+
+% [val,bp,wf]=gaussint(a,b,n,fun) integrates
+% a function from a to b using an n-point
+% Gauss rule which is exact for a polynomial
+% of degree 2*n-1. Concepts on page 93 of
+% 'Methods of Numerical Integration' by
+% Philip Davis and Philip Rabinowitz yield
+% the base points and weight factors.
+%
+% fun - function being integrated
+% a,b - integration limits
+% n - order of formula
+% val - value of the integral
+% bp,wf - Gauss base points and weight factors on [-1,1]
+
+fun=vectorize(fun);
+
+u=(1:n-1)./sqrt((2*(1:n-1)).^2-1);
+[vc,bp]=eig(diag(u,-1)+diag(u,1));
+[bp,k]=sort(diag(bp));
+wf=2*vc(1,k)'.^2;
+x=(a+b)/2+(b-a)/2*bp;
+f=feval(fun,x)*(b-a)/2;
+val=wf(:)'*f(:);
@@ -0,0 +1,37 @@
+function [integral,difference,ratio]=gausstable(a,b,n0,f)
+%
+% function [integral,difference,ratio]=gausstable(a,b,n0,f)
+%
+% This calls GAUSSINT implementation of nth order Gaussian quadrature
+% rule to integrate the function f over the interval [a,b]. The values
+% of n used are
+% n = n0,2*n0,4*n0,...,256*n0
+% The corresponding numerical integrals are returned in the
+% vector integral. The differences of successive numerical
+% integrals are returned in the vector difference:
+% difference(i) = integral(i)-integral(i-1), i=2,...,9
+% The entries in ratio give the rate of decrease in these
+% differences.
+%
+% Ex:
+% [integral,difference,ratio]=gausstable(0,2,2,'exp(-x^2)')
+
+% Convert string to inline function
+f=inline(f);
+
+% Initialize output vectors.
+integral = zeros(9,1);
+difference = zeros(9,1);
+ratio = zeros(9,1);
+
+for i=1:9
+ n = n0*2^(i-1);
+ integral(i) = gaussint(a,b,n,f);
+end
+
+% Calculate the differences of the successive
+% gaussian quadrature rule integrals and the ratio
+% of decrease in these differences.
+difference(2:9)=abs(integral(2:9)-integral(1:8));
+ratio(3:9)=difference(2:8)./difference(3:9);
+
View
@@ -0,0 +1,58 @@
+function [integral,difference,ratio]=simpson(a,b,n0,f)
+%
+% function [integral,difference,ratio]=simpson(a,b,n0,f)
+%
+% This uses Simpson's rule with n subdivisions to integrate the function
+% f over the interval [a,b]. The values of n used are
+% n = n0,2*n0,4*n0,...,256*n0
+% The value of n0 MUST be nonzero and even.
+% The corresponding numerical integrals are returned in the vector
+% integral. The differences of successive numerical integrals are
+% returned in the vector difference:
+% difference(i) = integral(i)-integral(i-1), i=2,...,9
+% The entries in ratio give the rate of decrease in these
+% differences.
+%
+% Ex:
+% [integral,difference,ratio]=simpson(0,2,2,'exp(-x^2)')
+
+% Convert string to inline function
+f=inline(f);
+
+% Initialize output vectors.
+integral = zeros(9,1);
+difference = zeros(9,1);
+ratio = zeros(9,1);
+
+% Initialize for Simpson integration.
+sumend = feval(f,a) +feval(f,b);
+sumodd = 0;
+sumeven = 0;
+
+% Initialize for case of n0 > 2.
+if(n0 > 2)
+ h = (b-a)/n0;
+ for i=2:2:n0-2
+ sumeven = sumeven + feval(f,a+i*h);
+ end
+end
+
+% Calculate the numerical integrals, doing each
+% by appropriately modifying the preceding case.
+for i=1:9
+ n = (n0)*(2^(i-1));
+ h = (b-a)/n;
+ sumeven = sumeven + sumodd;
+ sumodd = 0;
+ for k=1:2:n-1
+ sumodd = sumodd + feval(f,a+k*h);
+ end
+ integral(i) = h*(sumend + 4*sumodd + 2*sumeven)/3;
+end
+
+% Calculate the differences of the successive
+% Simpson rule integrals and the ratio
+% of decrease in these differences.
+difference(2:9)=abs(integral(2:9)-integral(1:8));
+ratio(3:9)=difference(2:8)./difference(3:9);
+
@@ -0,0 +1,56 @@
+function [integral,difference,ratio]=trapezoidal(a,b,n0,f)
+%
+% function [integral,difference,ratio]=trapezoidal(a,b,n0,f)
+%
+% This uses the trapezoidal rule with n subdivisions to
+% integrate the function f over the interval [a,b]. The
+% values of n used are
+% n = n0,2*n0,4*n0,...,256*n0
+% The corresponding numerical integrals are returned in the
+% vector integral. The differences of successive numerical
+% integrals are returned in the vector difference:
+% difference(i) = |integral(i)-integral(i-1)|, i=2,...,9
+% The entries in ratio give the rate of decrease in these
+% differences.
+%
+% Ex:
+% [integral,difference,ratio]=trapezoidal(0,2,2,'exp(-x^2)')
+
+% Convert string to inline function
+f=inline(f);
+
+% Initialize output vectors.
+integral = zeros(9,1);
+difference = zeros(9,1);
+ratio = zeros(9,1);
+
+% Initialize for trapezoidal rule.
+sumend = (feval(f,a) +feval(f,b))/2;
+sum = 0;
+
+% Initialize for case of n0 > 2.
+if(n0 > 2)
+ h = (b-a)/n0;
+ for i=2:2:n0-2
+ sum = sum + feval(f,a+i*h);
+ end
+end
+
+% Calculate the numerical integrals, doing each
+% by appropriately modifying the preceding case.
+for i=1:9
+ n = n0*2^(i-1);
+ h = (b-a)/n;
+ for k=1:2:n-1
+ sum = sum + feval(f,a+k*h);
+ end
+ integral(i) = h*(sumend + sum);
+end
+
+% Calculate the differences of the successive
+% trapezoidal rule integrals and the ratio
+% of decrease in these differences.
+difference(2:9)=abs(integral(2:9)-integral(1:8));
+ratio(3:9)=difference(2:8)./difference(3:9);
+
+
View
View
Oops, something went wrong.

0 comments on commit caaded6

Please sign in to comment.