Monday, May 27, 2013

Simpson composite integration algorithm in Matlab

Simpson composite integration algorithm in Matlab :

function simps(a, b, n)


%simps(a, b, n) approximates the integral of a function f(x) in the 
%interval [a;b] by the composite simpson rule
%n is the number of subintervals
%the user needs to specify the function f(x) at the bottom

%Author: Alain G. Kapitho
%Date  : Jan 2006

h = (b-a)/n;

sum_even = 0;

for i = 1:n/2-1
   x(i) = a + 2*i*h;
   sum_even = sum_even + f(x(i));
end

sum_odd = 0;

for i = 1:n/2
   x(i) = a + (2*i-1)*h;
   sum_odd = sum_odd + f(x(i));
end

integral = h*(f(a)+ 2*sum_even + 4*sum_odd +f(b))/3


%this needs to be changed accordingly with the specific problem you have at
%hand, before proceeding to the command line
function y = f(x)
y = exp(x);

Composite trapezoidal rule algorithm in Matlab

Composite trapezoidal rule algorithm in Matlab :

function trapez(a, b, n)

%trapez(a, b, n) approximates the integral of a function f(x) in the 
%interval [a;b], by the composite trapezoidal rule
%n is the number of subintervals
%the user needs to specify the function f(x) at the bottom

%Author: Alain G. Kapitho
%Date  : Jan 2006

h = (b-a)/n;

sum = 0;

for i = 1:n-1
   x(i) = a + i*h;
   sum = sum + f(x(i));
end

integral = h*(f(a) + 2*sum + f(b))/2

%this needs to be changed accordingly with the specific problem you have at
%hand, before proceeding to the command line
function y = f(x)
y = 2 + sin(2*sqrt(x));

Newton interpolation method algorithm in Matlab

Newton interpolation method algorithm in Matlab :

function [yy,E,d,n_iter]=interpola_newton(x,y,xx,error)
% Syntaxis: [yy,d,n_iter]=interpola_newton(x,y,xx,error)
% ----------------------------------------------------------------------
% Devuelve el valor interpolado yy en un determinado punto xx solicitado 
% mediante el metodo de las diferencia divididas de Newton, empleando (n-1)
% nodos para la interpolación y un nodo para determinar el error de la 
% interpolacion.
%
% Entrada:
% x     --> Vector que contiene las abcisas de los nodos de interpolacion 
% y     --> Vector que contiene las ordenadas de los nodos de interpolacion
% xx    --> Abcisa del nodo que se desea incorporar
% error --> Tolerancia de la interpolacion {default | 0.05}
%         
%
% Salida: 
% yy     --> Valor interpolado para el nodo con abcisa xx
% E      --> Error real de la interpolacion
% d      --> Tabla de Diferencias Divididas
% n_iter --> Numero de iteraciones necesitadas para verificar la condicion 
%            de parada
%
%
% Ejemplo:
% 
% x=[1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30];
% y=[1.00000, 0.97350, 0.95135, 0.93304, 0.91817, 0.90640, 0.89747];
% xx=1.17; error=0.00005;
% [yy,E,d,n_iter] = interpola_newton(x,y,xx,error)
%
% yy = 0.9267
%  
% E = 2.0944e-05
%  
% d =
% 
%    -0.5300    0.8700   -0.6800    0.7333   -1.3333    6.2222
%    -0.4430    0.7680   -0.5333    0.4000    0.5333         0
%    -0.3662    0.6880   -0.4533    0.5333         0         0
%    -0.2974    0.6200   -0.3467         0         0         0
%    -0.2354    0.5680         0         0         0         0
%    -0.1786         0         0         0         0         0
%  
% n_iter = 4
%
%
% Ver tambien: interpola_lagrange
%
% ----------------------------------------------------------------------
%
% Elaborado por:
%
% Msc. Alexeis Companioni Guerra
% Lic. Vianka Orovio Cobo
%
% Revision: 1.0 
% Fecha: 30/09/2007


%% Validacion de las entradas
if nargin<3
    error('Verifique los datos de entrada.');
elseif nargin<4
    error = 0.05;
end

if min(size(x))>1 || min(size(y))>1
    error('Tanto x como y deben ser vectores.');
elseif max(size(x)) ~= max(size(y))
    error('Verifique la dimension de los vectores.');
