Math 365 / Comp 365: Homework 1

Due: Thursday, February 4, 2016, beginning of class

Please submit a stapled hard copy of your work

Reminder: you are allowed — and in fact encouraged — to work on and discuss homework together. You should, however, write up your own assignments in your own words. If you work with another person or if you get significant help from a classmate or external resource, you should give that person or resource credit in your solution (at no penalty to you)

Anna Grim

Problem 1

The point of this problem is twofold: (i) to illustrate what can happen if you accumulate many truncations of numbers and (ii) to give you practice writing programs.

In the 1999 movie Office Space, a character creates a program that takes fractions of cents that are truncated in a bank’s transactions and deposits them into his own account. This is not a new idea, and hackers who have actually attempted it have been arrested. In this exercise you will simulate the program to determine how long it would take to become a millionaire this way.

Assume the following details:

  • You have access to 10,000 bank accounts

  • Initially, the bank accounts have values that are uniformly distributed between $100 and $100,000

  • The nominal annual interest rate on the accounts is 5%

  • The daily interest rate is thus .05/365

  • Interest is compounded each day and added to the accounts, except that fractions of a cent are truncated

  • The truncated fractions are deposited into an illegal account that initially has a balance of $0

  • The illegal account can hold fractional values and it also accrues daily interest

Your job is to write an R script that simulates this situation and finds how long it takes for the illegal account to reach a million dollars.

Here is some R help:

The following code generates the initial accounts:

accounts=runif(10000,100,100000)
accounts = floor(accounts*100)/100

The first line sets up 10,000 accounts with values uniformly between 100 and 100,000 The second line removes the fractions of cents (look at the data before and after that line is applied) To calculate interest for one day:

interest = accounts*(.05/365)

Depending on how you do it, you might want to use an if-then statement. For example, you might use something like

if (illegal > 1000000) break

The break command breaks out of a loop. Or, perhaps more elegantly, you might use a while loop

while (illegal < 1000000) { do stuff here }

You can find help on syntax in the Help Menu under “The R Language Definition.” That’s where I went to remind myself about the syntax for an if-then statement and a while loop.

accounts = runif(10000,100,100000); #Create accounts
accounts = floor(accounts*100)/100; 
illegal  = 0;
cnt      = 0;

while(illegal<1000000){
  cnt      = cnt+1;
  interest = accounts*(.05/365);
  deposit  = floor(interest*100)/100;
  illegal  = illegal+sum(interest-deposit); 
  accounts = accounts+deposit;
}
print(cnt/365);
## [1] 54.77808

After running the program several times, the approximate average number of years to accumulate a million dollars is 54 years.

Problem 2

The function

f(x)=(cos(x))^2+sin(x)-1

has roots at both \(\pi\) and \(\pi/2\).

  1. Using the convergence properties we discussed in class, estimate how many iterations of the bisection method are necessary to find the accurate solution to eight correct decimal places, when we start with a bracketing interval of [2.5,3.5].

On the i-th iteration, the error is approximately (b-a)/2^i. Therefore, if we would like to have an estimate correct to the eighth decimal place with the starting interval [2.5,3.5], then solving the equation (3.5-2.5)/2^i = 10^-8 will give the number of iterations. Solving this equation, we find that

        (3.5-2.5)/2^i = 10^-8
         log_{2}(2^i) = log_{2}(10^8)
                    i = log_{2}(10^8)  
                    i = 26.57

We need the ceiling of this answer because i must an integer, so we must perform 27 iterations.

  1. Use your bisect function to see how many iterations it actually takes.

Using the bisection method, only 20 iterations are necessary to find a solution that is correct to eight decimal places.

  1. By hand, compute a formula for Newton’s iteration for this particular choice of \(f(x)\); i.e., write a formula for \(x_{i+1}\) in terms of \(x_{i}\).

\(x_{i+1} = x_{i}-(cos(x_{i})^2+sin(x_{i})-1)/(-2cos(x_{i})sin(x_{i})+cos(x_{i}))\)

  1. Set \(x_0\) to 4.2, run 50 iterations of your Newton’s formula from c) to compute \(\{x_i\}_{i=0,1,\ldots,50}\), and compute \(e_i=|x_i-\pi|\) for all \(i\). Now plot \(e_i / e_{i-1}\) over all \(i\). What is the convergence rate to \(\pi\)?
  est = matrix(0,nrow=20,ncol=4);
  cnt = 0;
  x   = 4.2

