ENGR30003辅导、讲解C/C++程序、讲解Matlab编程、C/Matlab讲解

ENGR30003: Numerical Programming for Engineers

Semester 2, 2018 / Assignment 2

Marks : This assignment is worth 35% of the overall assessment for this course.

Due Date : Monday, 15th October 2018, 5:00pm, via ① submission. You will lose

penalty marks at the rate of 10% per day or part day late, and the assignment will no 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, performing least-squares approximations and interpolations, regressing solutions and

solving a differential equation using different integration and differentiation schemes.

2 Rootfinding [9/35 marks]

M

q<qmax

b

q>qmax

M

Figure 1: Difference in shock wave type for different wedge angles.

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.

The equation relating the wedge half-angle θ to the shock oblique angle, β, and the Mach

number, M is given by

tan (θ) = 2cot (β)

M2

sin2

(β)?1

M2(γ+cos (2β))+2

, (1)

where γ = 1.4 is the ratio of specific heats. For a given M, as you keep increasing θ, there is a

critical angle, θmax where the shock becomes detached. We can also recast Eq. (1) into the

following form

f (β) = 2cot (β)

M2

sin2

(β)?1

M2(γ+cos (2β))+2

tan (θ) . (2)

Your tasks are the following:

2018 The University of Melbourne 1

3

5 10 15 20 25 30 35 40 45

-

0 10 20 30 40 50 60 70 80

90 Figure 2: θ ?β? M diagram for M = 5.0.

2.1 Analytical solution for θ = 0

[1 mark]

For θ = 0, i.e. the object would be a flat plate, show that two possible solutions for β are

βL = arcsinμ

1M?, βU = 90?

(3)

βU and βL are usually called the strong shock and weak shock solutions, respectively. Even

though we can mathematically obtain two possible solutions, in reality, or physically, only the

weak shock solution occurs. Note also that in order for the solution to be physically realisable,

θ < β < 90?

2.2 Graphical solution [2 mark]

Plot f (β) vs β for

a. M = 1.5 and θ = 5

, 10? and 15?

b. M = 5.0 and θ = 20?

, 30? and 45?

Indicate how βU and βL change with θ and M. Can you identify from your plots the approximate

value of θmax?

PLEASE REMEMBER to change the angle from degrees to radians when you use sinusoidal

functions in your computer program.

2.3 C program to solve shock-wave equation [6 marks]

Your task is to 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 = 5.0 and θ = 20?

, you should find that

βL = 29.80092...

, and βU = 84.55625...

(b) Extend your C program to find βU and βL values for M = 5.0 and 0 ≤ θ ≤ θmax. Remember

that for θ = 0, βL = arcsin?1M￠

and βU = 90?

. Plot values of θ on the horizontal axis and

corresponding values of β on the vertical axis. Your solution to this part of the assignment

should look like Fig. 2. Note that you can plot your results obtained from your C code with

your program of choice, e.g. Matlab, Excel, etc.

2018 The University of Melbourne 2

(c) Use your computer program to solve f (β) = 0 for M = 2.0, 3.0, 4.0, 6.0, 7.0, 8.0. Plot β vs θ

for all the M’s. 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.

The implementation for this task is to be done inwhere input parameters are to

be read in from. The input file consists of three parts, as shown below:

M, theta , be ta_ l , beta_u ,gamma

5. 0 , 2 0. 0 , 0. 0 , 0. 0 , 1. 4

M

5.0

M

2.0

3.0

4.0

6.0

7.0

8.0

The first two lines corresponds to the part (a) where for M = 5.0 and θ = 200

, you need to compute

the values of βL and βU . You will notice the initial values set here are 00 apiece so that you can

modify them to find the right initial guess to compute the angles. The next line corresponds to

part (b) where you will for M = 5 evaluate the values of βL and βU for different values of θ

(increments of 10

from 00 up to θmax). The next set of lines corresponds to part (c) where you will

evaluate for different M, the values of βL and βU for different values of θ (increments of 10

from

0 up to θmax). Outputting the results of this task are to be done only for part (c) into in the format as shown below. This example only shows part of the results for

θ = 0

for the set of Mach numbers chosen (you can use this to validate your code). The output

of M, βU ,βL is to be done up to 6 decimal places while θ as an integer:

M, theta , beta_ lower , beta_upper

3.000000 ,0 ,19.471221 ,90.000000

3.000000 ,1 ,20.157561 ,89.649912

4.000000 ,0 ,14.477513 ,90.000000

4.000000 ,1 ,15.131265 ,89.719941

You must dynamically allocate space for the Mach numbers so that your implementation may

work for a different set of M. For each Mach number, upon increasing θ by 10

, 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.

