本篇博客介绍了如何使用NumPy进行科学计算，包括通过矩阵表示计算斐波那契数列，绘制利萨茹曲线以及使用傅里叶级数绘制方波、锯齿波和三角波。此外，还探讨了位操作函数和比较函数在整数处理中的应用。
第五章 矩阵和通用函数
5.13 斐波那契数列
斐波那契（Fibonacci）数列是基于递推关系生成的。直接用NumPy代码来解释递推关系是比较麻烦的，不过我们可以用矩阵的形式或者黄金分割公式来解释它。因此，我们将介绍matrix和rint函数。使用matrix函数创建矩阵，rint函数对浮点数取整，但结果仍为浮点数类型。
5.14 动手实践：计算斐波那契数列
斐波那契数列的递推关系可以用矩阵来表示。斐波那契数列的计算等价于矩阵的连乘。
- (1) 创建斐波那契矩阵：
import numpy as np
F = np.matrix([[1, 1], [1, 0]])
print( "F", F)
创建的斐波那契矩阵如下所示：
F [[1 1]
[1 0]]
- (2) 计算斐波那契数列中的第8个数，即矩阵的幂为8减去1。计算出的斐波那契数位于矩阵的对角线上：
print( "8th Fibonacci", (F ** 7)[0, 0])
输出的第8个斐波那契数为：
8th Fibonacci 21
- (3) 利用黄金分割公式或通常所说的比奈公式（Binet’ s Formula），加上取整函数，就可以直接计算斐波那契数。计算前8个斐波那契数：
n = np.arange(1, 9)
sqrt5 = np.sqrt(5)
phi = (1 + sqrt5)/2
fibonacci = np.rint((phi**n - (-1/phi)**n)/sqrt5)
print( "Fibonacci", fibonacci)
输出的斐波那契数为：
Fibonacci [ 1. 1. 2. 3. 5. 8. 13. 21.]
案例完整代码如下：
import numpy as np
F = np.matrix([[1, 1], [1, 0]])
print( "F", F)
print( "8th Fibonacci", (F ** 7)[0, 0])
n = np.arange(1, 9)
sqrt5 = np.sqrt(5)
phi = (1 + sqrt5)/2
fibonacci = np.rint((phi**n - (-1/phi)**n)/sqrt5)
print( "Fibonacci", fibonacci)
5.15 利萨茹曲线
在NumPy中，所有的标准三角函数如sin、 cos、 tan等均有对应的通用函数。利萨茹曲线（Lissajous curve）是一种很有趣的使用三角函数的方式。我至今仍记得在物理实验室的示波器上显示出利萨茹曲线时的情景。利萨茹曲线由以下参数方程定义：
x = A sin(at + n/2)
y = B sin(bt)
5.16 动手实践：绘制利萨茹曲线
利萨茹曲线的参数包括A、 B、 a和b。为简单起见，我们令A和B均为1。
- (1) 使用
linspace函数初始化变量t，即从-pi到pi上均匀分布的201个点：
a = 9
b = 8
t = np.linspace(-np.pi, np.pi, 201)
- (2) 使用
sin函数和NumPy常量pi计算变量x：
x = np.sin(a * t + np.pi/2)
- (3) 使用
sin函数计算变量y：
y = np.sin(b * t)
(4) 绘制的曲线如下所示。
案例完整代码如下：
import numpy as np
import matplotlib.pyplot as plt
a = 9
b = 8
t = np.linspace(-np.pi, np.pi, 201)
x = np.sin(a * t + np.pi/2)
y = np.sin(b * t)
plt.plot(x, y)
plt.show()
5.17 方波
方波也是一种可以在示波器上显示的波形。方波可以近似表示为多个正弦波的叠加。事实上，任意一个方波信号都可以用无穷傅里叶级数来表示。
傅里叶级数（Fourier series）是以正弦函数和余弦函数为基函数的无穷级数，以著名的法国数学家Jean-Baptiste Fourier命名。
方波可以表示为如下的傅里叶级数：
5.18 动手实践：绘制方波
与前面的教程中一样，我们仍将以相同的方式初始化t和k。我们需要累加很多项级数，且级数越多结果越精确，这里取k=99以保证足够的精度。绘制方波的步骤如下。
- (1) 我们从初始化
t和k开始，并将函数值初始化为0：
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(-np.pi, np.pi, 201)
k = np.arange(1, 99)
k = 2 * k - 1
f = np.zeros_like(t)
- (2) 接下来，直接使用
sin和sum函数进行计算：
for i in range(len(t)):
f[i] = np.sum(np.sin(k * t[i])/k)
f = (4 / np.pi) * f
- (3) 绘制波形：
plt.plot(t, f)
plt.show()
案例完整代码如下：
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(-np.pi, np.pi, 201)
k = np.arange(1, 99)
k = 2 * k - 1
f = np.zeros_like(t)
for i in range(len(t)):
f[i] = np.sum(np.sin(k * t[i])/k)
f = (4 / np.pi) * f
plt.plot(t, f)
plt.show()
5.19 锯齿波和三角波
在示波器上，锯齿波和三角波也是常见的波形。和方波类似，我们也可以将它们表示成无穷傅里叶级数。对锯齿波取绝对值即可得到三角波。锯齿波的无穷级数表达式如下：
5.20 动手实践：绘制锯齿波和三角波
与前面的教程中一样，我们仍将以相同的方式初始化t和k。同样，取k=99以保证足够的精度。绘制锯齿波和三角波的步骤如下。
- (1) 将函数值初始化为
0：
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(-np.pi, np.pi, 201)
k = np.arange(1, 99)
f = np.zeros_like(t)
- (2) 同样，直接使用
sin和sum函数进行计算：
for i in range(len(t)):
f[i] = np.sum(np.sin(2 * np.pi * k * t[i])/k)
f = (-2 / np.pi) * f
- (3) 同时绘制锯齿波和三角波并不难，因为三角波函数的取值恰好是锯齿波函数值的绝对值。使用如下代码绘制波形：
plt.plot(t, f, lw=1.0)
plt.plot(t, np.abs(f), lw=2.0)
plt.show()
在下图中，较粗的曲线为三角波。
案例完整代码如下：
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(-np.pi, np.pi, 201)
k = np.arange(1, 99)
f = np.zeros_like(t)
for i in range(len(t)):
f[i] = np.sum(np.sin(2 * np.pi * k * t[i])/k)
f = (-2 / np.pi) * f
plt.plot(t, f, lw=1.0)
plt.plot(t, np.abs(f), lw=2.0)
plt.show()
5.21 位操作函数和比较函数
位操作函数可以在整数或整数数组的位上进行操作，它们都是通用函数。 ^、 &、 |、 <<、 >>等位操作符在NumPy中也有对应的部分， <、>、 ==等比较运算符也是如此。有了这些操作符，你可以在代码中玩一些高级技巧以提升代码的性能。不过，它们会使代码变得难以理解，因此需谨慎使用。
5.22 动手实践：玩转二进制位
我们将学习三个小技巧——检查两个整数的符号是否一致，检查一个数是否为2的幂数，以及计算一个数被2的幂数整除后的余数。我们会分别展示两种方法，即使用位操作符和使用相应的NumPy函数。
- (1) 第一个小技巧需要用
XOR或者^操作符。XOR操作符又被称为不等运算符，因此当两个操作数的符号不一致时，XOR操作的结果为负数。在NumPy中，^操作符对应于bitwise_xor函数，<操作符对应于less函数。
import numpy as np
x = np.arange(-9, 9)
y = -x
print( "Sign different?", (x ^ y) < 0)
print( "Sign different?", np.less(np.bitwise_xor(x, y), 0))
不出所料，除了都等于0的情况，所有整数对的符号均相异。
- (2) 在二进制数中，
2的幂数表示为一个1后面跟一串0的形式，例如10、100、1000等。而比2的幂数小1的数表示为一串二进制的1，例如11、111、1111（即十进制里的3、7、15）等。如果我们在2的幂数以及比它小1的数之间执行位与操作AND，那么应该得到0。在NumPy中，&操作符对应于bitwise_and函数，==操作符对应于equal函数。
print( "Power of 2?\n", x, "\n", (x & (x - 1)) == 0)
print( "Power of 2?\n", x, "\n", np.equal(np.bitwise_and(x, (x - 1)), 0))
结果如下：
Power of 2?
[-9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8]
[False False False False False False False False False True True True
False True False False False True]
Power of 2?
[-9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8]
[False False False False False False False False False True True True
False True False False False True]
- (3) 计算余数的技巧实际上只在模为
2的幂数（如4、8、16等）时有效。二进制的位左移一位，则数值翻倍。在前一个小技巧中我们看到，将2的幂数减去1可以得到一串1组成的二进制数，如11、111、1111等。这为我们提供了掩码（mask），与这样的掩码做位与操作AND即可得到以2的幂数作为模的余数。在NumPy中，<<操作符对应于left_shift函数。
print( "Modulus 4\n", x, "\n", x & ((1 << 2) - 1))
print( "Modulus 4\n", x, "\n", np.bitwise_and(x, np.left_shift(1, 2) - 1))
结果如下：
Modulus 4
[-9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8]
[3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0]
Modulus 4
[-9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8]
[3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0]
案例完整代码如下：
import numpy as np
x = np.arange(-9, 9)
y = -x
print( "Sign different?", (x ^ y) < 0)
print( "Sign different?", np.less(np.bitwise_xor(x, y), 0))
print( "Power of 2?\n", x, "\n", (x & (x - 1)) == 0)
print( "Power of 2?\n", x, "\n", np.equal(np.bitwise_and(x, (x - 1)), 0))
print( "Modulus 4\n", x, "\n", x & ((1 << 2) - 1))
print( "Modulus 4\n", x, "\n", np.bitwise_and(x, np.left_shift(1, 2) - 1))
扫一扫
4365
到【灌水乐园】发言
