Numerical Analysis

Homework 2

Due March 11th before class

As we have learn in class, floating point arithmetic can be tricky sometimes. Bad implementation can

introduce numerical instability. We have expericed one important class of such instability - catastrophic

cancellation. In this homework, you will explore other types of numerical instability. More tips and tricks

are given in the textbook in chapter 2.7 - you are encouraged to read it through.

Problem 1

(5 points) Consider the following integral:

In =

∫ 1 0

xn

x + 5 dx, n = 0, 1, 2, · · · ,

(i) Prove the recurrence relation   I0 = ln 1.2

In = −5In−1 + 1n, n = 1, 2, · · · ;

(ii) Use the above recurrence to compute I1, I2, · · · , I30. (The best way to do it is to write a very short

program - but it is not meant to be a programming problem, you can use a calculator to compute

30 numbers if you prefer.) Print out your answers to see if they look reasonable. If not, give an

explanation.

Problem 2

(5 points) Assume we are doing arithmetics in a floating point system with base 10 and 8 digits in the

mantissa. (i.e. numbers are represented as: x̃ = 0.a1, · · · , a8 · 10e.) We do not specify the lower and upper

bound of the exponent, but assume they are sufficiently small and large respectively.

Floating point addition is carried out in the following fashion: consider the addition of x1 = m1 · 10e1

and x2 = m2 · 10e2 , where e1 > e2. We first rewrite x2 as: x2 = (m2 · 10e2−e1 ) · 10e1 := m̃2 · 10e1 , where

m̃2 = m2 · 10e2−e1 is the new mantissa. Then we perform the addition, rounding the result to the closest

floating point number.

1

With the information given above, consider the addition of the three floating point numbers:

x = 0.23371258 × 10−4

y = 0.33678429 × 102

z = −0.33677811 × 102

(i) What are the results of (x+y)+z and x+(y+z) respectively?

(ii) What is the exact result (in the real number system)? Compare the results and make a conclusion.

(Does the associative law still hold in the floating number system? What’s the lesson learned?)

Problem 3

(0 point - You do not have to submit your work on this one. But if you have time, it is a fun problem

to think over! If you do submit your answer to this one, we will take a look and provide feedback.) Bessel

functions (of the first kind), denoted as Jn(x), is an important class of special functions. They are formally

defined to be solutions to the following ordinary differential equation:

x2y ′′

+ xy ′ + (x2 − n2)y = 0,

where n is a parameter. (You do not need to know ordinary differential equations to solve this problem!)

They satisfy the following three-term recurrence relation:

2x

n Jn(x) = Jn−1(x) + Jn+1(x).

It provides us a tool to evaluate Jn(x) for all n given two initial values. Now let’s fix x = 1.0 and try it out.

(i) Given J0(x) = 0.765197686557967, J1(x) = 0.440050585744934. Use the recurrence to compute

J0(x), · · · , J20(x). Plot (n, Jn(x)) using the semilogy plot.

(ii) Given J20(x) = 3.873503008524651e − 25, J19(x) = 1.548478441211653e − 23. Use the recurrence to

compute J0(x), · · · , J20(x). Plot (n, Jn(x)) using the semilogy plot.

(iii) The exact values of the Bessel functions can be obtained by calling the function besselj in MATLAB.

plot the exact (n, Jn(x)), compare with the results obtained above, make a conclusion.

(iv) Other interesting facts and algorithms about Bessel functions can be found here. Take a look when you

have time! Enjoy!

2