2018 The University of Melbourne 3

3 Regression [4/35 marks]

Here, an alternative way of solving a Least Squares Problem is considered. Recall from the

lecture that??

(4)

is the system of equations for linear regression y? = ax+ b.

Show that:

b = y? ax, a =PN

i=1

(xi ? x)(yi ? y)

PN

i=1

(xi ? x)2

(5)

Hint: x =PN

i=1

xi/N is the mean of xi.

? Hint: You know the inverse of a 2×2 matrix!

By considering Eq. 5, when will linear regression fail to find a solution?

4 Linear Algebraic Systems [5/35 marks]

Consider the following tri-diagonal system:?

a1 b1 0 0 ... 0

c2 a2 b2 0 ... 0

0 c3 a3 b3 ...

0 ... 0 cN?1 aN?1 bN?1

0 ... 0 0 cN aN?????????

(6)

Using Gauss elimination, show that the above matrix can be rewritten as、

(7)

where

for i = 2,3,..., N

2018 The University of Melbourne 4

Thus the solution to the original tri-diagonal matrix can be written as

xi =Qi/ai

for i = N

(Qi bixi+1)/ai

for i = N 1, N ?2,...,1

The algorithm outlined above is called the Thomas algorithm. It is a very efficient method for

solving linear tri-diagonal systems.

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

, bi

, ci and Qi

. The output from your C code should be

the solution vector xi.

Use your code to solve the following tri-diagonal system?

1 ?2 0 0 0 0

2 4 5 0 0 0

0 8 ?9 2 0 0

0 0 6 3 4 0

0 0 0 6 3 2

0 0 0 0 1 3

(8)

You will implement the code for this task in the function , where you will read in

as an input the vectors ai

, bi

, ci and Qi from the file②. The output from your

implementation should be the solution vector xi

, written out to(up to 6

decimal places). Your implementation should allocate dynamically the space for the values of ai

,

bi

, ci and Qi such that your implementation would work for different problem sizes as well.

5 Interpolation [4/35 marks]

Given the following data

x0 = 0 f (x0) = 0

x1 = 1 f (x1) = 11

x2 = 4 f (x2) = 19

x3 = 9 f (x3) = 19

Write C code that uses both second-order Lagrange interpolating polynomials and cubic splines

in order to estimate f (x). Estimate the value of f (x = 5) using both methods. You will implement

both methods in , reading in from the file the values of x and f (x) and

outputting the values of the interpolated value to in the following format (up to

6 decimal places). The values shown are just an example and do not represent the correct values.

lagrange

11.145601

cub ic

23.098532

You must dynamically allocate space for x, f (x) and write your code such that it can work for a

different set of input and different xo. As an example, for the cubic splines, your code should be

capable of identifying the interval where xo lies and compute f (xo), where xo is 5 for this task but

will be different during submission. Plot the interpolated function from both methods using

either MATLAB or Excel. What conclusions can you draw from this plot?

2018 The University of Melbourne 5

6 Differentiation, differential equations [13/35 marks]

Write a C program that solves the modified wave equation

(9)

on the interval x = [0;1], and for times 0 ≤ t ≤ 10, using the phase speed c = 1.

The wave equation is to be solved on the spatial interval 0 ≤ x ≤ 1, that is discretized using

Nx +1 points xi

, with i = 0,1,2,..., Nx. The equidistant grid spacing therefore is x = 1/Nx.

To solve the equations numerically, the modified wave equation is discretized using the

following nomenclature. n is the time index so that fn

i

is the function value at time level n at

point xi and fn+1

i

is the function value at time level n+1 at point xi

. ?t is the numerical

timestep for the time integration.

You have to consider the two different approaches:

a. Euler explicit with upwinding:

This means that for the spatial derivative of f ,

, an upwinding scheme needs to be codedup

for i = 1,..., Nx

For i = 0, use the boundary stencil

which implies a so-called ‘periodic’ boundary condition. This means that when the

propagating wave reaches the outlet boundary of the domain, at x = 1, it enters at the

inflow again, at x = 0.

b. Lax–Wendroff Method:

(11)

This method requires evaluation of both the first and second derivatives of

Both derivatives are approximated with central difference schemes, with the first

derivative calculated as

for i = 1,..., Nx ?1

Again, periodic boundary conditions are to be used, so that for i = 0 the boundary stencil is

2018 The University of Melbourne 6

and for i = Nx the boundary stencil is

The second derivative is computed as 2

for i = 1,..., Nx 1

Again, periodic boundary conditions are to be used, so that for i = 0 the boundary stencil is

and for i = Nx the boundary stencil is ?2

fx2≈fn

02fn

Nx+ fn

Nx1

x2

