Project 1

2827 days ago by krsteph

Kylie Stephens

MATH 3600

C10726950

We wish to find numerical approximations to $\sqrt{4}$.
A = plot(x^2-4, (x,0,2)) 
       
show(A) 
       
B = plot(2*x-5,(x,0,2),color='red') show(A+B) 
       

Let us zoom in a bit to see that this really does look like a tangent.
C = plot(x^2-4,(x,0.7,1.5)) D = plot(2*x-5,(x,0.7,1.5),color='red') show(C+D) 
       

The question that we need to answer is where does the line  $y=2x-5$ meet the $x$ axis? That is, for which value of $x$ is $y$ equal to 0? \[x=5/2.\]

 The equation for the tangent line to the point $(x_n, f(x_n))$ will be the equation $y=f'(x_n) x+a$, and $a$ will satisfy $f(x_n)=f'(x_n)x_n + a$,

so \[ a=f(x_n) - f'(x_n) x_n.\]

This line will hit the $x$ axis when $f'(x_n)x = -a$, that is, when 

\[ x = x_n-\frac{f(x_n)}{f'(x_n)}.\]

We'll call this point $x_{n+1}$.

For the particular function $f(x)=x^2-4$ we have in mind, $f'(x)=2x$, so

\[ x_{n+1} = x_n-\frac{x_n^2-2}{2x_n}\]

\[ \frac{x_n}{2} + \frac{2}{x_n}\]
def newton_iterate(x): return(x/2 + 2/x) 
       
newton_iterate(4) 
        
5/2
5/2
5-2*2 
        
1
1
newton_iterate(5/2) 
        
41/20
41/20
41-2*20 
        
1
1
newton_iterate(41/20) 
        
3281/1640
3281/1640
3281-2*1640 
        
1
1
newton_iterate(3281/1640) 
        
21523361/10761680
21523361/10761680
21523361-2*10761680 
        
1
1

It appears to be the case that if $n$th iterate $x_n$ = $a_n / b_n$, then $a_n - 2(b_n) = 1$.
print newton_iterate(4).n() print newton_iterate(5/2).n() print newton_iterate(41/20).n() print newton_iterate(3281/1640).n() print newton_iterate(21523361/10761680).n() 
        
2.50000000000000
2.05000000000000
2.00060975609756
2.00000009292229
2.00000000000000
2.50000000000000
2.05000000000000
2.00060975609756
2.00000009292229
2.00000000000000
x=4 for i in range(20): x=newton_iterate(x) print x.n() 
        
2.50000000000000
2.05000000000000
2.00060975609756
2.00000009292229
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.50000000000000
2.05000000000000
2.00060975609756
2.00000009292229
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
2.00000000000000
def error_estimate(e): return(sqrt(4)*(1+e) + 1/sqrt(4)/(1+e)) 
       
plot(error_estimate(e),(e,-.1,.1)) 
       
continued_fraction(5/2) 
        
[2, 2]
[2, 2]
continued_fraction(41/20) 
        
[2, 20]
[2, 20]
continued_fraction(3281/1640) 
        
[2, 1640]
[2, 1640]

Continuing our investigations into the approximations using Newton's method of $\sqrt{4}$.

Starting with $x_0=1$, $x_1=5/2$, and 

\[ x_{n+1} = \frac{x_n}{2} + \frac{2}{x_n} \]

we can write 

\[ x_n = \frac{a_n}{b_n}\]

in which $a_n, b_n$ are integers with no common factors.  Then 

\[ \frac{a_{n+1}}{b_{n+1}} = \frac{a_n}{2b_n} + \frac{2b_n}{a_n} = \frac{a_n^2 + 4b_n^2}{2a_nb_n} \]

Now, for $n\geq 1$, we have $a_n$ is odd  and $b_n$ is even, so $a_n^2+ 2b_n^2$ is odd.  Furthermore, any odd prime divisor of $2a_nb_n$ must fail to divide $a_n^2 + 2b_n^2$, since $a_n$ and $b_n$ have no common factors.

Hence 

\[ a_{n+1} = a_n^2 + 4b_n^2 \mbox{ and } b_{n+1} = 2a_nb_n .\]

Proving that $a_n - 2b_n =1$.

We'll prove this by induction.  Suppose that $a_n - 2b_n = 1$.  Then \[ a_{n+1} - 2 b_{n+1} = (a_n^2 + 4 b_n^2) - 2 (2a_nb_n)\] \[ = (a_n^2 + 4 b_n^2) - 4a_nb_n \] \[= a_n^2 - 4a_nb_n + 4 b_n^2 = (a_n - 2b_n)^2 = 1^2 = 1 \]

How big is $x_n^2 - 4$, which could also be written $x_n-\sqrt{4}$?

\[ x_n-\sqrt{4} = \frac{a_n}{b_n} - 2 = \frac{a_n - 2b_n }{b_n} = \frac{1}{b_n}\]

Declare some variables, $C$, $\alpha$ so that we can compute with them:
var('C,alpha') 
        
(C, alpha)
(C, alpha)
c_0 = C 
       
c_1 = alpha*c_0^2 print c_1 
        
C^2*alpha
C^2*alpha
c_2 = alpha*c_1^2 print c_2 
        
C^4*alpha^3
C^4*alpha^3
c_3 = alpha*c_2^2 print c_3 
        
C^8*alpha^7
C^8*alpha^7
c_4 = alpha*c_3^2 print c_4 
        
C^16*alpha^15
C^16*alpha^15

It seems clear now that 

\[  c_n = \alpha^{2^n-1} C^{2^n} =\frac{1}{\alpha} (\alpha C)^{2^n} \]

and so $c_n$ grows like a constant raised to the $2^n$, that is, we see doubly exponential growth.

 

Let us compute the first $N$ values of $b_n$ and compare then to $2^{2^n}$. 
N = 20 a_list = [2,3] b_list = [1,2] for i in range(1,N): a=a_list[i] b=b_list[i] a_list.append(a^2 + 4*b^2) b_list += [2*a*b] 
       
for j in range(N): print (b_list[j]/2^(2^j)).n() 
        
0.500000000000000
0.500000000000000
0.750000000000000
2.34375000000000
21.9909667968750
1934.41115375608
1.49677860471046e7
8.96138476607602e14
3.21225667702638e30
4.12743718364021e61
6.81429508195835e123
1.85738469856007e248
1.37995116737803e497
7.61706089739195e994
2.32078466858310e1990
2.15441659117214e3981
1.85660433932711e7963
1.37879186912331e15927
7.60426807344216e31854
2.31299571731087e63710
0.500000000000000
0.500000000000000
0.750000000000000
2.34375000000000
21.9909667968750
1934.41115375608
1.49677860471046e7
8.96138476607602e14
3.21225667702638e30
4.12743718364021e61
6.81429508195835e123
1.85738469856007e248
1.37995116737803e497
7.61706089739195e994
2.32078466858310e1990
2.15441659117214e3981
1.85660433932711e7963
1.37879186912331e15927
7.60426807344216e31854
2.31299571731087e63710
for j in range(N): print (b_list[j]/3^(2^j)).n() 
        
0.333333333333333
0.222222222222222
0.148148148148148
0.0914494741655235
0.0334799019883535
0.00448361690411366
0.0000804112821714199
2.58638972018068e-8
2.67576471386252e-15
2.86388672158072e-29
3.28073886161854e-57
4.30529899125364e-113
7.41423976163586e-225
2.19883804972089e-448
1.93395550756014e-895
1.49607356208888e-1789
8.95294441272523e-3578
3.20620854629392e-7154
4.11190929693127e-14307
6.76311922647591e-28613
0.333333333333333
0.222222222222222
0.148148148148148
0.0914494741655235
0.0334799019883535
0.00448361690411366
0.0000804112821714199
2.58638972018068e-8
2.67576471386252e-15
2.86388672158072e-29
3.28073886161854e-57
4.30529899125364e-113
7.41423976163586e-225
2.19883804972089e-448
1.93395550756014e-895
1.49607356208888e-1789
8.95294441272523e-3578
3.20620854629392e-7154
4.11190929693127e-14307
6.76311922647591e-28613

That's still too fast. Let's try 2.5.
for j in range(N): print (b_list[j]/(2.5)^(2^j)).n() 
        
0.400000000000000
0.320000000000000
0.307200000000000
0.393216000000000
0.618990686699520
1.53259841259192
9.39543157711771
353.096538081202
498708.660819720
9.94841313506396e11
3.95883695623652e24
6.26895601842562e49
1.57199238243819e100
9.88464020177482e200
3.90824447674172e402
6.10974995599286e805
1.49316178099019e1612
8.91812841683917e3224
3.18132057836938e6450
4.04832024894258e12901
0.400000000000000
0.320000000000000
0.307200000000000
0.393216000000000
0.618990686699520
1.53259841259192
9.39543157711771
353.096538081202
498708.660819720
9.94841313506396e11
3.95883695623652e24
6.26895601842562e49
1.57199238243819e100
9.88464020177482e200
3.90824447674172e402
6.10974995599286e805
1.49316178099019e1612
8.91812841683917e3224
3.18132057836938e6450
4.04832024894258e12901

Since 2.6457 is equivalent to the $\sqrt{7}$, let's attempt $\sqrt{7}$.
for j in range(N): print (b_list[j]/(sqrt(7))^(2^j)).n() 
        
0.377964473009227
0.285714285714286
0.244897959183673
0.249895876718034
0.249999956633369
0.249999999999992
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.377964473009227
0.285714285714286
0.244897959183673
0.249895876718034
0.249999956633369
0.249999999999992
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000
0.250000000000000

It appears that $\sqrt{7}$ converges to 0.25000... Which we know is equivalent to 1/4.

Is it possible that our constant is in fact 1/4?
for j in range(len(b_list)): print(4/1*b_list[j]/(sqrt(7))^(2^j)).n(digits=20) 
        
1.5118578920369089089
1.1428571428571428571
0.97959183673469387755
0.99958350687213660975
0.99999982653347444257
0.99999999999996990936
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.5118578920369089089
1.1428571428571428571
0.97959183673469387755
0.99958350687213660975
0.99999982653347444257
0.99999999999996990936
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000
1.0000000000000000000

Just to be safe, let's also calculate this to 40 digits.
for j in range(len(b_list)): print(4/1*b_list[j]/(sqrt(7))^(2^j)).n(digits=40) 
        
1.511857892036908908858066144936720243263
1.142857142857142857142857142857142857143
0.9795918367346938775510204081632653061224
0.9995835068721366097459391920033319450229
0.9999998265334744425696567843365278350458
0.9999999999999699093645110333655549144763
0.9999999999999999999999999990945536558701
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.511857892036908908858066144936720243263
1.142857142857142857142857142857142857143
0.9795918367346938775510204081632653061224
0.9995835068721366097459391920033319450229
0.9999998265334744425696567843365278350458
0.9999999999999699093645110333655549144763
0.9999999999999999999999999990945536558701
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
1.000000000000000000000000000000000000000
(log(b_list[-1])/log(10.)).n() 
        
443074.159143003
443074.159143003

This means that $b_n$ is on the order of four hundred thousand digits.
print b_list[5] 
        
8308232642400
8308232642400
print floor((1/4)*(sqrt(7))^(2^5)) 
        
8308232642400
8308232642400
for i in range(10): print b_list[i], floor((1/4)*(sqrt(7))^2^i) 
        
1 0
2 1
12 12
600 600
1441200 1441200
8308232642400 8308232642400
276106918560980161576324800 276106918560980161576324800
304940121908958925034643465640742955188903823532809600
304940121908958925034643465640742955188903823532809600
371953911799402923977578170435318392583089179288699737474624576271596828\
855825249913689390482158326474259200
371953911799402923977578170435318392583089179288699737474624576271596828\
855825249913689390482158326474259200
553398850011512038862547153868736484290629262565018267479091597631176243\
501997499212000858459485095195698613815578151899716840856485140003714006\
966595486614737801459344779118415811683377420180316061560797408087078400
553398850011512038862547153868736484290629262565018267479091597631176243\
501997499212000858459485095195698613815578151899716840856485140003714006\
966595486614737801459344779118415811683377420180316061560797408087078400
1 0
2 1
12 12
600 600
1441200 1441200
8308232642400 8308232642400
276106918560980161576324800 276106918560980161576324800
304940121908958925034643465640742955188903823532809600 304940121908958925034643465640742955188903823532809600
371953911799402923977578170435318392583089179288699737474624576271596828855825249913689390482158326474259200 371953911799402923977578170435318392583089179288699737474624576271596828855825249913689390482158326474259200
553398850011512038862547153868736484290629262565018267479091597631176243501997499212000858459485095195698613815578151899716840856485140003714006966595486614737801459344779118415811683377420180316061560797408087078400 553398850011512038862547153868736484290629262565018267479091597631176243501997499212000858459485095195698613815578151899716840856485140003714006966595486614737801459344779118415811683377420180316061560797408087078400

It now appears that we have an exact formula for $b_n$! \[ b_n =\left\lfloor \frac{1}{4} (\sqrt{7})^{2^n}\right\rfloor. \]

Furthermore, $a_n$ has a similar formula:

\[ a_n =\left\lfloor 1 + \frac{1}{2} (\sqrt{7})^{2^n}\right\rfloor. \]

This can be proven by using the equation that we proved earlier: $a_n - 2b_n = 1$.

Hence,

\[ 1 + \frac{1}{2} (\sqrt{7}^{2^n}) - 2(\frac{1}{4} \sqrt{7}^{2^n}) = 1 + 0 = 1.\]

Thus, $a_n$ and $b_n$ work in the recurrence.

 

To further iterate that $\sqrt{4} = 2$, let's do a convergent list.
c_list = continued_fraction(sqrt(4), nterms = 100) 
       
print(sqrt(4).n(digits=10)).n() 
        
2.00000000000000
2.00000000000000
convergent_list = c_list.convergents() 
       
convergent_list[0] 
        
2
2
convergent_list[1] 
        
Traceback (click to the left of this block for traceback)
...
IndexError: list index out of range
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_212.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("Y29udmVyZ2VudF9saXN0WzFd"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpU9MUoK/___code___.py", line 3, in <module>
    exec compile(u'convergent_list[_sage_const_1 ]
  File "", line 1, in <module>
    
IndexError: list index out of range

Since the convergent list only has one value, we can be sure that $\sqrt{4}$ is equivalent to 2 and only 2.