Homework 3

INSERT YOUR NAME HERE

Directions: Add work to this notebook to solve the problems below.

Be careful to name the objects you create as described in the problem. Some problems may include some tests you can run to test your code.

Before you submit your work, make sure it runs correctly without errors. Click “Kernel > Restart & Run All” to restart and run everything.

Some problems are derived from problems in these books:

  • LL: Programming for Computations - Python by Svein Linge and Hans Petter Langtangen, 2nd edition.
  • L: A Primer on Scientific Programming with Python by Hans Petter Langtangen, 2nd edition.
  • TAK: Applied Scientific Computing With Python by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh.

Standard imports

In [ ]:
import math as m

1. Partial sums

It is often natural to implement a sequence as a function. For example, we might consider the sequence $$a_k = \frac{1}{k+1}$$ This can be implemented in Python as:

def a(k):
    return 1/(k+1)

Write a function partial_sum(f, n) which takes as input a function $f$ and an integer $n \geq 0$ and returns the sum $$\sum_{k=0}^{n-1} f(k).$$ If $n=0$, the function should return zero. The function $f$ represents a sequence and should be assumed to take as input an integer and return a number.

In [ ]:
 

Tests for your code:

Here we test that you are computing $\sum_{k=0}^{n-1} 1$ correctly:

In [ ]:
def b(k):
    return 1
ans = partial_sum(b, 100)
if ans ==100:
    print('That is correct, the answer is 100.')
else:
    print(f'That is incorrect, the answer is 100 not {ans}.')

Here we test the result with the a(k) function from the problem.

In [ ]:
def a(k):
    return 1/(k+1)
ans = partial_sum(a, 1000)
if abs(ans - 7.4854) < 10**-4:
    print('That is correct, the answer is approximately 7.4854.')
else:
    print(f'That is incorrect, the answer is approximately 7.4854 not {ans}.')

2. Sum until small

As in the previous problem Partial sums, we will sum sequences given by functions.

Write a function sum_until_small(f, epsilon) which returns $$\sum_{k=0}^{n-1} f(k)$$ where $n$ is the smallest non-negative integer such that $|f(n)|<\epsilon$. The empty sum is zero, so if $|f(0)| \geq \epsilon$ then the function should return zero.

Warning: The evaluating the function will enter an infinite loop if $f(n)$ is not eventually less than zero. If your program ever ends up in an infinite loop, you can interupt the kernel (Python) by clicking Kernel > Interupt (or clicking the black square button on the toolbar). This should break the infinite loop and allow you to fix things. If this doesn't work you can click Kernel > Restart (or the restart button on the toolbar). This is more drastic. It will shut down the running Python instance and start a new one. This means you will loose access to all previously constructed variables. You will need to reexecute your code cells to get access to those variables.

In [ ]:
 

Tests for your code:

If $c(k)=\frac{1}{2^k}$ then $\sum_{k=0}^n = 1 - \frac{1}{2^{n+1}}$. From this (and the fact that computers can do exact arithmetic with some powers of two) sum_until_small(c, 1/2**i) should return 2-1/2**i for each integer $i \geq 0$.

In [ ]:
def c(k):
    return 1/2**k
for i in range(10):
    assert sum_until_small(c, 1/2**i) == 2 - 1/2**i, \
        f'Failure when i={i}'

Also sum_until_small(c, 2) should return zero since the first term is smaller than $2$.

In [ ]:
sum_until_small(c, 2)

3. Geometric series

A geometric series is a series of the form $$c + cr + cr^2 + cr^3+\ldots=\sum_{k=0}^\infty c r^k.$$

Use closures (discussed near the end of the Plotting Graphs notebook from class) to define a function geometric_sequence(c, r) which returns a function f such that f(n) is the term of the sequence $\{c, cr, cr^2, cr^3, \ldots.\}$ with index $n$. (Observe $f(0)=c$, $f(1)=cr$, $f(2)=cr^2$ and so on.)

In [ ]:
 

Tests for your code:

In [ ]:
f = geometric_sequence(3, 1/4) # The sequence f(n) = 3/4^n
assert f(0) == 3
assert f(1) == 3/4
assert f(2) == 3/16

If you have already done the partial sum problem, you can check the formula $$\sum_{k=0}^\infty c r^k=\frac{c}{1-r}.$$ In $f$ as above, we have $c=3$ and $r=1/4$, so the sum should converge to four.

The following should print True:

In [ ]:
abs(partial_sum(f, 100) - 4) < 10**-15

