Script for Gauss Elimination method
This matlab script can solve a system of linear equations by Gauss elimination method with partial pivoting. Fun fact is, this script can show the calculation steps. So, A lot less effort is needed while preparing assignments.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Title: Gauss Elimination with partial pivoting %
% Author: Dhiman Roy, Dept. of ME (BUET) %
% Date: August 23, 2021 %
% Licensed under Creative Commons %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
global showSteps;
showSteps = true; %can be true or false. true shows calculation steps
A = input('Please Enter the Co-efficient Matrix, A:\n' );
b = input('Please Enter the Constants(Must be a column vector, b:\n');
solByGEPP = gaussElimPartPiv(A, b)
function [x] = gaussElimPartPiv(A, b)
global showSteps;
[m, n] = size(A);
if m~=n
error("'A' must be a square matrix.");
else
aug = [A b];
if showSteps == 1
fprintf("Gauss elimination with partial pivoting.\nAugmented Matrix:\n");
disp(aug);
fprintf("\n");
end
for k = 1: n-1
%partial pivoting
for p = k+1:n
if abs(A(k,k) < A(p,k))
aug([k p],:) = aug([p k],:);
end
end
%forward elimination
for i = k+1: n
if showSteps == 1
fprintf("Row%d = Row%d - (%0.4f)/(%0.4f) * Row%d\n", i, i, aug(i,k), aug(k,k), k);
end
factor = aug(i,k)/aug(k,k);
aug(i,k:n+1) = aug(i,k:n+1)-factor*aug(k,k:n+1);
end
%display augmented matrix after each step
if showSteps == 1
fprintf("\n");
disp(aug)
fprintf("\n");
end
end
%back substitution
x = zeros(n,1);
x(n) = aug(n,n+1)/aug(n,n);
for i = n-1:-1:1
x(i) = (aug(i,n+1)-aug(i,i+1:n)*x(i+1:n))/aug(i,i);
end
end
end
Script for Gauss Jordan method
This matlab script can solve a system of linear equations by Gauss Jordan method with partial pivoting. Co-efficeint matrix and contant vector are be prompted for user input. Fun fact is, this script can show the calculation steps.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Title: Gauss Jordan with partial pivoting %
% Author: Dhiman Roy, Dept. of ME (BUET) %
% Date: August 23, 2021 %
% Licensed under Creative Commons %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
global showSteps;
showSteps = true; %can be true or false. true shows calculation steps
A = input('Please Enter the Co-efficient Matrix, A:\n' );
b = input('Please Enter the Constants(Must be a column vector, b:\n');
solByGJPP = gaussJordanPartPiv(A, b)
function [x] = gaussJordanPartPiv(A, b)
global showSteps;
[m, n] = size(A);
if m~=n
error("'A' must be a square matrix.");
else
aug = [A b];
if showSteps == 1
fprintf("Gauss Jordan with partial pivoting.\nAugmented Matrix:\n");
disp(aug);
fprintf("\n");
end
for k = 1: n-1
%partial pivoting
for p = k+1:n
if abs(A(k,k) < A(p,k))
aug([k p],:) = aug([p k],:);
end
end
%forward elimination
for i = k+1: n
if showSteps == 1
fprintf("Row%d = Row%d - (%0.4f)/(%0.4f) * Row%d\n", i, i, aug(i,k), aug(k,k), k);
end
factor = aug(i,k)/aug(k,k);
aug(i,k:n+1) = aug(i,k:n+1) - factor*aug(k,k:n+1);
end
%display augmented matrix after each step
if showSteps == 1
fprintf("\n");
disp(aug);
end
end
%elimination from upper triangular matrix
if showSteps == 1
fprintf("Elimination from upper triangular matrix:\n\n");
end
for k = n:-1:2
for i = k-1:-1:1
if showSteps == 1
fprintf("Row%d = Row%d - (%0.4f)/(%0.4f) * Row%d\n", i, i, aug(i,k), aug(k,k), k);
end
factor = aug(i,k)/aug(k,k);
aug(i,n+1:-1:i) = aug(i,n+1:-1:i) - factor*aug(k,n+1:-1:i);
end
if showSteps == 1
fprintf("\n");
disp(aug)
end
end
%converting to Identity Matrix
for i = 1:n
aug(i,:) = aug(i,:)/aug(i,i);
end
if showSteps == 1
fprintf("\nFinal form of the augmented matrix:\n\n");
disp(aug);
fprintf("\n");
end
end
x = aug(:,n+1);
end
Script for Gauss Seidel method
This matlab script can solve a system of linear equations by Gauss Seidel method, an iterative method. Co-efficeint matrix, contant vector, initial guesses and tolerence are be prompted for user input.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Title: Gauss Seidel method %
% Author: Dhiman Roy, Dept. of ME (BUET) %
% Date: August 23, 2021 %
% Licensed under Creative Commons %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
global showSteps;
showSteps = false; %true shows calculation steps
A = input('Please Enter the Co-efficient Matrix, A:\n' );
b = input('Please Enter the Constants(Must be a column vector, b:\n');
solByGaussSiedel = gaussSeidel(A, b)
function [x] = gaussSeidel(A, b)
global showSteps;
xIni = input("Enter initial guesses (Must be a column vector:\n");
eLimit = input("Enter tolerance:\n");
[m,n] = size(A);
for k = 1: n-1
for p = k+1:n
if abs(A(k,k) < A(p,k))
b([k p]) = b([p k]);
A([k p],:) = A([p k],:);
end
end
end
check = ones(n,1);
for i = 1:n
if 2*abs(A(i,i)) < sum(abs(A(i,:)))
check(i) = 0;
end
end
if all(check) ~= 1
error("Partial Pivoting does not make the coefficient matrix a diagonally dominant matrix.")
else
disp("Partial Pivoting makes the coefficient matrix a diagonally dominant matrix.")
end
x = xIni;
xOld = x;
for k = 1:100
err = 0;
for i = 1:n
x(i) = (b(i)-A(i,[1:i-1,i+1:n])*x([1:i-1,i+1:n]))/A(i,i);
err = err + abs((x(i)-xOld(i))/x(i));
xOld(i) = x(i);
end
if showSteps == 1
fprintf("Iteration %d\n",k);
disp(x);
disp(err);
end
if err < eLimit
break
end
end
end
Script for Jacobi method
Another iterative method for solving a system of linear equations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Title: Jacobi method %
% Author: Dhiman Roy, Dept. of ME (BUET) %
% Date: August 23, 2021 %
% Licensed under Creative Commons %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
global showSteps;
showSteps = false; %true shows calculation steps
A = input('Please Enter the Co-efficient Matrix, A:\n' );
b = input('Please Enter the Constants(Must be a column vector, b:\n');
solByJacobi = jacobi(A, b)
function [x] = jacobi(A, b)
global showSteps;
xIni = input("Enter initial guesses (Must be a column vector:\n");
eLimit = input("Enter tolerance:\n");
[m,n] = size(A);
for k = 1: n-1
for p = k+1:n
if abs(A(k,k) < A(p,k))
b([k p]) = b([p k]);
A([k p],:) = A([p k],:);
end
end
end
check = ones(n,1);
for i = 1:n
if 2*abs(A(i,i)) < sum(abs(A(i,:)))
check(i) = 0;
end
end
if all(check) ~= 1
error("Partial pivoting does not make the coefficient matrix a diagonally dominant matrix.")
else
disp("Partial pivoting makes the coefficient matrix a diagonally dominant matrix.")
end
xOld = xIni;
xNew = zeros(n,1);
for k = 1:100
err = 0;
for i = 1:n
xNew(i) = b(i)/A(i,i)-A(i,[1:i-1,i+1:n])*(xOld([1:i-1,i+1:n])/A(i,i));
err = err + abs((xNew(i)-xOld(i))/xNew(i));
end
if showSteps == 1
fprintf("Iteration %d\n",k);
disp(xNew);
abs((xNew-xOld)/xNew);
end
if err < eLimit
break;
end
xOld = xNew;
end
x = xNew;
end
Script for SOR method
Solves a system of linear equations by successive over relaxation, prompts for co-efficient matrix, contant vector, initial guesses and relaxation factor.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Title: SOR method %
% Author: Dhiman Roy, Dept. of ME (BUET) %
% Date: August 23, 2021 %
% Licensed under Creative Commons %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
global showSteps;
showSteps = false; %true shows calculation steps
A = input('Please Enter the Co-efficient Matrix, A:\n' );
b = input('Please Enter the Constants(Must be a column vector, b:\n');
solBySOR = sor(A, b)
function [x] = sor(A, b)
global showSteps;
xIni = input("Enter initial guesses (Must be a column vector:\n");
eLimit = input("Enter tolerance:\n");
w = input("Enter Relaxation Factor:\n");
[m,n] = size(A);
for k = 1: n-1
for p = k+1:n
if abs(A(k,k) < A(p,k))
b([k p]) = b([p k]);
A([k p],:) = A([p k],:);
end
end
end
check = ones(n,1);
for i = 1:n
if 2*abs(A(i,i)) < sum(abs(A(i,:)))
check(i) = 0;
end
end
if all(check) ~= 1
error("Partial Pivoting does not make the coefficient matrix a diagonally dominant matrix.")
else
disp("Partial Pivoting makes the coefficient matrix a diagonally dominant matrix.")
end
x = xIni;
xOld = x;
for k = 1:100
err = 0;
for i = 1:n
x(i) = (1-w)*x(i) + w*(b(i)-A(i,[1:i-1,i+1:n])*x([1:i-1,i+1:n]))/A(i,i);
abs((x(i)-xOld(i))/x(i))
err = err + abs((x(i)-xOld(i))/x(i));
xOld(i) = x(i);
end
if showSteps == 1
fprintf("Iteration %d\n",k);
disp(x);
disp(err);
end
if err < eLimit
break
end
end
end
Scheme code for extacting xy plot data
;==========================================
; Title: Data Extract for XY plots
; Author: Dhiman Roy, ME-14, BUET
; Date: February 08, 2021
; Licensed under Creative Commons.
;==========================================
(newline)
(display "Generating VERTICAL LINES...")(newline)
(newline)
;(define avg "Area-Weighted Average")
(define ramp_x-cor (list 0.060 0.070))
(define ramp_y_1-cor (list -0.018 -0.018))
(define ramp_y_2-cor (list -0.033 -0.033))
(define lip_x-cor (list 0.060 0.070))
(define lip_y_1-cor (list -0.018 -0.018))
(define lip_y_2-cor (list -0.033 -0.033))
(
Do ((i 1 (+ i 1))) ((> i (length ramp_x-cor)))
(Ti-menu-load-string (format #f "surface/line-surface line_ramp_~a ~a ~a ~a ~a" i (list-ref ramp_x-cor (- i 1)) (list-ref ramp_y_1-cor (- i 1)) (list-ref ramp_x-cor (- i 1)) (list-ref ramp_y_2-cor (- i 1))))
(newline)
(Ti-menu-load-string (format #f "plot/plot yes ramp_x-vel_~a.dat no yes 0 1 no no x-velocity line_ramp_~a ()" (* 10000 (list-ref ramp_x-cor (- i 1))) i))
(Ti-menu-load-string (format #f "plot/plot yes ramp_temp_~a.dat no yes 0 1 no no temperature line_ramp_~a ()" (* 10000 (list-ref ramp_x-cor (- i 1))) i))
)
(newline)
(newline)
(
Do ((i 1 (+ i 1))) ((> i (length lip_x-cor)))
(Ti-menu-load-string (format #f "surface/line-surface line_lip_~a ~a ~a ~a ~a" i (list-ref lip_x-cor (- i 1)) (list-ref lip_y_1-cor (- i 1)) (list-ref lip_x-cor (- i 1)) (list-ref lip_y_2-cor (- i 1))))
(newline)
(Ti-menu-load-string (format #f "plot/plot yes lip_x-vel_~a.dat no yes 0 1 no no x-velocity line_lip_~a ()" (* 10000 (list-ref lip_x-cor (- i 1))) i))
)
(newline)
(newline)
(
Do ((i 1 (+ i 1))) ((> i 51))
(display "*")
)
(newline)
(display "Vertical lines are generated.")(newline)
;(display "Monitors are created..")(newline)
(display "Files are saved in current working directory...")(newline)
(
Do ((i 1 (+ i 1))) ((> i 51))
(display "*")
)
(newline)
(newline)