(Answer)PMS1A-Sequence.ipynb

Kernel: Python 3 (Anaconda 5)

数列与 $\pi$ 的探究

课程介绍

目標

本节的 python + Math ，带大家来算算数列与用数列来探索 $\pi$。

数列就是体现程式“重复运算”的特性，要使用数列就需要理解抽象出“规律形式”的能力。

本节先从 Python 的 list 结构来探索数列，最后来用以下数列观察数列如何逼近 $\pi$。

$$\frac{\pi}{4} = \frac{1}{1} - \frac{1}{3} + \frac{1}{5}-\frac{1}{7} + \frac{1}{9} \cdots$$$$\frac{\pi^2}{6} =1+ \frac{1}{2^2} + \frac{1}{3^2} + \cdots +\frac{1}{n^2}+\cdots$$$$\frac{\pi}{2} =\frac{2}{1} \cdot \frac{2}{3} \cdot \frac{4}{3}\cdot \frac{4}{5} \cdot \frac{6}{5} \cdots$$

課程架構

0. 前言与故事： $e = 1+ \frac{1}{1!} + \frac{1}{2!}+ \frac{1}{3!}+\cdots$

1. list 与数列：[15,5,16,16,28,32,51,38,26]

2. range 与等差数列： $1^3+ 2^3 + 3^3 + 4^3 =(1+2+3+4)^2$

3. for-loop 与三角垛数列： $1 + (1+2) + (1+2+3) + (1+2+3+4)$

4. import math 与 pi 的探讨 ：$\frac{\pi}{2} =\frac{2}{1} \cdot \frac{2}{3} \cdot \frac{4}{3}\cdot \frac{4}{5} \cdot \frac{6}{5} \cdots$

核心知識

• 語法

  • 熟練 list 的用法： [2i+1 for i in range(1,n)]

  • 熟練對 list 求和： sum([i for i in range(1,101) ])

  • 了解格式化輸出： "1/{0:5} = {1:9.7f}".format(i,1.0/i)

  • 引入函式庫：import math,

• 數學

  • 能製造等差數列並求其和

  • 能了解無窮數列有收斂與發散的區別

  • pi 的數列近似表達

• 思維與能力

  • 能觀察數列的關係，並寫出其表示式

  • 能分辨遞迴與迴圈的分別並了解其使用時機

課程內容

0 前言与故事

0.a Lotus 123

本次的 Python ＋ Math 要来谈数列。操作数列是我早期接触电脑的回忆。在 14 岁时(1995)时，有了一台黑白的二手电脑，可能是 386 或 486 的电脑，电脑有个试算表软体叫 Lotus123 ，那时电脑也没什么游戏，我只是用试算表来做些计算。

0.b 神秘有趣的数学

那时还蛮喜欢跑书店，就是员林车站对面的金石堂书店，我还记得二楼有个书架，摆放些科普书，关于数学的书也没几本，但其中有本很吸引国二时的注意，就是这本孙文先编译的《神秘有趣的数学》。

里面列举很多数学的各方面的知识，也没什么证明，但就提些比较特别的结论。有章节提到了数列，谈到

$$1+ \frac{1}{2} + \frac{1}{2^2} + \cdots +\frac{1}{2^n}+\cdots =2$$

又提到

$$1+ \frac{1}{2} + \frac{1}{3} + \cdots +\frac{1}{n}+\cdots \to \infty$$

0.C 用 Loutus 初次探索数列

如同这书名，这对 14 岁的我就是感到神秘有趣，那时让我印象深刻的还有

$$\frac{\pi^2}{6} =1+ \frac{1}{2^2} + \frac{1}{3^2} + \cdots +\frac{1}{n^2}+\cdots$$

竟然左式的无理数会等于右式的有理数之和！我当时拿到 Lotus 时，就想到可以去验证他，那时看到荧幕上出现结果与 $\frac{\pi^2}{6}$ 很靠近时，也真的很感动。有些事，虽然你理智上已经觉得他是正确，但实际操作见证时，就还是有不同层次的感觉。

算完平方的调和数列后，我就开始探索以下这数列：

$$\frac{1}{1!}+ \frac{1}{2!} + \frac{1}{3!} + \cdots +\frac{1}{n!}+\cdots$$

竟然发现这串数列会靠近 $1.718281828...$ 而不是发散的情况，就感到很兴奋，我还记得隔天是周日早上，就拿着这自以为是新发现的成果去找数学老师。 现在想想也蛮好笑的，事隔 24 年这段天真的过程，还烙印在我脑海中。

1. list 与数列：[15,5,16,16,28,32,51,38,26]

1.a list 为一堆有顺序的资料：[1,2,3,4,5]

在 Python 中，list 是一串有顺序的资料，要制造 list 的方法很简单，就是把一串资料用 [ ] 包住，例如 [1,2,3,4,5]。

而 list 的结构，有包含些基本运算，例如 sum 求总和。

• 请执行以下程式码，并观察其结果

a = [1,2,3,4,5]
sum(a)

• 执行完后，更新数列为 1,4,9,16,25 并重新执行观察其结果。

a = [1,4,9,16,25]
sum(a)

55

• 以下为 1984 年后，中国队每届奥运金牌数量，请计算其总和。

