Numerical
Analysis
Algorithms in
MATLAB
Sat, Dec 9
Numerical Analysis Algorithms in MATLAB
Done by
ﺿﯿﺎء أﺣﻤﺪ ﻋﺒﺪة ﻟﻄﻒ
Eng. Supervised
ﻣﺤﻤﺪ ﻣﻘﻨﻢ
Page 1 of 15
CONTENTS
1. PREFACE
a. Introduction
b. MATLAB
c. MATLAB as a Numerical Analysis Tool
2. Basic Mathematical Operations
3. Matrices Operations
4. Solving Linear Equations
5. Taylor Series
6. Maclaurin Series
7. Solving Linear Equations using Cramer’s Method
8. Solving Linear Equations using Gauss Method
9. Finding the Inverse Matrix using Gauss-Jordan Method
10. Non-linear Equations
11. Bisection Method
12. Fixed Point Method
13. Fixed Point Method for a System of Equations
14. Finding Iterations Using Jacobian Method
15. Finding Iterations Using Gauss-Seidel Method
Page 2 of 15
PREFACE
Introduction
Numerical analysis is a branch of mathematics that solves continuous problems using
numeric approximation. It involves designing methods that give approximate but
accurate numeric solutions, which is useful in cases where the exact solution is
impossible to calculate. For example, numerical analysis can be used to compute the
trajectory of a spacecraft, the optimal price of a product, or the area under a curve.
Numerical analysis is also closely related to computer science, as it relies on algorithms
and software to implement the numerical methods.
MATLAB
MATLAB is a programming and numeric computing platform used to analyze data,
develop algorithms, and create models. MATLAB stands for Matrix Laboratory, and it
allows matrix manipulations, plotting of functions and data, implementation of
algorithms, creation of user interfaces, and interfacing with programs written in other
languages. MATLAB also has many toolboxes and apps that provide specialized
functionality for various applications and domains, such as control systems, deep
learning, image processing, signal processing, and robotics. MATLAB can also run on
the cloud, on embedded devices, and integrate with model-based design.
MATLAB as a Numerical Analysis Tool
MATLAB is widely used for applied numerical analysis in engineering, computational
finance, and computational biology. It provides a range of numerical methods for
solving problems such as interpolation, differentiation, integration, linear systems of
equations, eigenvalues, ordinary and partial differential equations, and more.
Page 3 of 15
Basic Mathematical Operations
• Addition
If x = 2 and y = 3, then x + y returns 5.
• Subtraction
If x = 2 and y = 3, then x - y returns -1.
• Multiplication
If x = 4 and y = 5, then x * y = 20.
• Power
If x = 6, then x^2 = 36.
• Division
If x = 5 and y = 2, then x / y = 2.5.
• Floor Division
If x = 5 and y = 2, then x / y = 2.
• Complex Number Operations
If x = 1 + 2i and y = 2 – 4i, then
x + y = 3 – 2i.
x – y = -1 + 6i.
x * y = 10.
X / y = -0.3 + 0.4i.
Page 4 of 15
Matrices Operations
a) Addition
If A = [1 2; 3 4] and B = [5 6; 7 8], then A + B returns [6 8; 10 12].
b) Subtraction
If A = [1 2; 3 4] and B = [5 6; 7 8], then A - B returns [-4 -4; -4 -4].
c) Multiplication
If A = [1 2; 3 4] and B = [5 6; 7 8], then A * B returns [19 22; 43 50].
d) Element-wise multiplication
If A = [1 2; 3 4] and B = [5 6; 7 8], then A .* B returns [5 12; 21 32].
e) Right array division
If A = [1 2; 3 4] and B = [5 6; 7 8], then A ./ B returns [0.2 0.3333; 0.4286 0.5].
f) Left array division
If A = [1 2; 3 4] and B = [5 6; 7 8], then B ./ A returns [5 3; 2.3333 2].
g) Power
If A = [1 2; 3 4], then A ^ 2 returns [7 10; 15 22].
h) Element-wise power
If A = [1 2; 3 4] and B = [2 3; 4 5], then A .^ B returns [1 8; 81 1024].
i) Transpose
If A = [1 2; 3 4], then A .' returns [1 3; 2 4].
j) Complex conjugate transpose
If A = [1+2i 3-4i; 5+6i 7-8i], then A' returns [1-2i 5-6i; 3+4i 7+8i].
Page 5 of 15
Solving Linear Equations
There are many ways to solve linear equations in MATLAB, some of them are as
follows.
If x + y = 5 and x – y = 1, then
Elimination
Code
A = [1 1; 1 -1];
b = [5; 1];
x = A\b;
x = x(1)
y = x(2)
Output
x=3
y=2
Substitution
Code
syms x y
eqn1 = x + y == 5;
eqn2 = x - y == 1;
sol = solve([eqn1, eqn2], [x, y]);
x = sol.x
y = sol.y
Output
x=3
y =2
Page 6 of 15
Taylor Series
Code
syms x
f = sin(x);
T = taylor(f,x,'Order',8)
T = x - x^3/6 + x^5/120 - x^7/5040
Output
T = - x^7/5040 + x^5/120 - x^3/6 + x
Maclaurin Series
Code
syms x
f = log(1+x);
T = taylor(f,x,'Order',5)
T = x - x^2/2 + x^3/3 - x^4/4
Output
T = - x^4/4 + x^3/3 - x^2/2 + x
Solving Linear Equations using Cramer’s Method
Code
A = [1 2 3; 4 5 6; 7 8 10];
b = [9; 15; 23];
function x = cramer(A,b)
n = size(A,1);
if det(A) == 0
error('The matrix is singular')
end
x = zeros(n,1);
d = det(A);
Page 7 of 15
for i = 1:n
Ai = A;
Ai(:,i) = b;
di = det(Ai);
x(i) = di/d;
end
end
disp('The solution is:')
disp(x)
end
Output
x , y , z =1.5 , 2 , 1.5
Solving Linear Equations using Gauss Method
Code
A = [1 2 3; 4 5 6; 7 8 10];
b = [9; 15; 23];
function x = gauss_elimination(A,b)
n = size(A,1);
if det(A) == 0
error('The matrix is singular')
end
for k = 1:n-1
for i = k+1:n
m = A(i,k)/A(k,k);
A(i,:) = A(i,:) - m*A(k,:);
b(i) = b(i) - m*b(k);
end
end
x = zeros(n,1);
x(n) = b(n)/A(n,n);
Page 8 of 15
for i = n-1:-1:1
s = 0;
for j = i+1:n
s = s + A(i,j)*x(j);
end
x(i) = (b(i) - s)/A(i,i);
end
end
fprintf('The solution is:\n')
disp(x)
end
Output
The solution is:
1.0000
-1.0000
2.0000
Finding the Inverse Matrix using Gauss-Jordan Method
Code
A = [1 2 3; 4 5 6; 7 8 10];
function invA = gauss_jordan(A)
n = size(A,1);
if det(A) == 0
error('The matrix is singular')
end
A = [A eye(n)];
for k = 1:n
A(k,:) = A(k,:)/A(k,k);
for i = [1:k-1 k+1:n]
A(i,:) = A(i,:) - A(i,k)*A(k,:);
end
end
Page 9 of 15
invA = A(:,n+1:end);
end
disp('The output matrix invA is:')
disp(invA)
end
Output
0.1472 −0.1444 0.0639
[−0.0611 0.0222 0.1056 ]
−0.0194 0.1889 −0.1028
Non-linear Equations
Code
syms x
f = x^3 + x - 1;
fun = matlabFunction(f);
x0 = 0;
x = fsolve(fun,x0)
Output
x = 0.6823
Bisection Method
Code
f = @(x) x^2 - 3;
a = 1;
b = 2;
eps = 1e-4;
function [root, iter] = bisection(f, a, b, tol, maxiter)
[root, iter] = untitled4(f, a, b, eps);
if f(a) * f(b) > 0
error('f(a) and f(b) must have opposite signs')
end
Page 10 of 15
iter = 0;
error = b - a;
while error > tol && iter < maxiter
iter = iter + 1;
c = (a + b) / 2;
fc = f(c);
if fc == 0 || error < eps
root = c;
break
end
if f(a) * fc < 0
b = c;
else
a = c;
end
error = b - a;
end
root = c;
end
fprintf('Root: %.4f\n', root);
fprintf('Iterations: %d\n', iter);
end
Output
Root: 1.7321
Iterations: 5
Fixed Point Method
Code
g = @(x) cos(x);
x0 = 0.5;
tol = 1e-4;
N = 100;
function [root, iter] = fixed_point(f, x0, tol, maxiter)
[root, iter] = fixedPointMethod(g, x0, tol, N);
iter = 0;
error = inf;
Page 11 of 15
while error > tol && iter < maxiter
iter = iter + 1;
x1 = f(x0);
error = abs(x1 - x0);
x0 = x1;
end
root = x1;
end
fprintf('Root: %.4f\n', root);
fprintf('Iterations: %d\n', iter);
end
Output
Root: 0.7391
Iterations: 7
Fixed Point Method for a System of Equations
Code
F = @(x) [sqrt(25 - x(2)^2); x(1) - 4];
x0 = [3; -1];
tol = 0.001;
function [root, iter] = fixed_point_system(F, x0, tol, maxiter)
[root, iter] = fixed_point_system(F, x0, tol, 100);
iter = 0;
error = inf;
while error > tol && iter < maxiter
iter = iter + 1;
x1 = F(x0);
error = norm(x1 - x0)
x0 = x1;
end
root = x1;
end
fprintf('The root is (%f, %f), found in %d iterations\n', root(1), root(2), iter);
end
Page 12 of 15
Output
The root is (4.7434, 0.7434), found in 7 iterations
Finding Iterations Using Jacobian Method
Code
F = @(x) [x(1)^2 + x(2)^2 - 1; x(1) - x(2) + 0.5];
J = @(x) [2*x(1), 2*x(2); 1, -1];
x0 = [0; 0];
tol = 0.001;
maxiter = 100;
function [x, iter] = jacobian_method(F, J, x0, tol, maxiter)
iter = 0;
error = inf;
while error > tol && iter < maxiter
iter = iter + 1;
x = x0 - J(x0) \ F(x0);
error = norm(x - x0);
x0 = x;
end
end
disp('The approximate root is:')
disp(x)
disp('The number of iterations is:')
disp(iter)
end
Output
𝑥𝑥 0.7862
[𝑦𝑦] = [ ]
0.2862
Page 13 of 15
Finding Iterations Using Gauss-Seidel Method
Code
A = [3 -0.1 -0.2; 0.1 7 -0.3; 0.3 -0.2 10];
b = [7.85; -19.3; 71.4];
x0 = [0; 0; 0];
tol = 0.001;
maxiter = 100;
function [x, iter] = gauss_seidel(A, b, x0, tol, maxiter)
n = size(A,1);
if det(A) == 0
error('The matrix is singular')
end
iter = 0;
error = inf;
while error > tol && iter < maxiter
iter = iter + 1;
x = zeros(n,1);
for i = 1:n
s = 0;
for j = 1:n
if j ~= i
s = s + A(i,j) * x(j);
end
end
x(i) = (b(i) - s) / A(i,i);
end
error = norm(x - x0);
x0 = x;
end
end
Page 14 of 15
disp('The solution is:')
disp(x)
disp('The number of iterations is:')
disp(iter)
end
Output
𝑥𝑥 3.0000
[𝑦𝑦] = [−2.5000]
𝑧𝑧 7.0000
Page 15 of 15
