diff options
Diffstat (limited to 'usth/MATH1.4/cons-sequences/README.md')
| -rw-r--r-- | usth/MATH1.4/cons-sequences/README.md | 296 |
1 files changed, 296 insertions, 0 deletions
diff --git a/usth/MATH1.4/cons-sequences/README.md b/usth/MATH1.4/cons-sequences/README.md new file mode 100644 index 0000000..493ae32 --- /dev/null +++ b/usth/MATH1.4/cons-sequences/README.md @@ -0,0 +1,296 @@ +# Infinite Sequences: A Case Study in Functional Python + +In this article, we will only consider sequences defined by a function whose +domain is a subset of the set of all integers. Such sequences will be +*visualized*, i.e. we will try to evaluate the first few (thousand) elements, +using functional programming paradigm, where functions are more similar to the +ones in math (in contrast to imperative style with side effects confusing to +inexperenced coders). The idea is taken from +[subsection 3.5.2 of SICP](https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-24.html#%_sec_3.5.2) +and adapted to Python, which, compare to Scheme, is significantly more popular: +Python is pre-installed on almost every modern Unix-like system, namely macOS, +GNU/Linux and the \*BSDs; and even at MIT, the new 6.01 in Python has recently +replaced the legendary 6.001 (SICP). + +One notable advantage of using Python is its huge **standard** library. For +example the "identity sequence" (sequence defined by the identity function) can +be imported directly from `itertools`: + +```python +>>> from itertools import count +>>> positive_integers = count(start=1) +>>> next(positive_integers) +1 +>>> next(positive_integers) +2 +>>> for _ in range(4): next(positive_integers) +... +3 +4 +5 +6 +``` + +To open a Python emulator, simply lauch your terminal and run `python`. If that +is somehow still too struggling, navigate to +[the interactive shell on Python.org](https://www.python.org/shell/). + +*Let's get it started* with somethings everyone hates: recursively defined +sequences, e.g. the famous Fibonacci (*F*<sub>*n*</sub> = *F*<sub>*n*-1</sub> + +*F*<sub>*n*-2</sub>, *F*<sub>1</sub> = 1, *F*<sub>0</sub> = 0). Since +[Python doesn't support](http://neopythonic.blogspot.com/2009/04/final-words-on-tail-calls.html) +[tail recursion](https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-11.html#call_footnote_Temp_48), +it's generally **not** a good idea to define anything recursively (which is, +ironically, the only trivial *functional* solution in this case) but since we +will only evaluate the first few terms (use the **Tab** key to indent the line +when needed): + +```python +>>> def fibonacci(n, a=0, b=1): +... # To avoid making the code look complicated, n < 0 is not handled here. +... return a if n == 0 else fibonacci(n - 1, b, a + b) +... +>>> fibo_seq = (fibonacci(n) for n in count(start=0)) +>>> for _ in range(7): next(fibo_seq) +... +0 +1 +1 +2 +3 +5 +8 +``` + +<details><summary>Note (click to expand)</summary> +The <code>fibo_seq</code> above is just to demonstrate how +<code>itertools.count</code> can be use to create an infinite sequence defined +by a function. For better performance, this should be used instead + +```python +def fibonacci_sequence(a=0, b=1): + yield a + yield from fibonacci_sequence(b, a + b) +``` +</details> + +It is noticable that the elements having been iterated through (using `next`) +will disappear forever in the void (oh no!), but that is the cost we are +willing to pay to save some memory, especially when we need to evaluate a +member of (arbitrarily) large index to estimate the sequence's limit. One case +in point is estimating a definite integral using +[left Riemann sum](https://en.wikipedia.org/wiki/Riemann_sum#Left_Riemann_sum): + +```python +>>> def integral(f, a, b): +... def leftRiemannSum(n): +... dx = (b-a) / n +... def x(i): return a + i*dx +... return sum(f(x(i)) for i in range(n)) * dx +... return leftRiemannSum +... +``` + +The function `integral(f, a, b)` as defined above returns a function taking *n* +as an argument. As *n* approaches ∞, its result approaches the value of the +integral of *f* on [*a, b*]. For example, we are going to estimate π as the +area of a semicircle whose radius is sqrt(2): + +```python +>>> from math import sqrt +>>> def semicircle(x): return sqrt(abs(2 - x*x)) +... +>>> pi = integral(semicircle, -sqrt(2), sqrt(2)) +>>> pi_seq = (pi(n) for n in count(start=2)) +>>> for _ in range(3): next(pi_seq) # the first few aren't quite close +... +2.000000029802323 +2.514157464087051 +2.7320508224700384 +``` + +At index around 1000, the result is somewhat acceptable: + + 3.1414873191059525 + 3.1414874770617427 + 3.1414876346231577 + +Since we are comfortable with sequence of sums, let's move on to sums of a +sequence, which are called series. For estimation, again, we are going to make +use of infinite sequences of partial sums, which are implemented as +`itertools.accumulate` by thoughtful Python developers. +[Geometric](https://en.wikipedia.org/wiki/Geometric_series) and +[*p*-series](https://math.oregonstate.edu/home/programs/undergrad/CalculusQuestStudyGuides/SandS/SeriesTests/p-series.html) +can be defined as follow: + +```python +>>> from itertools import accumulate as partial_sums +>>> def geometric_series(r, a=1): return partial_sums(a*r**n for n in count(0)) +... +>>> def p_series(p): return partial_sums(1 / n**p for n in count(1)) +... +``` + +We can then use these to determine whether the series is convergent or +divergent. For instance, the fact that *p*-series with *p* = 2 converges to +π<sup>2</sup>/6 ≈ 1.6449340668482264 can be verified via + +```python +>>> s = p_series(p=2) +>>> for _ in range(11): next(s) +... +1.0 +1.25 +1.3611111111111112 +1.4236111111111112 +1.4636111111111112 +1.4913888888888889 +1.511797052154195 +1.527422052154195 +1.5397677311665408 +1.5497677311665408 +1.558032193976458 +``` + +We can observe that it takes quite a lot of steps to get the precision we would +generally expect (*s*<sub>11</sub> is only precise to the first decimal place; +second decimal places: *s*<sub>101</sub>; third: *s*<sub>2304</sub>). Luckily, +many techniques for series acceleration are available. +[Shanks transformation](https://en.wikipedia.org/wiki/Shanks_transformation), +for instance, can be implemented as follow: + +```python +>>> from itertools import islice, tee +>>> def shanks(seq): +... return map(lambda x, y, z: (x*z - y*y) / (x + z - y*2), +... *(islice(t, i, None) for i, t in enumerate(tee(seq, 3)))) +... +``` + +In the code above, `lambda x, y, z: (x*z - y*y) / (x + z - y*2)` denotes the +anonymous function (*x*, *y*, *z*) ↦ (*xz*-*y*<sup>2</sup>)/(*x*+*z*-2*y*) and +`map` is a higher order function applying that function to respective elements +of subsequences starting from index 1, 2, 3 of `seq`. On Python 2, one should +import `imap` from `itertools` to get the same +[lazy](https://en.wikipedia.org/wiki/Lazy_evaluation) behavior of `map` on +Python 3. + +```python +>>> s = shanks(p_series(2)) +>>> for _ in range(10): next(s) +... +1.4500000000000002 +1.503968253968257 +1.53472222222223 +1.5545202020202133 +1.5683119658120213 +1.57846371882088 +1.5862455815659202 +1.5923993101138652 +1.5973867787856946 +1.6015104548459742 +``` + +The result was quite satisfying, yet we can do one step futher by continuously +applying the transformation to the sequence: + +```python +>>> def compose(transform, seq): +... yield next(seq) +... yield from compose(transform, transform(seq)) +... +>>> s = compose(shanks, p_series(2)) +>>> for _ in range(10): next(s) +... +1.0 +1.503968253968257 +1.5999812811165188 +1.6284732442271674 +1.6384666832276524 +1.642311342667821 +1.6425249569252578 +1.640277484549416 +1.6415443295058203 +1.642038043478661 +``` + +Shanks transformation works on every sequence (not just sequences of partial +sums). Back to previous example of using left Riemann sum to compute definite integral: + +```python +>>> pi_seq = compose(shanks, map(pi, count(2))) +>>> for _ in range(10): next(pi_seq) +... +2.000000029802323 +2.978391111182236 +3.105916845397819 +3.1323116570377185 +3.1389379264270736 +3.140788413965646 +3.140921512857936 +3.1400282163913436 +3.1400874774021816 +3.1407097229603256 +>>> next(islice(pi_seq, 300, None)) +3.1415061302492413 +``` + +Now having series defined, let's see if we can learn anything +about power series. Sequence of partial sums of power series +∑*c*<sub>*n*</sub>(*x*-*a*)<sup>*n*</sup> can be defined as + +```python +>>> from operator import mul +>>> def power_series(c, start=0, a=0): +>>> return lambda x: partial_sums(map(mul, c, (x**n for n in count(start)))) +... +``` + +We can use this to compute functions that can be written as +[Taylor series](https://en.wikipedia.org/wiki/Taylor_series): + +```python +>>> from math import factorial +>>> def exp(x): return power_series(1/factorial(n) for n in count(0))(x) +... +>>> def cos(x): +... c = ((1 - n%2) * (1 - n%4) / factorial(n) for n in count(0)) +... return power_series(c)(x) +... +>>> def sin(x): +... c = (n%2 * (2 - n%4) / factorial(n) for n in count(1)) +... return power_series(c, start=1)(x) +... +``` + +Amazing! Let's test 'em! + +```python +>>> e = compose(shanks, exp(1)) # this should converges to 2.718281828459045 +>>> for _ in range(4): next(e) +... +1.0 +2.749999999999996 +2.718276515152136 +2.718281825486623 +``` + +Impressive, huh? For sine and cosine, series acceleration is not even necessary: + +```python +>>> from math import pi as PI +>>> s = sin(PI/6) +>>> for _ in range(5): next(s) +... +0.5235987755982988 +0.5235987755982988 +0.49967417939436376 +0.49967417939436376 +0.5000021325887924 +>>> next(islice(cos(PI/3), 8, None)) +0.500000433432915 +``` + +[](http://creativecommons.org/licenses/by-sa/4.0/) +This work is licensed under a +[Creative Commons Attribution-ShareAlike 4.0 International License](http://creativecommons.org/licenses/by-sa/4.0/). |