MATH 541 NUMERICAL ANALYSIS

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)

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

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

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

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

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

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

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)

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;

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

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

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

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

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

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

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

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

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)

NATURAL CUBIC SPLINE

%Natural cubic spline based on Cubic Hermite Formula

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

%this is just the Cubic Hermite Graphing Code

for i=1:n-1,

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

LINEAR APPROXIMATION

%Example of a Linear Approximation

clear,clf reset

%DATA

n=length(xx);

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

xma=xx-xave;

%Solution of Diagonal Normal Equation

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

%Plot the Line and the Data

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)

EXPONENTIAL APPROXIMATION

This program requires the function subprogram FAapproxEx.m

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

clear,clf reset

%DATA

n=length(x);

%LINEARIZING TRANSFORMATION

xx=x; yy=log(y);

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

xma=xx-xave;

%Solution of Diagonal Normal Equation

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

%Plot the Line and the Transformed Data

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

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

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)

y=DATA(2,:);

%L2 Norm function

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

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

end

US OIL PRODUCTION DATA

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

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.

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)
% Quadrature using Simpson`s 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)/(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;

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