else
    x = x(:); y = y(:);
    n = length(x);
end

if sum(size(xx))~=2
   error('La abcisa del nodo de entrada debe ser un escalar.'); 
end

if sum(size(error))~=2
   error('El error maximo de la interpolacion debe ser un escalar.'); 
end

%% Conformacion de las Diferencias Divididas
d = diff(y)./diff(x);
for j=2:n-1
   for k=1:n-j
        d(k,j) = ( d(k+1,j-1) - d(k,j-1) ) / ( x(k+j) - x(k) );
   end
end

%% Seleccion de la primera fila de la matriz de Diferencias Divididas
p = [y(1),d(1,:)];

%% Conformacion iterativa de la solucion 
yy = p(1);
producto = xx-x(1); 
E = p(2)*producto;
i = 2;
while abs(E)>error && i<n
    yy = yy + E;
    producto = producto*(xx-x(i));
    i = i + 1;
    E = p(i)*producto;
end
n_iter = i - 1;

Lagrange interpolation algorithm in Matlab

Lagrange interpolation algorithm in Matlab :

function [yy]=interpola_lagrange(x,y,xx)
% Syntaxis: [yy]=interpola_lagrange(x,y,xx)
% ----------------------------------------------------------------------
% Devuelve el valor interpolado yy en un determinado punto xx solicitado 
% mediante el metodo de interpolacion de Lagrange.
%
% Entrada:
% x  --> Vector que contiene las abcisas de los nodos de interpolacion 
% y  --> Vector que contiene las ordenadas de los nodos de interpolacion
% xx --> Abcisa del nodo que se desea incorporar 
%         
%
% Salida: 
% yy --> Valor interpolado para el nodo con abcisa xx
%
%
% Ejemplo:
% 
% x=[1; 2; 4]; y=[2; 3; 1]; xx=2.5;
% [yy]=interpola_lagrange(x,y,xx)
%
% yy =
%      3
%
%
% Ver tambien: interpola_newton
%
% ----------------------------------------------------------------------
%
% Elaborado por:
%
% Msc. Alexeis Comanioni Guerra
% Lic. Vianka Orovio Cobo
%
% Revision: 1.0 
% Fecha: 30/09/2007


%% Validación de las entradas
if nargin<3
    error('Verifique los datos de entrada.');
end

if min(size(x))>1 || min(size(y))>1
    error('Tanto x como y deben ser vectores.');
else
    x=x(:); y=y(:);
    n=length(x);
end

if sum(size(xx))~=2
   error('La abcisa del nodo de entrada debe ser un escalar.'); 
end

%% Creación de los polinomios de Lagrange
U=ones(n-1,1); 
for i=1:n
    x_temp=x; 
    x_temp(find(x_temp==x_temp(i)))=[];
    L(i)=prod( xx*U-x_temp ) / prod( x(i)*U-x_temp );
end

%% Obtencion del valor interpolado
yy=dot(y,L);

Seidel method algorithm in Matlab

Seidel method algorithm in Matlab :

function [Sol,E,iter,csuf]=seidel(M,b,x0,error,n_iter)
% Syntaxis: [Sol,E,iter,csuf]=seidel(M,b,x0,error,n_iter)
% ----------------------------------------------------------------------
% Esta funcion resuelve un sistema de ecuaciones lineales (SEL) de la
% forma Ax=b por el metodo de Seidel segun
%
% x0 in V (arbitrario)
% x(k+1)=( I-inv(D-E)*A )*x(k) + inv(D-E)*b
%
% Entrada:
% M   --> Matriz del sitema
% b   --> Vector independiente (debe ser un vector columna)
% x0  --> Vector de aproximaciones iniciales
% error  --> Tolerancia del calculo {default | 0.00001}
% n_iter --> Numero maximo de iteraciones {default | 1000} 
%
% Salida: 
% Sol  --> Vector solucion del SEL
% E    --> Error de la aproximacion
% csuf --> Cumplimiento de la condición suficiente de convergencia "alpha<1"
%
%
% Ejemplo:
% 
% M=[9 -1 2; 1 8 2; 1 -1 11]; b=[9; 19; 10];
% [Sol,E,iter,csuf]=seidel(M,b)
%
% Sol =
%      1.0000
%      2.0000
%      1.0000
%
% E = 1.2249e-06
%
% iter = 6
% 
% csuf = 1
%
%
% Ver tambien: jacobi
%
% ----------------------------------------------------------------------
%
% Elaborado por:
%
% Msc. Alexeis Companioni Guerra
% Lic. Vianka Orovio Cobo
%
% Revision: 1.0 
% Fecha: 30/09/2007


