II ‫البرمجة المتقدمة‬

Numerical Methods with Arrays

2nd Class-2nd Semester

Created by
Dr. Saib A. Yousif
ADDITION AND SUBTRACTION
A ± B = C, a i j ± bi j = ci j ------ (1)
• Matrices must be of same (size).
Example:

Then the matrix C that is obtained by adding A

and B :

C=

2
Example:

3
Note: K +A = B, k + aij=bij --- (2)
Where, k=scalar
Examples:

4
ARRAY MULTIPLICATION
A B = C, ci j = 𝑚
𝑘=1 𝑎 𝑖 𝑘 ∗ 𝑏𝑘 𝑗
----- (2)

Example:

5
Examples:

6
7
8
Identity Matrix (I)
• It is a square matrix in which the diagonal elements
are 1s, and the rest of the elements are 0s.
• When the identity matrix multiplies another matrix
(or vector), that matrix (or vector) is unchanged.

• If a matrix A is square, it can be multiplied by the

identity matrix, I, from the left or from the right.
A I = I A = A ---- (3)

9
Determinant of A matrix
Det A = ai1 Ci1+ ai2 Ci2 + -----+ ain Cin ----- (4)
where Cij is called the cofactor and itself is given as
Cij = (-1)i+j det Mij ----- (5)
where Mij is called the minor and is formed by deleting
row i and column j of the array A.
Example: find det A where
Solution:

= (-1)(-7) +(1)(-13) + (2) (8) = 10

Note: In Matlab, = det(A)
10
Inverse of a matrix:
𝟏
𝑨−𝟏 = [𝑪𝒊𝒋 ]’ ------ (6)
𝒅𝒆𝒕 𝑨
Where Cij = cofactor and [𝑪𝒊𝒋 ]’ = transpose
A A-1 = I ----- (7)
Where , A = square matrix, A-1 = inverse of A
Example: Find A-1 of
Solution:

Example: find the inverse of the matrix

Solution:

11
 −7 − 13 8
C= 2 −2 2 −7 2 3
𝐶 ′ = −13 − 2 7
3 7 −2
8 2 −2

Eq. (4): det A = (-1)(-7) +(1)(-13) + (2) (8) = 10

𝐶′
𝐴−1 =
det 𝐴

Note: In MATLAB, inverse of A = inv(A)

12
Example:

13
Solution of system of linear algebraic equations:
The system of equations can be expressed in
matrix form: A X = b ---- (8)
Where:

And the augmented matrix is

14
1. Cramer’s rule
Let the equations:

𝐴1 𝐴2 𝐴3
𝑥1 = , 𝑥2 = , 𝑥3 = ----- (9)
𝐴 𝐴 𝐴

𝐴 , |𝐴1| ,, |𝐴2|

𝐴3
15
Example:
Solve: x1 + 0.306 x3 = 101.48
x2 + 0.702 x3 = 225.78
-2x1+ x2 = 0, by Cramer’s rule using MATLAB program.
Solution:
1 0 0.306 101.48 0 0.306 1 101.48 0.306
A= 0 1 0.702 , A1 = 225.78 1 0.702 A2 = 0 225.78 0.702
−2 1 0
0 1 0 −2 0 0
1 0 101.48
A3 = 0 1 225.78
−2 1 0
𝐴 = 1[(1)(0) – (0.702)(1)] -0[(0)(0)-(0.702)(-2)] + 0.306[(0)(1) –(1)(-
2)] = - 0.702 + 0.306 (2) = - 0.09
𝐴1 = 101.48 [(1)(0) – (0.702)(1)] -0[(225.78)(0)-(0.702)(0)] +
0.306[(225.78)(1) –(1)(0)] = 101.48(-0.702) + 0.306 (225.78) = -2.1503
16
X1 = -2.1503/(-0.09) = 23.89
𝐴2 = 1[(225.78)(0) – (0.702)(0)] -101.48[(0)(0)-(0.702)(-2)] + 0.306[(0)(0) –
(225.78)(-2)] = -101.48 (1.404) + 0.306 (225.78) (2) = -4.3006

𝐴3 = 1[(1)(0) – (225.78)(1)] -0[(0)(0)-(225.78)(-2)] + 101.48[(0)(1) –(1)(-2)] =

- 0.702 + 0.306 (2)= -225.78+101.48(2) = -22.82

X3 = -22.82/(-0.09) = 253.5556

17
Cramer’s rule using MATLAB

18
2. Gaussian Elimination
The matrix A reduces to lower triangular matrix which is a special type of
square matrix where the elements below the diagonal are zero
Example: Solve

Solution:

Step 1. Elimination of x1

Step 2. Elimination of x2
19
Step 3: Back substitution
x3 = -190/-95 = 2

4
2

20
 Example: Solve by Gaussian elimination
6 − 1 1 13 Exchange R2 by R1
1 1 1 9
10 1 − 1 19
Step 1. Elimination of x1
1 1 1 9
6 −1 1 13 R2 – 6 R1 1 1 1 9
10 1 − 1 19 R3 – 10 R1 6−6 −1−6 1−6 13 − 54
10 − 10 1 − 10 − 1 − 10 19 − 90
Step 2. Elimination of x2
1 1 1 9
0 −7 −5 − 41
0 −9 − 11 − 71 R3 – (9/7) R2

1 1 1 9 1 1 1 9
0 − 7 −5 − 41 0 −7 −5 − 41
9 9 9 32
0 −9− −7 − 11 − −5 − 71 − ( )(−41)
7 7 7 0 0 − − 128/7
7

21
Back substitution
-32/7 x3= -128/7 X3 = 4

-7 x2 – 5 x3= -41 -7 x2 -5(4) =-41 X2 =3

X1 +x2 +x3 = 3
X1 +3+4 =9
X1 = 2

22
 Gauss Elimination using Matlab-1

% Elimination of x1
d21= a(2,1)/a(1,1); % Back Substitution
a(2,1)= a(2,1)-d21*a(1,1); X3=b(3)/a(3,3)
a(2,2)= a(2,2)-d21*a(1,2); X2=(b(2)-X3*a(2,3))/a(2,2)
a(2,3)= a(2,3)-d21*a(1,3); X1=(b(1)-X3*a(1,3)-X2*a(1,2))/a(1,1)
b(2)=b(2)-d21*b(1);
d31= a(3,1)/a(1,1);
a(3,1)= a(3,1)-d31*a(1,1);
a(3,2)= a(3,2)-d31*a(1,2);
a(3,3)= a(3,3)-d31*a(1,3);
b(3)=b(3)-d31*b(1);
% Elimination of x2
d32= a(3,2)/a(2,2);
a(3,2) = a(3,2)-d32*a(2,2);
a(3,3) = a(3,3)-d32*a(2,3);
b(3)=b(3)-d32*b(2);

23
Gauss elimination using MATLAB-2
% Matlab Program to solve (nxn) system equation
% by using Gaussian Elimination method
clear ; clc ; close all
n = input('Please Enter the size of the equation system n = ') ;
C = input('Please Enter the elements of the Matrix C ' ) ;
b = input('Please Enter the elements of the Matrix b ' ) ;
dett = det(C)
if dett == 0
print('This system unsolvable because det(C) = 0 ')
else
b = b'
A=[C b]
for j = 1:(n-1)
for i= (j+1) : n
mult = A(i,j)/A(j,j) ;
for k= j:n+1
A(i,k) = A(i,k) - mult*A(j,k) ;
A
end
end
end
for p = n:-1:1
for r = p+1:n
x(p) = A(p,r)/A(p,r-1)
end
end
24
end
3. Gauss - Jordan
The elimination step results in an identity matrix rather
than a triangular matrix and the solution is obtained in
the right-hand-side vector.
Eq.8: A X = b ( multiply by A-1 )

I X = A-1 B
X = A-1 B ----- (9)
Example: Solve by Guass-Jourdan the following equations
system

25
 Divide R2
by (-7)

1 1 1 9 R1 – R2 1−0 1−1 1−
5
9 − 41/7
5 7
0 1 41/7 5
7 0 1
7
41/7
0 −9 − 11 − 71 5
R3+9 R2 0−0 −9+9 − 11 + 9 ∗ − 71 + 9 ∗ 41/7
7

2 Divide R3 by 1 0
2
22/7 R1 – (2/7)R3
1 0 22/7 7
7 (-32/7) 5
5 0 1 41/7 R2-(5/7)R3
0 1 41/7 7
7 0 0 1 4
32
0 0 − − 128/7
26 7
 2 2 22 2
1 0 − − ∗4 1 0 0 2
7 7 7 7
0 1 0 3
5 5 41 5
0 1 − − ∗4 0 0 1 4
7 7 7 7
0 0 1 4

X3 = 4, x2 = 3, x1 = 2

27
Example: Gauss-Jordan using MATLAB ,Solve

28
 4. Left division,

Example

29
5. Right division, /
The right division is used to solve the matrix equation, X C = D .
In this equation X and D are row vectors. This equation can be
solved by multiplying, on the right, both sides by the inverse of
C: X C = D --- (10)

X = D/C ---- (11)

Example:

30
6. Jacobi method
It is an iterative method.
1. The Initial values of x1, x2, and x3 are guessed (0,0,0).
2. The guessed values of x1, x2 and x3 are substituted in Eqs.(12-
14), to find new values of x1, x2, and x3.

3. The entire procedure is repeated until our solution converges

closely enough to the true values. Convergence can be
checked using the criterion.

31
Example:
Solve by Jacobi method with initial guess (x1,x2,x3)=(0,0,0).

Solution: The above equations are rearranged to

Now, make the initial guess x1 = 0, x2 = 0, x3 = 0.

We can continue this iterations for the values k = 0, 1, 2,3,….

32
 Jacobi method using MATLAB

33
6. Gauss -Seidel
It is an iterative method.
1. The Initial values of x1, x2, and x3 are guessed (0,0,0).
2. The guessed values of x2 and x3 are substituted in Eq.(12) to
find x1.
3. The guessed value of x3 and refined value of x1 are
substituted in Eq.(13) to find x2.
4. The refined values of x1 and x2 are substituted in Eq.(13) to
find x3.
5. The entire procedure is repeated until our solution
converges closely enough to the true values. Convergence
can be checked using the criterion Eq.(15).

34
Example:
Solve by Gauss-Seidel method with initial guess
(x1,x2,x3)=(0,0,0). The true solution is x1 = 3,
x2 = −2.5, and x3 = 7.
Solution: The above equations are rearranged to

From above equations

For the second iteration, the same process is repeated to compute

35
• Additional iterations could be applied to
improve the answers.

36
Gauss Seidel implemented in MATLAB

37
ELEMENT-BY-ELEMENT OPERATIONS

If two vectors a and b are

then element-by-element multiplication, division, and
exponentiation of the two vectors gives:

38
If two matrices A and B are

then element-by-element multiplication and division of the two

matrices give:

Element-by-element exponentiation of matrix A gives:

39
Examples using matlab

40
41
42
BUILT-IN FUNCTIONS FOR ANALYZING ARRAYS

43
44
Problems
1. For the following distillation column calculate the values of F1, F3 and F4?

Ans: F1 = 1000, F2= 500, F3 = 1000

2. Calculate Reynold Number at a diameter D=0.2 m, using different velocity's
as u=0.1,0.2,0.3 …1 m/s knowing that: Re= (ρud)/µ , µ=0.001 , ρ =1000
Ans: Re =
1.0e+005 *
0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 1.8000 2.0000

45
3. Estimate the average density of a Water/Ethanol mixture at
different water compositions knowing that, Water density =1000
kg/m3, Ethanol density =780 kg/m3. Mixture density= Xwater ×
Water density + Xethanol × Ethanol density .
Ans: Pav = 780 802 824 846 868 890 912 934 956 978 1000

4. convert the temperature in Celsius into °F, TF = 1.8 Tc + 32 and

then into °R, TR = TF+459.69 for every temperature from 0
increasing 15 to 100°C. Combine the three results into one
matrix.
Ans:

46
5. a. Find the physical properties for water in the range of
temperatures from 273 to 323 K.

b.
Define the 5 x 4 matrix and find the content of the following matrices and check
your results for content using Matlab.

6. Calculate the values of XA,XB,YA,YB, L and V for the vapor

liquid separator shown in fig.

47
Ans: xa = 0.3077, xb = 0.6923, ya = 0.5846, yb =0.4154, L = 66.6667, V =
33.3333
7. The variation of vapor pressure p (in units of mm Hg) of benzene
with temperature in the range of can be modeled with
the equation

where a = 34172 and b = 7.9622 are material constants and T is

absolute temperature (K). calculates the pressure for various
temperatures. Create a vector of temperatures from 00 C to 45 0C with
increments of 2 degrees,

48
and display a two - column table p and T, where the first column
temperatures in 0C, and the second column the corresponding
pressures in mm Hg.
Ans:

8. An ice cream container shaped as a frustum of a cone with is

designed to have a volume of 1,000 cm3. Determine R1, R2, and
the surface area, S, of the paper for containers with heights h of
8, 10, 12, 14, and 16 cm. Display the results in a table. The
volume of the container, V, and the surface area of the paper are
given by:

Ans:

49
9. An aluminum sphere ( cm) is dropped in a glass cylinder filled
with glycerin. The velocity of the sphere as a function of time can
be modeled by the equation

where V is the volume of the sphere, g = 9.81 m/s2 is the

gravitational acceleration, k = 0.0018 is a constant, and ρal =
2700 kg/m3 and ρgl = 1260 kg/m3 are the density of aluminum
and glycerin, respectively. Determine the velocity of the sphere
for t = 0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, and 0.35 s.
Ans:

50
10.

Ans:

11.

Ans: 42
51
12. Solve by Gauss elimination, Gauss-Jordan, Jacobi and Gauss –
Seidel methods:
a. Ans.
d.

b. Ans.

c. Ans.

52