Matlab implementation of the algorithms presented in the lecture notes of Numerical methods for large algebraic systems

• Chapter 1 Direct solution methods

1. Direct solution methods

 

poisson
k = 1;
x=a(:,1);
alpha = gauss(x,k);
aprod = gaussxma ( k, alpha, a );

 

1. poisson Program to construct a test matrix
2. gauss Computes a Gauss vector
3. gaussxma Computes the product of a Gauss vector with a matrix A
4. ludecom Computes an LU decomposition
5. lsol Computes a solution of a lower triangular system
6. usol Computes a solution of an upper triangular system
7. ucondCondition estimation
8. choldecom Cholesky decomposition

 
% This program fills a matrix with the discretized Poisson equation
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-02-2005

clear a f
disp(' ')
disp(' ')
disp('This program is used to fill the matrix and righthandside')
disp(' ')
disp('This is a discretization of the following equation:')
disp(' ')
disp(' -c  - c   +u1*c  + u2*c  = f  x,y in [0,1]x[0,1],')
disp('   xx   yy      x       y')
disp(' ')
disp(' with Dirichlet boundary conditions everywhere.')
disp(' ')
disp(' ')
n1 = 3;
n2 = 3;
u1 = 0;
u2 = 0;
upwind = 0;
if norm([u1 u2]) > 0
   upwind = input(' Do you want upwind (choose 1) or not (choose 0)
:');
  disp(' ')
   disp(' ')
end
m = n1*n2;
h1 = 1/(n1+1);
h2 = 1/(n2+1);
if upwind ==0
for i = 1:m
   a(i,i) = 2/h1^2+2/h2^2;
end
for i = 1:m-n1
   a(i,i+n1) = -1/h2^2+u2/(2*h2);
   a(i+n1,i) = -1/h2^2-u2/(2*h2);
end
for i = 1:m-1
   a(i,i+1) = -1/h1^2+u1/(2*h1);
   a(i+1,i) = -1/h1^2-u1/(2*h1);
end
end
%
% end while loop
%
if upwind ==1
for i = 1:m
   a(i,i) = 2/h1^2+2/h2^2+abs(u1)/h1+abs(u2)/h2;
end
for i = 1:m-n1
 if u2 > 0
   a(i,i+n1) = -1/h2^2;
   a(i+n1,i) = -1/h2^2-u2/(h2);
 else
   a(i,i+n1) = -1/h2^2+u2/(h2);
   a(i+n1,i) = -1/h2^2;
 end
end
for i = 1:m-1
 if u1 > 0
   a(i+1,i) = -1/h1^2-u1/(h1);
   a(i,i+1) = -1/h1^2;
 else
   a(i,i+1) = -1/h1^2+u1/(h1);
   a(i+1,i) = -1/h1^2;
 end
end
end
%
% end while loop
%

for i = 1:n2-1
   a(i*n1,i*n1+1) = 0;
   a(i*n1+1,i*n1) = 0;
end
for i = 1:n1
   for j = 1:n2
      f(i+(j-1)*n1) = 1+sin(pi*i*h1)*sin(pi*j*h2);
   end
end
%a = sparse(a);
 

 
% This function computes a gauss vector when an original vector
% is given and k is given such that the resulting vector should
% be zero from k+1 to the end  
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-02-2005

function alpha = gauss(x,k);

n = max(size(x));

if x(k) ==0
   disp(' The original vector contains a zero element at position k')
else
   for i = k+1: n
      alpha(i,1) = x(i)/x(k);
   end
end
 

 
% This function computes the product of a gauss vector with a matrix A
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-02-2005

function aprod = gaussxma ( k, alpha, a );

[n,q] = size(a);

for j = 1: q
   for i = 1 : k
      aprod(i,j) = a(i,j);
   end
   for i = k+1: n
      aprod(i,j) = a(i,j) - alpha(i)*a(k,j);
   end
end
 

 
% This function computes an LU decomposition.
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-02-2005

function [l,u] = ludecom(a)

clear l u

n=max(size(a));
u = a;
for k = 1: n-1
   for i = k+1:n
      l(i,k) = u(i,k)/u(k,k);
      for j = k+1:n
         u(i,j) = u(i,j)-l(i,k)*u(k,j);
      end
   end
end
l(n,n) = 0;
l = l+eye(n);
for k = 1: n-1
   for i = k+1:n
      u(i,k) = 0;
   end
end
 

 
% This function computes a solution of a lower triangular system
% of equations.  
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-02-2005

function y = lsol(l,b);

n = max(size(l));
nulpiv = 1;

for i = 1 : n
   if l(i,i) == 0
      disp(' The main diagonal of the lower triangular matrix')
      disp(' contains a zero element')
      nulpiv = 0;
   else
      if nulpiv == 1
         y(i,1) = b(i);
         for j = 1: i-1
            y(i,1) = y(i,1)-l(i,j)*y(j,1);
         end
         y(i,1) = y(i,1)/l(i,i);
      end
   end
end
 

 
% This function computes a solution of an upper triangular system
% of equations. 
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-02-2005

function y = usol(u,b);

n = max(size(u));
nulpiv = 1;

for i = n : -1 : 1
   if u(i,i) == 0
      disp(' The main diagonal of the upper triangular matrix')
      disp(' contains a zero element')
      nulpiv = 0;
   else
      if nulpiv == 1
         y(i,1) = b(i);
         for j = i+1: n
            y(i,1) = y(i,1)-u(i,j)*y(j,1);
         end
         y(i,1) = y(i,1)/u(i,i);
      end
   end
end
 

 
% This function computes an estimate of the condition number of an
% upper triangular matrix
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-02-2005

function con = ucond(u);

n = max(size(u));
nulpiv = 1;

for i = 1 : n
   p(i,1) = 0;
end

for i = n : -1 : 1
   if u(i,i) == 0
      disp(' The main diagonal of the upper triangular matrix')
      disp(' contains a zero element')
      nulpiv = 0;
   else
      if nulpiv == 1
         if p(i) > 0
            y(i,1) = ( -1-p(i,1))/u(i,i);
         else
            y(i,1) = ( 1-p(i,1))/u(i,i);
         end
         for j = 1: i-1
            p(j,1) = p(j,1)+u(j,i)*y(i,1);
         end
      end
   end
end

con = norm(y,inf)*norm(u,inf);
 

 
% This function computes a Choleski decomposition of a symmetric
% postive definite matrix 
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-02-2005

function l = choldecom(a);

n = max(size(a));
nulpiv = 1;

for k = 1 : n
   for p = 1 : k
      l(k,p) = a(k,p);
   end
end

for k = 1 : n
   for p = 1 : k-1
      l(k,k) = l(k,k) - l(k,p)*l(k,p);
   end
   if l(k,k) <= 0
      disp(' The matrix is not postive definite')
      nulpiv = 0;
   else
      if nulpiv == 1
         l(k,k) = sqrt(l(k,k));
         for i = k+1: n
            for p = 1 : k-1
               l(i,k) = l(i,k) - l(i,p)*l(k,p);
            end
            l(i,k) = l(i,k)/l(k,k);
         end
      end
   end
end
 

Contact information:

Back to Numerical methods for large algebraic systems or home page or educational page of Kees Vuik