%% Validacion de las entradas
if nargin<2
    error('Verifique los datos de entrada.');
else
    if det(M) == 0
        error('M debe ser no singular.');
    end
    
    if diff(size(M))~=0
        error('La matriz de entrada debe ser cuadrada.');
    end
    
    if min(size(b))>1
        error('b debe ser un vector.');
    end
    
    if max(size(M)) ~= max(size(b))
        error('Verifique la dimension del vector independiente respecto a la matriz del sistema.');
    end
end

b = b(:);
n = length(b);

if nargin<3
    x0 = zeros(n,1);
else
    if min(size(x0))>1
        error('x0 debe ser un vector.');
    end
    
    if max(size(M)) ~= max(size(x0))
        error('Verifique la dimension del vector x0 respecto a la matriz del sistema.');
    end
end

if nargin<4
    error = 0.00001;
else
    if sum(size(error))~=2
        error('La tolerancia debe ser un escalar.'); 
    end
end

if nargin<5
    n_iter = 1000;
else
    if sum(size(n_iter))~=2
        error('El numero de iteraciones debe ser un escalar.'); 
    end
    
    if round(n_iter)-n_iter~=0
       error('El numero de iteraciones debe ser un entero.');
    end
end
   
%% Determinacion del factor de convergencia "beta"
Jac = eye(n)-( diag( diag(M).^(-1) )*M ); 
Jac_up = triu(Jac,1); 
Jac_lo = tril(Jac,-1);

beta = max( sum(abs(Jac_up),2) ./ ( ones(n,1)-sum(abs(Jac_lo),2) ) );
if beta>=1
    csuf = 0;
else
    csuf = 1;
end

%% Ciclo iterativo
iter = 0; Sol = zeros(n,1); E = inf; 
s_inv = inv(diag(diag(M)) - tril(M,-1).*(-1));
S = eye(n)-s_inv*M;        % Matriz de Seidel 

while (E>=error) && (iter<n_iter) 
    Sol = S*x0 + s_inv*b;
    E = abs(beta/(1-beta))*norm(Sol-x0,inf);
    x0 = Sol; 
    iter = iter + 1;
end

Jacobi iterative method algorithm in Matlab

Jacobi iterative method algorithm in Matlab :

function jacobi(A, b, N)

%Jacobi(A, b, N) solve iteratively a system of linear equations whereby
%A is the coefficient matrix, and b is the right-hand side column vector.
%N is the maximum number of iterations.
%The method implemented is the Jacobi iterative. 
%The starting vector is the null vector, but can be adjusted to one's needs.
%The iterative form is based on the Jacobi transition/iteration matrix
%Tj = inv(D)*(L+U) and the constant vector cj = inv(D)*b.
%The output is the solution vector x. 

%This file follows the algorithmic guidelines given in the book 
%Numerical Analysis, 7th Ed, by Burden & Faires

%Author: Alain G. Kapitho
%Date  : Dec. 2005
%Rev.  : Aug. 2007

n = size(A,1);
%splitting matrix A into the three matrices L, U and D
D = diag(diag(A));
L = tril(-A,-1);
U = triu(-A,1);

%transition matrix and constant vector used for iterations
Tj = inv(D)*(L+U);
cj = inv(D)*b;

tol = 1e-05;
k = 1;
x = zeros(n,1);   %starting vector

while k <= N
   x(:,k+1) = Tj*x(:,k) + cj;
   if norm(x(:,k+1)-x(:,k)) < tol
      disp('The procedure was successful')
      disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations')
      disp(k); disp('x = ');disp(x(:,k+1));
      break
   end
   k = k+1;
end

