COSC2500/COSC7500—Semester 2 2024Python

Java Python COSC2500/COSC7500—Semester 2, 2024

MODULE 2. SOLVING EQUATIONS

Exercises—due 3:00 pm Friday 20th September

Required

The grade of 1–5 awarded for these exercises is based on the required exercises. If all of the required exercises are completed correctly, a grade of 5 will be obtained.

R2.1 Implement the bisection method, fixed-point iteration, and Newton’s method. Compare these for the following tasks:

a) Find the first two positive roots of f(x) = cos(x).

b) Find both roots of f(x) = (x − 1)3 (x + 2).

You can compare the methods by (a) the number of iterations required to find the roots to some target accuracy and also (b) by the difficulty of choos-ing starting points or intervals that give convergence.

Note that Sauer provides code for bisection and fixed point iteration (Pro-gram 1.1 and Program 1.2). You can use these, or write your own code, or use other code.

Checklist

Did you search for the roots using the required methods?

Did you show the results for each method?

Did you explain how you did it?

Did you compare the methods?

Did you discuss any difficulties?

Did you discuss adjusting fixed-point iteration when needed?

Did you discuss any unusual results?

Did you include code where appropriate?

Did you reference the source of any code you didn’t write?

R2.2 Implement Golden Section Search, Successive Parabolic Interpolation, and Newton’s method for finding minima. Compare these for the following task: find all maxima and minima of

You can compare the methods by (a) the number of iterations required to find the maxima and minima to some target accuracy and also (b) by the difficulty of choosing starting points or intervals that give convergence.

Note that Sauer provides code for golden section search and successive parabolic interpolation iteration (Program 13.1 and Program 13.2). You can use these, or write your own code, or use other code.

Checklist

Did you search for the maxima and minima using the required methods?

Did you show the results for each method?

Did you explain how you did it?

Did you compare the methods?

Did you discuss any difficulties or unusual results?

Did you include code where appropriate?

Did you reference the source of any code you didn’t write?

R2.3 From Sauer computer problems 2.3, slightly modified (pg 94 in 2E, 101 in 1E):

a) For the n × n Hilbert matrix, with Hij = 1/(i + j − 1), and x being a vector of all ones, calculate b = Ax. Use this b and your A to find xc, the double-precision computed solution, using either naive Gaussian elimination and Matlab’s backslash operation. Find the infinity norm of the forward error and the error magnification factor of the problem Ax = b (see below for definitions and details), and compare with the condition number of A. Do this for n = 2, 5 and 10. A naive Gaussian elimination code (based on the code in Sauer) is available on Black-board. You might want to use the Matlab function hilb() to generate the Hilbert matrix.

c) Repeat for Ai j = |i − j| + 1, for n = 100, 200, 300, 400, and 500.

e) For what values of n do the solutions in the above problems have no correct significant digits?

The infinity norm of a vector is the largest absolute value of the elements. In Matlab, you can use either norm(x,Inf) or max(abs(x)) to find the infinity norm of a vector x, ||x||∞. The forward error is the infinity norm of the dif-ference between the actual solution x and your computed solution xc. The backward error is the error in the constant vector b; we will assume that the relative backward error in the machine epsilon ϵ (Sauer uses a different def-inition for the backward error). The error magnification factor is the ratio of the relative forward error to the relative backward error, so

The error magnification factor tells us how errors in our constant vector b will affect our computed solution xc.

Checklist

Did you show the results for each case?

Did you explain what you did?

Did you compare your errors with the condition number?

If you couldn’t reach a size which gives zero correct digits, did you state the largest size you tried?

Did you discuss any difficulties or unusual results?

Did you include code where appropriate?

Did you reference the source of any code you didn’t write?

R2.4 Discuss the code find_eq.m or find_eq.py available on Blackboard. (Dis-cuss either the Matlab version or the Python version; you don’t have to discuss both. You are not expected to run this code.)

Note: What kind of things can you write about for a question like this? You might like to answer some of the following questions: What does the code do? What method does it use? Does it appear to be correct? Is there anything unusual or interesting about it?

In this particular case, why might there be code for both Newton’s method and bisection?

Checklist—Your discussion should include:

The general function of the code.

Methods used by the code.

Unusual features of the code.

Does the code appear to be correct? Why?

Are there any problems with the code?

Can you suggest any improvements to the code?

Additional

Attempts at these exercises can earn additional marks, but will not count towards the grade of 1–5 for the exercises. Completing all of these exercises does not mean that 6 marks will be obtained—the marks depend on the quality of the answers. It is possible to earn all 4 marks without completing all of these additional exercises.

A2.1 How are methods for root-finding affected by error? Try adding random errors such as you used in the Module 1 exercises.

A2.2 For fixed-point iteration to converge to a fixed point, the absolute value of the slope need to be less than 1. If we are using fixed-point iteration to find a root (rather th COSC2500/COSC7500—Semester 2, 2024Python an being interested in the fixed-point as such), what restrictions does this place on the original function? Is it possible to trans-form. a function for which the required conditions are not satisfied into one where they are satisfied? Is this possible for only some functions, or for all continuous functions?