4. Guess the limit

Print some partial sums for the series $$\sum_{n=0}^\infty \frac{1}{(n+1)(n+2)(n+3)}.$$ You may use the functions you have already written in this assignment.

Store your guess for the limit in the variable guess_for_limit as a float. Your grade for this problem will depend on getting this guess correct.

In [ ]:
 

5. Guess the remainders

This problem also concerns the series $$\sum_{k=0}^\infty \frac{1}{(k+1)(k+2)(k+3)}.$$ Let $L$ be the limit of the series. The remainders are the differences with the partial sums $$R_n = L-S_n = L - \sum_{n=0}^{n-1} \frac{1}{(k+1)(k+2)(k+3)}.$$ These quantities are computed and stored in a list below and printed out, assuming your partial_sum code and guess_for_limit was correct.

In [ ]:
def a(k):
    return 1/((k+1)*(k+2)*(k+3))
R = 10*[None]
for n in range(10):
    R[n] = guess_the_limit - partial_sum(a, n)
    print(f'R_{n} = {R[n]}')

The actual remainders are fractions. Guess a formula for the numerator and denominator. Write a function remainder_guess(n) which returns a pair of integers (p, q) such that $R_n=\frac{p}{q}$.

Hint: Consider the recipricals of the remainders. Try to look for patterns.

In [ ]:
 

Tests for your code:

In [ ]:
# Test if remainder guess works
for n in range(10):
    numerator,denominator = remainder_guess(n)
    assert type(numerator) == int, 'The numerator should be an int'
    assert type(denominator) == int, 'The denominator should be an int'
    if numerator/denominator - R[n] > 10**-8:
        print(f'Your guess seems to be incorrect when n={n}')

6. Inductive proof

In this problem you will inductively prove your conjectures. If $L$ is the proposed limit (guess_for_limit) and $p_n/q_n$ is your guess for the remainder, then you are guessing that the partial sums $$S_n = \sum_{k=0}^{n-1} \frac{1}{(k+1)(k+2)(k+3)}$$ are given by $L-\frac{p_n}{q_n}$. Prove this by induction.

It follows that the series sums to $L$ once you observe that $\lim_{n \to \infty} \frac{p_n}{q_n}=0$.

Remarks: You may want to review the inductive proof provided in the Series notebook from class to get the form correct. Include your proof in a Markdown cell in the notebook. (Including a picture of your hand-written proof is okay. You should be able to insert an image in Jupyter by clicking Edit > Insert Image.)

In [ ]:
 

7. Effective comparison test

Suppose $a_k$ and $b_k$ are sequences satisfying $|a_k| \leq b_k$. The Comparison test tells us that if $B=\sum_{k=0}^\infty b_k$ exists as a finite number, then so does $A=\sum_{k=0}^\infty a_k$. Furthermore, this test tells us that $|A| \leq B$. We will see that if the limit of the series $\sum_{k=0}^\infty b_k$ is known, then the limit of $\sum_{k=0}^\infty a_k$ can be estimated.

Suppose $A = \sum_{k=0}^\infty a_k$ and $B = \sum_{k=0}^\infty b_k$. We're assuming $B$ is known and we want to estimate $A$. Use $S^a$ and $S^b$ to denote the partial sums $$S^a_n = \sum_{k=0}^{n-1} a_k \quad \text{and} \quad S^b_n = \sum_{k=0}^{n-1} b_k.$$ We'll similarly denote the remainders by $$R^a_n = A - S^a_n = \sum_{k=n}^\infty a_k \quad \text{and} \quad R^b_n = B - S^b_n = \sum_{k=n}^\infty b_k.$$ We can apply the Comparison test to these remainder formulas to see that $|R^a_n| \leq R^b_n$. Therefore, for each $n$ we can confine $A$ to a closed interval: $$A = S^a_n + R^a_n \in [S^a_n - R^b_n, S^a_n + R^b_n].$$ The remainder $R^b_n=B-S^a_n$ so, we have $$ \star \qquad A \in [S^a_n - B + S^b_n, S^a_n + B - S^b_n]. \qquad $$

Write a function comparison_limit(a, b, B, epsilon) which returns the partial $S^a_n$ where $n$ is the smallest number such that the interval described in $\star$ has width less than $\epsilon$. Here a and b are functions defining sequences. The variable B stores the actual limit of the series $\sum_{k=0}^\infty b_k$ as above. The variable epsilon stores a positive number.

Remark: This approach gives guaranteed control of the truncation error!