15,5,16,16,28,32,51,38,26

gold = [15,5,16,16,28,32,51,38,26]
sum(gold)

227

1.b 计算数列平均

要求得数列的项数可以用 len 指令

• 请用 len(gold) 求上述数列共几项

len(gold)

9

在计算完总和后,要求平均,可以用 sum(a)/len(a).

其中 len 表示数列的长度(length), 也就是项数。

• 请计算 1984 年后，中国对奥运金牌的平均获奖数。

gold = [15,5,16,16,28,32,51,38,26]
sum(gold)/len(gold)

25.22222222222222

除了 sum, len 外，也有 max, min

max(gold)

51

也有个 sorted 函数可以排序

sorted(gold)

[5, 15, 16, 16, 26, 28, 32, 38, 51]

• 请计算 1,4,9,16,25 的平均

• 计算前请估算平均是否大於中位数 9

seq = [1,4,9,16,25]
print(sum(seq)/len(seq))

11.0

1.c 调整输出格式

当有多资料要输出时可利用 print

请分别输出以下两个数列的总和：

c1: 1, 2, 3, 4, 5

c2: 1, 8,27,64,125

c1 = [1,2,3,4,5]
c2 = [1,8,27,64,125]
print(sum(c1))
print(sum(c2))

15
225

可以在一个 print 内，输出多个资料，并用 , 分隔。 例如：

c1 = [1,2,3,4,5]
print(c1,sum(c1))

c1 = [1,2,3,4,5]
c2 = [1,8,27,64,125]
print(c1,sum(c1),sum(c1)/len(c1))
print(c2,sum(c2),sum(c2)/len(c2))

[1, 2, 3, 4, 5] 15 3.0
[1, 8, 27, 64, 125] 225 45.0

要在上述的输出再加入文字，可利用 format 的格式

c1 = [1,2,3,4,5]
c2 = [1,8,27,64,125]
print("数列 {} 的总和为 {} 平均为 {}".format(c1,sum(c1),sum(c1)/len(c1)))
print("数列 {} 的总和为 {} 平均为 {}".format(c2,sum(c2),sum(c2)/len(c2)))

数列 [1, 2, 3, 4, 5] 的总和为 15 平均为 3.0
数列 [1, 8, 27, 64, 125] 的总和为 225 平均为 45.0

1.d 取得数列的个别项与子数列

利用 index 可以取得某一元素在数列的位置(index)。

注意 list 的第一项的 index 为 0。

• 请找出 gold 的最大值在 list 的位置

• 请利用 1984 为第一项对应的年份，则最大值发生在那一年。（奥运四年举办一次）

gold = [15,5,16,16,28,32,51,38,26]
print(gold.index(max(gold)))
print(gold.index(max(gold))*4+1984)

6
2008

若要取得数列 xxx 的第 k 项可用 xxx[k]

• 请执行以下程式，观察其 index 与这数在数列的位置

fib = [1,1,2,3,5,8,13,21,34,55,89]
print(fib[0],fib[5],fib[10])

fib = [1,1,2,3,5,8,13,21,34,55,89]
print(fib[0],fib[5])
print(fib[10])

1 8
89

利用 [-1] 则可以由后往前数

若 index 中放入不在 index 的范围内的数 则会有错误讯息。

fib = [1,1,2,3,5,8,13,21,34,55,89]
print(fib[-1],fib[-2],fib[-10])
print(fib[1000])

89 55 1

【习题1a】数列的建立

以下哪些是 Python 内可以正确建立 list 的方式。

sa = {1, 3, 6, 10, 15, 21}
sa = {1; 3; 6; 10; 15; 21}
sa = [1, 3, 6, 10, 15, 21]
sa = [1 3 6 10 15 21]

sa = [1, 3, 6, 10, 15, 21]
print(sa)

[1, 3, 6, 10, 15, 21]

【习题1b】计算平均

请完成以下程式来计算 [1,2,4,8,16,32,64] 的数列的总和与平均。

sb = [1,2,4,8,16,32,64]
print("{} 的总和为 {} 的平均为 {}".format(sb, sum(sb), sum(sb)/len(sb)))

[1, 2, 4, 8, 16, 32, 64] 的总和为 127 的平均为 18.142857142857142

【习题1c】判断 index

请先想想以下程式码执行的结果。再执行一次程式来验证

spi = [3,1,4,1,5,9,2,6,5,3,5]
print(spi[2],"=", 4)
print(spi[-2],"=", 3)
print(len(spi),"=",11)

4 = 4
3 = 3
11 = 11

【习题1d】奥运会金牌

以下为美国在 1984年后奥运会获得金牌的数量

83,36,37,44,37,36,36,46,46

请计算其平均获奖数，并列出其在 2008 年的获奖数

goldA =[83,36,37,44,37,36,36,46,46]
print(sum(goldA)/len(goldA))
print(len(goldA))

44.55555555555556
9

2 利用 for 製造有规律数列

2.a range 的使用

当我们想要运算 1+2+3+...+100 时，若直接输入 [1,2,3,4,...,100] 也太麻烦。

在 python 就要善用 range 这指令，你尝试输入以下指令来观看其结果

