ENGR30003 – Numerical Programming for Engineers – Semester 2, 2020
1 Root finding [9/35 Marks]
Imagine a wedge-shaped object flying through air at supersonic speeds (see Fig. 1). In compressible flow
theory (covered in MCEN90008 Fluid Dynamics), it is known that an oblique shock wave forms at the
front tip of this object under certain conditions.
C program to solve shock-wave equation
Write a C code that solves Eq. (2) using the Newton–Raphson method to find the root of f(β), regarding
θ and M as parameters and solving for β.
(a) Write your C program such that it uses the Newton–Raphson method to solve f(β) = 0. What
values of βL and βU do you think might be appropriate to use as an initial guess? For M = 4.5 and
θ = 26◦
, you should find that βL = 38.17378…
◦ and βU = 81.9852…
(b) Use your computer program to solve f(β) = 0 (i.e. find βU and βL values) for M = 2.5, 3.5, 4.5,
5.5, 6.5, 7.5. Remember that for θ = 0, βL = arcsin
and βU = 90◦
. Plot β vs θ for all the
M’s. That is, plot values of θ on the horizontal axis and corresponding values of β on the vertical
axis. This plot is called the θ − β − M diagram and you will find it very useful if you decide to
do MCEN90008 Fluid Dynamics in the future. Your solution to this part of the assignment should
look like Fig. 2 which shows the solution for M = 5. You can plot your results obtained from your
C code with your program of choice, e.g. Matlab, Excel, etc.
The implementation for this task is to be done in the function named shockwave(). The input parameters
are to be read in from the input file in_shock.csv.
The first two lines correspond to the part (a) where for M = 4.5 and θ = 260
, you need to compute the
values of βL and βU . You will notice the initial values set here are 0
0 apiece so that you can modify
them to find the right initial guess to compute the angles. The next set of lines corresponds to part (b)
where you will evaluate for different M, the values of βL and βU for different values of θ (ranging from
0 up to θmax at an increment of 1
). Outputting the results of this task are to be done only for part
(b) into a file named out_shock.csv in the format as shown below. This example only shows part of the
results for θ = 00
for the set of Mach numbers chosen (you can use this to validate your code).
You must dynamically allocate space for the Mach numbers so that your implementation may work for
a different set of M values (we will evaluate your code for a different set of M values!). For each Mach
number, upon increasing θ by 1
, you will reach the maximum θ beyond which the solutions of β are not
physically relevant. You will write to file only up to this maximum θ per Mach number.
2 Regression [4/35 marks]
Here, an alternative way of solving a Least Squares Problem is considered.
3 Linear Algebraic Systems [5/35 marks]
Consider the following tri-diagonal system:
(a) Using Gauss elimination, show that the above matrix can be rewritten as
The algorithm outlined above is called the Thomas algorithm. When applicable, it is a very efficient
method for solving linear tri-diagonal systems.
(b) Write a C code using the Thomas algorithm to solve the tri-diagonal system shown in Eq. 6. Since
the tri-diagonal system is a banded matrix, you need not store all the zeros! Instead, your code
should take as an input the vectors ai
, ci and ri