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