A MATLAB Algorithm for a Pentadiagonal Matrix

Problem Statement: Develop a computer algorithm for the pentadiagonal matrix with band of b and solve following equations as well.

MATLAB Program to Solve a Pentadiagonal Matrix with a band b:

% A function to solve equations which coefficients have form of penta

% diagonal matrix with band b

 

function x=penta_diagonal(A,b)

 

[M,N]=size(A);

% Dimension checking

 

if M~=N

    

    error('Matrix must be square');

    

end

 

if length(b)~=M

    

    error('Matrix and vector must have the same number of rows');

    

end

 

x=zeros(N,1);

    

 

if A==A'    % Matrix symmetry checking

    

    % Band Extractions

    

    d=diag(A);

    f=diag(A,1);

    e=diag(A,2);

    

    disp(d);

    disp(f);

    disp(e);

        

    alpha=zeros(N,1);

    gamma=zeros(N-1,1);

    delta=zeros(N-2,1);

    c=zeros(N,1);

    z=zeros(N,1);

    

    % Factor A=LDL'

    

    alpha(1)=d(1);

    gamma(1)=f(1)/alpha(1);

    delta(1)=e(1)/alpha(1);

    

    alpha(2)=d(2)-f(1)*gamma(1);

    gamma(2)=(f(2)-e(1)*gamma(1))/alpha(2);

    delta(2)=e(2)/alpha(2);

    

    for k=3:N-2

        

        alpha(k)=d(k)-e(k-2)*delta(k-2)-alpha(k-1)*gamma(k-1)^2;

        gamma(k)=(f(k)-e(k-1)*gamma(k-1))/alpha(k);

        delta(k)=e(k)/alpha(k);

    end

    

    alpha(N-1)=d(N-1)-e(N-3)*delta(N-3)-alpha(N-2)*gamma(N-2)^2;

    gamma(N-1)=(f(N-1)-e(N-2)*gamma(N-2))/alpha(N-1);

    alpha(N)=d(N)-e(N-2)*delta(N-2)-alpha(N-1)*gamma(N-1)^2;

    

    % Updating Lx=b, Dc=z

    

    z(1)=b(1);

    z(2)=b(2)-gamma(1)*z(1);

    

    for k=3:N

        z(k)=b(k)-gamma(k-1)*z(k-1)-delta(k-2)*z(k-2);

    end

    

    c=z./alpha;

        

    % Back Substitution L'x=c

    

    x(N)=c(N);

    x(N-1)=c(N-1)-gamma(N-1)*x(N);

    

    for k=N-2:-1:1

        

        x(k)=c(k)-gamma(k)*x(k+1)-delta(k)*x(k+2);

    end  

   

end

 

x;

end

Program Outputs:

ans =

   0.673740053050398

  -1.021220159151194

   1.389920424403183

  -1.148541114058356

   2.055702917771884

  -2.018567639257295

Application of Gaussian Elimination to Solve System of Linear Algebraic Equations

Problem: In this example, we are going to solve the following system of linear algebraic equations by the widely-used Gaussian elimination method. Here, we are going to develop a MATLAB program that implements the Gaussian elimination to solve the following linear algebraic equations.

The following program implements Gaussian elimination method with partial pivoting and scaling to solve system of linear algebraic equations. The explanations of the codes are mentioned just right of each line.

MATLAB Program for Gauss Elimination Method

% Gaussian elimination with partial pivoting and scaling

 

function x = gausselimination(a,b)

 

% a - (nxn) matrix

% b - column vector of length n

 

format long;

 

m=size(a,1); % get number of rows in matrix a

 

n=length(b); % get length of b

 

if (m ~= n)

    

    error('a and b do not have the same number of rows')

    

end

 

a(:,n+1)=b;  % Forming (n,n+1) augmented matrix

 

for c=1:n

              

for r=c:m

     

% Scaling of the input matrix "a"

       

    Z(r,c) = a(r,c)/max(a(r,c:n));  % Scaling prior to pivoting 

 

    K(r:r,c:c)=Z(r,c); % Normalized coefficient parameters saved temporality 

                       % to another matrix "K"     

end

 

disp(K);

           

[~,f]=max(abs(K(1:r,c:c))); % Finding maximum absolute value from "K" matrix 

                            % for switching rows in matrix "a"      

disp(f);  

 

a([c,f],1:n+1)=a([f,c],1:n+1); % Switching between two rows 

                                                                                    

disp(a);

                                        

for i=c

% Making diagonal element of matrix "a" into 1.0

    

    a(i,i:n+1) = a(i,i:n+1)/a(i,i);   

 

for j=i+1:n

    

    % Making all elements below the diagonal element into zero

 

    a(j,i:n+1) = a(j,i:n+1) - a(j,i)*a(i,i:n+1);

    

end

 

end

 

disp(a);

 

end

 

% Process of back substitution

 

for j=n-1:-1:1

 

    a(j,n+1) = a(j,n+1) - a(j,j+1:n)*a(j+1:n,n+1);

 

end

    

% Returning solution

 

x=a(:,n+1);

 

end

Program Outputs:

ans =

 

  -1.761817043997860

   0.896228033874012

   4.051931404116157

  -1.617130802539541

   2.041913538501913

   0.151832487155935

Therefore, the solutions of the six linear algebraic equations are -1.76; 0.89; 4.05; -1.62; 2.04 and 0.15.