a1 = [k for k in range(1,10)]
a2 = [k**2 for k in range(1,10)]
print("数列",a1,"总和为",sum(a1))
# print("数列",???,"总和为",???)

---

• 注意这数列的首项、末项为多少？

• 这数列共有几项？

a1 = [k for k in range(1,10)]
a2 = [k**2 for k in range(1,10)]
print("数列",a1,"总和为",sum(a1))
print("数列",a2,"总和为",sum(a2))

数列 [1, 2, 3, 4, 5, 6, 7, 8, 9] 总和为 45
数列 [1, 4, 9, 16, 25, 36, 49, 64, 81] 总和为 285

在 for 前面的式子，就是數學上的一般項。

当要输出偶数数列时就用 2*k 来表示 要输出 2的次方数列时就用 2**k 来表示

a3 = [ 2*k for k in range(1,11)] 
a4 = [ 2**k for k in range(1,11)] 

a3 = [ 2*k for k in range(1,11)] 
a4 = [ 2**k for k in range(1,11)] 
print(a3)
print(a4)

[2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
[2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]

2.b Range 的参数省略

若要产生奇数数列该如何制造呢？

以下哪些可以制造出 1,3,5,7,9呢？

b1 =[ 3*k for k in range(1,5)] 
b2 =[ 2*k-1 for k in range(1,6)] 
b3 =[ 2*k+1 for k in range(1,6)] 
b4 =[ 2*k+1 for k in range(5)] 

• 请执行程式前先想想

b1 =[ 3*k for k in range(1,5)] 
b2 =[ 2*k-1 for k in range(1,6)] 
b3 =[ 2*k+1 for k in range(1,6)] 
b4 =[ 2*k+1 for k in range(5)] 
print(b1)
print(b2)
print(b3)
print(b4)

[3, 6, 9, 12]
[1, 3, 5, 7, 9]
[3, 5, 7, 9, 11]
[1, 3, 5, 7, 9]

在上述数列中， 使用 range(5) 其实与 range(0,5) 相同。

也就是当 range 只有一个参数时，只规定上界，而初始值为 0

• 已知西元 1900 年為鼠年，请列出 1900 年以後，生肖屬鼠的出生年。

sx = [1900+12*k for k in range(10)]
print(sx)

[1900, 1912, 1924, 1936, 1948, 1960, 1972, 1984, 1996, 2008]

2.c Range 的第 3 个参数

对於先前的 1,3,5,7,9 初了使用 2k+1, 2k-1 来表示外，

还可使用 Range 的第三个参数来调整。

• 执行前先想想，以下哪些可产生 1,3,5,7,9

c1 = [k for k in range(1,9,2)]
c2 = [k for k in range(1,10,2)]
c3 = [k for k in range(1,11,2)]
c4 = [k+1 for k in range(0,10,2)]

c1 = [k for k in range(1,9,2)]
c2 = [k for k in range(1,10,2)]
c3 = [k for k in range(1,11,2)]
c4 = [k+1 for k in range(0,10,2)]
print(c1)
print(c2)
print(c3)
print(c4)

[1, 3, 5, 7]
[1, 3, 5, 7, 9]
[1, 3, 5, 7, 9]
[1, 3, 5, 7, 9]

利用 range 的第三个参数，请输出一个 1900 到 2050 的鼠年。

并统计共几个鼠年。

sx2 = [k for k in range(1900,2050,12)]
print(sx2)
print(len(sx2))

[1900, 1912, 1924, 1936, 1948, 1960, 1972, 1984, 1996, 2008, 2020, 2032, 2044]
13

2.d 計算 1+2+...+10000 会花很久吗？

• 計算 1~100 個連續整數和

d1 = [k for k in range(1,101)]
print(sum(d1))

d1 = [k for k in range(1,101)]
print(sum(d1))

5050

在上述求和中，其实计算机是逐项相加。

虽然计算机的计算速度很快，但项数很多时，还是要占很多空间，与时间。

• 试着将项数改为 $10^4, 10^5,\cdots, 10^8$ 来观察其计算时间。

d2 = [k for k in range(1,10**6+1)]
print(sum(d2))

500000500000

此时数学就扮演着重要角色

等差数列的和可用公式 $\frac{(1+n) n }{2}$ 来得到，

其推导的过程为『倒写相加』

$$\begin{align} S =& 1 +2 + 3 +\cdots + 100 \\ S =& 100+ 99 + 98 + \cdots + 1 \\ 2S = & 101 + 101 + 101 + \cdots + 101\\ S = & \dfrac{101\times 100}{2} = 5050 \end{align}$$

对电脑而言，利用公式来计算，会比直接生成数列快很多。

所以多些数学背景是可以大大提升程式的执行效率。

【试验】请调整不同的 n 值来计算以下等差的求和

n= 100000000
print(n*(n+1)/2)

n= 100000000
print(n*(n+1)/2)

5000000050000000.0

【补充】电影《丈量世界》中，有高斯计算 1+2+...+100 的课堂片段

[YouTube] Gauss 計算 1+2+...+100

【习题2a】等差数列

请试着用 list 的 for 来制造数列 1,4,7,10,...,31 并算其总和与平均

e2a = [1+3*k for k in range(11)]
print("数列 {} 的和為 {} 平均为 {}".format(e2a,sum(e2a),sum(e2a)/len(e2a)))

数列 [1, 4, 7, 10, 13, 16, 19, 22, 25, 28, 31] 的和為 176 平均为 16.0

【习题2b】龙年

请试着写出一个数列来列出出 1900~2050 的龙年？（注意，1900年为鼠年，那哪年为龙年呢？）

• 解答中给了三种表示法，你可比较哪种表示法适合此题的叙述。

sx1 = [k for k in range(1904,2051,12)]
sx2 = [8+12*k for k in range(158,170)]
sx3 = [1904+12*k for k in range((2051-1904)//12+1)]
print(sx1)
print(sx2)
print(sx3)

[1904, 1916, 1928, 1940, 1952, 1964, 1976, 1988, 2000, 2012, 2024, 2036, 2048]
[1904, 1916, 1928, 1940, 1952, 1964, 1976, 1988, 2000, 2012, 2024, 2036, 2048]
[1904, 1916, 1928, 1940, 1952, 1964, 1976, 1988, 2000, 2012, 2024, 2036, 2048]

【习题2c】哈雷彗星

哈雷彗星的轨道周期约为 76 年。此慧星在 1608 年时被观测过，则接下来 10 次被观测的年份约为哪些年。

halley1 = [1608+76*k for k in range(10)]
halley2 = [k for k in range(1608,1608+76*10,76)]
print(halley1)
print(halley2)

[1608, 1684, 1760, 1836, 1912, 1988, 2064, 2140, 2216, 2292]
[1608, 1684, 1760, 1836, 1912, 1988, 2064, 2140, 2216, 2292]

【习题2d】西洋棋与米粒

传说中，西洋棋的发明人向国王要了如下的奖赏

在第 $1$ 格放入 $1$粒米、在第 $2$ 格放入 $2$ 粒米、在第 $3$ 格放入 $2^2$ 粒米、 在第 $n$ 格放入 $2^n$ 粒米、

请计算所有米粒的数量和。

print(sum([2**k for k in range(64)]))

18446744073709551615

承上，假设每碗饭有 2000 颗米、每人一天吃了 3 碗饭，则这些米可供地球上 70 亿人吃多久？

riceNum = sum([2**k for k in range(64)])
riceDay = riceNum/(riceNum)
riceYear = riceDay/(2000*3*7*10**9)
print(riceDay)
print(riceYear)

1.0
2.380952380952381e-14

3 for-loop 与三角垛数列 1 + (1+2) + (1+2+3) + (1+2+3+4)

3.a for-loop 来观察 1+3+5+...+(2n-1)

在前一节，要观察一个数列在不同项数时的变化，需重覆输入并修正程式码

odd5 = [2*k+1 for k in range(5)]
print(odd5, sum(odd5))
odd6 = [2*k+1 for k in range(6)]
print(odd6, sum(odd6)
odd7 = [2*k+1 for k in range(7)] 
print(odd7, sum(odd7)

odd5 = [2*k+1 for k in range(5)]
print(odd5, sum(odd5))
odd6 = [2*k+1 for k in range(6)]
print(odd6, sum(odd6))
odd7 = [2*k+1 for k in range(7)] 
print(odd7, sum(odd7))

[1, 3, 5, 7, 9] 25
[1, 3, 5, 7, 9, 11] 36
[1, 3, 5, 7, 9, 11, 13] 49

上述程式码有很多重覆的结构，可用 for-loop 来简化。

请输入以下程式码

for n in [4,5,6,10]:
  odd = [2*k+1 for k in range(n)] 
  print(n, odd, sum(odd))

for n in [5,6,7]:
  odds = [2*k+1 for k in range(n)] 
  print(n, odds, sum(odds)) 

5 [1, 3, 5, 7, 9] 25
6 [1, 3, 5, 7, 9, 11] 36
7 [1, 3, 5, 7, 9, 11, 13] 49

3.b 三角垛

利用 for 回圈，可以来计算出以下三角垛的数量。

第一层是 1 个、

第二层 1+2 个、两层共 1+(1+2)= 4 个。

第三层 1+2+3 个、三层共 1+(1+2)+(1+2+3) = 10 个。

第四层 1+2+3+4 个、四层共 1+(1+2)+(1+2+3)+(1+2+3+4) = 20 个。

请问累积到第五层共多少个？

ans = 0
n = 5
for i in range(1,n+1):
    print(i, sum([k for k in range(1,i+1)]))
    ans += sum([k for k in range(1,i+1)])
print(ans)

1 1
2 3
3 6
4 10
5 15
35

上式中的 sum([k for k in range(1,i+1)])

其實就是等差數列 $1+2+3+...+ i = \frac{i\cdot(i+1)}{2}$

因此程式码可改写如下

ans = 0
n = 5
for i in range(1,n+1):
    print(i, n*(n+1)/2)
    ans +=  n*(n+1)/2
print(ans)

1 15.0
2 15.0
3 15.0
4 15.0
5 15.0
75.0

对于上述的 for-loop 也可用 list 的方式呈现

n = 5
triDuo = [i*(i+1)/2 for i in range(1,n+1)]
print(triDuo,sum(triDuo))

[1.0, 3.0, 6.0, 10.0, 15.0] 35.0

3.c 三角垛 与 for-loop

可对 triDuo 再用一层 for 来比较不同的 n 值

N = 7
triDuo = 0 
for n in range(2,N+1):
    triDuo = [i*(i+1)/2 for i in range(1,n+1)]
    print(triDuo,sum(triDuo))

[1.0, 3.0] 4.0
[1.0, 3.0, 6.0] 10.0
[1.0, 3.0, 6.0, 10.0] 20.0
[1.0, 3.0, 6.0, 10.0, 15.0] 35.0
[1.0, 3.0, 6.0, 10.0, 15.0, 21.0] 56.0
[1.0, 3.0, 6.0, 10.0, 15.0, 21.0, 28.0] 84.0

对上述的总和其实也可用公式 $\dfrac{n(n+1)(n+2)}{3\times 2\times 1}$ 来计算

• 请写个数列来显示 $\dfrac{n(n+1)(n+2)}{3\times 2\times 1}$

triDuoM = [int(i*(i+1)*(i+2)/(2*3))for i in range(1,n+1)]
print(triDuoM,sum(triDuoM))

[1, 4, 10, 20, 35, 56, 84] 210

上列的总和竟然也可写为 $\dfrac{n(n+1)(n+2)(n+3)}{1\times 2\times 3\times 4}$ !!!

用數學的 $\sum$ 表達有如下關係:

$$\begin{array}{rcl} \sum\limits_{k=1}^n k &=& \dfrac{n(n+1)}{2} $.2cm] \sum\limits_{k=1}^n\sum\limits_{i=1}^k i &=& \dfrac{n(n+1)(n+2)}{1\cdot 2\cdot 3} $.2cm] \sum\limits_{k=1}^n\sum\limits_{i=1}^k \sum\limits_{j=1}^i j &=& \dfrac{n(n+1)(n+2)(n=3)}{1\cdot 2\cdot 3 \cdot 4} \end{array}$$

3.d 一尺之棰,日取其半,万世不竭

请计算以下数列的总和：

$1 + \frac{1}{2} + \frac{1}{2^2} + \frac{1}{2^3} +\cdots +\frac{1}{2^{n-1}}$

n = 10
seqG2 = [float(1.0/(2**i)) for i in range(0,n)]
print(seqG2,sum(seqG2))

[1.0, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125, 0.00390625, 0.001953125] 1.998046875

N = 10
for n in range(1,N+1):
    seqG2 = [float(float(1.0/(2**i))) for i in range(0,n)]
    print("{} 项的和为 {}".format(n,sum(seqG2)))

1 项的和为 1.0
2 项的和为 1.5
3 项的和为 1.75
4 项的和为 1.875
5 项的和为 1.9375
6 项的和为 1.96875
7 项的和为 1.984375
8 项的和为 1.9921875
9 项的和为 1.99609375
10 项的和为 1.998046875

可利用 {:2d}、{:10.8f} 来调整输出格式

N = 15
for n in range(1,N+1):
    seqG2 = [1/2**i for i in range(0,n)]
    print("{:2d} 项的和为 {:10.8f}".format(n,sum(seqG2)))

 1 项的和为 1.00000000
 2 项的和为 1.50000000
 3 项的和为 1.75000000
 4 项的和为 1.87500000
 5 项的和为 1.93750000
 6 项的和为 1.96875000
 7 项的和为 1.98437500
 8 项的和为 1.99218750
 9 项的和为 1.99609375
10 项的和为 1.99804688
11 项的和为 1.99902344
12 项的和为 1.99951172
13 项的和为 1.99975586
14 项的和为 1.99987793
15 项的和为 1.99993896

【习题3a】立方和数列

请写个回圈比较以下两个数列的关系

$$\begin{align} \mbox{seq}1(n) &= 1+2+3+4+\cdots+n \\ \mbox{seq}3(n) &= 1^3+2^3+3^3+4^3+\cdots+n^3 \\ \end{align}$$

N = 10
print(" n | 连续整数和 连续立方和")
for i in range(2,N+1):
    seq1 = sum([k for k in range(i)])
    seq3 = sum([k**3 for k in range(i)])
    print("{:2d} | {:8d} {:8d}".format(i,seq1,seq3))

 n | 连续整数和 连续立方和
 2 |        1        1
 3 |        3        9
 4 |        6       36
 5 |       10      100
 6 |       15      225
 7 |       21      441
 8 |       28      784
 9 |       36     1296
10 |       45     2025

【习题3b】n角数

(1) 下图为三角数，其圆圈的数量可表示为公差为 1 的数列和。

$1 + 2 + 3 + \cdots + n$

(2) 下图为四角数，其圆圈的数量可表示为公差为 3 的数列和。

$1 + 3 + 5 + \cdots + 2n-1$

(3) 下图为五角数，其圆圈的数量可表示为公差为 2 的数列和。

$1 + 4 + 7 + \cdots + 3n-2$

先观察第 2 层比第 1 层多 4 个，共 1 + 4 = 5 个。

先观察第 3 层比第 2 层多 7 个，共 1 + 4 + 7 = 12 个。

先观察第 4 层比第 3 层多 10 个，共 1 + 4 + 7 + 10 = 22 个。

(4) 下图为六角数，请列出每边有 $3 \sim 9$ 个数的六角数。

$1 + 5 + 9 + \cdots + 4n-3$

先观察第 2 层比第 1 层多 5 个，共 1 + 5 = 6 个。

先观察第 3 层比第 2 层多 9 个，共 1 + 5 + 9 = 15 个。

先观察第 4 层比第 3 层多 13 个，共 1 + 5 + 9 + 13 = 28 个。

• 請列出一个表格显示 數量為 $3 \sim 10$ 的 三角數、四角數、五角數、六角數。

print(" i | 三角数 四角数 五角数 六角数")
print("---|------------------------")
for i in range(3,11):
    S3 = sum([k for k in range(1,i+1)])
    S4 = sum([2*k-1 for k in range(1,i+1)])
    S5 = sum([3*k-1 for k in range(1,i+1)])
    S6 = sum([4*k-1 for k in range(1,i+1)])
    print("{:2d} |{:6d}{:6d}{:6d}{:6d}".format(i,S3,S4,S5,S6))

 i | 三角数 四角数 五角数 六角数
---|------------------------
 3 |     6     9    15    21
 4 |    10    16    26    36
 5 |    15    25    40    55
 6 |    21    36    57    78
 7 |    28    49    77   105
 8 |    36    64   100   136
 9 |    45    81   126   171
10 |    55   100   155   210

【习题3c】四角垛

而 n 阶四角垛的总个数为 $1 +2^2+ 3^2 + 4^2 + \cdots + n^2$ ，也可用公式来求 $\frac{n(n+1)(2n+1)}{6}$。

• 请分别用数列的和与数学公式来验证上述公式。

【补充】上述图形也可表示

$$\begin{align} & 1 + 2^2 + 3^2 + 4^2 + 5^2 \\ = & 1 + (1+2) + (1+2+3) + (1+2+3+4) + (1+2+3+4+5)\\ & \quad + 1 + (1+2) + (1+2+3) + (1+2+3+4) \\ \end{align}$$

利用公式来代换也可得到

$\frac{n(n+1)(2n+1)}{6} = \frac{n(n+1)(n+2)}{6} + \frac{(n-1)n(n+1)}{6}$

N = 10
for i in range(1,N+1):
    seq = [k**2 for k in range(1,i)]
    print("{:2d} {:6.0f} {:6.0f}".format(i,sum(seq),i*(i+1)*(2*i+1)/(6)))

 1      0      1
 2      1      5
 3      5     14
 4     14     30
 5     30     55
 6     55     91
 7     91    140
 8    140    204
 9    204    285
10    285    385

【习题3d】等比数列

• 请计算以下数列的总和

• 输出时利用 {:10.8f} 来控制小数点以下有 8 位。 $$1 + \frac{1}{3} + \frac{1}{3^2} + \frac{1}{3^3} + \cdots +\frac{1}{3^{n-1}}$$

【数学】当 $r<1$ 时， $1+r+r^2+r^3+\cdots + r^n + \cdots$ 的极限为 $\frac{1}{1-r}$.

N = 15
r= 1.0/3
for n in range(1,N+1):
    seq = [r**k for k in range(0,n)]
    print('{:2d} 项的和为 {:10.8f} 与极限差 {:10.8f}'.format(n, sum(seq),1.0/(1-r)-sum(seq)))

 1 项的和为 1.00000000 与极限差 0.50000000
 2 项的和为 1.33333333 与极限差 0.16666667
 3 项的和为 1.44444444 与极限差 0.05555556
 4 项的和为 1.48148148 与极限差 0.01851852
 5 项的和为 1.49382716 与极限差 0.00617284
 6 项的和为 1.49794239 与极限差 0.00205761
 7 项的和为 1.49931413 与极限差 0.00068587
 8 项的和为 1.49977138 与极限差 0.00022862
 9 项的和为 1.49992379 与极限差 0.00007621
10 项的和为 1.49997460 与极限差 0.00002540
11 项的和为 1.49999153 与极限差 0.00000847
12 项的和为 1.49999718 与极限差 0.00000282
13 项的和为 1.49999906 与极限差 0.00000094
14 项的和为 1.49999969 与极限差 0.00000031
15 项的和为 1.49999990 与极限差 0.00000010

4 import math 与 pi 的探讨 ：$\frac{\pi}{2} =\frac{2}{1} \cdot \frac{2}{3} \cdot \frac{4}{3}\cdot \frac{4}{5} \cdot \frac{6}{5} \cdots$

4.a 调和数列

请计算下列算式到第 n 项的总和

$1+ \dfrac{1}{2} + \dfrac{1}{3} + \cdots + \dfrac{1}{n}$

N = 10
seqH = [float(1.0/k) for k in range(1,N+1)] 
print(seqH,sum(seqH))

[1.0, 0.5, 0.3333333333333333, 0.25, 0.2, 0.16666666666666666, 0.14285714285714285, 0.125, 0.1111111111111111, 0.1] 2.9289682539682538

请利用 for 上述数列在 n = 10, 100, 1000, 10000 项的总和变化

for n in [10,100,1000,10000]:
    seqH = [1/k for k in range(1,n+1)] 
    print(n,sum(seqH)) 

10 2.9289682539682538
100 5.187377517639621
1000 7.485470860550343
10000 9.787606036044348

将上述输出结果格式化

for n in [10,100,1000,10**4,10**5]:
    seqH = [1/k for k in range(1,n+1)] 
    print("{:6d} {:10.7f}".format(n,sum(seqH)))

    10  2.9289683
   100  5.1873775
  1000  7.4854709
 10000  9.7876060
100000 12.0901461

这个调和数列是发散数列。

但这数列的近似值为 $\log_e(n) + \frac{1}{2}$

请比较这两个数值的差。（要用 log 需要 import math）

from math import log
for n in [10,100,1000,10**4,10**5]:
    seqH = [1/k for k in range(1,n+1)] 
    print("{:6d} {:10.7f} {:10.7f} {:10.7f}".format(n, sum(seqH), log(n)+0.5 ,sum(seqH)-log(n)-0.5 ))

    10  2.9289683  2.8025851  0.1263832
   100  5.1873775  5.1051702  0.0822073
  1000  7.4854709  7.4077553  0.0777156
 10000  9.7876060  9.7103404  0.0772657
100000 12.0901461 12.0129255  0.0772207

4.b 调和平方数列

请计算比较下列算式到第 n 项的总和 $(n=10,100,1000,10000,100000)$

$1+ \dfrac{1}{2^2} + \dfrac{1}{3^2} + \cdots + \dfrac{1}{n^2}$

for n in [10,100,1000,10000,100000]:
    seqSH = [1/k**2 for k in range(1,n+1)] 
    print(n,sum(seqSH)) 

10 1.5497677311665408
100 1.6349839001848923
1000 1.6439345666815615
10000 1.6448340718480652
100000 1.6449240668982423

请比较此数列的总和与 $\dfrac{\pi^2}{6}$ 的差距

import math
for n in [10,100,1000,10000,100000]:
    seqSH = [1/k**2 for k in range(1,n+1)] 
    print("{:7d} {:10.7f} {:10.7f}".format(n,sum(seqSH),math.pi**2/6-sum(seqSH)))

     10  1.5497677  0.0951663
    100  1.6349839  0.0099502
   1000  1.6439346  0.0009995
  10000  1.6448341  0.0001000
 100000  1.6449241  0.0000100

4.c 交错数列

下列 $\pi$ 的近似数列公式为 格雷果理在 1671 年发现，而莱布尼兹在 1674 年发现。

$$\dfrac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9}+\cdots$$

请先输出 1,-1,1,-1,1,-1,1 的数列

n = 10
[(-1)**(k+1) for k in range(1,n+1)]

[1, -1, 1, -1, 1, -1, 1, -1, 1, -1]

请输出 $\dfrac{1}{1},\dfrac{1}{3},\dfrac{1}{5},\dfrac{1}{7},\cdots$ 的数列

n = 10
[1/(2*k-1) for k in range(1,n+1)]

[1.0,
 0.3333333333333333,
 0.2,
 0.14285714285714285,
 0.1111111111111111,
 0.09090909090909091,
 0.07692307692307693,
 0.06666666666666667,
 0.058823529411764705,
 0.05263157894736842]

请输出 $\dfrac{1}{1},\dfrac{-1}{3},\dfrac{1}{5},\dfrac{-1}{7},\cdots$ 的数列

n = 10
[(-1)**(k+1)/(2*k-1) for k in range(1,n+1)]

[1.0,
 -0.3333333333333333,
 0.2,
 -0.14285714285714285,
 0.1111111111111111,
 -0.09090909090909091,
 0.07692307692307693,
 -0.06666666666666667,
 0.058823529411764705,
 -0.05263157894736842]

import math
for n in [10,100,1000,10000,100000]:
    seqAH = [(-1)**(k+1)/(2*k-1) for k in range(1,n+1)] 
    print("{:6d} {:10.7f} {:10.7f}".format(n,sum(seqAH),math.pi/4-sum(seqAH)))

    10  0.7604599  0.0249383
   100  0.7828982  0.0024999
  1000  0.7851482  0.0002500
 10000  0.7853732  0.0000250
100000  0.7853957  0.0000025

4.d 用乘法来估算 pi

约翰沃利斯在 1655 年发现了沃利斯乘积，是欧洲第二个无穷项圆周率公式。

$$\frac{\pi}{2} =\frac{2}{1} \cdot \frac{2}{3} \cdot \frac{4}{3}\cdot \frac{4}{5} \cdot \frac{6}{5} \cdots$$

请验算计算以下数列，先试着列出 2,2,4,4,6,6 与 1,3,5,7,9,13

N = 10
evens = [2*((i+1)//2) for i in range(1,N+1)]
odds = [2*((i)//2)+1 for i in range(1,N+1)]
print(evens)
print(odds)

[2, 2, 4, 4, 6, 6, 8, 8, 10, 10]
[1, 3, 3, 5, 5, 7, 7, 9, 9, 11]

使用 for-loop 将上述数逐项连乘再相除

N = 20
wallis = [1]*(N+1)
evens = [2*int((i+1)//2) for i in range(1,N+1)]
odds = [2*int((i)//2)+1 for i in range(1,N+1)]
for i in range(1,N-1):
    wallis[i] = wallis[i-1]*evens[i-1]/odds[i-1]
    print(i, wallis[i],evens[i],odds[i])

1 2.0 2 3
2 1.3333333333333333 4 3
3 1.7777777777777777 4 5
4 1.422222222222222 6 5
5 1.7066666666666663 6 7
6 1.4628571428571426 8 7
7 1.6718367346938774 8 9
8 1.4860770975056687 10 9
9 1.6511967750062988 10 11
10 1.5010879772784533 12 11
11 1.6375505206674035 12 13
12 1.51158509600068 14 13
13 1.627860872616117 14 15
14 1.5193368144417092 16 15
15 1.6206259354044898 16 17
16 1.5252949980277553 18 17
17 1.6150182332058585 18 19
18 1.530017273563445 20 19

利用 pyplot 来绘制 pyplot

N = 100
wallis = [1]*(N)
evens = [2*int((i+1)//2) for i in range(1,N+1)]
odds = [2*int((i)//2)+1 for i in range(1,N+1)]
for i in range(1,N):
    wallis[i] = 1.0*wallis[i-1]*evens[i-1]/odds[i-1]

import matplotlib.pyplot as plt
plt.plot(wallis)

[<matplotlib.lines.Line2D at 0x7fbbfc279160>]

【习题4a】Telescope 数列

请计算以下数列到 100, 1000, 10000 项的总和

$$\mbox{ts}(n) = \frac{1}{1\times 2} + \frac{1}{2\times 3}+ \frac{1}{3\times 4}+\cdots+\frac{1}{n\times (n+1)}$$

请判断，此数列是否收收敛或发散？

for i in [100,1000,10000]:
    seq = [1/(k*(k+1)) for k in range(1,i+1)]
    print(i, sum(seq))

100 0.9900990099009898
1000 0.9990009990009997
10000 0.9999000099990007

【习题4b】交错数列与 ln

请写个函式来观察 $\ln(2)$ 与 $1-\frac{1}{2} + \frac{1}{3} - \frac{1}{4} +\cdots$ 的差距。

并用 matplotlib.pyplot 来绘制其收敛情况

• Python 中， $\ln(2)$ 要用 math.log(2) 来输出

import math
def seqln2(n):
  return float(sum([1/k for k in range(1,n+1)]))

for k in [10,100, 1000, 10000]:
  print("{0:5} {1:9.7f} {2:9.7f} {3:9.7f}".format(k,seqln2(k), math.log(2), seqln2(k)-math.log(2)))

   10 2.9289683 0.6931472 2.2358211
  100 5.1873775 0.6931472 4.4942303
 1000 7.4854709 0.6931472 6.7923237
10000 9.7876060 0.6931472 9.0944589

【习题4c】Viete 的 pi 的逼近式

以下為 1593 年，由 Francois Viete 發現可用以下無窮數列來逼近 pi 的值的。請利用迴圈來計算其求 $n$ 項的值。

$\frac{2}{\pi} = \frac{\sqrt{2}}{2}\frac{\sqrt{2+\sqrt{2}}}{2}\frac{\sqrt{2+\sqrt{2+\sqrt{2}}}}{2}\cdots$

import math

for n in [10,100,1000,10000]:
    result = 1
    seq = [0]
    for i in range(1,n+1):
      seq.append(math.sqrt(2+seq[i-1]))
      result *= seq[i]/2
    print("{:5d}{:13.9f}{:13.9f}".format(n,2/result,2/result-math.pi))

   10  3.141591422 -0.000001232
  100  3.141592654  0.000000000
 1000  3.141592654  0.000000000
10000  3.141592654  0.000000000

【习题4d】多重根式与黄金比例

请用回圈计算以下数列 $a_n$ ，并比较此数列与 $\phi = \frac{1+\sqrt{5}}{2}$ 的关系。

$$\begin{align} a_1&= \sqrt{1}\\ a_2&= \sqrt{1+\sqrt{1}}\\ a_3&= \sqrt{1+\sqrt{1+\sqrt{1}}}\\ a_4&= \sqrt{1+\sqrt{1+\sqrt{1+\sqrt{1}}}}\\ \end{align}$$

N = 10
fib = [1]
phi = (1+5**0.5)/2

n = 1

for i in range(1,N+1):
  fib.append((fib[-1]+1)**0.5)
  if (abs(fib[-1]-phi)>0.001):
    n = i+1
  
print("a_n 與 phi 在 n >= {} 時，差距小於 0.001".format(n))


import matplotlib.pyplot as plt
print(fib)
plt.plot(fib,'-o')

a_n 與 phi 在 n >= 6 時，差距小於 0.001
[1, 1.4142135623730951, 1.5537739740300374, 1.5980531824786175, 1.6118477541252516, 1.616121206508117, 1.6174427985273905, 1.617851290609675, 1.6179775309347393, 1.6180165422314876, 1.6180285974702324]

[<matplotlib.lines.Line2D at 0x7fbbf9df9d68>]