In [ ]:
 

Tests for your code:

The Fibonacci sequence is defined inductively by $a_0=1$, $a_1=1$ and $a_{k+2}=a_{k+1}+a_k$ for $k \geq 0$. It can be implemented in Python as below:

In [ ]:
from functools import lru_cache

@lru_cache
def fibonacci(k):
    if k <= 1:
        return 1
    else:
        return fibonacci(k-1) + fibonacci(k-2)

(Remark: The @lru_cache decorator is not necessary but stores the results in memory, making the recusive code quicker. You can read more about about this Python feature on the Geeks for Geeks website.)

The Fibonacci sequence grows exponentially at a rate asymptotic to the golden mean $\varphi = \frac{1 + \sqrt{5}}{2}$. It can be shown that the Fibonnacci sequence satisfies that $$\varphi^{n-1} \leq a_n$$ for all $n$.

The inverse Fibonacci series is $\sum_{k=0}^\infty 1/a_k$. We can see it is bounded above by a geometric series: $$\frac{1}{a_n} \leq c r^n \quad \text{where $c=\varphi$ and $r=1/\varphi$}.$$ From our formula for the limit of the geometric series, we can see the limit of this geometric series is $$L =\sum_{n=0}^\infty c r^k = \frac{c}{1-r} = \frac{\varphi^2}{\varphi-1}.$$

The golden mean, this limit L, and the inverse Fibonacci sequence are coded below.

In [ ]:
golden_mean = (1+m.sqrt(5))/2

L = golden_mean**2 / (golden_mean-1)

def inverse_fibonacci(k):
    return 1/fibonacci(k)

Given your work above, we should be able to estimate the value taken by the inverse Fibonacci series to an acuracy of $10^{-j}$ using the following code:

In [ ]:
for j in range(9):
    limit = comparison_limit(inverse_fibonacci, 
                             geometric_sequence(golden_mean, 1/golden_mean),
                             L,
                             10**-j)
    print(f'Acurate to {j+1} digits: {limit}')

My code returns:

Acurate to 1 digits: 3.0333333333333337
Acurate to 2 digits: 3.3304690407631585
Acurate to 3 digits: 3.3572331499135393
Acurate to 4 digits: 3.3594986693502227
Acurate to 5 digits: 3.359850770755157
Acurate to 6 digits: 3.3598825197189863
Acurate to 7 digits: 3.359885382521269
Acurate to 8 digits: 3.3598856248487095
Acurate to 9 digits: 3.3598856625106417

If your code returns a similar printout, then it should be correct.

There is a Wikipedia page devoted to the limit.

8. Zeta of 3

For $s > 1$, the Reimann zeta function is given by the convergent series $$\zeta(s)=\sum_{k=1}^\infty \frac{1}{k^s}=\frac{1}{1^s}+\frac{1}{2^s}+\frac{1}{3^s}\ldots.$$ In the Series notebook, we calculated that $\zeta(2)=\frac{\pi^2}{6}$. Now we will approximate $\zeta(3)$.

Observe that $$\frac{1}{(k+3)^3} \leq \frac{1}{(k+1)(k+2)(k+3)}$$ for all $k \geq 0$. We already know the value of the series $$ \sum_{n=0}^\infty \frac{1}{(k+1)(k+2)(k+3)} $$ from work above. (The limit of this series can be rigorously proven

Therefore we can use our function comparison_limit to estimate the infinite sum $$\frac{1}{3^3} + \frac{1}{4^3} +\ldots = \sum_{n=0}^\infty \frac{1}{(n+3)^3}.$$ Note that $\zeta(3)$ is this limit plus the missing $9/8$ because of the two missing terms.

Write a function zeta_3(epsilon) which takes as input an $\epsilon>0$. The function should return a number which is within $\frac{\epsilon}{2}$ of $\zeta(3)$ (so that $\zeta(3)$ is confined to an interval of length $\epsilon$, $[x-\frac{\epsilon}{2}, x+\frac{\epsilon}{2}]$ where $x$ is the returned value).

You may use your work on prior problems. You should make use of comparison_limit for example.

In [ ]:
 

Tests for your code:

The following code should print $\zeta(3)$ to increasing accuracy.

In [ ]:
for i in range(1,10):
    ans = zeta_3(10**(1-i))
    ans_str = f'{ans:.{i-1}f}'
    print(f'zeta(3) to {i:2} digits is {ans_str}')

When run, more and more digits of the same number should appear above.