Skip to main content

辅导案例-ENGR20005-Assignment 1

By May 15, 2020No Comments

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−200×3+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 Bisec tion.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 th ings 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 (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 ([email protected]) 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


Author admin

More posts by admin