while(cnt<20)
{
  cnt = cnt+1;
  x           = x-(cos(x)^2+sin(x)-1)/(-2*cos(x)*sin(x)+cos(x));
  est[cnt,1]  = x;
  est[cnt,2]  = abs(x-pi);
  
  if(cnt>1)
  {
    est[cnt,3] = est[cnt,2]/est[cnt-1,2];
    est[cnt,4] = est[cnt,2]/est[cnt-1,2]^2;
  }
}

print(est);
##           [,1]         [,2]         [,3]      [,4]
##  [1,] 2.987070 1.545222e-01 0.000000e+00 0.0000000
##  [2,] 3.177469 3.587602e-02 2.321739e-01 1.5025276
##  [3,] 3.142778 1.185795e-03 3.305259e-02 0.9213005
##  [4,] 3.141594 1.402228e-06 1.182521e-03 0.9972390
##  [5,] 3.141593 1.966427e-12 1.402359e-06 1.0000932
##  [6,] 3.141593 0.000000e+00 0.000000e+00 0.0000000
##  [7,] 3.141593 4.440892e-16          Inf       Inf
##  [8,] 3.141593 0.000000e+00 0.000000e+00 0.0000000
##  [9,] 3.141593 4.440892e-16          Inf       Inf
## [10,] 3.141593 0.000000e+00 0.000000e+00 0.0000000
## [11,] 3.141593 4.440892e-16          Inf       Inf
## [12,] 3.141593 0.000000e+00 0.000000e+00 0.0000000
## [13,] 3.141593 4.440892e-16          Inf       Inf
## [14,] 3.141593 0.000000e+00 0.000000e+00 0.0000000
## [15,] 3.141593 4.440892e-16          Inf       Inf
## [16,] 3.141593 0.000000e+00 0.000000e+00 0.0000000
## [17,] 3.141593 4.440892e-16          Inf       Inf
## [18,] 3.141593 0.000000e+00 0.000000e+00 0.0000000
## [19,] 3.141593 4.440892e-16          Inf       Inf
## [20,] 3.141593 0.000000e+00 0.000000e+00 0.0000000
plot(est[,3], type="o", col="blue")

plot(est[,4], type="o", col="blue")

The rate of convergence to pi is quadratic because lim \(e_{i}/e_{i-1}^2\)->c where c is approximately one. This can be seen graphically in the second plot where the limit of this ratio of errors tends to one.

  1. Repeat d) with \(x_0\)=1.5. What is the convergence rate to \(\pi/2\)? If your answer is different from d), explain why.
  est = matrix(0,nrow=20,ncol=3);
  cnt = 0;
  x   = 1.5;

while(cnt<20)
{
  cnt = cnt+1;
  x           = x-(cos(x)^2+sin(x)-1)/(-2*cos(x)*sin(x)+cos(x));
  est[cnt,1]  = x;
  est[cnt,2]  = abs(x-(pi/2));
  
  if(cnt>1)
  {
    est[cnt,3] = est[cnt,2]/est[cnt-1,2];
  }
}

print(est);
##           [,1]         [,2]      [,3]
##  [1,] 1.535502 3.529421e-02 0.0000000
##  [2,] 1.553162 1.763427e-02 0.4996363
##  [3,] 1.561981 8.815535e-03 0.4999093
##  [4,] 1.566389 4.407568e-03 0.4999773
##  [5,] 1.568593 2.203759e-03 0.4999943
##  [6,] 1.569694 1.101876e-03 0.4999986
##  [7,] 1.570245 5.509378e-04 0.4999996
##  [8,] 1.570521 2.754688e-04 0.4999999
##  [9,] 1.570659 1.377344e-04 0.5000000
## [10,] 1.570727 6.886721e-05 0.5000000
## [11,] 1.570762 3.443360e-05 0.5000000
## [12,] 1.570779 1.721680e-05 0.5000000
## [13,] 1.570788 8.608397e-06 0.4999998
## [14,] 1.570792 4.304186e-06 0.4999985
## [15,] 1.570794 2.152087e-06 0.4999986
## [16,] 1.570795 1.076060e-06 0.5000077
## [17,] 1.570796 5.381064e-07 0.5000710
## [18,] 1.570796 2.690647e-07 0.5000213
## [19,] 1.570796 1.345496e-07 0.5000640
## [20,] 1.570796 6.688806e-08 0.4971257
plot(est[,3], type="o", col="blue")