A2.3 How are methods for finding minima affected by error? Try adding random errors such as you used in the Module 1 exercises.

A2.4 Find the (absolute, not local) maximum of the function calculated by spikey fn() (provided in spikey fn.m).

A2.5 Implement one of the methods for root-finding or optimisation from this module in a compiled language of your choice. Compare the performance of your compiled code with your Matlab (or Python) version.

A2.6 Write code that automatically chooses starting points or intervals for one or more of the methods for root-finding or optimisation in R2.1 or R2.2 that will successfully find at least one root or minimum for each of the cases in those questions. Find a function that it will fail for.

A2.7 Implement a steepest descent method to find minima of a function of two variables (i.e., z = f(x, y)).

A2.8 Consider the carbon monoxide levels in a restaurant, which consists of a main room (with a volume of 200 m3 ), with two smaller rooms of volume 100 m3 . One small room is a smoker’s section, and the other is a chil-dren’s section. The small rooms are connected to the main room, one at each end, but not to each other. Ventilator fans move air at 200 m3/hr into the smoker’s room, 50 m3/hr into the children’s room, 150 m3/hr out of the children’s room, and 100 m3/hr out of the children’s room end of the main room. The small rooms exchange 25 m3/hr with their ends of the main room. We can divide the main room into two parts, one at each end, with air exchanged at 50 m3/hr. The outside air has a carbon monoxide concentration of 2 mg/m3 . The smokers add carbon monoxide at the rate of 1000 mg/hr. A faulty grill at the smoker’s end of the main room adds 2000 mg/hr.

a) What is the equilibrium carbon monoxide concentration in each room?

b) What would the equilibrium concentration be if the faulty grill was not producing carbon monoxide?

c) Explore the possibilities of changing the ventilation system (you might like to write your code in the first place to make this easy to do). Would you make any recommendations about the ventilation?

(From S. C. Chapra, R. P. Canale, Numerical methods for engineers, 3rd ed, McGraw-Hill 1998, pp 323–324.)

A2.9 On Blackboard, you will find C code (from the book by M. Pal) that im-plements Gaussian elimination. Compare the execution time of the C im-plementation with the naive Gaussian elimination Matlab code from Black-board. You can also compare with other methods (e.g., combining this ques-tion with A2.10).

There are several things (steps along the way) that you will need to do to complete this exercise (often the case when solving computational prob-lems).

Some hints:

• You will need to make a few simple modifications to the C code. For example, it will currently only handle systems of size up to n = 10.

• The program is set up to read in the matrix and constants (right-hand vector) from the user at runtime. You probably don’t want to type in large matrices(!). Perhaps you can create the matrices in matlab and cut-and-paste them with the mouse into your C code (declaring the variables), and then remove the code that reads from the user?

• To measure computation time for your C program, you will need to use functions provided in the library time.h (specifically clock() and time())

A2.10 The time required to solve a linear system depends on the size of the ma-trix. It will also depend on the method of solution, and whether or not the matrix is sparse. Find the dependence of the time on the size of the ma-trix by solving linear systems of different sizes (you can use tic and toc to time each case). Compare naive Gaussian elimination, the backslash oper-ator and inv(), and two or more iterative methods. Do this for both sparse and non-sparse matrices. You could use the systems from Sauer computer problems 2.5 1 & 2, and the Matlab code from the chapter. Note that it’s easy to convert a sparse matrix to a full matrix, using full(). Plot the times on a log–log graph. What do the slopes of the lines tell you? For what size full matrix do we want to use iterative methods? For what size sparse matrix do we need to store it as a sparse matrix? Do we need to use iterative methods for sparse matrices? Are these choices due to memory or execution time?

Note that different methods and different types of matrices will be able to reach different sizes before becoming too slow. This makes it difficult to have a single loop over sizes covering all cases. It is easier to use separate loops, to different maximum sizes, for different cases.

A2.11 a) Solve the following sparse system, using Matlab’s backslash and one or more iterative solvers (from Matlab or Sauer or elsewhere), and compare the errors (as in R3.2) and time taken.

The correct solution is xi = 1. Solve for sizes n = 100 and n = 100000 (and other sizes, if you wish).

b) Do the same for the sparse system with (1, 2, 1) replacing the (−1, 3, 1) in (a), and the constant vector with [1, 0, .., 0, −1]. The solution is [1, −1, 1, −1, 1, −1, ...].

These matrices are from Sauer computer problems 2.5 (pg 121–122 in 3E, 116 in 2E, 122 in 1E). The preceding sections is Sauer will be useful if you use the iterative solvers from Sauer.

A2.12 Consider a system described by

If we have a set ym of M measurements of y at times tm, we have a linear system that, if M ≥ N, can be used to find an if the values of bn are known.

a) Generate and solve a test system of this type, with a random error to simulate measurement errors in ym.

b) How does the solution compare with the actual values of an used to generate the system? How does the accuracy depend on M/N?

c) What happens if 2 values of bn are similar?

A2.13 Make a table or list comparing methods for solving systems of non-linear equations. Test one or more of these methods