{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Homework 3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"INSERT YOUR NAME HERE"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Directions:__ Add work to this notebook to solve the problems below.\n",
"\n",
"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.\n",
"\n",
"Before you submit your work, make sure it runs correctly without errors. Click “Kernel > Restart & Run All” to restart and run everything.\n",
"\n",
"Some problems are derived from problems in these books:\n",
"* __LL__: *Programming for Computations - Python* by Svein Linge and Hans Petter Langtangen, 2nd edition.\n",
"* __L__: *A Primer on Scientific Programming with Python* by Hans Petter Langtangen, 2nd edition.\n",
"* __TAK__: *Applied Scientific Computing With Python* by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Standard imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import math as m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Partial sums"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is often natural to implement a sequence as a function. For example, we might consider the sequence\n",
"$$a_k = \\frac{1}{k+1}$$\n",
"This can be implemented in Python as:\n",
"```python\n",
"def a(k):\n",
" return 1/(k+1)\n",
"```\n",
"Write a function `partial_sum(f, n)` which takes as input a function $f$ and an integer $n \\geq 0$ and returns the sum\n",
"$$\\sum_{k=0}^{n-1} f(k).$$ \n",
"If $n=0$, the function should return zero.\n",
"The function $f$ represents a sequence and should be assumed to take as input an integer and return a number."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we test that you are computing $\\sum_{k=0}^{n-1} 1$ correctly:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def b(k):\n",
" return 1\n",
"ans = partial_sum(b, 100)\n",
"if ans ==100:\n",
" print('That is correct, the answer is 100.')\n",
"else:\n",
" print(f'That is incorrect, the answer is 100 not {ans}.')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we test the result with the `a(k)` function from the problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def a(k):\n",
" return 1/(k+1)\n",
"ans = partial_sum(a, 1000)\n",
"if abs(ans - 7.4854) < 10**-4:\n",
" print('That is correct, the answer is approximately 7.4854.')\n",
"else:\n",
" print(f'That is incorrect, the answer is approximately 7.4854 not {ans}.')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Sum until small"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As in the previous problem *Partial sums*, we will sum sequences given by functions.\n",
"\n",
"Write a function `sum_until_small(f, epsilon)` which returns\n",
"$$\\sum_{k=0}^{n-1} f(k)$$\n",
"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.\n",
"\n",
"**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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def c(k):\n",
" return 1/2**k\n",
"for i in range(10):\n",
" assert sum_until_small(c, 1/2**i) == 2 - 1/2**i, \\\n",
" f'Failure when i={i}'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Also `sum_until_small(c, 2)` should return zero since the first term is smaller than $2$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sum_until_small(c, 2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Geometric series"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A geometric series is a series of the form\n",
"$$c + cr + cr^2 + cr^3+\\ldots=\\sum_{k=0}^\\infty c r^k.$$\n",
"\n",
"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.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f = geometric_sequence(3, 1/4) # The sequence f(n) = 3/4^n\n",
"assert f(0) == 3\n",
"assert f(1) == 3/4\n",
"assert f(2) == 3/16"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you have already done the partial sum problem, you can check the formula\n",
"$$\\sum_{k=0}^\\infty c r^k=\\frac{c}{1-r}.$$\n",
"In $f$ as above, we have $c=3$ and $r=1/4$, so the sum should converge to four.\n",
"\n",
"The following should print `True`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"abs(partial_sum(f, 100) - 4) < 10**-15"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. Guess the limit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Print some partial sums for the series\n",
"$$\\sum_{n=0}^\\infty \\frac{1}{(n+1)(n+2)(n+3)}.$$\n",
"You may use the functions you have already written in this assignment.\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5. Guess the remainders"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This problem also concerns the series\n",
"$$\\sum_{k=0}^\\infty \\frac{1}{(k+1)(k+2)(k+3)}.$$\n",
"Let $L$ be the limit of the series. The remainders are the differences with the partial sums\n",
"$$R_n = L-S_n = L - \\sum_{n=0}^{n-1} \\frac{1}{(k+1)(k+2)(k+3)}.$$\n",
"These quantities are computed and stored in a list below and printed out, assuming your `partial_sum` code and `guess_for_limit` was correct."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def a(k):\n",
" return 1/((k+1)*(k+2)*(k+3))\n",
"R = 10*[None]\n",
"for n in range(10):\n",
" R[n] = guess_the_limit - partial_sum(a, n)\n",
" print(f'R_{n} = {R[n]}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Hint*: Consider the recipricals of the remainders. Try to look for patterns."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Test if remainder guess works\n",
"for n in range(10):\n",
" numerator,denominator = remainder_guess(n)\n",
" assert type(numerator) == int, 'The numerator should be an int'\n",
" assert type(denominator) == int, 'The denominator should be an int'\n",
" if numerator/denominator - R[n] > 10**-8:\n",
" print(f'Your guess seems to be incorrect when n={n}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6. Inductive proof"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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 \n",
"$$S_n = \\sum_{k=0}^{n-1} \\frac{1}{(k+1)(k+2)(k+3)}$$\n",
"are given by $L-\\frac{p_n}{q_n}$. Prove this by induction.\n",
"\n",
"It follows that the series sums to $L$ once you observe that $\\lim_{n \\to \\infty} \\frac{p_n}{q_n}=0$.\n",
"\n",
"*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`.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 7. Effective comparison test"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"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\n",
"$$S^a_n = \\sum_{k=0}^{n-1} a_k \\quad \\text{and} \\quad S^b_n = \\sum_{k=0}^{n-1} b_k.$$\n",
"We'll similarly denote the remainders by \n",
"$$R^a_n = A - S^a_n = \\sum_{k=n}^\\infty a_k \n",
"\\quad \\text{and} \\quad \n",
"R^b_n = B - S^b_n = \\sum_{k=n}^\\infty b_k.$$\n",
"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:\n",
"$$A = S^a_n + R^a_n \\in [S^a_n - R^b_n, S^a_n + R^b_n].$$\n",
"The remainder $R^b_n=B-S^a_n$ so, we have\n",
"$$\n",
"\\star \\qquad A \\in [S^a_n - B + S^b_n, S^a_n + B - S^b_n]. \\qquad\n",
"$$\n",
"\n",
"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.\n",
"\n",
"*Remark:* This approach gives guaranteed control of the truncation error!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from functools import lru_cache\n",
"\n",
"@lru_cache\n",
"def fibonacci(k):\n",
" if k <= 1:\n",
" return 1\n",
" else:\n",
" return fibonacci(k-1) + fibonacci(k-2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(*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](https://www.geeksforgeeks.org/python-functools-lru_cache/).)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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\n",
"$$\\varphi^{n-1} \\leq a_n$$ for all $n$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The inverse Fibonacci series is $\\sum_{k=0}^\\infty 1/a_k$. We can see it is bounded above by a geometric series:\n",
"$$\\frac{1}{a_n} \\leq c r^n \\quad \\text{where $c=\\varphi$ and $r=1/\\varphi$}.$$\n",
"From our formula for the limit of the geometric series, we can see the limit of this geometric series is\n",
"$$L =\\sum_{n=0}^\\infty c r^k = \\frac{c}{1-r} = \\frac{\\varphi^2}{\\varphi-1}.$$\n",
"\n",
"The golden mean, this limit L, and the inverse Fibonacci sequence are coded below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"golden_mean = (1+m.sqrt(5))/2\n",
"\n",
"L = golden_mean**2 / (golden_mean-1)\n",
"\n",
"def inverse_fibonacci(k):\n",
" return 1/fibonacci(k)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for j in range(9):\n",
" limit = comparison_limit(inverse_fibonacci, \n",
" geometric_sequence(golden_mean, 1/golden_mean),\n",
" L,\n",
" 10**-j)\n",
" print(f'Acurate to {j+1} digits: {limit}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"My code returns:\n",
"```\n",
"Acurate to 1 digits: 3.0333333333333337\n",
"Acurate to 2 digits: 3.3304690407631585\n",
"Acurate to 3 digits: 3.3572331499135393\n",
"Acurate to 4 digits: 3.3594986693502227\n",
"Acurate to 5 digits: 3.359850770755157\n",
"Acurate to 6 digits: 3.3598825197189863\n",
"Acurate to 7 digits: 3.359885382521269\n",
"Acurate to 8 digits: 3.3598856248487095\n",
"Acurate to 9 digits: 3.3598856625106417\n",
"```\n",
"If your code returns a similar printout, then it should be correct.\n",
"\n",
"There is a [Wikipedia page devoted to the limit](https://en.wikipedia.org/wiki/Reciprocal_Fibonacci_constant)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 8. Zeta of 3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $s > 1$, the Reimann zeta function is given by the convergent series\n",
"$$\\zeta(s)=\\sum_{k=1}^\\infty \\frac{1}{k^s}=\\frac{1}{1^s}+\\frac{1}{2^s}+\\frac{1}{3^s}\\ldots.$$\n",
"In the `Series` notebook, we calculated that $\\zeta(2)=\\frac{\\pi^2}{6}$. Now we will approximate $\\zeta(3)$.\n",
"\n",
"Observe that\n",
"$$\\frac{1}{(k+3)^3} \\leq \\frac{1}{(k+1)(k+2)(k+3)}$$\n",
"for all $k \\geq 0$. We already know the value of the series\n",
"$$\n",
"\\sum_{n=0}^\\infty \\frac{1}{(k+1)(k+2)(k+3)}\n",
"$$\n",
"from work above. (The limit of this series can be rigorously proven\n",
"\n",
"Therefore we can use our function `comparison_limit` to estimate the infinite sum\n",
"$$\\frac{1}{3^3} + \\frac{1}{4^3} +\\ldots = \\sum_{n=0}^\\infty \\frac{1}{(n+3)^3}.$$\n",
"Note that $\\zeta(3)$ is this limit plus the missing $9/8$ because of the two missing terms.\n",
"\n",
"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).\n",
"\n",
"You may use your work on prior problems. You should make use of `comparison_limit` for example."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code should print $\\zeta(3)$ to increasing accuracy."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for i in range(1,10):\n",
" ans = zeta_3(10**(1-i))\n",
" ans_str = f'{ans:.{i-1}f}'\n",
" print(f'zeta(3) to {i:2} digits is {ans_str}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When run, more and more digits of the same number should appear above."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}