数学建模学习（四）

Python的numpy库使用

标准安装的Python中用列表（list）保存的一组值，可以用来当作数组使用，但是由于列表的元素可以是任意对象，因此列表中所保存的是对象的指针。这样为了保存一个简单的[1, 2, 3]，需要有3个指针和3个整数对象。对于数值运算来说，这种结构显然比较浪费内存和CPU的计算时间。

此外，Python还提供了一个array模块。array对象和列表不同，它直接保存数值，和C语言的一维数组比较类似。但是由于它不支持多维，也没有各种运算函数，因此也不适合做数值运算。下面将以具体的代码来学习numpy库。

数组生成示例一

import numpy as np
a1 = np.array([1, 2, 3, 4])   #生成整型数组
a2 = a1.astype(float)
a3 = np.array([1, 2, 3, 4], dtype=float)   #浮点数
print(a1.dtype); print(a2.dtype); print(a3.dtype)
b = np.array([[1, 2, 3], [4, 5, 6]])
c = np.arange(1,5)        #生成数组[1, 2, 3, 4]
d = np.linspace(1, 4, 4)  #生成数组[1, 2, 3, 4]
e = np.logspace(1, 3, 3, base=2)  #生成数组[2, 4, 8]

 数组生成示例二

import numpy as np
a = np.ones(4, dtype=int)     #输出[1, 1, 1, 1]
b = np.ones((4,), dtype=int)  #同a
c= np.ones((4,1))             #输出4行1列的数组
d = np.zeros(4)               #输出[0, 0, 0, 0]
e = np.empty(3)               #生成3个元素的空数组行向量
f = np.eye(3)                 #生成3阶单位阵
g = np.eye(3, k=1)  #生成第k对角线的元素为1，其他元素为0的3阶方阵
h = np.zeros_like(a)          #生成与a同维数的全0数组

数组元素索引示例

import numpy as np
a = np.arange(16).reshape(4,4)  #生成4行4列的数组
b = a[1][2]   #输出6
c = a[1, 2]   #同b
d = a[1:2, 2:3]  #输出[[6]]
x = np.array([0, 1, 2, 1])
print(a[x==1])  #输出a的第2、4行元素

矩阵合并示例

import numpy as np
a = np.arange(16).reshape(4,4)  #生成4行4列的数组
b = np.floor(5*np.random.random((2, 4)))
c = np.ceil(6*np.random.random((4, 2)))
d = np.vstack([a, b])  #上下合并矩阵
e = np.hstack([a, c])  #左右合并矩阵

矩阵分割示例

import numpy as np
a = np.arange(16).reshape(4,4)  #生成4行4列的数组
b = np.vsplit(a, 2)             #行分割
print('行分割：\n', b[0], '\n', b[1])
c = np.hsplit(a, 4)             #列分割
print('列分割：\n', c[0], '\n', c[1], '\n', c[2], '\n', c[3])

数组元素的索引示例

import numpy as np
a = np.arange(16).reshape(4,4)  #生成4行4列的数组
b = a[1][2]   #输出6
c = a[1, 2]   #同b
d = a[1:2, 2:3]  #输出[[6]]
x = np.array([0, 1, 2, 1])
print(a[x==1])  #输出a的第2、4行元素

矩阵元素求和示例

import numpy as np
a = np.array([[0, 3, 4], [1, 6, 4]])
b = a.sum()   #使用方法，求矩阵所有元素的和
c1 = sum(a)   #使用内置函数，求矩阵逐列元素的和
c2 = np.sum(a, axis=0) #使用函数，求矩阵逐列元素的和
c3 = np.sum(a, axis=0, keepdims=True)   #逐列求和
print(c2.shape, c3.shape)  #c2是(3,)数组，c3是(1,3)数组

逐个元素运算示例

import numpy as np
a = np.array([[0, 3, 4], [1, 6, 4]])
b = np.array([[1, 2, 3], [2, 1, 4]])
c = a / b   #两个矩阵对应元素相除
d = np.array([2, 3, 2])
e = a * d   #d先广播成与a同维数的矩阵，再逐个元素相乘
f = np.array([[3],[2]])
g = a * f   #f先广播成与a同维数的矩阵，再逐个元素相乘
h = a ** (1/2)  #a矩阵逐个元素的1/2次幂。

矩阵运算与线性代数

1.范数计算

import numpy as np
a = np.array([[0, 3, 4], [1, 6, 4]])
b = np.linalg.norm(a, axis=1)  #求行向量2范数
c = np.linalg.norm(a, axis=0)  #求列向量2范数
d = np.linalg.norm(a)   #求矩阵2范数
print('行向量2范数为：', np.round(b, 4))
print('列向量2范数为：', np.round(c, 4))
print('矩阵2范数为：', round(d, 4))

 2.求解线性方程组的唯一解

 

import numpy as np
a = np.array([[3, 1], [1, 2]])
b = np.array([9, 8])
x1 = np.linalg.inv(a) @ b  #第一种解法
#上面语句中@表示矩阵乘法
x2 = np.linalg.solve(a, b) #第二种解法
print(x1); print(x2)
求得x=2，y=3。

 3.求超定线性方程组的最小二乘解

import numpy as np
a = np.array([[3, 1], [1, 2], [1, 1]])
b = np.array([9, 8, 6])
x = np.linalg.pinv(a) @ b  
print(np.round(x, 4))
    求得的最小二乘解为x=2，y=3.1667。

4.求特征值和特征向量

import numpy as np
a = np.eye(4)
b = np.rot90(a)
c, d = np.linalg.eig(b)
print('特征值为：', c)
print('特征向量为：\n', d)