The initial condition for f is

0 ≤ x < 0.125 → f (x,t = 0) = 0

0.125 ≤ x ≤ 0.375 → f (x,t = 0) = 0.5[1?cos{8π(x?0.125)}]

0.375 < x ≤ 1 → f (x,t = 0) = 0.

6.1 Suggested Code Structure

In the following possible steps for coding up the methods are suggested.

a. Set up main program that does the following

Allocate arrays for the functions fni

, fn+1

, etc, where i = 0,1,2,..., Nx (keep Nx a

parameter so that you can change it).

Write the initial condition into the array f

, and write to a file for later visualization.

Set up loop over the number of timesteps you want to run (each with ?t) - within this

loop you want to call your time integration routine.

Write intermediate solutions to file for later visualization at several time instances.

b. Code up your time integration routine:

if using the Euler upwind method, call a routine that computes ?f?x

based on the

function values at time level n.

Then evaluate the solution at the next time level n+1 by using Eq. 10.

if using the Lax–Wendroff method, call a routine that computes both ?f

and 2

fx

2 based

on the function values at time level n.

Then evaluate the solution at the next time level n+1 by using Eq. 11.

c. Code up your routines for computation of the spatial derivatives fx

and 2

fx2

2018 The University of Melbourne 7

You will implement your routines in the function ， reading in the values of c, Nx and

CFL from . This will allow you to compile your code just once and run it for

different values of c, Nx and CFL by just changing the infile. For code assessment, you will write

out the solutions from your Euler upwind and Lax–Wendroff versions to

andrespectively, the k

th time the time loop executes (this depends on

what Nx,CFL are) in the following format (up to 6 decimal places).

x , f ( x )

0.000000 ,0.000000

0.010000 ,0.000394

0.020000 ,0.000907

The k

th time will be read in from the infile from the parameter ??t??t?r.

6.2 Tasks

Run your code until a time of t = 10 for both numerical methods (Euler upwind and

Lax–Wendroff) for the following cases:

Chose two different resolutions ?x, by setting Nx to 80 or 200.

For each resolution, chose the timestep ?t from the CFL number, defined as ct

x

), using

CFL=1, 0.75 and 0.25.

Output your results for the time-levels t = 2,4,6,8,10 and plot the exact solution as well, which is

the same as the initial condition. How well does the solution at t = 10 agree with the exact

solution Discuss:

How does the agreement of the numerical prediction with the exact solution depend on the

grid resolution 

How does the agreement of the numerical prediction with the exact solution depend on the

CFL number?

What happens if you chose CFL>1 for either method, for example CFL=1.002?

7 Submission

This assignment, unlike assignment 1, consists of two parts

a. A project report, detailing any derivations and solutions and displaying the required graphs

b. C programs developed to solve some of the problems

You need to submit your programs and report for assessment; Submission of the report and

the code will be done via①. You may discuss your assignment work during your

workshop, and with others in the class, but what gets typed into your programs and the report

must be individual work, not copied from anyone else. So, do not give a hard or soft copy of

your work to anyone; do not “lend” your “Uni backup” memory stick to others for any reason at

all; and do not ask others to give you their programs “just so that I can take a look and get some

2018 The University of Melbourne 8

ideas, I won’t copy, honest”. The best way to help your friends in this regard is to say a very firm

“no” when they ask for a copy of, or to see, your programs or report, pointing out that your “no”,

and their acceptance of that decision, is the only thing that will preserve your friendship. A

sophisticated program that undertakes deep structural analysis of C code identifying regions of

similarity will be run over all submissions in “compare every pair” mode. Students whose

programs or reports are so identified will be referred to the Student Center. See②for more information.

7.1 Project Report

Your project report need not be a full technical report but should state all approximations made

and use figures of results to back up any conclusions. Be sure to include enough detail (using

appendices as necessary) so that your results could be reproduced by another researcher (or

yourself at a future date!) wishing to check or extend your findings. Your report will be primarily

assessed on the completeness of the results, and the visual/logical effectiveness of the manner in

which they are presented. Please type your report - scanned handwritten notes are sometimes

too difficult to read and mark and therefore will not be accepted.

7.2 C programs

The C codes are to be submitted on① where they would be tested on different inputs from

the ones you were provided to test your implementations.

7.2.1 Provided Code

where the parsing of data from command line is to be done and timing for each

task implemented.

where you will implement four functions (Question 2),

②(Question 4),(Question 5) and (Question 6), one for

each task., which acts as a header file to link the C file to the main file. You may edit this to

add any strt’s or the input arguments to the functions.

Input files (Question 2),(Question 4),

(Question 5) and (Question 6) which you must use as input during

execution of your program. During your submission, your code will be tested on different

