Introduction to Computational Mathematics Assignment - Matlab and Finite Precision Arithmetic

1. Matlab: prime numbers.

Write a function sol=myprime(n), which returns 1 if natural number n is prime and 0 otherwise. Note that, by convention, 0 and 1 are not considered prime numbers. Use the command mod to check whether a number divides n. Test your code with n=0, n=1, n=5 and n=25 and report the results.

2. Matlab: comparing elements of two matrices.

Write a function num=howmanybigger(A,B) that does the following. If the dimensions of the matrices A and B do not match, it uses the error command to break with an appropriate error message. Otherwise it returns the number of elements in A which are bigger than the elements in B at the corresponding positions.

Submit a script test howmanybigger.m that tests your code with

and report the results.

3. Floating point number system.

On a computer with a binary number system, the distance between 7 and the next larger floating point number is 2-12. What is the (exact) distance between 71 and the next largest floating point number?

4. Avoiding catastrophic cancellation.

Consider the expression - z_A = 1/(√(1 + x²) -√(1 - x²)).

(a) Explain why the formula for z_A is susceptible to catastrophic cancellation errors for x close to 0.

(b) Use reformulation to find an alternative expression z_B for expression z_A in a way that avoids catastrophic cancellation for x close to 0. (Hint: p² - q² = (p-q)(p+q).)

(c) Illustrate how algorithm B that uses the reformulated expression z_B you propose improves the accuracy for calculations in the floating point system F[b = 10, m = 4, e = 1] with rounding, compared to algorithm A based on the original expression z_A. Use x = 1.234 · 10^-3.

(c1) Find a highly accurate approximate value for z_A(x) = z_B(x) (for example by calculating the expression in Matlab; for our purposes it does not matter much whether you use expression z_A or z_B to compute this highly accurate value).

(c2) Compute the machine epsilon for the floating point system (with rounding!).

(c3) Find the value of z¯_A, obtained by evaluating the original expression z_A in the floating point system.

(c4) Find the value of z¯_B, obtained by evaluating the reformulated expression z_B in the floating point system.

(c5) Compute the relative errors in the approximations of z_A and z_B relative to the highly accurate value for z_A = z_B,

|δz_A| = |(z_A - z¯_A)/zA|

and

|δz_B| = |(z_B - z¯_B)/z_B|.

Compare the two relative errors to each other, and compare them to the machine epsilon. Explain how this illustrates that the alternative expression succeeds in avoiding the catastrophic cancellation step, showing that algorithm B is a numerically stable algorithm to compute z_A = z_B, while algorithm A based on the original expression is numerically unstable.

5. Matlab: recursion

It can be convenient to solve a problem recursively. A recursive function is a function that calls itself. The following function sums up the numbers from m to n in an unusual (recursive) way.

function rs=recursivesum(m,n)

if (n>m)

rs=m+recursivesum(m+1,n)

elseif (n==m)

rs=m;

else

rs=0;

end

end

(a) Copy this code into function recursivesum.m. Report and explain what gets printed to the screen when you invoke recursivesum(1,5), and why it is printed in that order.

(b) Next, write a function nf=recursivefactorial(n), which computes the factorial of n in a recursive way.

(c) Test your code with n=0, n=1 and n=5 and report the results.