College of Engineering
Department of Mechanical & Industrial Engineering
Numerical Methods for Engineers (MEIE4183)
Course Project (Fall)
Submission NOTE:
1- Write the problem and code script, its final results/solution/iteration table on the Word
file for each problem as per the provided format.
2- Prepare a program for each problem as a separate Matlab script file.
3- Zip all the program files, Excel, and this Word file, and upload them to Moodle as a
onezipped folder.
Problem-1:
Write the given data in an Excel file and save it as Read_Excel.xlsx. Write a code to read the
data from the Excel file, and store it in matrix variables. Store the header in one array variable
and data in another array variable. Rewrite the stored array variables in another Excel file with
the name Write_Excel.xlsx, and in the same format as in the read file.
Xi Yi XiYi Xi2 Yi-Ymean Predicted Y Yi-Predicted Y
1 0.5 0.5 1 8.5762 0.9107 0.1686
2 2.5 5 4 0.8622 1.75 0.5625
3 2 6 9 2.0408 2.5892 0.3472
4 4 16 16 0.3265 3.4285 0.3265
5 3.5 17.5 25 0.0051 4.2678 0.5896
6 6 36 36 6.6128 5.1071 0.7971
7 5.5 38.5 49 4.2908 5.9464 0.1992
Solution 1:
clear
clc
% Read data from the Excel file
filename = 'Read_Excel.xlsx'; % this is the input file
[num, txt, raw] = xlsread(filename);% identifying the types of input
% Store the header in one variable and data in another variable
headerArray = txt(1, :); % First row contains headers
dataArray = num; % Numeric data
% Write the stored arrays to Write_Excel.xlsx
% Prepare data to write, combining header and data
writeData = [headerArray; num2cell(dataArray)];
% Step 4: Write to the new Excel file
write_filename = 'Write_Excel.xlsx';
xlswrite(write_filename, writeData);
Note: The program will create an excel file and in case a file with the same
name exists it will overwrite the old file, and we can check by noticing the
time of save of the write file.
Snapshots:
Problem-2:
Write a Matlab program to find the root(s) of the following function:
𝑓(𝑥)=4𝑥3−11𝑥2+15𝑥−5 Using:
(a) Bisection method (filename: bisection.m)
(b) False position method (filename: FPM.m)
(c) Fixed point iteration method (filename: FPIM.m)
(d) Secant method (filename: SM.m)
NOTE: take a suitable initial guess (s) wherever necessary.
Code:
1. Bisection method
% Define the function
f = @(x) 4*x^3 - 11*x^2 + 15*x - 5;
% Initial interval [xl, xu]
xl = 0; % Lower bound
xu = 2; % Upper bound
% Check if the initial guesses bracket a root
if f(xl) * f(xu) > 0
error('The function does not change sign over the interval. Please
choose different initial guesses.');
end
% Tolerance for stopping criterion
tolerance = 1e-6;
% Maximum number of iterations to prevent infinite loop
max_iterations = 1000;
% Bisection method loop
for iteration = 1:max_iterations
% Calculate midpoint
xm = (xl + xu) / 2;
% Evaluate function at midpoint
fm = f(xm);
% Check if the midpoint is a root or if the interval is within the
desired tolerance
if abs(fm) < tolerance || (xu - xl) / 2 < tolerance
fprintf('The approximate root is: %.6f\n', xm);
fprintf('Number of iterations: %d\n', iteration);
return;
end
% Determine the subinterval that contains the root
if f(xl) * fm < 0
xu = xm; % Root is in the lower subinterval
else
xl = xm; % Root is in the upper subinterval
end
end
% If the loop completes without finding a root
error('The method did not converge after %d iterations.', max_iterations);
2. False position method
% Define the function
f = @(x) 4*x^3 - 11*x^2 + 15*x - 5;
% Initial guesses
xl = 0; % Lower bound
xu = 2; % Upper bound
% Check if the initial interval brackets a root
if f(xl) * f(xu) > 0
error('The function does not change sign over the interval. Please
choose different initial guesses.');
end
% Tolerance for stopping criterion
tolerance = 1e-6;
% Maximum number of iterations to prevent infinite loop
max_iterations = 1000;
% False Position Method loop
for iteration = 1:max_iterations
% Calculate the point where the secant line intersects the x-axis
xm = xu - (f(xu) * (xl - xu)) / (f(xl) - f(xu));
% Evaluate function at xm
fm = f(xm);
% Check if the approximation is within the desired tolerance
if abs(fm) < tolerance
fprintf('The approximate root is: %.6f\n', xm);
fprintf('Number of iterations: %d\n', iteration);
return;
end
% Update the bounds
if f(xl) * fm < 0
xu = xm; % Root is in the lower subinterval
else
xl = xm; % Root is in the upper subinterval
end
end
% If the loop completes without finding a root
error('The method did not converge after %d iterations.', max_iterations);
3. fixed-point method
% Define the function f(x)
f = @(x) 4*x^3 - 11*x^2 + 15*x - 5;
% Rearrange f(x) = 0 to x = g(x)
g = @(x) sqrt((15*x - 5) / 4);
% Initial guess
x0 = 1.5;
% Tolerance for stopping criterion
tolerance = 1e-6;
% Maximum number of iterations to prevent infinite loop
max_iterations = 1000;
% Initialize variables
x_prev = x0; % Previous value of x
iteration = 0; % Counter for iterations
% Fixed-Point Iteration loop
while iteration < max_iterations
% Compute the next approximation
x_next = g(x_prev);
% Check if the current approximation is within the desired tolerance
if abs(x_next - x_prev) < tolerance
fprintf('The approximate root is: %.6f\n', x_next);
fprintf('Number of iterations: %d\n', iteration);
break;
end
% Update for the next iteration
x_prev = x_next;
iteration = iteration + 1;
end
% If the loop completes without convergence
if iteration == max_iterations
fprintf('The method did not converge after %d iterations.\n',
max_iterations);
end
4.Secant method
% Define the function f(x)
f = @(x) 4*x.^3 - 11*x.^2 + 15*x - 5;
% Set the tolerance for stopping criterion
tolerance = 1e-6;
% Initialize error to a large value
error = 100;
% Initial guesses (x0 and x1)
x0 = 0; % First initial guess
x1 = 1; % Second initial guess
% Secant method loop
while error > tolerance
% Compute the next approximation using the Secant formula
x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0));
% Calculate the relative error as a percentage
error = abs((x2 - x1) / x2) * 100;
% Update x0 and x1 for the next iteration
x0 = x1;
x1 = x2;
end
% Display the root found using the Secant method
fprintf('The root using Secant method is: %.6f\n', x2);
Execution Snapshots /snapshots of result:
1. Bisection method
2. False position method
3. Fixed-point method
4.Secant method
Problem-3:
Write a Matlab program to solve the following set of equations
x + y - 2z = -3
6x + 2y + 2z =3
-3x + 5y + z = 1
Using:
(a) Gauss-elimination with partial pivoting (filename: GEWPP.m)
(b) Lower upper decomposition method (filename: LUDM.m)
(c) Gauss-Seidel Method (filename: GSM.m)
Code:
1. Gauss-elimination
A = [1, 1, -2; 6, 2, 2; -3, 5, 1];
b = [-3; 3; 1];
Ab = [A b];
n = length(b);
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
f21 = Ab(2, 1) / Ab(1, 1);
Ab(2, :) = Ab(2, :) - Ab(1, :) * f21;
f31 = Ab(3, 1) / Ab(1, 1);
Ab(3, :) = Ab(3, :) - Ab(1, :) * f31;
f32 = Ab(3, 2) / Ab(2, 2);
Ab(3, :) = Ab(3, :) - Ab(2, :) * f32;
a(3) = Ab(3, 4) / Ab(3, 3);
a(2) = (Ab(2, 4) - Ab(2, 3) * a(3)) / Ab(2, 2);
a(1) = (Ab(1, 4) - Ab(1, 3) * a(3) - Ab(1, 2) * a(2)) / Ab(1, 1);
transpose(a)
2. Lower upper
% L-U Decomposition Example 10.2
clear;
clc;
% Define matrix A
A = [1, 1, -2; 6, 2, 2; -3, 5, 1];
% Define matrix b
b = [-3; 3; 1];
% Initialize U with A
U = A;
% Compute L-U decomposition
f21 = U(2,1) / U(1,1);
U(2,:) = U(2,:) - f21 * U(1,:);
f31 = U(3,1) / U(1,1);
U(3,:) = U(3,:) - f31 * U(1,:);
f32 = U(3,2) / U(2,2);
U(3,:) = U(3,:) - f32 * U(2,:);
% Display the Lower (L) and Upper (U) matrices
L = eye(3); % Identity matrix
L(2,1) = f21;
L(3,1) = f31;
L(3,2) = f32;
disp('Lower Triangular Matrix (L):');
disp(L);
disp('Upper Triangular Matrix (U):');
disp(U);
% Solve for x using L-U decomposition
y = L \ b;
x = U \ y;
disp('Solution (x):');
disp(x);
3. Gauss seidel
% Clear workspace and command window
clear;
clc;
% Define the coefficient matrix A and the constant vector B
A = [1, 1, -2; 6, 2, 2; -3, 5, 1];
B = [-3; 3; 1];
% Define the initial guess vector P
P = [1; 0; 1]; % Initial guess for the solution
% Get the size of the system (number of equations)
N = length(B);
% Initialize the solution vector X with zeros
X = zeros(N, 1);
% Perform iterations for a maximum of 100 steps
for j = 1:100
% Iterate over each equation
for i = 1:N
% Compute the i-th component of the solution vector X
X(i) = (B(i) / A(i, i)) - ...
(A(i, [1:i-1, i+1:N]) * P([1:i-1, i+1:N])) / A(i, i);
end
% Display the current iteration number and solution vector
fprintf('Iteration no %d\n', j);
disp(X);
% Update the guess vector for the next iteration
P = X;
end
Execution Snapshots /snapshots of result:
1. Gauss-elimination
2. Lower upper
3. Gauss seidel
Problem-4:
Write a program to find the inverse of the coefficient matrix using lower-upper decomposition
method for the following set of equations: (Matlab filename: LUDINV.m)
10x + 2y - z = 29
-3x - 6y + 2z = -61
x + y + 5z = -21
Note: Matlab has a direct function to find the inverse of a matrix, you should not use that rather
write a program as per the method taught in Chapter 10 using the concept of lower-upper
decomposition.
Code:
clear
clc
A=[10 2 -1; -3 -6 2; 1 1 5]
b=[29; -61; -21]
f21=A(2,1)/A(1,1);
A(2,:)=A(2,:)-A(1,:)*f21;
f31=A(3,1)/A(1,1);
A(3,:)=A(3,:)-A(1,:)*f31;
f32=A(3,2)/A(2,2);
A(3,:)=A(3,:)-A(2,:)*f32;
U=A;
L=eye(3);
L(2,1)=f21;
L(3,1)=f31;
L(3,2)=f32;
I=eye(3);
d1=L\I(:,1);
d2=L\I(:,2);
d3=L\I(:,3);
x(:,1)=U\d1;
x(:,2)=U\d2;
x(:,3)=U\d3;
x
snapshoots:
Results:
Problem-5:
Write a Matlab program to solve the following multiple regression problem to fit the given data
using multiple linear regression. (Matlab filename: MLRM.m)
X1 X2 Y
0.1 0.1 5
2.1 1 10
2.5 2 9
1.1 3.5 0.1
4 6.3 3
7 2.1 27
Code:
clc
clear
x1=[0.1 2.1 2.5 1.1 4 7];
x2=[0.1 1 2 3.5 6.3 2.1];
y=[5 10 9 0.1 3 27];
n=length(x1);
sum(y);
sum(x1);
sum(x2);
sum(x1.^2);
sum(x2.^2);
sum(x1.*x2);
sum(x1.*y);
sum(x2.*y);
A=[n sum(x1) sum(x2); sum(x1) sum(x1.^2) sum(x1.*x2); sum(x2)
sum(x1.*x2) sum(x2.^2)]
B=[sum(y); sum(x1.*y); sum(x2.*y)]
coeff_A= A\B
snapshoots:
Results:
Problem-6:
Write a Matlab program to solve the following problem:
Determine the resistance at 33oC using the Lagrange interpolation method of third-order
polynomial. (Matlab filename: LIM3ORD.m)
R (ohm) T (oC)
1111 25.11
910 30.13
636 40.12
455 50.13
Code:
clc
clear
T=[25.11 30.13 40.12 50.13];
R=[1111 910 636 455];
Tgiven=33;
n=length(T);
sum=0;
for i=1:n
a=1;
for j=1:n
if j~=i
a=a*(Tgiven-T(j))/(T(i)-T(j));
end
end
sum=sum+a*R(i);
end
fprintf('The resistance at temperature 35 is:%.2f\n',sum);
Snapshots of results:
Problem-7:
Write a Matlab program to evaluate the integral of the following function using Simpson’s 3/8
rule: (Matlab filename: SIMP1BY3.m)
Code:
clc
clear
f=@(x) 1-x-4*x.^3+2*x.^5;
a=-2; b=4;
n=3;
h=(b-a)/n;
x=[a,a+h,a+2*h,b];
y=f(x);
I=(3*h/8)*(y(1)+3*y(2)+3*y(3)+y(4));
fprintf('The integration by using Simpsons 3/8 rule is %.0f',I);
Snapshots of results:
