Numerical Methods

By Nguyen Khai Minh & Lan Anh

from USTH Learning Support
May 2025

Contents
1 Introduction 2

2 Errors and uncertainty (Lecture 1) 3

2.1 Approximations and Round-Off errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Significant figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.3 Accuracy and precision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.4 Error definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3 MATLAB 5
3.1 Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2 Mathematical operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.3 Vectors and Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.3.1 Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.3.2 Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.3.3 Matrix operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.4 The colon (:) operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.5 Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.6 Some functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4 Roots of Non-linear equations (Lecture 2) 13

4.1 Bisection method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.2 Newton-Raphson method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

5 Roots of systems of non-linear equations : Newton - Raphson method 17

5.1 Some example exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

6 Gauss elimination 19
6.1 Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
6.2 General steps of Gauss Elimination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
6.3 Naive Gauss Elimination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
6.4 Gauss elimination : Pivoting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

7 Iterative methods 26
7.1 Gauss - Seidel method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
7.2 Jacobi method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
7.3 Improvement of convergence using relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

1
1 Introduction
Part in test: Theory (7 MCQs, 4 points. Include: motivation; roles, applications, formulae and insights of different
types of numerical methods

What is numerical methods:

Numerical methods are techniques that use basic arithmetic and computers to find approximate solutions to mathe-
matical problems that are too hard to solve for exact solution.

Motivation of learning numerical methods:

Real-life problems are usually too complex to solve by hand - or analytically. We apply numerical methods to these
problems, using computers, to find good enough (approximation) and practical answers. And it should be done
quickly. Numerical methods are used widely in engineering and science.

Examples of applications:
• Solving polynomial equations: Simple linear or quadratic equations, for example: x − 3 = 0 and x2 − 2x −
3 = 0, can be solved analytically. But as the degree increases, usually for degree higher than 5, analytical
solutions become impossible. → Numerical methods: Newton-Raphson, bisection methods to find the roots
approximately.
• Fluid flow - Navier-Stokes Equations: Navier–Stokes momentum equation (conservative form):

∂u 1 1
+ u · ∇u = − ∇p + ν∇2 u + ν∇(∇ · u) + g
∂t ρ 3
This partial differential equations (PDEs) describes how fluid velocity evolves. They are too complex to solve
analytically. → Numerical methods: CDF - Computational Fluid Dynamics, techniques like finite volume
or finite difference methods to discretize and solve flow fields in space and time.
• Heat Transfer - Fourier’s Law: Steady state equation:

⃗q = −k∇T

Fourier’s law models heat conduction. Analytical methods fail in complex geometries or with nonuniform materi-
als. → Numerical methods: finite difference or finite element methods to simulate temperature distributions
in systems.

• Structural and mechanical engineering: To predict mechanical response of structures, we can apply New-
ton’s 2nd law: X
F⃗ = m⃗a
Since real-world mechanical systems, turbines, or machines are very complex, analytical methods are insufficient.
→ Numerical methods: finite element analysis (FEA) to evaluate stress, strain, and deformation in complex
geometries under various loading conditions.

Characteristics of Numerical methods:

• The procedure of numerical methods is iterative, meaning it is repeated many times, to get closer and closer
to true solution (enhance accuracy after each iteration).
• The solution of numerical methods is an approximation to the true solution.
• The algorithm is simple and can be easily programmed.
• There is a possibility that numerical method’s procedure diverge from the true solution, which is not what we
want. (we want the solution to converge to the true solution)

Note: Various applications are presented in the slides and you should definitely read them thoroughly.

2
2 Errors and uncertainty (Lecture 1)
2.1 Approximations and Round-Off errors
• Numerical methods will always give errors since the solutions are approximated. These are approximation
and round-off errors.
• In cases that we can use analytical method, we can compute the exact value of error of numerical method.
Otherwise, we would calculate estimates of the errors.

2.2 Significant figures

• Significant figures (sf ) are digits of a number that are necessary for the number to be ”accurate enough”. We
can use these significant figures with confidence (in a numerical result).
• Some rules and examples:

– Non-zero digits are always significant. Ex: 3.141 has 4 sf.

– Zeros between non-zero digits are significant. Ex: 1005 has 4 sf.
– Leading zeros are NOT significant. Ex: 0.0001845 has 4 sf; 0.01845 and 0.1845 also have 4 sf.
– Trailing zeros after decimal point are significant. Ex: 1.00 has 3 sf.
– Trailing zeros in whole number are significant only if they have a decimal point. Ex: 100. has 3 sf; 100 has
1 sf.
– Exact numbers, and irrationally defined numbers like Euler’s number (e) and pi (π), have an infinite number
of significant figures.
– A × 10x has the same number of sf as A. Ex: 4.53 × 104 has 3 sf, 4.530 × 104 has 4 sf and 4.5300 × 104 has
5 sf.
• Positions of significant figures:

2.3 Accuracy and precision

• Accuracy measures how close the computed result is to the true value.
• Precision measures how close the results are with each other.

The concept is illustrated in Figure 1, with the center of the circles represent true value. (true solution)

3
 Figure 1: Accuracy and Precision

There are also 2 types of measurement errors:

• Inaccuracy (also called bias): a systematic fault that lead to deviation from the truth.
• Imprecision (also called uncertainty): magnitude of the scatter.

In Figure 1, (a) is inaccurate and imprecise, (b) is accurate but has low precision (high uncertainty), (c) is inaccurate
but precise and (d) is accurate and precise.

2.4 Error definitions

There are two types of error arise from approximations:

• Truncation errors: happen when simplifying mathematical processes to make computations more manageable.
• Round-off errors: happen when numbers having limited significant figures.
For both type, we have the following relationship:

True value = approximation + error

=⇒ True error: Et = true value − approximation

But using just the absolute error ignores the size of the quantity being measured. For instance, a 1 cm error is more
serious when measuring a small object (like a rivet) than a large one (like a bridge). To address this, we use the true
fractional relative error:
true error
True fractional relative error =
true value
Since the true error is defined as:

True error = true value − approximation

we can express the true percent relative error as:
true error
εt = × 100%
true value

4
 where εt represents the true percent relative error.
εt is for analytical solutions, for numerical approximation, we have:
approximation error
εa = × 100%
approximation
curent approximation − previous approximation
=⇒ εa = × 100%
current approximation

Another quantity that we have to consider is εs , the pre-specified percent tolerance, which is the limit of the error of
our computation. That means, our computation will be repeated until:

|εa | < εs

Example: Calculation of errors

Suppose that you have the task of measuring the lengths of a bridge and a rivet and come up with 9999 and 9 cm,
respectively. If the true values are 10,000 and 10 cm, respectively, compute (a) the true error and (b) the true percent
relative error for each case. (Example 3.1 in Textbook: Numerical Methods for Engineer)
Solution.
The error for measuring the bridge is:
Et = 10, 000 − 9999 = 1 cm
and for the rivet it is
Et = 10 − 9 = 1 cm
The percent relative error for the bridge is:
1
εt = × 100% = 0.01%
10, 000
and for the rivet it is
1
× 100% = 10%
εt =
10
Thus, although both measurements have an error of 1 cm, the relative error for the rivet is much greater. We would
conclude that we have done an adequate job of measuring the bridge, whereas our estimate for the rivet needs to be
much more accurate.

3 MATLAB
3.1 Variables
To assign a (new) variable, simply put:

>> x = 5.71;

If you type without the semicolon (;) the output will be printed immediately in the Command Window after Enter.
So to view variable content, simply type the variable and Enter. You can also see all the variables you are working
with in the Workspace Window.

If you do not explicitly create a variable, the result will be assigned to the default variable named ans, for example:

>> sin(1)

ans =
0.8415

5
 Operator Description Example Result

^ Exponentiation 2^3 8
- Negation -5 -5
*, / Multiplication and Division 6/2 3
\ Left Division 2\6 3
+, - Addition and Subtraction 3+4 7

3.2 Mathematical operations

You got this!

Some relational operators:

Operator Description Example Result

== Equal to 2==3 0
=
~ Not equal to 2~=3 1
> Greater than 6>2 1
< Less than 6<2 0
>= Greater than or equal to 3>=4 0
<= Less than or equal to 3<=4 1

Some logical operators:

Operator Description Example Result

& AND A&B 1 (true) only where both matrices contain nonzero values.
| OR A|B 1 (true) where either matrix contains a nonzero value.
~ NOT ~A 1 (true) where A is zero and 0 (false) where A is nonzero.

3.3 Vectors and Matrix

3.3.1 Vectors
To create a row vector (dimension: 1 × 5 or 1 row, 5 columns), there are 2 ways. Note that there is no semicolon (;)
between elements.

>> a = [1 2 3 4 5];
>> a = 1:5

To create a column vector, separate the elements with a semicolon (;):

>> a = [1; 2; 3; 4; 5];

>> a = [1 2 3 4 5]’

The 2nd way, using ′ is transpose of a vector. In particular, we transposed a row vector to create a column vector.

6
3.3.2 Matrix
To create a matrix of 3 rows and 3 columns, we type each row then separate the rows by the semicolon ; as the
following:

>> A = [1 2 3; 4 5 6; 7 8 9];

Special matrices: matrices of all zeros or ones. You type the dimension of the matrix inside the bracket in the order
(# of row, # of column)

>> z = zeros(2, 4);

>> x = ones(2, 3);

To access an element of a matrix, use round bracket and the position of the element in the order (row index, column index)
For example, to access the element in the first row and 2nd column:

>> A = [1 2 3; 4 5 6; 7 8 9];
>> A(1,2)

ans=
2

Note that indexing in MATLAB starts from 1.

3.3.3 Matrix operations

A’ - Transpose

>> A = [3 2 1; 5 1 0; 2 1 7];
>> A’

ans=
3 5 2
2 1 1
1 0 7

B*A - Multiplication

>> B = [1 3 1; 4 9 5; 2 7 2];
>> B*A

ans=
20 6 8
67 22 39
45 13 16

Remind: this is the normal matrix multiplication and usually B ∗ A ̸= A ∗ B.

B.*A - Multiplication by each element (With the . you can perform the following operation on each element of the
matrix)

>> B.*A

7
 ans=
3 6 1
20 9 0
4 7 14

For this case, B. ∗ A = A. ∗ B.

B/A - Matrix division

>> B/A

ans=
2.0652 -0.9783 -0.1522
5.9130 -2.6957 -0.1304
4.8913 -2.3696 -0.4130

We also have B./A - Matrix division by elements. (B./A ̸= A./B just like B/A ̸= A/B)

[B A] - Joining matrices horizontally

>> [B A]

ans =
1 3 1 3 2 1
4 9 5 5 1 0
2 7 2 2 1 7

[B; A] - Joining matrices vertically:

>> [B; A]

ans =
1 3 1
4 9 5
2 7 2
3 2 1
5 1 0
2 1 7

3.4 The colon (:) operator

Create a vector from 1 to 10:

>> x = 1:10

x=
1 2 3 4 5 6 7 8 9 10

Create a vector from 1 to 10 with specified step of 2 (difference between 2 elements):

>> x = 1:2:10

x=
1 3 5 7 9

8
We can also create a decreasing vector with step like this:

>> x = 10:-2:0

x=
10 8 6 4 2 0

Indexing using colon operator:

• A(:, n) - Extract the nth column of matrix A:

>> A = [1 2 3; 4 5 6; 7 8 9];
>> A(:,2) % extract elements in every row & in 2nd column

ans=
2
5
8

• A(m, :) - Extract the mth row of matrix A:

>> A = [1 2 3; 4 5 6; 7 8 9];
>> A(1,:) % extract elements in the 1st row & all column

ans=
1 2 3

• A(:) - Extract all elements of matrix A: (note that the result has dimension 6 × 1)

>> A = [1 2 3; 4 5 6];
>> A(:) % extract all elements

ans=
1
2
3
4
5
6

• A(n, j:k) - Extract elements in nth row and from column jth to kth:

>> A = [4 5 6; 7 8 9];
>> A(1, 2:3) % 1st row, column 2 to 3

ans=
5 6

9
 • A(:) - Extract all elements of A and put the result in form of a column: (the matrix is read by column)

>> A = [4 5 6; 7 8 9];
>> A(:) % all elements

ans=
4
7
5
8
6
9

An example from Labwork:

One solution:

Note: the . is to perform the operation on each element of a matrix.

3.5 Visualization
The 4 basic lines to make a plot in MATLAB:

>> figure
>> x = 1:10;
>> y = 2*x + 1;
>> plot(x, y)

Note that for the plot() function, we specify x and y. Style of the plot (type of plot, color, etc.) can be also be
added. More details: Link.
To make titles and labels for x and y-axis:

10
 title(’2-D Line Plot’)
xlabel(’x’)
ylabel(’cos(5x)’)

One way to plot multiple functions in one plot can be:

>> x = 0:pi/10:2*pi;
>> y1 = sin(x);
>> y2 = sin(x-0.25);
>> y3 = sin(x-0.5);

>> figure
>> plot(x,y1,’g’,x,y2,’b--o’,x,y3,’c*’)

Figure 2: Result.

3.6 Some functions

sum: Calculates the sum of elements in an array or matrix.

a = [1 2 3 4 5];
s = sum(a);
% s = 15

syms: Declares symbolic variables.

syms x;
expr = x^2 + 3*x + 2;

symsum: Computes symbolic summation.

syms k;
s = symsum(1/k^2, k, 1, Inf);

11
diag: Creates a diagonal matrix or extracts the diagonal of a square matrix.

A = diag([1 2 3])

A=
1 0 0
0 2 0
0 0 3

diag(A)

ans=
1
2
3

linspace: Generates a linearly spaced vector. Example: create a vector from 0 to 1, have 5 elements equally spaced
(step = (1-0)/5 = 0.25)

x = linspace(0, 1, 5)

x=
0 0.25 0.5 0.75 1

f zero: Finds a root of a nonlinear function. (note: different initial points may give different results)

f = @(x) x^2 - 4; % function

x0 = 1; % initial point
root = fzero(f, x0)

root=
2

f = @(x) x^2 - 4; % function

x0 = -1; % initial point
root = fzero(f, x0)

root=
-2

solve: Solves equations symbolically.

syms x;
sol = solve(x^2 - 4 == 0, x)

sol=
-2
2

sin: Computes sine of input in radians. (similar for cos, tan, cot)

sin(pi/2)

12
 ans=
1
f or loop:
Example (from Labwork): Write a script to calculate: 12 + 22 + ... + 10002

result = 0;
for i = 1:1000
result = result + i^2;
end
fprintf(’Result: %d’, result)

Result: 333833500
while loop:
Example (from Labwork): Write code with while loop that makes the computer count only the odd number from 1 to
N (N: integer), for example from 1 to 5, there are 3 odd number.

N = 5;
i = 1;
odd_count = 0;
while i <= N
if mod(i, 2) == 1
odd_count = odd_count + 1;
end
i = i + 1;
end
fprintf(’Count of odd numbers from 1 to %d: %d.’, N, odd_count)

Count of odd numbers from 1 to 5: 3

4 Roots of Non-linear equations (Lecture 2)

Non-linear equations are equations with exponents greater than one, for example, ax2 + bx + c = 0, y 2 + dy/dx = 0, ...

4.1 Bisection method

Method
Bisection method is one method of the type: Bracketing methods. The result is bracketed by two initial guesses (2
endpoints).

Figure 3: Illustration of bisectional method

13
Prerequisite: consider the function on an interval [a; b] given that f (a)f (b) < 0. This guarantees that the function =
0 at some point.
Steps:
1. Calculate the midpoint of the interval: c = 0.5 ∗ (a + b)

2. Calculate f (c)
3. Check if |f (c)| < ε or |b − a| < ε, stop and return c as the root. (ε is a prespecified stopping criterion)
4. Examine the sign of f(c) and change the interval to [a, c] or [c, b] if f (a)f (c) < 0 or f (c)f (b) < 0. Back to step 1.

We have estimated error:

xnew
r − xold
r
|εa | = new
× 100%
xr

• xnew
r : the root for the present iteration

• xold
r : the root from the previous iteration

When |εa | ≤ εs , the computation is terminated.

εs : prespecified stopping criterion (the prespecified acceptable level)

Limitations
Limitations of Bisection method:
• The bisection method only finds roots where the function crosses the x axis. It cannot find roots where the
function is tangent to the x axis.
There would be no sign change in this case.

• The bisection method can be fooled by singularities in the function.

A singularity is a point where the function is not defined or goes to infinity. The sign change could happen
around this singularity, for example if the function jumps from very large positive value to a very large negative
value.
• The bisection method cannot find complex roots of polynomials.
Because complex roots do not lie on the real number line.

Example
Find the root of f (x) = x3 − x − 2
Choose a and b that f (a)f (b) < 0: a = 1, b = 2.
Set stopping criterion: relative error εs = 5% = 0.05

Solution:
Step 1: Compute midpoint.
The midpoint is:
1+2
xr = = 1.5
2
Step 2 : Evaluate function

f (1.5) = 1.53 − 1.5 − 2 = 3.375 − 1.5 − 2 = −0.125

(this value is not sufficiently small)

Step 3 : Compute approximate relative error
Previous root: take either a or b, take a = 1, xcurrent
r = 1.5

1.5 − 1 0.5
|εa | = × 100 = × 100 ≈ 33.33%
1.5 1.5

14
 εa > εs =⇒ continue

Step 4 : Update interval

Sincef (1) = −2, f (1.5) = −0.124 and f (2) = 4, the root lies in [1.5, 2]. New interval is:

a = 1.5, b=2

=⇒ End of iteration 1. Iteration 2: New root: midpoint of [1.5, 2]:

1.5 + 2
xr = = 1.75
2
Evaluate:
f (1.75) = 1.609375
(this value is not sufficiently small) Approximate error:

1.75 − 1.5 0.25

|εa | = × 100 = × 100 ≈ 14.29%
1.75 1.75

εa > εs =⇒ continue

Update interval:
f (1.5) = −0.125
f (1.75) = 1.609375
f (2) = 4
=⇒ new interval: [1.5, 1.75]
Continue with another iteration, until εa <= εs
From iteration 3, we have a table

Iteration a b c εa
1 1 2 1.5 100%
2 1.5 2 1.75 14.29%
3 1.5 1.75 1.625 7.69%
4 1.5 1.625 1.5625 4%

εa <= εs =⇒ iterations end

4.2 Newton-Raphson method

Method
This method only requires a starting point a near the root. At this point, we draw a the tangent line to the function
f (x), which has a slope equal to the derivative f ′ (a). This tangent line intersects the x-axis at a new point a′ , which
is given by:
f (a)
a′ = a − ′ (the Newton-Raphson formula)
f (a)
Good news is this point is usually closer to the root of f than a, so we use it as a better approximation. By repeating
this step, we can get closer and closer to the real root.

15
 Figure 4: Illustration of Newton-Raphson method.

Steps:
1. Evaluate the derivative f ′ (x) symbolically.

2. Have an initial guess xi , estimate the new value xi+1 by the Newton-Raphson formula:

f (xi )
xi+1 = xi −
f ′ (xi )

3. Find the relative approximate error:

xi+1 − xi
|εa | = × 100%
xi+1
Terminate if |εa | ≤ εs , else start another iteration (from step 2).

Limitations
• This algorithm can have slow convergence due to the nature of the function or a poor first initial guess.
• Newton-Raphson can oscillate around a local maximum or local minimum, as near-zero slope is reached in this
case. In general, this algorithm has problems with near-zero slopes, sometimes, it can point us to a really far
away region from the real root. Look at Figure 5, notice the sequence of x1 , x2 , ...

Quoted from textbook: “Thus, there is no general convergence criterion for Newton-Raphson. Its convergence
depends on the nature of the function and on the accuracy of the initial guess. The only remedy is to have an initial
guess that is “sufficiently” close to the root. And for some functions, no guess will work! Good guesses are usually
predicated on knowledge of the physical problem setting or on devices such as graphs that provide insight into the
behavior of the solution. The lack of a general convergence criterion also suggests that good computer software should
be designed to recognize slow convergence or divergence.”

16
 Figure 5: Illustration of problems with Newton Raphson Method

Example
Example: Estimate the root of f (x) = e−x − x, employing an initial guess of x0 = 0.

Solution:

Set stopping criterion: εs = 5%

Step 1: Find f ′ (x) symbolically:

f ′ (x) = −e−x − 1
Step 2: Estimate xi+1
f (xi ) e−xi − xi
xi+1 = xi − = x i −
f ′ (xi ) −e−xi − 1
For the first iteration, x0 = 0 and:
e−0 − 0
x1 = 0 − = 0.5
−e−0 − 1
Step 3: Estimate the relative error:
0.5 − 0
|εa | = × 100% = 100%
0.5
|εa | > εs =⇒ End of iteration 1, start iteration 2.
Iteration 2:
e−0.5 − 0.5
x2 = 0.5 − = 0.5663
−e−0.5 − 1
0.5663 − 0.5
|εa | = × 100% = 11.7%
0.5663
|εa | > εs =⇒ End of iteration 2, start iteration 3.
Continue until |εa | < εs . It should take 1 more iteration. (Check this yourself, i might be wrong.) We have a table of
all iterations:
i xi εa
1 0.5 100%
2 0.5663 11.7%
3 0.5671 0.15%
We can see that |εa | <= εs =⇒ iterations end

5 Roots of systems of non-linear equations : Newton - Raphson method

• Target : Find the root of a system of nonlinear equations.

• Recall that the Newton-Raphson method was predicated on employing the derivative (that is, the slope) of a
function to estimate its intercept with the axis of the independent variable—that is, the root.
• This estimate based on a first - order Taylor series expansion :

f (xi+1 ) = f (xi ) + (xi+1 − xi )f ′ (xi )

17
 – xi : the initial guess
– xi+1 : the point where the slope intercepts the x axis → f (xi+1 ) = 0 → we have the single - equation form
of the Newton - Raphson method :
f (xi )
xi+1 = xi − ′
f (xi )
→ The multi - equation form is derived in an identical fashion. However, a multivariable Taylor series must be
used.
Suppose that we have a systems of 2 nonlinear equations, with 2 unknown variable x, y :
(
u(x, y) = 0
v(x, y) = 0

We can write a first - order Taylor for each nonlinear equation :

∂ui ∂ui
ui+1 = ui + (xi+1 − xi ) + (yi+1 − yi )
∂x ∂y
∂vi ∂vi
vi+1 = vi + (xi+1 − xi ) + (yi+1 − yi )
∂x ∂y
Just as for the single - equation form, ui+1 and vi+1 will be 0
→ 2 above equations can be rewritten as :
∂ui ∂ui ∂ui ∂ui
xi+1 + yi+1 = −ui + xi + yi
∂x ∂yi ∂x ∂y
∂vi ∂vi ∂vi ∂vi
xi+1 + yi+1 = −vi + xi + yi
∂x ∂y ∂x ∂y
Because every value subscripted with i′ s are known from the latest guess → the only unknowns are xi+1 and yi+1 :

ui ∂v ∂ui
∂y − vi ∂y
i

xi+1 = xi − ∂ui ∂vi ∂ui ∂vi

∂x ∂y − ∂y ∂x

vi ∂u
∂x
i
− ui ∂v
∂x
i

yi+1 = yi − ∂ui ∂vi ∂ui ∂vi

∂x ∂y − ∂y ∂x

The dominator of each of these equations is the determinant of Jacobian of the system.

5.1 Some example exercises

Exercise 1 : Use the multiple - equation Newton - Raphson method to determine roots of the following equations.
Note that a correct pair of roots is x = 2 and y = 3. Initiate the computation with guesses of x = 1.5 and y = 3.5
(
u(x, y) = x2 + xy − 10 = 0
v(x, y) = y + 3xy 2 − 57 = 0

Solution : First, we compute the partial derivatives and evaluate them at the initial guesses of x and y :
∂u0 ∂u0
= 2x + y = 2(1.5) + 3.5 = 6.5 = x = 1.5
∂x ∂y
∂v0 ∂v0
= 3y 2 = 3(3.5)2 = 36.75 = 1 + 6xy = 1 + 6(1.5)(3.5) = 32.5
∂x ∂y
Thus, the determinant of the Jacobian for the first iteration is
6.5(32.5) − 1.5(36.75) = 156.125
The values of the functions can be evaluated at the initial guesses as :
u0 = (1.5)2 + 1.5(3.5) − 10 = −2.5
v0 = 3.5 + 3(1.5)(3.5)2 − 57 = 1.625

18
These values can be substituted into Newton - Raphson method for mutiple - equation to give :

−2.5(32.5) − 1.625(1.5)
x = 1.5 − = 2.03603
156.125
1.625(6.5) − (−2.5)(36.75)
y = 3.5 − = 2.84388
156.125
Thus, the results are converging to the true values of x = 2 and y = 3. The computation can be repeated until an
acceptable accuracy is obtained.

Example 2 : Use the multiple-equation Newton-Raphson method to determine the roots of equations by using 3
iterations with initial guess : x = 0.8, y = 1.8. Round to 5 decimal places.
(
x−y+1=0
x2 + y 2 − 4 = 0

Solution : We will calculate by each iteration

Iteration 1 :

First, we have to calculate the values of partial derivatives at those initial guess
∂u0 ∂u0
=1 = −1
∂x ∂y
∂v0 ∂v0
= 2x = 1.6 = 2y = 3.6
∂x ∂y
Then, we compute the determinant of the Jacobian :

1(3.6) − (−1)(1.6) = 5.2

The values of the functions can be evaluated at the initial guess:

u0 = 0.8 − 1.8 + 1 = 0
v0 = 0.82 + 1.82 − 4 = −0.12

Finally, we substitute into Newton - Raphson equation to get

0(1) − (−0.12)(−1)
x = 0.8 − = 0.82308
5.2
(−0.12)(1) − (0)(1.6)
y = 1.8 − = 1.82308
5.2

Iteration 2 and Iteration 3 : The steps are similar to Iteration 1, so do it yourself.

Answer : x = 0.82288, y = 1.82288

6 Gauss elimination
6.1 Matrix
In this part, I will introduce some matrix notations, special types of matrix, and how to represent linear algebraic
equations in matrix form.

19
• [A] : a shorthand notation for m by n matrix.
• aij an individual element of the matrix.
• a11 , a22 , a33 : the principal or main diagonal of the matrix.
• aij = aji : symmetric matrix.
• m = n : square matrices.

• Diagonal matrix  
a11 0 0
[A] =  0 a22 0 
0 0 a33

• Upper triangular matrix  

a11 a12 a13
[A] =  0 a22 a23 
0 0 a33

• Lower triangular matrix  

a11 0 0
[A] = a21 a22 0 
a31 a32 a33

• Identity matrix  
1 0 0
[I] = 0 1 0
0 0 1

• Banded matrix  
a11 a12 0 0
a21 a22 a23 0 
[A] = 
 0

a32 a33 a34 
0 0 a43 a44

• Matrices with row dimension n = 1 ⇒ row vectors

 
[B] = b1 b2 · · · bm

• Matrices with column dimension n = 1 ⇒ column vectors

 
c1
 c2 
[C] =  . 
 
 .. 
cn

20
 • A 3 × 3 set of linear equations:
a11 x1 + a12 x2 + a13 x3 = b1
a21 x1 + a22 x2 + a23 x3 = b2
a31 x1 + a32 x2 + a33 x3 = b3

• Can be expressed in matrix form:

     
a11 a12 a13 x1 b1
[A]{x} = {b} where [A] = a21 a22 a23  , {x} = x2  , {b} = b2 
a31 a32 a33 x3 b3

• Solution:
{x} = [A]−1 {b}

6.2 General steps of Gauss Elimination

The procedure consisted of two steps:
• Forward elimination of unknown : We will eliminate one of the unknowns from the equations. The result is
that we have 1 equation with 1 variable
• Back substitution: The equation above can be solved directly and the result is back substituted into one of the
original equations to solve for the remaining variables.

6.3 Naive Gauss Elimination

The method still consist of 2 steps :
• Forward elimination of unknowns
– Reduce the coefficient matrix [A] to an upper triangular system
– Eliminate x1 from the 2nd to the nth equation.
– Eliminate x2 from the 3rd to the nth equation
– Continue process until the nth equation has only 1 non-zero coefficient
• Back - substituted

21
 – From the last equation, we can find xn
– Substitute to (n − 1)th equation → find xn−1
– Continue until we find all the unknowns.
Example 1 : Solve the following system of linear equations by using Naive Gauss elimination


 4x1 + x2 − x3 = 3
2x1 + 7x2 + x3 = 19

x1 + −3x2 + 12x3 = 31


Solution :
Step 1 :
First, we transfer the system to a coefficient matrix :
 
4 1 −1 3
[A] =  2 7 1 19 
1 −3 12 31
Then, we apply the Naive Gauss elimination step by step :

Step 2 : Forward elimination.

We try to eliminate [A] to upper triangular matrix by using row operations.
−1
L2 + L1 → L2
2
−1
L3 + L1 → L3
4
1
L2 + L3 → L3
2
After step 2, [A] become :  
4 1 −1 3
[A] =  0 6.5 1.5 17.5 
0 0 13 39
As you can see, [A] ı́s in upper triangular form, and this matrix correspond to the system below :

4x1 + x2 − x3 = 3

6.5x2 + 1.5x3 = 17.5

13x3 = 39


→ x3 = 3
Step 3 : Back-substituted.
Substitute back x3 = 3 to the remaining equations :
 
4x1 + x2 = 6
 x 1 = 1

6.5x2 = 13 → x2 = 2
 
13x3 = 39 x3 = 3
 

You can practice naive Gauss elimination by solving some exercises :

Exercise 1 : Solve : 

 x1 + x2 − x3 = −3
6x1 + 2x2 + 2x3 = 2

−3x1 + 4x2 + x3 = 1


Answer : x1 = −0.25 x2 = −0.5 x3 = 2.25

Exercise 2 : Solve : 

 2x1 + x2 − x3 = 1
5x1 + 2x2 + 2x3 = −4

3x1 + x2 + x3 = 5


Answer : x1 = 14 x2 = −32 x3 = −5

22
6.4 Gauss elimination : Pivoting
• The reason why naive Gauss elimination is called ”naive” because during 2 main steps, division by zero can
occur. For example : 

 2x2 + 3x3 = 8
4x1 + 6x2 + 7x3 = −3

2x1 − 3x2 + 6x3 = 5


The first coefficient of the first equation is 0, so we cannot use it to turn other coefficients to 0 → method failed.
• In the process of calculation, you may have round off error when you divide by a number that is approximately
0.

• Determine the coefficient with the largest absolute value in the column below the pivot element
• The rows can then be switched so that the largest element is the pivot element. This is called partial pivoting

Example 1 : Solve the following system. 


 2x2 + 5x3 = 1
2x1 + x2 + x3 = 1

3x1 + x2 = 2


Step 1: We transfer to coefficient matrix :

 
0 2 5 1
[A] =  2 1 1 1 
3 1 0 2

Step 2 : We consider elements in the first column :

 
0
2
3

We can see that first element is 0, and the largest element is 3 at the last row → we need to exchange row 1 and row
3:
L3 ←→ L1
The matrix [A] become  
3 1 0 2
[A] =  2 1 1 1 
0 2 5 1
Step 3 : Now we use row operations as naive Gauss elimination
−2
L2 + L1 → L2
3
L3 + −6L2 → L3
A become upper triangular form :  
3 1 0 2
1 −1
[A] =  0 3 1.5 3

0 0 −1 3
Step 4 : The corresponding system of linear equations will be

3x1 + x2 = 2



1 −1

x2 + x3 =
3 3


−x3 = 3

23
→ x3 = −3
Step 5 : Substitute back x3 = −3 to the remaining equations:

3x1 + x2 = 2
 
x1 = −2

 
1 −1

x2 + −3 = → x2 = 8
 3 3 
x3 = −3

 
x3 = −3

Example 2 : Suppose that a team of three parachutists is connected by a weightless cord while free-falling at a velocity
of 5m/s. Calculate the tensions T and R in each section of cord and the acceleration a of the team (round to 4
decimal places), given the following:

Figure 6: Illustration of 3 parachutists

Solution : Formula of air resistance : Fdrag = c · v, c is drag coefficient, v is velocity.

24
 Figure 7: Forces acting on each parachutists

Because when the parachutists are falling down, there exists air resistance, so we have drag coefficients c in the
table below :

Parachutist Mass, kg Drag Coefficient, kg/s

1 70 10
2 60 14
3 40 17

Table 1: Mass and drag coefficients for each parachutist.

Now we apply Newton’s Second law to each parachutist, with positive direction is downward.
• Parachutist 1 : m1 g − T − c1 v = m1 a
• Parachutist 2 : m2 g + T − c2 v − R = m2 a
• Parachutist 3 : m3 g − c3 v + R = m3 a
Substitute m, c from Table 1, g = 9.8m/s2 , v = 5m/s, we have a system of linear equations
 

 70 · 9.8 − T − 10 · 5 = 70a 70a + T = 636

60 · 9.8 + T − 14 · 5 − R = 60a → 60a − T + R = 518
 
40 · 9.8 − 17 · 5 + R = 40a 40a − R = 307
 

The rest is up to you, use Gauss elimination Pivoting to solve this system of linear equations

Answer : a = 8.5941 T = 34.4118 R = 36.7647

You can also practice Gauss elimination by solving some exercises :

Exercise 1 : Solve 

 x2 + 2x3 = −1
3x1 + 2x2 + x3 = 1

x1 + x3 = 3


Answer : x1 = 2 x2 = −3 x3 = 1

Exercise 2 : Solve 

 2x2 + 1x3 = 4
x1 + x2 + x3 = 6

2x1 − x2 + x3 = 5


Answer : x1 = 5 x2 = 3 x3 = −2

25
Exercise 3 : Suppose that a team of three parachutists is connected by a weightless cord while free-falling at a
velocity v = 9m/s, g = 9.8m/s2 . Calculate the tensions T and R in each section of cord and the acceleration a of
the team (round to 4 decimal places), given the following:

Parachutist Mass, kg Drag Coefficient, kg/s

1 60 10
2 75 12
3 90 15

Table 2: Mass and drag coefficients for each parachutist.

Answer : a = 8.32 T = −1.2 R = 1.8

7 Iterative methods
7.1 Gauss - Seidel method
• The most commonly used iterative method for solving linear equations
• Employ initial guesses
• Iterates to obtain refined estimates of the solution
• Convergence can be checked using the criterion
→ round-off errors controlled by the number of iterations.

Principle
Suppose that we have a system of linear equations :

a11 x1 + a12 x2 + a13 x3 = b1

a21 x1 + a22 x2 + a23 x3 = b2

a31 x1 + a32 x2 + a33 x3 = b3


If the diagonal elements (a11 , a22 , a33 ) are all nonzero, then
• the first equation can be solved for x1
• the second for x2
• the third for x3
To be more specific :
b1 − a12 x2 − a13 x3


 x1 =



 a11
 b2 − a21 x1 − a23 x3
x2 =

 a22
b − a31 x1 − a32 x2


x 3 = 3


a33

Gauss-Seidel Iterative Formula

This method updates the values of x1 , x2 , x3 sequentially using:

(j) b1 − a12 xj−1

2 − a13 x3j−1
x1 =
a11
(j) b2 − a21 x1 − a23 xj−1
j
3
x2 =
a22
(j) b3 − a31 xj1 − a32 xj2
x3 =
a33

26
where:
• (j) represents the current iteration.
• (j − 1) represents the previous iteration.

A simple approach is to assume that the initial guesses are all zero. (x01 = x02 = x03 = 0). You can use other guesses
(if prior knowledge exists), but these are usually sufficient.

As each new x value is computed, it is immediately used in the next equation to determine another x value. For
example, xj1 is used immediately to calculate xj2 .

Stopping Criterion
Convergence can be checked for each iteration using the criterion:

xji − xj−1
i
εa,i = × 100% ≤ εs
xji

where εs is a prespecified tolerance for all unknowns (x1 , x2 , x3 )

Converge condition
If the following condition holds, Gauss-Seidel will converge:
n
X
|aii | > |aij |
j=1
j̸=i

Example: Consider the linear system: 

4x1 + x2 + 2x3 = 4

3x1 + 5x2 + x3 = 7

x1 + x2 + 3x3 = 3


The coefficient matrix is:

 
4 1 2
A = 3 5 1
1 1 3
To check whether the Gauss-Seidel method will converge, verify if the matrix is strictly diagonally dominant. That is,
for each row i, check:
n
X
|aii | > |aij |
j=1
j̸=i

Row 1: |4| > |1| + |2| (true)

Row 2: |5| > |3| + |1| (true)
Row 3: |3| > |1| + |1| (true)
Since all rows satisfy the condition, the Gauss-Seidel method will converge for this system.

Example
Use the Gauss-Seidel method to obtain the solution for

 3x1 − 0.1x2 − 0.2x3 = 7.85

0.1x1 + 7x2 − 0.3x3 = −19.3

0.3x1 − 0.2x2 + 10x3 = 71.4


27
Solution
Step 1: Solve each equations for its unknown on the diagonal:
7.85 + 0.1x2 + 0.2x3

x 1 =



 3
−19.3 − 0.1x1 + 0.3x3

x2 =

 7
71.4 − 0.3x 1 + 0.2x2


x 3 =

10
Step 2: Assume x2 = 0 and x3 = 0, we can calculate:

7.85 + 0.1(0) + 0.2(0)

x1 = = 2.616667
3
Step 3: Using calculated x1 = −2.616667 and x3 = 0, we can calculate:

−19.3 − 0.1(2.616667) + 0.3(0)

x2 = = −2.794524
7
Step 3: Using calculated x1 and x2 , we can calculate:

71.4 − 0.3(2.616667) + 0.2(−2.794524)

x3 = = 7.005610
10
=⇒ End of iteration 1. For this iteration, all approximate errors will be 100%. Hence we continue with another
iteration.
Iteration 2:
7.85 + 0.1(−2.794524) + 0.2(7.005610)
x1 = = 2.990557
3
−19.3 − 0.1(2.990557) + 0.3(7.005610)
x2 = = −2.499625
7
71.4 − 0.3(2.990557) + 0.2(−2.499625)
x3 = = 7.000291
10
Errors of this iteration:
2.990557 − 2.616667
|εa,1 | = × 100% = 12.5%
2.990557
−2.499625 + 2.794524
|εa,2 | = × 100% = 11.8%
−2.499625
7.000291 − 7.005610
|εa,3 | = × 100% = 0.076%
7.000291
Depending on εs , we can stop (if all |εa,i | < εs ) or continue with other iterations.

7.2 Jacobi method

• An algorithm for determining the solutions of a diagonally dominant system of linear equations

• Transforms a matrix to a diagonal matrix by eliminating off-diagonal terms in a systematic fashion

• Requires an infinite number of operations
• The approach is iterative in that it is repeated until the off-diagonal terms are ”sufficiently” small

• An alternative approach, with Gauss-Seidel, utilizes a somewhat different tactic

28
Jacobi vs Gauss Seidel

As you can see, Jacobi method will calculate all the values of x1 , x2 , x3 in the first iteration, then use those values
for the second iteration.

This is different from the Gauss-Seidel method using a new value of x immediately to find another x value.

Example
Solve the following system of linear equations by Jacobi method :
10x1 − x2 + 2x3

 =6

−x + 11x − x + 3x

= 25
1 2 3 4

 2x 1 − x 2 + 10x 3 − x 4 = −11

3x2 − x3 + 8x4

= 15
Solution :
Step 1 : Solve each equations for its unknown on the diagonal:
6 − 2x3 + x2


 x1 =


 10
25 − 3x 4 + x3 + x1


x 2 =


11
x 3 =
 −11 + x 4 + x2 − 2x1



 10
x4 = 15 − 3x2 + x3



8
Step 2 : We choose x1 = x2 = x3 = x4 = 0 as intial approximation, then the solution of the first iteration are :
6−0+0


 x1 = = 0.6


 10
25 − 0 + 0 + 0


x 2 =


11
−11 + 0+0−0
x3 =





 10
x4 = 15 − 0 + 0 = 1.875



8
Step 3 : The iteration will continue until the desired accuracy has been reached. Here are the approximated solutions
after five iterations

29
 Iteration x1 x2 x3 x4
1 0.6 2.27272 -1.1 1.87500
2 1.04727 1.71590 -0.80522 0.88522
3 0.93263 2.05330 -1.0493 1.13088
4 1.01519 1.95369 -0.9681 0.97384
5 0.98899 2.0114 -1.0102 1.02135

7.3 Improvement of convergence using relaxation

Formula
xnew
i = λxnew
i + (1 − λ)xold
i

Types of relaxation
Types I:
• If λ is between 0 and 1 → this type of modification is underrelaxation.
• It is used to make a nonconvergent system converge or to hasten convergence.
Types II :
• If λ is between 1 and 2 → this type of modification is overrelaxation.
• accelerate the convergence of an already convergent system
• It is also called as succesive or simultaneous overrelaxation.

Example
Solve the following system with Gauss Seidel using overrelaxation (λ = 1.2)
(
−3x1 + 12x2 = 9
10x1 − 2x2 = 8

Solution :
Step 1 : We rearrange the equations so that they are diagonally dominant :
(
10x1 − 2x2 = 8
−3x1 + 12x2 = 9

Step 2 : Solve the first equation for x1 , the second equation for x2 .
8 + 2x2

x 1 =
 = 0.8 + 0.2x2
10
x2 = 9 + 3x1 = 0.75 + 0.25x1

12
Step 3 : First iteration : using initial guesses x1 = x2 = 0, we can solve for x1
x1 = 0.8 + 0.2(0) = 0.8
Step 4 : Before solving x2 , we first apply relaxation for x1 :
x1,r = 1.2(0.8) − 0.2(0) = 0.96
r means it is ”relaxed” value.
Step 5 : This result is then used to compute x2 :
x2 = 0.75 + 0.25(0.96) = 0.99
Step 6 : Before moving to the next iteration, we then apply relaxation for x2 :
x2,r = 1.2(0.99) − 0.2(0) = 1.188
The procedure continue until it reach the desired accuracy.

30