本专栏内容为：数学建模原理 记录学习数学建模
💓博主csdn个人主页：小小unicorn
⏩专栏分类：数学建模
🚚代码仓库：小小unicorn的代码仓库🚚
🌹🌹🌹关注我带你学习编程知识
目录
- numpy库
- 导入
- 数组创建
- 示例1：
- 示例2：
- 数组元素索引：
- 示例：
- 矩阵运算与线性代数
- 求解线性方程组的唯一解
- 例1：求解线性方程组：
- 求超定方程组的最小二乘解
- 例2：求解线性方程组：
- 第一种：
- 第二种：
- 两种方法的比较：
- 什么是超定方程？
- 求矩阵的特征值特征向量：
- Pandas库
- 生成二维数组：
- 读写文件
- 数据预处理：
- SciPy库
- 求解非线性方程组
- 积分：
- 最小二乘解
- 求矩阵最大模特征值与特征向量
- SymPy库
- 例1：求符号解
- 例2：求符号解
- 例3：求矩阵的特征值与特征向量（符号解）
- Matplotlib库
- 子图：
- 三维绘图
- 例1：三维曲线
numpy库
导入
import numpy as
数组创建
示例1：
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]
print(b)
print(c)
print(d)
print(e)
示例2：
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行元素
矩阵运算与线性代数
线性代数只要使用numpy.linalg模块，其常用函数如下：
求解线性方程组的唯一解
例1：求解线性方程组：
代码：
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)
结果：
求超定方程组的最小二乘解
例2：求解线性方程组：
代码：
第一种：
import numpy as np
a = np.array([[3, 1], [1, 2], [1, 1]])
b = np.array([9, 8, 6])
#a*b
x = np.linalg.pinv(a) @ b
print("求得的最小二乘解为：")
print(np.round(x, 4))
第二种：
使用 lstsq 求解最小二乘问题:
- a: 系数矩阵。
- b: 右侧的结果向量。
- rcond: 用来裁剪奇异值的一个参数（可以先设置为 None，默认行为会使用合理的值）
lstsq 函数返回四个值：
- x: 最小二乘解，即满足 a @ x ≈ b 的解。
- residuals: 残差的平方和。如果方程是过定的（更多的方程比未知数多），这个值才有意义。
- rank: 矩阵 a 的秩。
- s: a 的奇异值。
import numpy as np
a = np.array([[3, 1], [1, 2], [1, 1]])
b = np.array([9, 8, 6])
x, residuals, rank, s = np.linalg.lstsq(a, b, rcond=None)
print("求得的最小二乘解为：")
#round对数组中的每个元素进行四舍五入。4就是保留四位
print(np.round(x, 4))
两种方法的比较：
-
pinv方法：先通过np.linalg.pinv(a)求解a的伪逆矩阵，再用矩阵乘法@乘以向量b来得到最小二乘解。适合直接计算伪逆矩阵的场合。 -
lstsq方法：直接求解最小二乘解，它比伪逆法更高效，因为它不需要显式计算伪逆矩阵，并且提供了残差、矩阵秩等额外信息。 -
在数值精度和计算效率上，
lstsq通常表现得更好，尤其是当矩阵a较大时
结果：
什么是超定方程？
超定线性方程组（overdetermined system of linear equations）指的是方程的数量多于未知数的线性方程组。这意味着方程的数量 mmm 大于未知数的数量 nnn（即 m>nm > nm>n）。
特点：
- 方程数量多于未知数：这类方程组往往没有精确的解，因为过多的方程对未知数施加了更多的约束条件，这些条件可能互相矛盾，无法同时满足。
- 无法精确求解：由于不可能满足所有方程，通常需要使用最小二乘法等数值方法来找到使方程“近似成立”的解。这种解不是精确解，而是尽可能让方程残差（误差）最小的解。
例子：
求矩阵的特征值特征向量：
例1：
代码：
import numpy as np
#它创建一个 4x4 的单位矩阵
#单位矩阵的定义是：对角线上（从左上到右下）所有元素为 1，其余元素为 0。
a = np.eye(4)
# a = [[1. 0. 0. 0.]
# [0. 1. 0. 0.]
# [0. 0. 1. 0.]
# [0. 0. 0. 1.]]
#将矩阵逆时针旋转 90 度
b = np.rot90(a)
# b = [[0. 0. 0. 1.]
# [0. 0. 1. 0.]
# [0. 1. 0. 0.]
# [1. 0. 0. 0.]]
#如果想顺时针旋转90度，加上一个参数k=-1；
# b=np.rot90(a,-1)
# b = [[0. 0. 0. 1.]
# [0. 0. 1. 0.]
# [0. 1. 0. 0.]
# [1. 0. 0. 0.]]
#顺时针旋转与逆时针旋转的关系：
#顺时针旋转 90 度相当于逆时针旋转 270 度。因此，你也可以通过 np.rot90(a, k=3) 达到相同的效果。
c, d = np.linalg.eig(b)
print('特征值为：', c)
print('特征向量为：\n', d)
结果：
Pandas库
生成二维数组：
例1：生成服从标准正态分布的24x4随机数矩阵，并保存为DataFrame数据结构
import pandas as pd
import numpy as np
# 创建数据
dates = pd.date_range(start='20191101', end='20191124', freq='D')
a1 = pd.DataFrame(np.random.randn(24, 4), index=dates, columns=list('ABCD'))
a2 = pd.DataFrame(np.random.rand(24, 4))
# 数据写入文件实例
a1.to_excel('data2_4_1.xlsx') # 保存到 Excel
a2.to_csv('data2_4_2.csv') # 保存到 CSV
#不包含行索引时
a1.to_excel('data2_4_1.xlsx', index=False) # 保存到 Excel，不包含行索引
a2.to_csv('data2_4_2.csv', index=False) # 保存到 CSV，不包含行索引
读写文件
例2：写入文件
import pandas as pd
import numpy as np
# 创建数据
dates = pd.date_range(start='20191101', end='20191124', freq='D')
a1 = pd.DataFrame(np.random.randn(24, 4), index=dates, columns=list('ABCD'))
a2 = pd.DataFrame(np.random.rand(24, 4))
# 数据写入文件实例
a1.to_excel('data2_4_1.xlsx') # 保存到 Excel
a2.to_csv('data2_4_2.csv') # 保存到 CSV
#不包含行索引时
a1.to_excel('data2_4_1.xlsx', index=False) # 保存到 Excel，不包含行索引
a2.to_csv('data2_4_2.csv', index=False) # 保存到 CSV，不包含行索引
# 使用 with 语句创建 ExcelWriter 对象
with pd.ExcelWriter('data2_4_3.xlsx') as writer:
a1.to_excel(writer, sheet_name="Sheet1") # 把a1写入Excel文件的Sheet1
a2.to_excel(writer, sheet_name="Sheet2") # 把a2写入另一个表单Sheet2
# 文件会自动保存和关闭，无需显式调用 save() 方法
结果：
如果写入数据时，不包含行索引：
import pandas as pd
import numpy as np
# 创建数据
dates = pd.date_range(start='20191101', end='20191124', freq='D')
a1 = pd.DataFrame(np.random.randn(24, 4), index=dates, columns=list('ABCD'))
a2 = pd.DataFrame(np.random.rand(24, 4))
#不包含行索引时
a1.to_excel('data2_4_4.xlsx', index=False) # 保存到 Excel，不包含行索引 header=False 不包含列索引
a2.to_csv('data2_4_5.csv', index=False) # 保存到 CSV，不包含行索引
# 使用 with 语句创建 ExcelWriter 对象
with pd.ExcelWriter('data2_4_6.xlsx') as writer:
a1.to_excel(writer, sheet_name="Sheet1", index=False) # 把a1写入Excel文件的Sheet1
a2.to_excel(writer, sheet_name="Sheet2",index=False) # 把a2写入另一个表单Sheet2
数据预处理：
例1：拆分，合并，分组计算
import pandas as pd
import numpy as np
d=pd.DataFrame(np.random.randint(1,6,(10,4)), columns=list("ABCD"))
d1=d[:4] #获取前4行数据
d2=d[4:] #获取第5行以后的数据
dd=pd.concat([d1,d2]) #数据行合并
s1=d.groupby('A').mean() #数据分组求均值
s2=d.groupby('A').apply(sum) #数据分组求和
print(d,'\n')
print(d1,'\n')
print(d2,'\n')
print(dd,'\n')
print(s1,'\n')
print(s2,'\n')
结果：
例2：DataFrame数据操作示例
import pandas as pd
import numpy as np
a = pd.DataFrame(np.random.randint(1,6,(5,3)),
index=['a', 'b', 'c', 'd', 'e'],
columns=['one', 'two', 'three'])
a.loc['a', 'one'] = np.nan #修改第1行第1列的数据
b = a.iloc[1:3, 0:2].values #提取第2、3行，第1、2列数据
a['four'] = 'bar' #增加第4列数据
a2 = a.reindex(['a', 'b', 'c', 'd', 'e', 'f'])
a3 = a2.dropna() #删除有不确定值的行
print(a,'\n')
print(b,'\n')
print(a2,'\n')
print(a3,'\n')
SciPy库
模块功能表：
求解非线性方程组
求方程：
代码实现：
from scipy.optimize import fsolve, root
fx = lambda x: x**980-5.01*x**979+7.398*x**978\
-3.388*x**977-x**3+5.01*x**2-7.398*x+3.388
x1 = fsolve(fx, 1.5, maxfev=4000) #函数调用4000次
x2 = root(fx, 1.5)
print(x1,'\n','-------------');print(x2)
注意：
例2：求下列方程组的一组数值解
from scipy.optimize import fsolve, root
fx = lambda x: [x[0]**2+x[1]**2-1, x[0]-x[1]]
#初始值为[1,1]
s1 = fsolve(fx, [1, 1])
s2 = root(fx, [1, 1])
print(s1,'\n','--------------'); print(s2)
积分：
例1：
from scipy.integrate import quad
def fun46(x, a, b):
return a*x**2+b*x
I1 = quad(fun46, 0, 1, args=(2, 1))
I2 = quad(fun46, 0, 1, args=(2, 10))
print(I1); print(I2)
结果：
最小二乘解
例1：
代码：
from scipy.optimize import least_squares
import numpy as np
a=np.loadtxt('Z:\Code\python\建模培训2024\code\data2_47.txt')
x0=a[0]; y0=a[1]; d=a[2]
fx=lambda x: np.sqrt((x0-x[0])**2+(y0-x[1])**2)-d
s=least_squares(fx, np.random.rand(2))
print(s, '\n', '------------', '\n', s.x)
求矩阵最大模特征值与特征向量
代码：
from scipy.sparse.linalg import eigs
import numpy as np
a = np.array([[1, 2, 3], [2, 1, 3], [3, 3, 6]], dtype=float) #必须加float,否则出错
b, c = np.linalg.eig(a)
d, e = eigs(a, 1)
print('最大模特征值为：', d)
print('对应的特征向量为：\n', e)
SymPy库
定义符号变量或符号函数命令如下：
例1：求符号解
import sympy as sp
a, b, c, x=sp.symbols('a,b,c,x')
x0=sp.solve(a*x**2+b*x+c, x)
print(x0)
例2：求符号解
代码实现：
import sympy as sp
# sp.var('x1,x2')
# s=sp.solve([x1**2+x2**2-1,x1-x2],[x1,x2])
# print(s)
x1,x2=sp.symbols('x1,x2')
s=sp.solve([x1**2+x2**2-1,x1-x2],[x1,x2])
print(s)
例3：求矩阵的特征值与特征向量（符号解）
import numpy as np
import sympy as sp
a = np.identity(4) #单位矩阵的另一种写法
b = np.rot90(a)
c = sp.Matrix(b)
print('特征值为：', c.eigenvals())
print('特征向量为：\n', c.eigenvects())
Matplotlib库
plot绘图常见的样式和颜色类型：
例1：
代码：
import pandas as pd
import pylab as plt
plt.rc('font',family='SimHei') #用来正常显示中文标签
plt.rc('font',size=16) #设置显示字体大小
a=pd.read_excel("Z:\Code\python\建模培训2024\code\data2_52.xlsx", header=None)
b=a.values #提取其中的数据
x=b[0]; y=b[1:]
#第三个参数。虚实线+形状+颜色
plt.plot(x,y[0],'-*b',label='钻石')
plt.plot(x,y[1],'--dr',label='铂金')
plt.xlabel('月份'); plt.ylabel('每月销量')
plt.legend(loc='upper left'); plt.grid(); plt.show()
例2：画出销售额的柱状图;
import pandas as pd
import pylab as plt
plt.rc('font',family='SimHei') #用来正常显示中文标签
plt.rc('font',size=16) #设置显示字体大小
a=pd.read_excel("Z:\Code\python\建模培训2024\code\data2_52.xlsx",header=None)
b=a.T
b.plot(kind='bar');
plt.legend(['钻石', '铂金'])
plt.xticks(range(6), b[0], rotation=0)
plt.ylabel('数量'); plt.show()
子图：
代码：
import pylab as plt
import numpy as np
plt.rc('font',family='SimHei') #用来正常显示中文标签
plt.rc('axes',unicode_minus=False)
#生成随机数，进行归一化处理
y1=np.random.randint(2, 5, 6)
#使其和为1
y1=y1/sum(y1)
#创建第一个子图，水平条形图
# 在一个 2x2 的子图网格中创建第一个子图，位于第一个位置。
plt.subplot(2, 2, 1)
str=['苹果', '葡萄', '桃子', '梨', '香蕉', '菠萝']
plt.barh(str,y1) #水平条形图
#创建第二个子图，饼图
plt.subplot(222)
plt.pie(y1, labels=str) #饼图
#创建第三个子图，折线图
plt.subplot(212)
#生成一个包含 100 个点的数组，这些点均匀分布在 0.01 到 10 之间。
x2=np.linspace(0.01, 10, 100)
y2=np.sin(10*x2)/x2
plt.plot(x2,y2); plt.xlabel('$x$')
plt.ylabel('$\\mathrm{sin}(10x)/x$'); plt.show()
结果：
三维绘图
例1：三维曲线
代码：
import pylab as plt
import numpy as np
ax=plt.axes(projection='3d') #设置三维图形模式
z=np.linspace(-50, 50, 1000)
x=z**2*np.sin(z); y=z**2*np.cos(z)
ax.plot(x, y, z, 'c')
plt.show()
例2：画出三维表面图：
import pylab as plt
import numpy as np
x=np.linspace(-4,4,100)
x,y=np.meshgrid(x,x)
z=50*np.sin(x+y)
ax=plt.axes(projection='3d')
ax.plot_surface(x, y, z, color='b')
plt.show()
例3：
代码：
import pylab as plt
import numpy as np
ax=plt.axes(projection='3d')
X = np.arange(-6, 6, 0.25)
Y = np.arange(-6, 6, 0.25)
X, Y = np.meshgrid(X, Y)
Z = np.sin(np.sqrt(X**2 + Y**2))
surf = ax.plot_surface(X, Y, Z, cmap='coolwarm')
#在图形旁边添加一个颜色条，用于表示不同颜色对应的 Z 值（即表面的高度）。
plt.colorbar(surf); plt.show()
