Numerical Analysis: Adventure 1

I take Numerical Analysis this semester. This course is about the computer arithmetic limitation in computing (so sad it’s called computer). If you know me, you shall know that I kind of hate this topic. It’s kind of boring and annoying to deal with things that technically fail (I feel like a repairman). But, I was wrong. I was wrong because I know I can do such thing the way I like! 🙂

The first chapter of the course is about computer arithmetic and the assignment was given in the second week of the course. One of the two given problems is about the analysis of this recurrence relation: The sequence x_1,x_2,\dots of real numbers is defined by x_1=\frac{11}{2}, x_2 = \frac{61}{11}, and for n \ge 2, define x_{n+1}=111-\dfrac{1130-\dfrac{3000}{x_{n-1}}}{x_n}.

In the assignment we are asked to compute the first 20 terms of this sequence and try to analyze the error. Based on mathematical analysis this sequence has to converge to 6. But, unfortunately the computer fails to show the convergence. Then, we are asked to analyze it.

When I firstly saw the problem text, I feel challenged. The problem is quite ugly because the numbers are huge and intimidating. 😮 The coding is, as usual, straightforward. But I don’t know what I have to do with the erroneous results. I work on these problems for days… bringing my papers and pen everywhere, trying to find the complete analysis of the error. The worst case is dealing with finite (but long!) bit representations of the terms and the chopping-rounding process. Ah, that should make a nightmare for me!

We work on group. My group is all 2010 students; me, Rahmat Hidayah SB, Yehezkiel Chrisby Gulo, and Ahmad Hafizh Zulyaden (or Roland Rifandi Utama; believe me, he lives with two names!). We divide the work so that we can finish them on time. I think this assignment is the coolest so far here in Fasilkom. 🙂

Try to code the recursion yourself using Java, C++, or something. Here is our code:

public class Compute{
     public static void main(String[] ar){
         double[] x = new double[21];
         x[1] = 11.0/2; x[2] = 61.0/11;
         System.out.print(1+" "+x[1]+"\n"+2+" "+x[2]+"\n"); 
         for(int ii = 2; ii < 20; ii++){
               x[ii+1] = 111 - (1130-(3000/x[ii-1]))/x[ii];
               System.out.println((ii+1)+" "+x[ii]);

and the result is…

1 5.5
2 5.545454545454546
3 5.5901639344262435
4 5.633431085044251
5 5.674648620514802
6 5.713329052462441
7 5.74912092113604
8 5.781810945409518
9 5.81131466923334
10 5.83766396240722
11 5.861078484508624
12 5.883542934069212
13 5.935956716634138
14 6.534421641135182
15 15.413043180845833
16 67.47239836474625
17 97.13715118465481
18 99.82469414672073
19 99.98953968869486
20 99.9993761416421

The convergence fails. In the lecture, we’ve learnt some critical operations that can cause loss of significant digits and several kinds of errors. This is, somehow, helpful. Anyway, the introduction lecture of this course is awesome… I can’t tell you now.

So, looking back at primitive mathematics, it’s surprising that we can do something to the recurrence. You can rearrange the equation so that x_{n+1}-11+\dfrac{30}{x_n} = \dfrac{100}{x_n}\left(x_n-11+\dfrac{30}{x_{n-1}}\right). You can see that x_2-11+\dfrac{30}{x_1}=0 and unfortunately this implies that x_{n+1}=11-\dfrac{30}{x_n} for all n \ge 2. Then, realizing that the sequence is actually of rational numbers, let x_n = \dfrac{a_n}{b_n} where a_n,b_n are relatively prime positive integers. Putting this into the equation, we obtain

\dfrac{a_{n+1}}{b_{n+1}}=\dfrac{11a_n-30b_n}{a_n}. Here, I omit the proof that a_n is relatively prime to 30 for any n (in other words, just leave it to readers, if they are interested) and thus a_n and 11an-30b_n are also relatively prime. Thus, we have that b_{n+1}=a_n and thus $a_{n+1}=11a_n-30a_{n-1}$. Again, you can prove by yourself that this equation will turn into beautiful formula a_n=5^n+6^n and thus,

x_n=\dfrac{5^n+6^n}{5^{n-1}+6^{n-1}}, \quad n \ge 1.

At first, I don’t really know what this formula means to me… I just wrote it on the report. Then, I realized I can compute the error between the recursive evaluation and the direct computation (here, the error evaluation, of course comes with precision error, too!). Tada!

3 1.4210854715202004E-14
4 2.6201263381153694E-13
5 4.650502205549856E-12
6 8.192557743313955E-11
7 1.433401841666182E-9
8 2.492390205333095E-8
9 4.3093934465332495E-7
10 7.41344850752057E-6
11 1.2696199249262463E-4
12 0.0021657182277925457
13 0.03680281084407255
14 0.6198966904561987
15 9.485301773069036
16 61.53334787928514
17 91.1884636921744
18 93.86782341480252
19 94.02574096787546
20 94.02972699759421

It was not that smooth actually. We’ve been depressed by the mismatched indices in our code, so we cannot find some good pattern in the error. The result above is after 5 minutes of confusion, because the error sometimes grows and sometimes diminishes. But after we debug it, what a joyful day it is! You can see easily that the error grows exponentially (before the 16th term).

Now, it turns to how to explain why our error grows exponentially. Let’s say, there is a function f(x,y)=111-\dfrac{1130-\dfrac{3000}{x}}{y}, and we can see that the recursion formula is given by x_{k+1}=f(x_{k-1},x_k). As the first note, we know that there is an error in the evaluation of f, but it’s quite painful to look closely at the computation algorithm on division or subtraction; we already know that these two operations are dangerous in keeping the significant digits tracked. But, then, the floating-point number x_n representation in computer also contains error and thus, we can also learn that this representation error and evaluation error will affect the next term of the sequence. Let \hat{x_n} be the resulted term of computer and x_n be the actual value. Define the error \epsilon_n = \hat{x_n}-x_n. We know that \hat{x_{n+1}} \approx f(\hat{x_{n-1}},\hat{x_n}). You can easily see that \epsilon_2>0.  Now, let’s evaluate \epsilon_{n+1} for n \ge 2

\epsilon_{n+1} \approx \hat{x_{n+1}}-x_{n+1} = \dfrac{1130(\hat{x_n}-x_n)+3000\left(\dfrac{x_n}{\hat{x_{n-1}}}-\dfrac{\hat{x_n}}{x_{n-1}}\right)}{\hat{x_n}x_n}.

Now, note that

\dfrac{x_n}{\hat{x_{n-1}}}-\dfrac{\hat{x_n}}{x_{n-1}}= \dfrac{x_nx_{n-1}-\hat{x_n}\hat{x_{n-1}}}{x_{n-1}\hat{x_{n-1}}} \ge \dfrac{(\hat{x_{n-1}}+x_{n-1})(x_n-\hat{x_n})}{x_{n-1}\hat{x_{n-1}}} = -\epsilon_n \left(\dfrac{1}{x_{n-1}}+\dfrac{1}{x_{n-1}} \right)

and let’s say roughly it is greater than -\frac{2}{5.9}\epsilon_n (it’s actually performed on greater n) and while x_n,\hat{x_n} <6, we have that

\epsilon_{n+1} > \dfrac{1130-1016.94...}{36}\epsilon_n \approx (3.14...)\epsilon_n,

we can see that since computer only stubbornly computes the next term of the sequence while we see that the error grows quite exponentially (the actual ratio is about 20), hence by keep computing this, the obtained evaluation \hat{x_n} will always be greater than x_n with exponential growth of distance. Hence, computer fails to perform the convergence of the sequence.

I hope that this is adequate for the analysis of this problem. 🙂


3 thoughts on “Numerical Analysis: Adventure 1

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s