MATH 541 NUMERICAL ANALYSIS

CHAPTER 1
 LOSS OF SIGNIFICANCE COSINE EXAMPLE
CHAPTER 3
 BISECTION METHOD FALSE POSITION SECANT  METHOD NEWTON'S METHOD RECURSION COMPUTATION OF PI
CHAPTER 4
 LAGRANGE INTERPOLATION NEWTON FORM DIVIDED DIFFERNCE TABLE PIECEWISE CUBIC HERMITE CHEBYSHEV NODES MOVIE NATURAL CUBIC SPLINE LINEAR APPROXIMATION EXPONENTIAL APPROXIMATION US OIL PRODUCTION DATA WORLD OIL PRODUCTION DATA
CHAPTER 5
CHAPTER 7
 LU FACTORIZATION FORWARD AND BACK SUBSTITUTION TRIDIAGONAL SYSTEM SOLVER JACOBI METHOD GAUSS-SEIDEL METHOD
HOME

LOSS OF SIGNIFICANCE EXAMPLE

%machine precision example
clear
x(1)=1-1/exp(1); s(1)=1; f(1)=1;
for i= 2:21
x(i)=1-(i-1)*x(i-1);
s(i)=1/i;
f(i)=(i-1)*f(i-1);
end
n=0:20;
z=[n;x;s;f];
fprintf(1,'\n\n n x(n) 1/(n+1) n!\n\n')
fprintf(1,'%2.0f %13.8f %10.8f %10.3g\n',z)

 TOP BACK

COSINE EXAMPLE

clear
format compact
x=-1:.001:4;
tic
p81=1-x.^2/2+x.^4/24-x.^6/(24*30)+x.^8/(24*30*56);
toc, tic
p81=1-x.*x/2+x.*x.*x.*x/24- x.*x.*x.*x.*x.*x/(24*30)+x.*x.*x.*x.*x.*x.*x.*x/(24*30*56);
toc, tic
p82=1-x.*x/2.*(1-x.*x/12.*(1-x.*x/30.*(1-x.*x/56)));
toc, tic
xs=x.*x;
p83=1-xs/2.*(1-xs/12.*(1-xs/30.*(1-xs/56)));
toc
p2=1-xs/2;
p4=1-xs/2.*(1-xs/12);
p6=1-xs/2.*(1-xs/12.*(1-xs/30));
p8=1-xs/2.*(1-xs/12.*(1-xs/30.*(1-xs/56)));
plot(x,p2,x,p4,x,p6,x,p8)
grid on
xlabel('x-axis'), ylabel('y-axis')
title('Taylor Series Approx to Cosine')

 TOP BACK

BISECTION METHOD (Note that this code needs to call the f.m file)