if norm(x(:,k+1)- x(:,k)) > tol || k > N
   disp('Maximum number of iterations reached without satisfying condition:')
   disp('||x^(k+1) - x^(k)|| < tol'); disp(tol);
   disp('Please, examine the sequence of iterates')
   disp('In case you observe convergence, then increase the maximum number of iterations')
   disp('In case of divergence, the matrix may not be diagonally dominant')
   disp(x');
end

Gauss-Seidel iterative method algorithm in Matlab

Gauss-Seidel iterative method algorithm in Matlab :

function gauss_seidel(A, b, N)

%Gauss_seidel(A, b, N) solve iteratively a system of linear equations 
%whereby A is the coefficient matrix, and b is the right-hand side column vector.
%N is the maximum number of iterations.
%The method implemented is the Gauss-Seidel iterative. 
%The starting vector is the null vector, but can be adjusted to one's needs.
%The iterative form is based on the Gauss-Seidel transition/iteration matrix
%Tg = inv(D-L)*U and the constant vector cg = inv(D-L)*b.
%The output is the solution vector x.

%This file follows the algorithmic guidelines given in the book 
%Numerical Analysis, 7th Ed, by Burden & Faires

%Author: Alain G. Kapitho
%Date  : Dec. 2005
%Rev.  : Aug. 2007


n = size(A,1);
%splitting matrix A into the three matrices L, U and D
D = diag(diag(A));
L = tril(-A,-1);
U = triu(-A,1);

%transition matrix and constant vector used for iterations
Tg = inv(D-L)*U; 
cg = inv(D-L)*b;

tol = 1e-05;
k = 1;
x = zeros(n,1);   %starting vector

while k <= N
   x(:,k+1) = Tg*x(:,k) + cg;
   if norm(x(:,k+1)-x(:,k)) < tol
      disp('The procedure was successful')
      disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations')
      disp(k); disp('x = ');disp(x(:,k+1));
      break
   end
   k = k+1;
end

if norm(x(:,k+1)- x(:,k)) > tol || k > N
   disp('Maximum number of iterations reached without satisfying condition:')
   disp('||x^(k+1) - x^(k)|| < tol'); disp(tol);
   disp('Please, examine the sequence of iterates')
   disp('In case you observe convergence, then increase the maximum number of iterations')
   disp('In case of divergence, the matrix may not be diagonally dominant')
   disp(x');
end

Gauss Jordan Elimination & Pivoting algorithm in Matlab

Gauss Jordan Elimination & Pivoting algorithm in Matlab :

%Gauss Jordan elimination with pivoting
%copy this file to the command window and press enter
%Bangladesh University of Engineering & Technology
%Nadim Chowdhury
%Department Of Electrical & Electronics Engineering
%if You have any problem in understanding this  feel free to e-mail me At
%nadim_eee_buet@yahoo.com
A=0;
x=0;
n=input('How many variables=');
disp('Enter the coefficients along with constants For instance if x+y+3z=-5 then enter 1 1 3 -5 each number followed by an enter not space');
for i=1:n
    for j=1:n+1
        A(i,j)=input('');
    end
end
%pivoting
for i=1:n-1
    for j=i+1:n
        if abs(A(j,i))>abs(A(i,i))
            T=A(j,:);
            A(j,:)=A(i,:);
            A(i,:)=T;
        end
    end
end
disp('After pivoting');
disp(A);
for k=1:n-1
    for i=k+1:n
        m=A(i,k)/A(k,k);
        for j=k:n+1
            A(i,j)=A(i,j)-m*A(k,j);
        end
    end
end
disp('Triangularize Form  ');
disp(A);
            
if A(n,n)==0
    disp('No unique solution');
end
    x(n)=A(n,n+1)/A(n,n);
    for j=n-1:-1:1
        sum=0;
    for i=1:n
        sprintf('x%.0f=%.10f',i,x(i))
    end
        

        for i=1:n-j
            sum=sum+A(j,n+1-i)*x(n+1-i);
        end
        x(j)=(A(j,n+1)-sum)/A(j,j);
    end    

Newton Method algorithm in Matlab

Newton Method algorithm in Matlab :

function [xvect,xdif,fx,nit]=mynewton(x0,nmax,fun,dfun,toll);
% MYNEWTON   Do the newton iteration to find the zeros of the given
% inline scalar function. 
%    [XVEC,XDIF,FX,NIT] = MYNEWTON(X0,NMAX,FUN,DFUN,TOLL)  
%       Input:  x0: starting value
%             nmax: maximum number of iteration
%             toll: tolerance, default is 1e-10
%              fun: given inline function
%             dfun: derivative of the function given as inline
%             function 
%
%    Output: xvect: stores values in all iterations (arg)
%             xdif: difference between two consequent values (arg) 
%               fx: stores function values  in all iterations 
%              nit: number of iteration required
%
%      Examples:
%         The first examples converge to unique zero
%            fun=inline('x^3+x^2-4');  x0=0.3; nmax=100;
%           dfun=inline('3*x^2+2*x');
%
%            fun=inline('x^4+x^3-5*x-12'); 
%           dfun=inline('4*x^3+3*x^2-5');
%
%            fun=inline('x^5+10*x^2-9*x+10');
%           dfun=inline('5*x^4+20*x-9');
%
%         Examples where different starting values converge to
%         different zeros
%            fun=inline('x^2-4*sin(x)');
%           dfun=inline('2*x-4*cos(x)');
%
%            fun=inline('x^3+x^2*(cos(x))^2-4');
%           dfun=inline('(3*x^2+2*x*cos(x)^2-2*x^2*cos(x)*sin(x))');
 
% Author: Bishnu Lamichhane, University of Stuttgart

if (nargin==4)  toll=1e-10; end

err=toll+1; nit=0; xvect=x0;x=x0;fx=feval(fun,x);xdif=[];
  while(nit<nmax & err>toll)
    nit=nit+1; x=xvect(nit);dfx=feval(dfun,x);
    if (dfx==0), err=toll*1e-10;
      disp('Stop for vanishing dfun');
    else,
      xn=x-fx(nit)/dfx;err=abs(xn-x);xdif=[xdif;err];
      x=xn;xvect=[xvect;x];fx=[fx;feval(fun,x)];
    end
  end
n=1:nit;
plot(n, xdif, '-*');
title(['Plot of error with respect to iteration, f(x)=',char(fun)]);
xc = get(gca,'XLim');
yc = get(gca,'YLim');
xc = (xc(1)+xc(2))/2;
text(xc,yc(2)*0.9,['x_{zero} = ',num2str(x),'; x_{start} =' , ...
          num2str(x0)], 'HorizontalAlignment', 'center'); 

Secant Method algorithm in Matlab

Secant Method algorithm in Matlab :

function [xvect,xdif,fx,nit] = secant(x1,x0,nmax,fun,toll);
% SECANT   Do the secant iteration to find the zeros of the given
% inline scalar function and its derivative. 
%    [XVEC,XDIF,FX,NIT] = SECANT(X1,X0,NMAX,FUN,TOLL)  
%       Input:x0 and x1 starting value
%             nmax: maximum number of iteration
%             toll: tolerance, default is 1e-10
%             fun: given inline function
%
%      Output: xvect: stores values in all iterations (arg)
%              xdif : difference between two successive values
%              (arg) 
%              fx: stores function values  in all iterations 
%              nit: number of iteration required
%
%     Examples:
%          fun=inline('x^3+x^2-4'); x0=0.3;x1=0.5; nmax=100;
%          fun=inline('x^5+10*x^2-9*x+10');

% Author: Bishnu Lamichhane, University of Stuttgart

if (nargin==4)  toll=1e-10; end
x=x1;
fx1=feval(fun,x);
xvect=[x];
 fx=[fx1];
 x=x0;
 fx0=feval(fun,x);
 xvect=[xvect;x];fx=[fx;fx0];err=toll+1;nit=0;xdif=[];
 while(nit<nmax & err>toll)
   nit=nit+1;
   if (abs(fx0-fx1)<eps)
     err=toll*1e-10;
     disp('Stop for vanishing dfun');
   else
     x=x0-fx0*(x0-x1)/(fx0-fx1);
     xvect=[xvect;x];
   fnew=feval(fun,x);
   fx=[fx;fnew];
   err=abs(x0-x);
   xdif=[xdif;err];
   x1=x0;fx1=fx0;x0=x;fx0=fnew;
   end;
 end;
n=1:nit;
plot(n, xdif, '-*');
title(['Plot of error with respect to iteration, f(x)=',char(fun)]);
xc = get(gca,'XLim');
yc = get(gca,'YLim');
xc = (xc(1)+xc(2))/2;
text(xc,yc(2)*0.9,['x_{zero} = ',num2str(x)], ...
          'HorizontalAlignment', 'center'); 

Regula-Falsi method algorithm in Matlab

Regula-Falsi method algorithm in Matlab :

function [c,E,fc]=regula(f,a,b,error,n_iter)
% syntaxis: [c,E,fc]=regula(f,a,b,error,n_iter)
% ---------------------------------------------------------------------
% Esta funcion permite determinar de manera aproximada el valor de la
% raiz de una ecuacion no lineal f(x)=0 mediante el metodo de Regula-Falsi.
% De manera adicional se proporciona una tabla con el valor de la 
% aproximacion y el error corespondiente en cada iteracion.
%
% Entrada:
% f      --> Funcion simbolica
% a, b   --> Extremos del intervalo
% error  --> Tolerancia del calculo {default | 0.00001}
% n_iter --> Numero maximo de iteraciones {default | 1000}
%
% Salida:
% c  --> Aproximacion encontrada
% E  --> Error de la aproximacion "c"
% fc --> Valor de la funcion en la aproximacion
%
%
% Ejemplo:
% 
% f=sym('x^2-1'); 
% [c,E,fc]=regula(f,0,3,0.005,20);
%
% --------|----------|------------|
% Iter.       Aprox.     Error.   
% --------|----------|------------|
% ini        0.3333      3.0000       
% 1          0.6000      0.2667       
% 2          0.7778      0.1778       
% 3          0.8824      0.1046       
% 4          0.9394      0.0570       
% 5          0.9692      0.0298       
% 6          0.9845      0.0153       
% 7          0.9922      0.0077       
% 8          0.9961      0.0039       
% --------|----------|------------|
%
% c = 0.9961
%
% E = 0.0039
%
% fc = -0.0078
%
%
% Ver tambien: biseccion
%
% ----------------------------------------------------------------------
%
% Elaborado por:
%
% Msc. Alexeis Comanioni Guerra
% Lic. Vianka Orovio Cobo
%
% Revision: 1.0 
% Fecha: 30/09/2007


%% Validacion de la entrada
if (nargin<3)  
    error('Revise los datos de entrada.');
else
    tipo = class(f);
    if all( (tipo(1:3)=='sym')>0 )~=1
        error('La funcion f debe ser simbolica.');
    end
    
    if sum(size(a))~=2 || sum(size(b))~=2
        error('Los extremos del intervalo deben ser escalares.'); 
    end  
end

if (nargin<4) 
    error=0.00001;
    n_iter=1000;
elseif (nargin<5) 
    n_iter=1000;
    if sum(size(error))~=2
        error('La tolerancia debe ser un escalar.'); 
    end
else
    if sum(size(error))~=2
        error('La tolerancia debe ser un escalar.'); 
    end
    
    if round(n_iter)-n_iter~=0
        error('El numero de iteraciones debe ser un entero.');
    end
end

%% Teorema de Bolzano
if subs(f,a)*subs(f,b)>0
 error('Imposible de aplicar el metodo: note que --> f(a)*f(b)>0');
end

%% Metodo de Regula-Falsi
x= a - (subs(f,a)*(b-a)/(subs(f,b)-subs(f,a))); 
E=b-a;
RES(1,1)=x; ERR(1,1)=E;
xa=x; i=1;
while (E>=error) && (i<n_iter)
    if subs(f,x)==0
        fprintf('La raiz es exactamente: %f',x);
        break;
 elseif subs(f,b)*subs(f,x)<0
     a=x; %b=b;
 else
     b=x; %a=a;
    end
    x= a - (subs(f,a)*(b-a)/(subs(f,b)-subs(f,a)));
    E= abs(x-xa);
    xa=x; i=i+1;
    RES(i,1)=x; ERR(i,1)=E;
end

%% Formateo de las salidas
disp('--------|----------|------------|');
fprintf('Iter.       Aprox.     Error.   \n');
disp('--------|----------|------------|');
fprintf('%-10s %-11.4f %-12.4f \n','ini',RES(1,1),ERR(1,1));
for i=2:max(size(RES))
    fprintf('%-10.0d %-11.4f %-12.4f \n',i-1,RES(i,1),ERR(i,1));
end
disp('--------|----------|------------|');

c=x; fc=subs(f,x);

The Fixed-point iteration algorithm in Matlab

The Fixed-point iteration algorithm in Matlab :

function fixed_point(p0, N)

%Fixed_Point(p0, N) approximates the root of the equation f(x) = 0
%rewritten in the form x = g(x), starting with an initial approximation p0. 
%The iterative technique is implemented N times.
%The user has to enter the function g(x)at the bottom


%Author: Alain G. Kapitho
%Date  : Jan. 2006

i = 1;
p(1) = p0;
tol = 1e-05;
while i <= N
   p(i+1) = g(p(i));
   if abs(p(i+1)-p(i)) < tol  %stopping criterion
      disp('The procedure was successful after k iterations')
      k = i
      disp('The root to the equation is')
      p(i+1)
      return
   end
   i = i+1;
end

if abs(p(i)-p(i-1)) > tol | i > N
   disp('The procedure was unsuccessful')
   disp('Condition |p(i+1)-p(i)| < tol was not sastified')
   tol
   disp('Please, examine the sequence of iterates')
   p = p'
   disp('In case you observe convergence, then increase the maximum number of iterations')
   disp('In case of divergence, try another initial approximation p0 or rewrite g(x)')
   disp('in such a way that |g''(x)|< 1 near the root')
end

%this part has to be changed accordingly with the specific function g(x)
function y = g(x)
%y = x - x.^3 - 4*x.^2 + 10;
%y = sqrt(10./x - 4*x);
y = x - (x.^3 + 4*x.^2 - 10)/(3*x.^2 + 8*x);

Algorithm of Bisection method in Mtalab

Algorithm of Bisection method in Mtalab :

function [c,E,fc]=biseccion(f,a,b,error,n_iter)
% syntaxis: [c,E,fc]=biseccion(f,a,b,E,n_iter)
% ---------------------------------------------------------------------
% Esta funcion permite determinar de manera aproximada el valor de la
% raiz de una ecuacion no lineal f(x)=0 mediante el metodo de Biseccion.
% De manera adicional se proporciona una tabla con el valor de la 
% aproximacion y el error corespondiente en cada iteracion.
%
% Entrada:
% f      --> Funcion simbolica
% a, b   --> Extremos del intervalo
% error  --> Tolerancia del calculo {default | 0.00001}
% n_iter --> Numero maximo de iteraciones {default | 1000}
%
% Salida:
% c  --> Aproximacion encontrada
% E  --> Error de la aproximacion "c"
% fc --> Valor de la funcion en la aproximacion
%
%
% Ejemplo:
% 
% f=sym('x^2-1'); 
% [c,E,fc]=biseccion(f,0,3,0.005,20);
%
% --------|----------|------------|
% Iter.       Aprox.     Error.   
% --------|----------|------------|
% ini        1.5000      1.5000       
% 1          0.7500      0.7500       
% 2          1.1250      0.3750       
% 3          0.9375      0.1875       
% 4          1.0313      0.0938       
% 5          0.9844      0.0469       
% 6          1.0078      0.0234       
% 7          0.9961      0.0117       
% 8          1.0020      0.0059       
% 9          0.9990      0.0029       
% --------|----------|------------|
%
% c = 0.9990
%
% E = 0.0029
%
% fc = -0.0020
%
%
% Ver tambien: regula
%
% ----------------------------------------------------------------------
%
% Elaborado por:
%
% Msc. Alexeis Comanioni Guerra
% Lic. Vianka Orovio Cobo
%
% Revision: 1.0 
% Fecha: 30/09/2007


%% Validacion de la entrada
if (nargin<3)  
    error('Revise los datos de entrada.');
else
    tipo = class(f);
    if all( (tipo(1:3)=='sym')>0 )~=1
        error('La funcion f debe ser simbolica.');
    end

    if sum(size(a))~=2 || sum(size(b))~=2
        error('Los extremos del intervalo deben ser escalares.'); 
    end  
end

if (nargin<4) 
    error=0.00001;
    n_iter=1000;
elseif (nargin<5) 
    n_iter=1000;
    if sum(size(error))~=2
        error('La tolerancia debe ser un escalar.'); 
    end
else
    if sum(size(error))~=2
        error('La tolerancia debe ser un escalar.'); 
    end
    
    if round(n_iter)-n_iter~=0
        error('El numero de iteraciones debe ser un entero.');
    end
end

%% Teorema de Bolzano
if subs(f,a)*subs(f,b)>0
    error('Imposible de aplicar el metodo: note que --> f(a)*f(b)>0');
end

%% Metodo de Biseccion
x=(a+b)/2; E=(b-a)/2;
RES(1,1)=x; ERR(1,1)=E;
i=1;
while (E>=error) && (i<=n_iter)
    if subs(f,x)==0
        fprintf('La raiz es exactamente: %f',x);
        break;
    elseif subs(f,b)*subs(f,x)<0
 a=x; %b=b;
    else
 b=x; %a=a;
    end
    x=(a+b)/2; E=(b-a)/2;
    i=i+1;
    RES(i,1)=x; ERR(i,1)=E;
end

%% Formateo de las salidas
disp('--------|----------|------------|');
fprintf('Iter.      Aprox.      Error.   \n');
disp('--------|----------|------------|');
fprintf('%-10s %-11.4f %-12.4f \n','ini',RES(1,1),ERR(1,1));
for i=2:max(size(RES))
    fprintf('%-10.0d %-11.4f %-12.4f \n',i-1,RES(i,1),ERR(i,1));
end
disp('--------|----------|------------|');

c=x; fc=subs(f,x);

Sunday, May 26, 2013

The LU Decomposition algorithm in Matlab

The LU Decomposition algorithm in Matlab :

function [L,U,P] = lu_dcmp(A)
%This gives LU decomposition of A with the permutation matrix P
% denoting the row switch(exchange) during factorization
NA = size(A,1);
AP = [A eye(NA)]; %augment with the permutation matrix.
for k = 1:NA - 1
%Partial Pivoting at AP(k,k)
[akx, kx] = max(abs(AP(k:NA,k)));
if akx < eps
error(’Singular matrix and No LU decomposition’)
end
mx = k+kx-1;
if kx > 1 % Row change if necessary
tmp_row = AP(k,:);
AP(k,:) = AP(mx,:);
AP(mx,:) = tmp_row;
end
% LU decomposition
for m = k + 1: NA
AP(m,k) = AP(m,k)/AP(k,k); %Eq.(2.4.8.2)
AP(m,k+1:NA) = AP(m,k + 1:NA)-AP(m,k)*AP(k,k + 1:NA); %Eq.(2.4.9)
end
end
P = AP(1:NA, NA + 1:NA + NA); %Permutation matrix
for m = 1:NA
for n = 1:NA
if m == n, L(m,m) = 1.; U(m,m) = AP(m,m);
elseif m > n, L(m,n) = AP(m,n); U(m,n) = 0.;
else L(m,n) = 0.; U(m,n) = AP(m,n);
end
end
end
if nargout == 0, disp(’L*U = P*A with’); L,U,P, end
%You can check if P’*L*U = A?

Gauss Elimination Algorithm to solve Ax = B in Matlab

Gauss Elimination Algorithm to solve Ax = B in Matlab :

function x = gauss(A,B)
%The sizes of matrices A,B are supposed to be NA x NA and NA x NB.
%This function solves Ax = B by Gauss elimination algorithm.
NA = size(A,2); [NB1,NB] = size(B);
if NB1 ~= NA, error(’A and B must have compatible dimensions’); end
N = NA + NB; AB = [A(1:NA,1:NA) B(1:NA,1:NB)]; % Augmented matrix
epss = eps*ones(NA,1);
for k = 1:NA
%Scaled Partial Pivoting at AB(k,k) by Eq.(2.2.20)
[akx,kx] = max(abs(AB(k:NA,k))./ ...
max(abs([AB(k:NA,k + 1:NA) epss(1:NA - k + 1)]’))’);
if akx < eps, error(’Singular matrix and No unique solution’); end
mx = k + kx - 1;
if kx > 1 % Row change if necessary
tmp_row = AB(k,k:N);
AB(k,k:N) = AB(mx,k:N);
AB(mx,k:N) = tmp_row;
end
% Gauss forward elimination
AB(k,k + 1:N) = AB(k,k+1:N)/AB(k,k);
AB(k,k) = 1; %make each diagonal element one
for m = k + 1: NA
AB(m,k+1:N) = AB(m,k+1:N) - AB(m,k)*AB(k,k+1:N); %Eq.(2.2.5)
AB(m,k) = 0;
end
end
%backward substitution for a upper-triangular matrix eqation
% having all the diagonal elements equal to one
x(NA,:) = AB(NA,NA+1:N);
for m = NA-1: -1:1
x(m,:) = AB(m,NA + 1:N)-AB(m,m + 1:NA)*x(m + 1:NA,:); %Eq.(2.2.7)
end