ENGR20005: Numerical Methods in Engineering
Semester 1, 2020 / Assignment 1
Marks : This assignment is worth 25% of the overall assessment for this course.
Due Date : Friday, 17th April 2020, 23:59, via Canvas submission (see section 5). You will lose
penalty marks at the rate of 10% per day or part day late, and the assignment will not be marked
if it is more than 5 days late.
1 Learning Outcomes
In this assignment, you will demonstrate your understanding of solving engineering problems
using numerical computations and assessing particular algorithms. The objectives of this
assignment are to program algorithms for root finding, solving systems of linear algebraic
equations and implementing an interpolation method.
2 Root finding [Σ 10 Marks]
The Newton Raphson method offers several advantages over bracketing methods such as the
Bisection method. It only requires a single initial value and it avoids the limitation of Bisection,
which cannot work when the function at the two initial values has equal sign. Additionally, it
converges faster than bracketing methods. However, there are two main disadvantages of using
the Newton Raphson method: the derivative of the function needs to be determined analytically
and it does not guarantee convergence.
The need for calculating the derivative of the function makes this method more difficult to be
automated in a programming code. It reduces its flexibility and there are many numerical
applications in which the analytical calculation of the derivative is not feasible. A way to remove
this limitation is to use a variation of the method, the Secant method, which can be thought of as
a finite-difference approximation of the Newton Raphson method. The Secant method requires
two initial values to start iterating, but in contrast to the Bisection method, these two initial
values do not need to bracket the root.
In summary, the Secant method can be considered an open method and is very flexible, as it does
not require any mathematical manipulation of the input. However, its convergence is not
guaranteed and it will diverge depending on the data analysed. The Bisection method is much
more robust, and as far as the two inserted values bracket the root, it will be stable, although the
convergence of the method is slower.
2.1 Function analysis with the Secant method [2 Marks]
The equation:
f (x)= x5−200x3+70 (1)
has three roots sitting at approximately x=−14.14, x= 0.70 and x= 14.14. The behaviour of this
function near the roots is shown in Figure 1.
a. Your task is to twice calculate by hand the first four iterations of a Secant method root
finding algorithm for this polynomial. Consider the following pairs of initial values [x1, x2]:
[-100, 100] and [-1000, 100].
© 2020 The University of Melbourne 1
Figure 1: Polynomial in the surroundings of the three roots
b. Discuss the limitations of the Secant method. Discuss how the selection of the two initial
values influences the convergence of the Secant method.
Note: The values of y increase rapidly due to the high order of the polynomial. You may round
the results in your report but do not round the results during your calculations. Use
MATLAB to perform the calculations if necessary.
Hint: It may be helpful (not mandatory) to plot the function and the secant lines calculated in
the consecutive iterations. You can use the following MATLAB commands:
• fplot(myfun, [xmin xmax]) to plot your function between two x-values.
• line([x0,x1],[y0,y1]) to plot the secant line between points (x0, y0) and (x1, y1).
• axis([xmin xmax ymin y max])if you want to crop the axes.
2.2 Development of Bisection method [1 Mark]
Develop a MATLAB function that takes as an input an anonymous function, two initial x-values
to estimate the root, a maximum number of iterations and a tolerance to determine when
convergence is reached. The user will define these inputs in a different script or in the command
window and then call the function in the following format:
• Bisection(x1, x2, MyFun, Maxiter, Tolerance);
The code should find the root using the Bisection method and must include the following
functionalities:
• If the two x-values inserted by the user result in an unsuccessful sign check, the code
should ask for a new pair of values until valid ones are inserted. You may need to consider
what happens if the user inserts directly an x-value corresponding to the one of the roots.
© 2020 The University of Melbourne 2
• The code will run until a root within the specified tolerance is found or the maximum
number of iterations is exhausted.
• The code will save the result of every iteration in a .txt file. Use the guidelines contained
in the .m skeleton file to create the text file.
• The code will not require any user interaction (with the exception of re-inserting the initial
x-values until a valid pair of boundaries is inserted).
2.3 Overview of Brent-Dekker method
The Brent-Dekker method improves the Secant method to make it more stable by combining it
with the Bisection method. This combination allows the Brent-Dekker method to converge for
any continuous function. For each iteration of the algorithm, the two methods are evaluated. If
the estimation of the root obtained by the Secant method lies between the last value calculated
by the algorithm and the new estimation obtained by the Bisection method, the value of the
Secant method is used as one of the boundaries for the next iteration. If not, it may indicate that
the Secant method is diverging and the Bisection iteration value is used as a boundary for the
next iteration. The Bisection method will need a sign check for the two initial values such that
the function values of both have the opposite sign, and Brent-Dekker algorithm needs to include
this check as well, to make sure that Bisection will be a valid method in the next iteration.
A representative work flow of the method is given as such; Consider that a previous iteration
ended with two estimations of the root being A and B. In the previous iteration, B was considered
a better estimation of the root than A because (abs(f(B)) < abs(f(A))). In the first iteration,
these values are defined by the user. The process to find the next estimated value can be
described as follows:
• Make a new estimation of the root using the Bisection method.
• Make a new estimation of the root using the Secant method.
• Compare if the Secant estimation lies between the new Bisection estimation and the
previous best estimation B. If so, the new estimation will be the one calculated by the
Secant method, and the Bisection’s will be discarded. If not, the new value will be the one
obtained by Bisection, and the Secant estimation will be discarded.
• Check for convergence at the new calculated value. If reached, end the algorithm.
• Check if f(new value) and f(B) have different signs. If so, they will be the two estimations
for the next iteration. If not, the values would be A and the newly calculated value.
• Refresh A and B values so that (abs(f(B)) < abs(f(A))).
• Start next iteration.
Figure 2 summarises the structure of the Brent-Dekker algorithm in a flow chart.
2.4 Development of Brent-Dekker method [5 Marks]
Develop a MATLAB function for root finding using the Brent-Dekker method. The call to the
function should be as follows:
• Dekker(x1, x2, MyFun, Maxiter, Tolerance);
© 2020 The University of Melbourne 3
Figure 2: Flow chart of a Brent-Dekker’s algorithm
The code should find the root using a combination of the Bisection and Secant method as
described in section 2.3 and must include the same functionalities as the Bisection code that you
have developed in section 2.2.
Note: Both implemented functions (Bisection and Dekker) will be marked using a random
function with different initial values, tolerances and maximum numbers of iterations. So it is a
good idea that you test and ensure robustness for different parameter inputs. The purpose of this
exercise is to give you an idea on how to develop codes with varying inputs. Please provide
separate function files for each method. Skeleton codes are provided with the files Bisection.m
© 2020 The University of Melbourne 4
and Dekker.m for both functions and therefore, the function name and the defined variables have
to match the convention of the skeleton codes.
2.5 Formatting of the output files
In the skeleton .m files you have some guidance in how to save the results of every iteration,
including the formatting of the floating numbers. You may play with the formatting of the
numbers while you are working, but remember to maintain the original formatting for the file
submission. Your output file should have a structure such as the following:
Iteration // Current estimation // Method used
0 100.0000000000 Initial estimation
1 1.0000000000 Secant method
2 1.0000000000 Secant method
3 0.4536997949 Secant method
4 -0.1268662164 Bisection method
5 0.0042974229 Secant method
6 -0.0000110026 Bisection method
7 0.0000000000 Secant method
7 0.0000000000 Final estimation
Convergence reached
Note: This is not the real solution.
2.6 Discussion of results [2 Marks]
Discuss your findings by comparing the Bisection method and the Brent-Dekker method. Reflect
on how different input parameters affect the convergence speed of both methods and compare
your solutions to the built-in MATLAB functions roots and fzero.
© 2020 The University of Melbourne 5
3 Linear Algebraic Systems [Σ 10 Marks]
In this task, we will be looking to solve the two dimensional Poisson equation using both direct
and iterative methods. A form of this equation can be found in the steady state heat equation,
from which the time derivative term is not taken into account. For this particular example, we
will consider the two dimensional steady state heat equation with a heat source represented by ψ
(inhomogeneous case). This equation is mathematically expressed as follows:
∂2φ
∂x2
+ ∂
2φ
∂y2
+ψ= 0 (2)
where ψ= 10 at all grid points. The domain is given by x ∈ [0,1] and y ∈ [0,1]. The boundary
conditions are given as follows:
φ(0, y)= 1
φ(1, y)= 1
φ(x,0)= 1
φ(x,1)= 1
(3)
A schematic representation of this problem is given as follows:
-0.2 0 0.2 0.4 0.6 0.8 1 1.2
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
φ
(0,
y)
= 1
φ
(1,
y)
= 1
φ(x,0) = 1
φ(x,1) = 1 Boundary Points
Figure 3: Illustration of the two dimensional grid
Equation (2) can be discretized using a second order central difference method on a uniformly
spaced grid (∆x=∆y=∆) as follows:
φi−1, j−2φi, j+φi+1, j
∆2
+ φi, j−1−2φi, j+φi, j+1
∆2
+ψi, j = 0. (4)
© 2020 The University of Melbourne 6
The resulting discretized equation can be written in the form of a system of equation:
[A][φ]= [C] (5)
Show that for a grid size of (Nx x Ny) of (5 x 5), [A] and [C] are given by [2 Marks]:
4 −1 0 −1 0 0 0 0 0
−1 4 −1 0 −1 0 0 0 0
0 −1 4 0 0 −1 0 0 0
−1 0 0 4 −1 0 −1 0 0
0 −1 0 −1 4 −1 0 −1 0
0 0 −1 0 −1 4 0 0 −1
0 0 0 −1 0 0 4 −1 0
0 0 0 0 −1 0 −1 4 −1
0 0 0 0 0 −1 0 −1 4
︸ ︷︷ ︸
[A]
φ2,2
φ3,2
φ4,2
φ2,3
φ3,3
φ4,3
φ2,4
φ3,4
φ4,4
=
∆2ψ2,2+φ2,1+φ1,2
∆2ψ3,2+φ3,1
∆2ψ4,2+φ5,2+φ4,1
∆2ψ2,3+φ1,3
∆2ψ3,3
∆2ψ4,3+φ5,3
∆2ψ2,4+φ1,4+φ2,5
∆2ψ3,4+φ3,5
∆2ψ4,4+φ5,4+φ4,5
︸ ︷︷ ︸
[C]
(6)
3.1 Direct Method [4 Marks]
Apply Crout’s LU decomposition method to solve the system of equation, [A][X] = [C] where [A] =
[L][U]. Write a MATLAB code to solve the questions below.
a. Find the [L] and [U] matrix. Check that [A] = [L][U].
b. Find [R] from [L][R] = [C] using forward substitution.
c. Find [X] from [U][X] = [R] using backward substitution.
d. Compare your solution with the built-in MATLAB matrix solver linsolve(A,C).
3.2 Iterative Methods [4 Marks]
Apply the Point Jacobi method and the Gauss Seidel method to solve the system of equations.
Plot the sum of error (|Err| =
√∑n
k=1 |Xk−X exact,k|2) of the iterated solution with respect to the
exact solution given by the direct method as a function of the number of iterations for the
different iterative methods. Discuss your results.
a. Implement the Point Jacobi Method in the given skeleton code (PJ_solve.m).
b. Implement the Gauss Seidel Method in the given skeleton code (GS_solve.m).
c. Compare the plot of error versus number of iterations required between the Point Jacobi
and Gauss Seidel Method. Set the tolerance to 1e-5.
Your implementation should work for different problem sizes as well. Please provide separate
function files for each method. A skeleton code is provided for each function and therefore the
function name and variables defined for the function have to be the same as the skeleton code.
The function is to be run on a main script (Question3_1.m, Question3_2.m) as provided. Submit
separate MATLAB files (.m) for questions 3.1 and 3.2.
© 2020 The University of Melbourne 7
4 Interpolation [Σ 10 Marks]
In this task, we are looking into the derivation of spline interpolation, implement a spline
interpolation function in MATLAB and apply the implemented function to plot the shapes of
compressor blades. In spline interpolation, polynomial functions Si(x) of a certain degree are
defined for each interval of two data points
S(x)= Si(x) for xi < x< xi+1 . (7)
The general formula for a cubic spline, which consists of third order polynomials, is
Si(x)= ai+bi(x− xi)+ ci(x− xi)2+di(x− xi)3 . (8)
4.1 Theory [2 Marks]
State in your report the six conditions to determine the natural cubic spline parameters ai, bi, ci
and di based on the function values f (xi) at the data points xi for i = 1 . . .N. Use those conditions
to derive the following equations (9) to (12), where the step size is simplified to hi = xi+1− xi.
Make sure that you display intermediate steps.
ai = f (xi) (9)
bi = 1hi (
ai+1−ai)− hi3 (2ci+ ci+1) (10)
h j−1c j−1+2
(
h j+h j−1
)
c j+h jc j+1 = 3h j
(
a j+1−a j
)+ 3
h j−1
(
a j−1−a j
)
(11)
di = ci+1− ci3hi
(12)
4.2 Implementation [4.5 Marks]
Use the skeleton file nat_cub_spline.m to implement natural cubic spline interpolation in
MATLAB. Your code should implement the function nat_cub_spline(x_itps, data_file),
which
• reads in a comma-separated value (.csv) file,
• determines the parameters ai, bi, ci and di,
• interpolates at the data points x_itps,
• sorts the interpolated values for each data point from small to large,
• returns the array res containing the interpolated values.
The required data type and format of the variables x_itps, data_file and res are defined in the
header of the skeleton file. Your code should have general validity and return interpolation
values at each data point that is specified in x_itps.
Note: One data point from parameter x_itps may occur in multiple spline interpolation
intervals.
© 2020 The University of Melbourne 8
4.3 Application [3.5 Marks]
The performance of an aircraft engine is significantly influenced by the design of its compressor
blades. The National Advisory Committee on Aeronautics (NACA) developed multiple airfoil
profile series from the 1930s to the 1950s, which can still be found on aircrafts in operation today.
Mean camber line
Chord line
Figure 4: Schematic of a four-digit NACA profile
The NACA four-digit series, e.g. NACA6315, describes the airfoil profile based on three
parameters:
• the maximum camber as percentage of the chord length (first digit, e.g. 6%),
• the distance of the maximum camber from the airfoil leading edge in tenth of the chord
length (second digit, e.g. 30%),
• the maximum airfoil thickness as percent of the chord length (last two digits, e.g. 15%).
The schematic of a four-digit NACA profile can be found in Figure 4, where the influence of the
colored parameters above is illustrated. The original data points of the NACA profiles
NACA6309, NACA6315 and NACA6321 can be found in the files naca6309.csv, naca6315.csv
and naca6321.csv, respectively.
a. Provide one plot for each NACA profile in your report. In each plot, the following things
should be considered: plot the original data points as red crosses. Use the developed
nat_cub_spline function to interpolate the NACA profile data at 100 sampling points in
each interval. Each sampling point is a value between two original data points. Plot the
interpolated data as a continuous curve in blue. Make sure that the axes are equally
spaced. How does the interpolated function look compared to the original data points?
b. To improve the accuracy of the NACA profile interpolation, we increase the number of
original data points and split the airfoil into three parts: the leading edge (LE) is the front
part of the airfoil, the suction side (SS) is the upper surface and the pressure side (PS) is the
lower surface. Your task is to create one plot of the NACA6309 profile, based on the
provided data points in the files naca6309_LE.csv, naca6309_SS.csv and
naca6309_PS.csv. Use your nat_cub_spline(x_itps, data_file) function to interpolate
between the original data points for the suction side and pressure side. To plot the leading
edge, the order of the spline interpolation is reduced to a quadratic spline. Use the
quad_spline(x_itps, data_file) function in the provided .p file for the leading edge.
Plot the airfoil in the same way as in part a. What differences to the NACA6309 plot from
the previous part can you observe?
c. Compare your nat_cub_spline function to the build-in MATLAB interp1 function. Load
the data of the suction side of the NACA6321 profile from the file naca6321_SS.csv. Create
three plots, in which you compare the original data points and the interpolated data from
© 2020 The University of Melbourne 9
the nat_cub_spline function to the options nearest, linear and spline of the build-in
interp1 function, respectively. Discuss the performance of the individual interp1methods.
Please provide the function file nat_cub_spline.m and separate main .m files which use the
function to generate the plots for parts a, b and c of task 4.3.
5 Instructions for submission
Please submit your assignment via file upload in Canvas. You should only submit one .zip file,
which is named ENGR20005_A1_.zip (replace with your personal
student ID). Your submission file should contain the following files:
• Report file ENGR20005_A1_Report_.pdf
• Function file Bisection.m
• Function file Dekker.m
• Result file Bisection.txt
• Result file Dekker.txt
• MATLAB .m file for question 2.6
• MATLAB file Question3_1.m
• Function file GS_solve.m
• Function file PJ_solve.m
• MATLAB file Question3_2.m
• Function file nat_cub_spline.m
• MATLAB .m file for question 4.3a
• MATLAB .m file for question 4.3b
• MATLAB .m file for question 4.3c
Note: All of the above files have to include your student ID.
6 Getting help
There are several ways for you to seek help with this assignment. Please have a look at the
MATLAB livescripts from the workshops to get an idea on how to develop the function files
needed for this assignment. Please do not post any source code on the Canvas discussion board.
You may ask questions during the workshops or send an email to Dr. Leon Chan
() directly.
Note: Students seeking extension for medical or other reasons outside of your control should
email Dr. Leon Chan as soon as possible.
© 2020 The University of Melbourne 10