%Method of Bisection for solving f(x)=0
%f is computed via the function f.m
clear
%initial guesses
a(1)=1.5; b(1)=2; t1=f(a(1)); t2=f(b(1));
if t1*t2 < 0
%number of iterations
n=20; nn=1:n;
for i=1:n-1,
c=(a(i)+b(i))/2; fc=f(c);
if t1*fc>0
a(i+1)=c; b(i+1)=b(i);
else
a(i+1)=a(i); b(i+1)=c;
end
end
else
'root not trapped'
end
format long
[nn',a',b']
%graph of f
u=-2:.1:2;
plot(u,f(u))
grid on
xlabel('x-axis')
ylabel('y-axis')

 TOP BACK

FALSE POSITION (Note that this code needs to call the file f.m)

%Method of False Position for solving f(x)=0
%f is computed via the function f.m
clear
%initial guesses
a(1)=1.5; b(1)=2; t1=f(a(1)); t2=f(b(1));
if t1*t2 < 0
%number of iterations
n=20; nn=1:n;
for i=1:n-1,
c=(a(i)*f(b(i))-b(i)*f(a(i)))/(f(b(i))-f(a(i))); fc=f(c);
if t1*fc>0
a(i+1)=c; b(i+1)=b(i);
else
a(i+1)=a(i); b(i+1)=c;
end
end
else
'root not trapped'
end
format long
[nn',a',b']
%graph of f
u=-2:.1:2;
plot(u,f(u))
grid on
xlabel('x-axis')
ylabel('y-axis')

 TOP BACK

SECANT METHOD (Note that this code needs to call the file f.m)

%Secant Method for solving f(x)=0
%f is computed via the function f.m
clear
%initial guesses
x(1)=-1.5; x(2)=-1;
%number of iterations
n=10; n=n+1;
for i=2:n,
x(i+1)= (x(i-1)*f(x(i))-x(i)*f(x(i-1)))/(f(x(i))-f(x(i-1)));
end
for i=2:n+1,
diff(i)=x(i)-x(i-1);
end
for i=3:n+1,
ratio(i)=diff(i)/diff(i-1);
qratio(i)=diff(i)/(diff(i-1))^2;
end
nn=1:n+1;
z=[nn;x; diff; ratio; qratio];
fprintf(1,'\n\n n x diff ratio quad ratio\n\n')
fprintf(1,'%2.0f %16.10f %16.10f %12.6f %12.6f\n',z)
%graph of f
u=-5:.01:5;
plot(u,f(u),'k')
axis([-2,2,-3,3])
grid on
xlabel('x-axis')
ylabel('y-axis')

 TOP BACK

NEWTON'S METHOD (Note that this code needs to call two files f.m and df.m)

Code for df.m

function y = df(x)
%computes the function df
y=4*x.^3-1;

%Newton's Method for solving f(x)=0
%f and df are computed via the functions f.m and df.m
clear
%initial guess
x(1)=1.4;
%number of iterations
n=14;
for i=1:n,
x(i+1)= x(i) - f(x(i))/df(x(i));
end
for i=2:n+1,
diff(i)=x(i)-x(i-1);
end
for i=3:n+1,
ratio(i)=diff(i)/diff(i-1);
qratio(i)=diff(i)/(diff(i-1))^2;
end
nn=1:n+1;
z=[nn;x; diff; ratio; qratio];
fprintf(1,'\n\n n x diff ratio quad ratio\n\n')
fprintf(1,'%2.0f %16.14f %16.14f %12.10f %12.10f\n',z)
%graph of f
u=-30:.1:30;
plot(u,f(u))
axis([-2,2,-3,3])
grid on
xlabel('x-axis')
ylabel('y-axis')

 TOP BACK

RECURSION (Note that this code needs the file f.m)

% Iterative Method for solving f(x)=x
%f is computed via the function f.m
clear
%initial guess
x(1)=.5;
%number of iterations
n=27;
diff(1)=0;ratio(1)=0; ratio(2)=0;
for i=1:n,
x(i+1)= f(x(i));
end
for i=1:n,
diff(i)=x(i+1)-x(i);
end
for i=1:n-1,
ratio(i)=diff(i+1)/diff(i);
xacc(i)=x(i) + diff(i)/(1-ratio(i));
end
diff(n+1)=0; ratio(n)=0;ratio(n+1)=0;xacc(n)=0;xacc(n+1)=0;
nn=1:n+1;
z=[nn;x; diff; ratio; xacc];
fprintf(1,'\n\n n x diff ratio x acc\n\n')
fprintf(1,'%2.0f %16.14f %12.10f %12.10f %16.14f\n',z)
%graph of f
u=-2:.1:2;
plot(u,f(u))
grid on
xlabel('x-axis')
ylabel('y-axis')

 TOP BACK

RECOMPUTES PI as in Homework #1. It also computes the differences ratios and accelerated values.

%Computation of Pi
clear
sc(1)=2; sc(2)=2*sqrt(2);
for i=2:29,
sc(i+1)=(sqrt(2)*sc(i))/sqrt(1+sqrt(1-(sc(i)/2^(i))^2));
end
for i=1:29,
diff(i)=sc(i+1)-sc(i);
end
for j=1:28,
ratio(j)=diff(j+1)/diff(j);
xacc(j)= sc(j)+ diff(j)/(1-ratio(j));
end
diff(30)=0; ratio(29)=0;ratio(30)=0;xacc(29)=0;xacc(30)=0;
format long
n=1:30; ns=2.^n;
z=[n;ns;sc; diff; ratio; xacc];
fprintf(1,'\n\n n # sides sc diff ratio sc acc\n\n')
fprintf(1,'%2.0f %10.0f %12.10f %12.10f %12.10f %12.10f\n',z)

 TOP BACK

LU FACTORIZATION FUNCTION

function [LU,row,det1] = lufact(A)
% [LU,row,det1] = lufact(A)
% A is an nxn matrix, input.
% LU is an nxn matrix, output.
% row is the row permutation information, output.
% det1 is the determinant of A, output.
% Partial pivoting is used.
[n n] = size(A); det1 = 1;
row = 1:n; % Initialize pointer vector.
for p = 1:n-1,
for k=p+1:n, % Find the pivot row.
if abs(A(row(k),p)) > abs(A(row(p),p)),
temp = row(p);
row(p) = row(k);
row(k) = temp;
det1 = - det1;
end
end
det1 = det1*A(row(p),p);
if det1 == 0, break, end
for k = p+1:n,
mult = A(row(k),p)/A(row(p),p);
A(row(k),p) = mult;
A(row(k),p+1:n) = A(row(k),p+1:n) - mult*A(row(p),p+1:n);
end
end
det1 = det1*A(row(n),n);
LU = A;

 TOP BACK

IMPLEMENTS FORWARD AND BACKWARD SUBSTITUTION

function [X,Y] = lusolv(A,B,row)
% [X,Y] = lusolv(A,B,row)
% A is the nxn matrix obtained from lufact, input.
% The matrix A must be in factored LU form,
% and row is the row permutation information, input.
% B is any nx1 vector, input.
% The solution to AX=B is X, output.
% The solution to LY=PB is Y, output.
% Forward substitution followed
% by back substitution is used.
[n n] = size(A); Y = zeros(n,1);
Y(1) = B(row(1));
for k = 2:n,
Y(k) = B(row(k)) - A(row(k),1:k-1)*Y(1:k-1);
end
X = zeros(n,1);
X(n) = Y(n)/A(row(n),n);
for k = n-1:-1:1,
X(k) = (Y(k) - A(row(k),k+1:n)*X(k+1:n))/A(row(k),k);
end

 TOP BACK

TRIDIAGONAL SYSTEM SOLVER

function X = tridiag(L,D,U,B)
% X = tridiag(L,D,U,B)
% Solution of a tridiagonal linear system (LDU)X = B.
% It is assumed that D and B have dimension n,
% and that L and U have dimension n-1;
% L sub diagonal, input.
% D diagonal vector, input.
% U super diagonal, input.
% B right hand side vector, input.
% X solution vector, output.
n = length(B);
for k = 2:n,
mult = L(k-1)/D(k-1);
D(k) = D(k) - mult*U(k-1);
B(k) = B(k) - mult*B(k-1); %forward substitution
end
X(n) = B(n)/D(n);
for k = (n-1):-1:1,
X(k) = (B(k) - U(k)*X(k+1))/D(k); %backward subtitution
end

 TOP BACK

JACOBI METHOD

function X = fjacobi(A,B,P)
% X = fjacobi(A,B,P)
% A is an nxn diagonally dominant matrix, input.
% B is a nx1 vector, input.
% P is an nX1 starting vector, input.
d = diag(A);
D = diag(d);
X = D\(B - (A-D)*P);

 TOP BACK

GAUSS-SEIDEL Method

function X = fgsiedel(A,B,P)
% X = fgsiedel(A,B,P)
% A is an nxn diagonally dominant matrix, input.
% B is a nx1 vector, input.
% P is an nX1 starting vector, input.
d = diag(A);
D = diag(d);
OFFD = A - D;
n = length(d);
X=P;
for j=1:n
X(j,1) = (B(j) - OFFD(j,:)*X)/d(j);
end

 TOP BACK

LAGRANGE INTERPOLATION

clear
hold off
DATA=[.25 1 1.5 2 2.4 5; 23.1 1.68 1 .84 .83 1.26];
n=length(DATA(1,:)); XDATA= DATA(1,:);
flops(0);
X=0:.1:6;
m=length(X); Z=zeros(1,m);
for i=1:n
Y=lagrange(X,XDATA,i);
Z=Z+DATA(2,i).*Y;
end
flops
plot(X,Z,'k',DATA(1,:),DATA(2,:),'r*')
axis([0,6,-5,30])
grid on
%Y=lagrange(X,XDATA,6);
%plot(X,Y,'b',DATA(1,:),zeros(1,n),'r*')
%axis([0,6,-2,2])
%grid on

LAGRANGE  BASIS FUNCTIONS

function Y = lagrange(X,XDATA,i)
%this function calculated the ith lagrange polynomial determined by the
%XDATA and return the y-values for the x-values in X
nn=length(XDATA); m=length(X);
Y=ones(1,m);
for k = 1:nn
if k ~= i
Y=Y.*(X - XDATA(k))./(XDATA(i)-XDATA(k));
end
end

 TOP BACK

NEWTON FORM

clear
hold off
DATA=[.25 1 1.5 2 2.4 5; 23.1 1.68 1 .84 .83 1.26];
n=length(DATA(1,:));
flops(0);
A=divdif(DATA');
X=0:.1:6;
m=length(X);
Z=A(n).*ones(1,m);
for i=n-1:-1:1
Z=A(i)+(X-DATA(1,i)).*Z;
end
flops
plot(X,Z,'k',DATA(1,:),DATA(2,:),'r*')
axis([0,6,-5,30])
grid on

 TOP BACK

DIVIDED DIFFERENCE TABLE

function A = divdif(DA)
%this function calculates the divided diff table for DA;
%DA is the n,2 input data;
%A is the vector of coef of the Newton Form;
[n,m] = size(DA);
D=zeros(n,n); D(:,1)= DA(:,2);
X=DA(:,1);
for j=2:n
for i=j:n
D(i,j)= (D(i,j-1)-D(i-1,j-1))/(X(i)-X(i-j+1));
end
end
%D
A=diag(D);

 TOP BACK

PIECEWISE CUBIC HERMITE

%piecewise hermite grapher
clear
x=[.25 1 1.5 2 2.4 5];
y=[23.1 1.68 1 .84 .83 1.26];
dy=[-50 -4 -0.5 0 0.1 0.5];
n=length(x);
t=[]; u=[];
for i=1:n-1,
dx=x(i+1)-x(i); r=dx/10; a=(y(i+1)-y(i))/dx;
b=(a-dy(i))/dx; c=(dy(i+1)-a)/dx; d=(c-b)/dx;
p=x(i):r:x(i+1);
s=y(i)+(p-x(i)).*(dy(i)+(p-x(i)).*(b+(p-x(i+1)).*d));
t=[t, p];
u=[u, s];
end
plot(t,u,x,y,'x')

 TOP BACK

CHEBYSHEV POLYNOMIAL MOVIE

clear all
n=6; p=1/2^(n-1);
NZ=[1 .1 1 0 1 -1]
CZ=2*n-1:-2:1;
CZ=cos(CZ*pi/(2*n))
t=0:.1:1;
M=moviein(11);
for j=1:11
Z=(1-t(j))*CZ+t(j)*NZ;
X=-1.5:.025:1.5;
Y=X-Z(1);
for i=2:n
Y=Y.*(X-Z(i));
end
plot(X,Y,'b',Z,zeros(1,n),'r*',[-1,1],[p,p],[-1,1],[-p,-p])
axis([-1.5,1.5,-1.5*p,1.5*p])
grid on
M(:,j)=getframe;
end
movie(M,-10,.8)

 TOP BACK

NATURAL CUBIC SPLINE

%Natural cubic spline based on Cubic Hermite Formula

clear
clf reset
x=[.25 1 1.5 2 2.4 5]
y=[23.1 1.68 1 .84 .83 .26]
n=length(x);

for i=1:n-1,

h(i)=x(i+1)-x(i); delta(i)= (y(i+1)-y(i))/h(i);

end

%Construct the Tridiagonal System

hp=h(2:n-1); hm=h(1:n-2);
deltap=delta(2:n-1); deltam=delta(1:n-2);
U=[1,hm]; L=[hp,1]; D=2*[1,hp+hm,1];
B=3*[delta(1), hp.*deltam+hm.*deltap,delta(n-1)];
dy = tridiag(L,D,U,B);

%this is just the Cubic Hermite Graphing Code

t=[];
u=[];

for i=1:n-1,

r=h(i)/10; a=delta(i);
b(i)=(a-dy(i))/h(i); c=(dy(i+1)-a)/h(i); d(i)=(c-b(i))/h(i);
p=x(i):r:x(i+1);
s=y(i)+(p-x(i)).*(dy(i)+(p-x(i)).*(b(i)+(p-x(i+1)).*d(i)));
t=[t, p];
u=[u, s];

end
plot(t,u,'b',x,y,'*r')
grid on

 TOP BACK

LINEAR APPROXIMATION

%Example of a Linear Approximation

clear,clf reset

%DATA

DATA=[1994:2002;

21375 22253 22818 23280 24097 24084 24139 26455 27113]

xx=DATA(1,:); yy=DATA(2,:);

n=length(xx);

%LINEAR FIT OF DATA Shifted by Ave M=SLOPE B=INTERCEPT

xave= sum(xx)/n, yave= sum(yy)/n;

xma=xx-xave;

%Solution of Diagonal Normal Equation

M=(xma*yy')/(xma*xma'),B=yave

%Plot the Line and the Data

t=linspace(xx(1)-1,xx(n)+1,20); s=B+M*(t-xave);

plot(t,s,'k',DATA(1,:),DATA(2,:),'b*','LineWidth',2,'Markersize', 8)

grid on, xlabel('x-axis'), ylabel('y-axis')

%Plot the equation on the graph

EQ=['y = ',num2str(B),' + ',num2str(M),' (x - ',num2str(xave),')'];

gtext(EQ,'FontSize',14)

 TOP BACK

EXPONENTIAL APPROXIMATION

This program requires the function subprogram FAapproxEx.m

%Example of an Approximation using Linearizing Transformation

%Want to fit y = C exp(a(x-1994))

clear,clf reset

%DATA

DATA=[1994:2002;

21375 22253 22818 23280 24097 24084 24139 26455 27113]

x=DATA(1,:)-1994; y=DATA(2,:);

n=length(x);

%LINEARIZING TRANSFORMATION

xx=x; yy=log(y);

%LINEAR FIT OF TRANS DATA M=SLOPE B=INTERCEPT

xave= sum(xx)/n, yave= sum(yy)/n;

xma=xx-xave;

%Solution of Diagonal Normal Equation

M=(xma*yy')/(xma*xma'),B=yave

%Plot the Line and the Transformed Data

t=linspace(xx(1)-1,xx(n)+1,20); s=B+M*(t-xave);

subplot(2,1,2)

plot(t,s,'k',xx,yy,'b*','LineWidth',2,'Markersize', 8)

grid on, xlabel('x-axis'), ylabel('y-axis')

%Plot the equation on the graph

EQ=['y = ',num2str(B),' + ',num2str(M),' (x - ',num2str(xave),')'];

gtext(EQ,'FontSize',14)

%Transform B to xx coords from x minus ave coords

B=B-M*xave

%Transform Back

C=exp(B), a=M

t=1993:.1:2003; s=C*exp(a*(t-1994));

hold on

subplot(2,1,1)

plot(t,s,'k',DATA(1,:),DATA(2,:), 'b*', 'linewidth',2,'markersize',8)

grid on, xlabel('x-axis'), ylabel('y-axis')

%Plot the equation of the graph

EQ=['y = ',num2str(C),' exp ( ',num2str(a), ' (x - 1994))'];

gtext(EQ, 'fontsize',14)

%Solve the Nonlinear Appro Problem

D=fminsearch('FApproxEx',[C,a])

hold on

subplot(2,1,1)

plot(t,D(1)*exp(D(2)*(t-1994)),'r','linewidth',2)

%Plot the equation on the graph

EQ=['y = ',num2str(D(1)),' exp ( ',num2str(D(2)), ' (x - 1994))'];

gtext(EQ,'color','r','fontsize',14)

FUNCTION SUBPROGRAM

function z=FApproxEx(D)

DATA=[1994:2002;

21375 22253 22818 23280 24097 24084 24139 26455 27113]

x=DATA(1,:)-1994;

y=DATA(2,:);

%L2 Norm function

z=(D(1)*exp(D(2)*x)-y)*(D(1)*exp(D(2)*x)-y)';

%L1 Norm function

%z=sum(abs(D(1)*exp(D(2)*x)-y));

%Linf Norm function

%z=max(abs(D(1)*exp(D(2)*x)-y));

end

 TOP BACK

US OIL PRODUCTION DATA

This is the total oil produced in millions of barrels within the US through the indicated year.

1919 , 5055 ; 1920 , 5497 ; 1921 , 5968 ; 1922 , 6526 ; 1923 , 7259 ; 1924 , 7972 ; 1925 , 8736 ; 1926 , 9507 ; 1927 , 10409 ; 1928 , 11310 ; 1929 , 12317 ; 1930 , 13213 ; 1931 , 14063 ; 1932 , 14847 ; 1933 , 15752 ; 1934 , 16658 ; 1935 , 17650 ; 1936 , 18737 ; 1937 , 20002 ; 1938 , 21204 ; 1939 , 22457 ; 1940 , 23792 ; 1941 , 25181 ; 1942 , 26555 ; 1943 , 28047 ; 1944 , 29716 ; 1945 , 31421 ; 1946 , 33149 ; 1947 , 34999 ; 1948 , 37001 ; 1949 , 38825 ; 1950 , 40776 ; 1951 , 42994 ; 1952 , 45250 ; 1953 , 47563 ; 1954 , 49877.988 ; 1955 , 52362.416 ; 1956 , 54979.699 ; 1957 , 57596.6 ; 1958 , 60045.587 ; 1959 , 62620.177 ; 1960 , 65195.11 ; 1961 , 67816.868 ; 1962 , 70493.057 ; 1963 , 73245.78 ; 1964 , 76032.602 ; 1965 , 78881.116 ; 1966 , 81908.879 ; 1967 , 85124.621 ; 1968 , 88294.207 ; 1969 , 91498.203 ; 1970 , 94848.869 ; 1971 , 98145.481 ; 1972 , 101438.88 ; 1973 , 104644.892 ; 1974 , 107701.828 ; 1975 , 110624.482 ; 1976 , 113470.437 ; 1977 , 116344.804 ; 1978 , 119393.533 ; 1979 , 122379.584 ; 1980 , 125384.615 ; 1981 , 128369.009 ; 1982 , 131384.3 ; 1983 , 134555.299 ; 1984 , 137804.995 ; 1985 , 141079.548 ; 1986 , 144247.8 ; 1987 , 147295.178 ; 1988 , 150274.301 ; 1989 , 153053.074 ; 1990 , 155737.761 ; 1991 , 158444.8 ; 1992 , 161069.432 ; 1993 , 163568.465 ; 1994 , 165999.941 ; 1995 , 168394.209 ; 1996 , 170760.226 ; 1997 , 173115.057 ; 1998 , 175396.976 ; 1999 , 177543.708 ; 2000 , 179674.415 ; 2001 , 181791.926 ; 2002 , 183915.105 ; 2003 , 186009.475 ;

 TOP BACK

WORLD OIL PRODUCTION DATA

The second column is the annual production and the third column is the cumulative production for the indicated year. All amounts are in millions of barrels.

1865 2.56        6.39
1875 13.14        84.86
1885 36.87        334.89
1895 103.66        1037.51
1905 214.99        2630.74
1915 432.16        5866.46
1925 1069.09        13372.69
1935 1654.91        26992.66
1945 2594.79        48241.14
1950 3832.50        64309.35
1955 5840.00        88490.60
1960 7661.35        122244
1961 8194.25        130438.25
1962 8887.75        139326
1963 9537.45        148863.45
1964 10285.7        159149.15
1965 11070.45       170219.6
1966 12030.4        182250
1967 12917.35       195167.35
1968 14099.95       209267.3
1969 15220.5        224487.8
1970 16749.85       241237.65
1971 17709.8        258947.45
1972 18666.1        277613.55
1973 20323.2        297936.75
1974 20337.8        318274.55
1975 19282.95       337557.5
1976 20929.1        358486.6
1977 21794.15       380280.75
1978 21958.4        402239.15
1979 22874.55       425113.7
1980 21754          446867.7
1981 20469.2        467336.9
1982 19520.2        486857.1
1983 19439.9        506297
1984 19888.85       526185.85
1985 19702.7        545888.55
1986 20523.95       566412.5
1987 20684.55       587097.05
1988 21440.1        608537.15
1989 21848.9        630386.05
1990 22106.59       652492.6
1991 21975.55       674468.2
1992 21977.75       696445.9
1993 21986.14       718432.1
1994 22261.72       740693.8
1995 22752.28       763446.1
1996 23254.52       786700.6
1997 23976.82       810677.4
1998 24426.34       835103.7
1999 24034.34       859138.1
2000 24945.54       884083.6
2001 24774.26       908857.9
2002 24376.22       933234.1
2003 25241.19       958475.3
2004 26458.01       984933.3
2005 26816.56      1011749.9

 TOP BACK

INTEGRATION M-FILES

MAIN

%computes the integral of f using Gaussian 2pt
%computes the accelerated value to determine error
%a = lower limit of integration, b = upper limit
%m = number of subintervals
a=-1; b=1; m=4; nn=[];
for i=1:7,
m=2*m; nn=[nn,m];
x(i)=gauss2(a,b,m);
end
diff(1)=0;ratio(1)=0; ratio(2)=0;xacc(1)=0;xacc(2)=0;
n=length(x);
for i=2:n,
diff(i)=x(i)-x(i-1);
end
for i=3:n,
ratio(i)=diff(i)/diff(i-1);
xacc(i)=x(i-2) + diff(i-1)/(1-ratio(i));
end
z=[nn;x; diff; ratio; xacc];
fprintf(1,'\n\n n x diff ratio x acc\n\n')
fprintf(1,'%2.0f %16.14f %12.10f %12.10f %16.14f\n',z)
s = (b - a)/300;  t = a:s:b;
y = f(t);
plot(t,y, 'k')

FUNCTION

function y = f(x)
%computes the integrand function f
y=sqrt(3*sin(x)+4);

TRAPEZOID

function s = trap(a,b,m)
% s = trap(a,b,m)
% Quadrature using the trapezoidal rule.
% f is the name of the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% m is the number of subintervals, input.
a; b; m;
h = (b - a)/m; s = 0;
for k=1:(m-1),
x = a + h*k;
s = s + f(x);
end
s = h*(f(a)+f(b))/2 + h*s;

SIMPSON'S

function s = simp(a,b,m)
% s = simp(a,b,m)
% f is the name of the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% m is the number of subintervals, input.
a; b; m; h = (b - a)/(2*m); s1 = 0; s2 = 0;
for k=1:m,
x = a + h*(2*k-1);
s1 = s1 + f(x);
end
for k=1:(m-1),
x = a + h*2*k;
s2 = s2 + f(x);
end
s = h*(f(a)+f(b)+4*s1+2*s2)/3;

GAUSSIAN 2PT

function s = gauss2(a,b,m)
% s = gauss2(a,b,m)
% Quadrature using the Gaussian 2pt rule.
% f is the name of the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% m is the number of subintervals, input.
a; b; m; h = (b - a)/m; s = 0;
x = a + h*(1-1/sqrt(3))/2;
y = a + h*(1+1/sqrt(3))/2;
for k=1:m,
s = s + f(x) + f(y);
x = x + h; y = y + h;
end
s = h*s/2;

GAUSSIAN 4PT

function s = gauss4(a,b,m)
% s = gauss4(a,b,m)
% Quadrature using the Gaussian 4pt rule.
% f is the name of the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% m is the number of subintervals, input.
a; b; m; h = (b - a)/m; s = 0;
x = a + h*(1-.861136311594053)/2;
y = a + h*(1-.339981043584856)/2;
z = a + h*(1+.339981043584856)/2;
t = a + h*(1+.861136311594053)/2;
for k=1:m,
s = s + .347854845137454*(f(x)+f(t)) + .652145154862546*(f(y)+f(z));
x = x + h; y = y + h; z = z + h; t = t + h;
end
s = h*s/2;

 TOP BACK

Revised 02/17/06
By Don Short, dshort@sciences.sdsu.edu