infiles from the ones you were provided. Please make sure any data structures used to read

in these infiles are dynamically allocated to avoid any errors during runtime.

Remember to fill in your details in each file you write, such as the student name and student

number. Key points about your code for this assignment you need to understand are as follows:

For the purposes of the report, you can output as many files or terminal outputs as you

need. These outputs can be used to generate graphs/plots and values for the different tasks.

Once you have all information you need for the report, your code must be made submission

worthy i.e. only output the outfiles described above (5 outfiles are expected). This means

your code must not expect user input once you execute it, all inputs would come from the

infiles or the terminal before execution.

You have to parse the command line arguments (all infiles and any command line values)

i.e. no hardcoding the names of the infiles or the value for interpolation task. This is

because we will be using our own infiles, with different filenames and different locations.

2018 The University of Melbourne 9

? Plan before you write your code. Cover all possibilities regarding different inputs.

Dynamically allocating structures for the infile contents is a must so that for infiles with

more or less entries don’t end up giving you errors during submission.

7.2.2 Submitting on Dimefox

All submissions must be made through the program. Assignments submitted through

any other method will not be marked. Transfer your files to your home drive on dimefox. Check

the transfer was successful by logging into dimefox and doing the command on the terminal.

Perform the following set of commands on the terminal from your home location on dimefox

(making the right folders and transfering the files in the right location):

Remember to check this folder contains only the .c or .h files (if you use multiple c files and h

files) you need for the assignment and the CSV data file. Then try compiling your code using

on the terminal from the folder, to see if it works. The following compilation procedure must

return no errors or warnings:

You can test your compiled code then using:

There are 5 command line arguments here, one for each of the 4 coding functions. The argument

number 4 (?) is part of the interpolation task: the value of xo at which the interpolated value is to

be outputted. Once your code works for the provided infiles, it would help you if you changed the

infiles and the 4th argument to different values and see if the code still works and outputs

acceptable results. If this works, your code is now submission worthy. Finally, you will submit

your files using the command as follows:

Do NOT submit your infiles as this would likely corrupt your submission and would take time to

fix. Wait for a few minutes and then carry out the verify and check steps:

Look through this text file to check (a) your program compiled (b) it executed without error.

Please read ① for confirmation about how to use, especially how to

② your submission. No special consideration will be given to any student who has not used

2018 The University of Melbourne 10

s properly.

You must also check that you have no memory leaks in your code as loss of memory from your

implementation will result in deductions. You can use for this as follows:

This should output a leak summary such as the best case shown below. It must be pointed out

that the runtime with this takes longer as compared with direct execution. Plan your submission

accordingly.

All submissions should be in C99, and use no functions outside of the C standard library and

maths library. Some key points to consider about your submission and verification are outlined

here:

Submissions that can not be compiled or run by s????t system will receive zero marks

for the programming part.

Submissions are also limited to a maximum runtime of 200 seconds and maximum file size

per task of 2 MB. It would help if you don’t write any additional files through your code. If

it is absolutely necessary, then make sure the files do not exceed 2 MB individually. This

would give errors during your submission.

will contain feedback for individual tasks as with the first assignment. For

question 6, you will receive feedback on both files individually. Since the contents of the

infiles used for your code testing by the markers will not be disclosed to you, the feedback

will only say whether your output was correct or not.

Because of the above point, you can work on each task individually and submit your code

with just one or two tasks implemented. The feedback will skip over tasks not implemented

and only look at the outfile of the tasks implemented. So please use s????t as frequently

as possible. Submitting your code for the first time a day before the deadline is not ideal for

you or us.

You must parse all command line arguments, even if you do submit just one or two tasks to

see if your implementation works. The arguments will be in the exact sequence as given in

the execution statement above. As an example, if you want to test only the interpolation

task, your code will have to parse the 3rd and 4th arguments as all 5 arguments will be

present during execution. The takeaway is that you should ensure you are associating the

right infile to the right task.

Only your last submission is stored in the system i.e. everytime you use the new

submission overwrites the previous version. Keep a backup of your previous versions

somewhere safe in case your latest submission works worse than the previous.

2018 The University of Melbourne 11

8 Getting Help

There are several ways for you to seek help with this assignment. First, check the Assignment 2

Frequently Asked Questions wiki in the LMS (subsection Assignments). It is likely that your

question has been answered here already. You may also discuss the assignment on the

“Assignment 2” discussion board. However, please do not post any source code on the discussion

board. You may also ask questions during the workshops or send me (Professor Sandberg,) an email directly.

Note: Students seeking extensions for medical or other “outside my control” reasons should

email Professor Sandberg, as soon as possible after those

circumstances arise.