The rate of convergence to the root \(pi/2\) is linear such that \(lim e_{i}/e_{i+1}->1\). The rates of convergence are different due to the shape of the graph of this function. Newton’s method converges quickly when the zero is pi because the derivative is nonzero at this point. In contrast, the derivative is zero at pi/2, which causes the slow convergence.

Problem 3

Computer Problem 7 from Section 1.1 of the book (involving a determinant)

The determinant is defined by the function d(x)=x4-202x2+1404x-2475. To find the x values where the determinant is 1000, set d(x)=1000 and use the bisection method over an interval containing a zero. One zero is at x=-17.188500 and this answer can be obtained with 34 iterations such that the initial interval is (-20,-15). The second zero is at x=9.708299, which can be obtained using the bisection method in 29 iterations with the initial interval (7,11). The code for the bisection method is included below.

  a  = -20;
  b  = -15;
  it = 0;
   f = function(x) {x^4-202*x^2+1404*x-3475};
  
  while(it<40){
    it = it+1;
    mid = (a+b)/2;
    
    if(f(a)*f(b)>0){
      print('Error!');
    }
    else{
      if (f(mid)*f(a)<0){
        b  = mid;
      }
      else{
        a  = mid;
      }
      
      if (abs(f(mid))<10^-6){
        break;}
    }
  }
print(mid);
## [1] -17.1885
print(it);
## [1] 34

Problem 4

Computer Problem 8 from Section 1.1 of the book (involving the Hilbert matrix)

The determinant of the Hilbert matrix is determined by the polynomial: f(x,y) = -(1/11113200000) + x/10668672000 + (766019 y)/53343360000 - ( 431 x y)/27783000 - (557483 y^2)/27783000 + (50593 x y^2)/ 2116800 + (930779 y^3)/2116800 - (248 x y^3)/315 + (248 y^4)/315 + x y^4 - y^5. In this polynomial, the value y is an eigenvalue and x is the value of the uppermost left entry. This polynomial was determined by denoting the uppermost left entry as x, subtracting y from each diagonal, and then computing the determinant. The value of x that defines the largest eigenvalue can be determined by the function g(x)=f(x,pi). By the Intermediate Value Theorem, there exists a zero in the neighborhood (2,4). Now using the bisection method with this initial interval, the x value that makes pi aN eigenvalue is x=2.948011. The code for the bisection method is included below.

  a  = 2;
  b  = 4;
  it = 0;
   f = function(x) {-(1/11113200000) + x/10668672000 + (766019*pi)/53343360000 - (
  431*x*pi)/27783000 - (557483*pi^2)/27783000 + (50593*x*pi^2)/
  2116800 + (930779*pi^3)/2116800 - (248*x*pi^3)/315 + (248*pi^4)/315 + 
  x*pi^4 - pi^5};
  
  while(it<40){
    it = it+1;
    mid = (a+b)/2;
    
    if(f(a)*f(b)>0){
      print('Error!');
    }
    else{
      if (f(mid)*f(a)<0){
        b  = mid;
      }
      else{
        a  = mid;
      }
      
      if (abs(f(mid))<10^-6){
        break;}
    }
    print(mid);
    print(it);
  }
## [1] 3
## [1] 1
## [1] 2.5
## [1] 2
## [1] 2.75
## [1] 3
## [1] 2.875
## [1] 4
## [1] 2.9375
## [1] 5
## [1] 2.96875
## [1] 6
## [1] 2.953125
## [1] 7
## [1] 2.945312
## [1] 8
## [1] 2.949219
## [1] 9
## [1] 2.947266
## [1] 10
## [1] 2.948242
## [1] 11
## [1] 2.947754
## [1] 12
## [1] 2.947998
## [1] 13
## [1] 2.94812
## [1] 14
## [1] 2.948059
## [1] 15
## [1] 2.948029
## [1] 16
## [1] 2.948013
## [1] 17
## [1] 2.948006
## [1] 18
## [1] 2.948009
## [1] 19
## [1] 2.948011
## [1] 20
## [1] 2.94801
## [1] 21
## [1] 2.948011
## [1] 22
## [1] 2.948011
## [1] 23
## [1] 2.948011
## [1] 24
## [1] 2.948011
## [1] 25

If you click on the question mark just to the left of the Knit HTML button, you can find a Markdown Quick Reference.