devbox
木道寻08
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

Python数学建模之线性代数模型_sp.simplify

作者：木道寻08 | 2024-08-06 17:46:54

赞

踩

sp.simplify
  • 本文只对已有数学模型进行Python代码的实现和解释，不分析模型的建立过程，模型的建立可查阅司守奎老师的《Python数学建模算法与应用》一书。

1.特征值与特征向量

1.差分方程

        解决此类问题需要使用sympy库函数，对于sympy库的详细介绍见这篇文章：

        http://t.csdnimg.cn/mOjYL

例：求斐波那契数列的通项公式。

解法一：运用特征值和特征向量求通项

  1. sp.var('k', positive=True, integer=True)    # 定义一个正整数变量k
  2. a = sp.Matrix([[0, 1], [1, 1]])     # 创建系数矩阵a，即A

 

  1. val = a.eigenvals()                 # 求A的特征值
  2. vec = a.eigenvects()                # 求A的特征向量

  1. P, D = a.diagonalize()              # 把A相似对角化，并将特征值矩阵和相似变换矩阵分别赋值给变量D和P
  2. ak = P @ (D**k) @ (P.inv())         # 求ak，即A^{k}。P.inv()表示矩阵P的逆矩阵

F = ak @ sp.Matrix([1, 1])          # 求F，即a_{k}

  1. s = sp.simplify(F[0])               # 对乘积F中的第一个元素，即F_{k}进行化简
  2. print(s)                            # 输出通项公式，即F_{k}

完整代码： 

  1. import sympy as sp
  2. sp.var('k', positive=True, integer=True)    # 定义一个正整数变量k
  3. a = sp.Matrix([[0, 1], [1, 1]])     # 创建系数矩阵a，即A
  4. val = a.eigenvals()                 # 求A的特征值
  5. vec = a.eigenvects()                # 求A的特征向量
  6. P, D = a.diagonalize()              # 把A相似对角化，并将特征值矩阵和相似变换矩阵分别赋值给变量D和P
  7. ak = P @ (D**k) @ (P.inv())         # 求ak，即A^{k}。P.inv()表示矩阵P的逆矩阵
  8. F = ak @ sp.Matrix([1, 1])          # 求F，即a_{k}
  9. s = sp.simplify(F[0])               # 对乘积F中的第一个元素，即F_{k}进行化简
  10. print(s)                            # 输出通项公式，即F_{k}
  11. """
  12. 将k逐个替换为0到19的整数，并计算s的数值，即输出斐波那契数列的前二十项。
  13. .subs(k, i)表示将符号表达式s中的符号变量k替换为变量i，得到一个新的符号表达式。
  14. .n()将这个新的符号表达式转换为数值，即进行数值计算。
  15. """
  16. sm = []
  17. for i in range(20):
  18.     sm.append(int(s.subs(k, i).n()))
  19. print(sm)

解法二：差分方程的特征根解法

  1. sp.var('t, c1, c2')             # 定义了三个符号变量t、c1和c2
  2. t0 = sp.solve(t**2-t-1)         # 求解特征方程t^{2}-t-1

  1. # 定义了两个方程eq1和eq2，分别表示c1+c2-1=0和c1*t0[0]+c2*t0[1]-1=0
  2. eq1 = c1+c2-1
  3. eq2 = c1*t0[0]+c2*t0[1]-1

  1. # 使用sp.solve()函数来解方程组，返回的是一个字典，包含方程组的解
  2. s = sp.solve([eq1, eq2])
  3. print('c1=', s[c1])
  4. print('c2=', s[c2])

 完整代码：

  1. import sympy as sp
  2. sp.var('t, c1, c2')             # 定义了三个符号变量t、c1和c2
  3. t0 = sp.solve(t**2-t-1)         # 求解特征方程t^{2}-t-1
  4. # 定义了两个方程eq1和eq2，分别表示c1+c2-1=0和c1*t0[0]+c2*t0[1]-1=0
  5. eq1 = c1+c2-1
  6. eq2 = c1*t0[0]+c2*t0[1]-1
  7. # 使用sp.solve()函数来解方程组，返回的是一个字典，包含方程组的解
  8. s = sp.solve([eq1, eq2])
  9. print('c1=', s[c1])
  10. print('c2=', s[c2])

        这种方法适用于不知道系数矩阵，只知道它的特征方程，否则可以直接用eigenvals()函数求解特征根。

解法三：直接利用Python软件求解

较为简单，不做解释。

  1. import sympy as sp
  2. sp.var('k')                                     # 定义了一个符号变量k
  3. y = sp.Function('y')                            # 定义了一个函数符号y，表示递归序列的通项公式
  4. f = y(k+2)-y(k+1)-y(k)                          # 定义了递归序列的递推关系式
  5. # 求解递归序列。sp.rsolve()函数接受三个参数：递推关系式、待求解的通项公式和初始条件，返回的结果为通项公式
  6. s = sp.rsolve(f, y(k), {y(0):1, y(1):1})
  7. print(s)

2.Leslie种群模型

  1. X0 = np.array([500, 1000, 500])
  2. L = np.array([[0, 4, 3], [0.5, 0, 0], [0, 0.25, 0]])

  1. X1 = L @ X0
  2. X2 = L @ X1
  3. X3 = L @ X2

  1. '''
  2. sp.var('lamda')                  # 定义符号变量
  3. p = Ls.charpoly(lamda)           # 计算矩阵Ls的特征多项式，并将结果赋值给变量p
  4. w1 = sp.roots(p)                 # 计算特征值
  5. '''
  6. w2 = Ls.eigenvals()              # 直接计算特征值
  7. v = Ls.eigenvects()              # 直接计算特征向量
  8. P, D = Ls.diagonalize()          # 相似对角化，将结果分别赋值给对角矩阵D和相似变换矩阵P

  1. Pinv = P.inv()                   # 求逆阵
  2. Pinv = sp.simplify(Pinv)         # 简化矩阵Pinv
  3. cc = Pinv @ X0
  4. print('c=', cc[0])

 完整代码：

  1. import numpy as np
  2. import sympy as sp
  3. X0 = np.array([500, 1000, 500])
  4. L = np.array([[0, 4, 3], [0.5, 0, 0], [0, 0.25, 0]])
  5. X1 = L @ X0
  6. X2 = L @ X1
  7. X3 = L @ X2
  8. # Rational()函数将有理数表示为分数的形式，以保留精确的数值
  9. Ls = sp.Matrix([[0, 4, 3],
  10.     [sp.Rational(1, 2), 0, 0],
  11.     [0, sp.Rational(1, 4), 0]])
  12. '''
  13. sp.var('lamda')                  # 定义符号变量
  14. p = Ls.charpoly(lamda)           # 计算矩阵Ls的特征多项式，并将结果赋值给变量p
  15. w1 = sp.roots(p)                 # 计算特征值
  16. '''
  17. w2 = Ls.eigenvals()              # 直接计算特征值
  18. v = Ls.eigenvects()              # 直接计算特征向量
  19. P, D = Ls.diagonalize()          # 相似对角化，将结果分别赋值给对角矩阵D和相似变换矩阵P
  20. Pinv = P.inv()                   # 求逆阵
  21. Pinv = sp.simplify(Pinv)         # 简化矩阵Pinv
  22. cc = Pinv @ X0
  23. print('c=', cc[0])

3.PageRank算法

1.基础的PageRank算法

  1. # 计算矩阵W
  2. L = [(1, 2), (2, 3), (2, 4), (3, 4), (3, 5),
  3.     (3, 6), (4, 1), (5, 6), (6, 1)]
  4. w = np.zeros((6, 6))                        # 邻接矩阵初始化
  5. for i in range(len(L)):
  6.     w[L[i][0]-1, L[i][1]-1] = 1
  7. # 计算矩阵P
  8. r = np.sum(w, axis=1, keepdims=True)
  9. P = w/r                                     # 这里利用矩阵广播

  1. """
  2. 用eigs函数来计算矩阵P的特征值和对应的特征向量
  3. P.T表示矩阵P的转置
  4. 参数1表示要计算的特征值和特征向量的数量,默认特征值按从大到小排列
  5. eigs函数的返回值是一个元组(val, vec)，val是一个包含特征值的一维数组，vec是一个包含特征向量的二维数组
  6. """
  7. val, vec = eigs(P.T, 1)
  8. V = vec.real                                # 将特征向量vec的实部提取出来，赋值给变量V
  9. V = V.flatten()                             # 将V展开成一维数组
  10. V = V/V.sum()
  11. print('V=', np.round(V, 4))

完整代码：

  1. import numpy as np
  2. from scipy.sparse.linalg import eigs
  3. import pylab as plt
  4. # 计算矩阵W
  5. L = [(1, 2), (2, 3), (2, 4), (3, 4), (3, 5),
  6.     (3, 6), (4, 1), (5, 6), (6, 1)]
  7. w = np.zeros((6, 6))                        # 邻接矩阵初始化
  8. for i in range(len(L)):
  9.     w[L[i][0]-1, L[i][1]-1] = 1
  10. # 计算矩阵P
  11. r = np.sum(w, axis=1, keepdims=True)
  12. P = w/r                                     # 这里利用矩阵广播
  13. """
  14. 用eigs函数来计算矩阵P的特征值和对应的特征向量
  15. P.T表示矩阵P的转置
  16. 参数1表示要计算的特征值和特征向量的数量,默认特征值按从大到小排列
  17. eigs函数的返回值是一个元组(val, vec)，val是一个包含特征值的一维数组，vec是一个包含特征向量的二维数组
  18. """
  19. val, vec = eigs(P.T, 1)
  20. V = vec.real                                # 将特征向量vec的实部提取出来，赋值给变量V
  21. V = V.flatten()                             # 将V展开成一维数组
  22. V = V/V.sum()
  23. print('V=', np.round(V, 4))
  24. # 用柱状图显示
  25. plt.bar(range(1, len(w)+1), V, width=0.6, color='b')
  26. plt.show()

 2.随机冲浪模型的PageRank值

 代码：

  1. import numpy as np
  2. from scipy.sparse.linalg import eigs
  3. import pylab as plt
  4. L = [(1, 2), (2, 3), (2, 4), (3, 4), (3, 5),
  5.     (3, 6), (4, 1), (5, 6), (6, 1)]
  6. w = np.zeros((6, 6))
  7. for i in range(len(L)):
  8.     w[L[i][0]-1, L[i][1]-1] = 1
  9. r = np.sum(w, axis=1, keepdims=True)
  10. P = (1-0.85)/w.shape[0]+0.85*w/r            # 这里利用矩阵广播
  11. val, vec = eigs(P.T, 1)
  12. V = vec.real
  13. V = V.flatten()                             # 展开成 (n, ) 形式的数组
  14. V = V/V.sum()
  15. print('V=', np.round(V, 4))
  16. plt.bar(range(1, len(w)+1), V, width=0.6, color='b')
  17. plt.show()

不做详细解释，使用时根据情况替换参数及可。 

2.矩阵的奇异值分解及应用

1.矩阵的奇异值分解

 求奇异值分解的代码如下，使用是将a替换为待求矩阵即可。

  1. import numpy as np 
  2. from numpy.linalg import svd
  3. a = np.array([[1, 0, 1], [0, 1, 1], [0, 0, 0]])
  4. u, s, vt = svd(a)       # a = u @ np.diag(s) @ vt
  5. print(u); print(s); print(vt)
'
运行

2.奇异值分解的应用

1.推荐系统的评分

txt文件的内容：

  1. 5	2	1	4	0	0	2	4	0	0	0
  2. 0	0	0	0	0	0	0	0	0	3	0
  3. 1	0	5	2	0	0	3	0	3	0	1
  4. 0	5	0	0	4	0	1	0	0	0	0
  5. 0	0	0	0	0	4	0	0	0	4	0
  6. 0	0	1	0	0	0	1	0	0	5	0
  7. 5	0	2	4	2	1	0	3	0	1	0
  8. 0	4	0	0	5	4	0	0	0	0	5
  9. 0	0	0	0	0	0	4	0	4	5	0
  10. 0	0	0	4	0	0	1	5	0	0	0
  11. 0	0	0	0	4	5	0	0	0	0	3
  12. 4	2	1	4	0	0	2	4	0	0	0
  13. 0	1	4	1	2	1	5	0	5	0	0
  14. 0	0	0	0	0	4	0	0	0	4	0
  15. 2	5	0	0	4	0	0	0	0	0	0
  16. 5	0	0	0	0	0	0	4	2	0	0
  17. 0	2	4	0	4	3	4	0	0	0	0
  18. 0	3	5	1	0	0	4	1	0	0	0

 解法一：非压缩数据的模型 

1）计算相似系数

2）评分估计

3）菜品推荐结果

  1. import numpy as np
  2. import pandas as pd
  3. # np.loadtxt()函数会从文本文件中加载数据，并将其存储为一个NumPy数组
  4. a = np.loadtxt('data3_6_1.txt')
  5. # corrcoef(a.T)根据数组a计算相关系数，得到相关系数矩阵
  6. # 再对结果进行归一化处理
  7. b = 0.5*np.corrcoef(a.T)+0.5
  8. # 使用Pandas将数组b转换为DataFrame，并将其保存为名为data3_6_2.xlsx的Excel文件，不包含索引
  9. c = pd.DataFrame(b)
  10. c.to_excel('data3_6_2.xlsx', index=False)
  11. print('请输入人员编号1-18')
  12. user = int(input())
  13. # 获取数组a的列数，即变量的个数，并将其保存在变量n中
  14. n = a.shape[1]
  15. # np.where返回值的格式为(array([ 0,  1,  2,  3,  4,  5,  6,  7,  8, 10], dtype=int64),)，取0号元素，即为所有未评分编号
  16. no = np.where(a[user-1,:]==0)[0]
  17. # set函数用于创建一个无序且不重复的集合对象，此处用集合操作得到已评分编号的集合
  18. yb = set(range(n))-set(no)
  19. # 转化为列表
  20. yb = list(yb)
  21. # 从数组a中获取第user-1行中已评分编号对应的分数，并将其保存在数组ys中
  22. ys = a[user-1, yb]
  23. """
  24. 创建一个长度为未评分编号的个数的零数组，并将其保存在数组sc中
  25. 对于每个未评分编号，计算其与已评分编号的相似度，并将计算结果保存在数组sc中
  26. 根据已评分的分数值估计出未评分的项目可能的分数
  27. """
  28. sc = np.zeros(len(no))                      
  29. for i in range(len(no)):
  30.     sim = b[no[i], yb]
  31.     sc[i] = ys @ sim/sum(sim)
  32. print('未评分项的编号为：', no+1)
  33. print('未评分项的分数为：', np.round(sc, 4))

解法二：基于奇异值分解压缩数据的模型

1）稀疏矩阵的降维处理

2）衡量菜品之间的相似性

3）评分估计

4）菜品推荐结果

  1. import numpy as np
  2. import pandas as pd
  3. a = np.loadtxt('data3_6_1.txt')
  4. # 使用numpy的linalg.svd()函数对矩阵a进行奇异值分解，分别得到左奇异向量矩阵u、奇异值向量sigma和右奇异向量矩阵vt
  5. u, sigma, vt = np.linalg.svd(a)
  6. print(sigma)
  7. # cumsum()函数的返回值是一个数组，第i个元素是sigma的0到i号元素的平方和
  8. cs = np.cumsum(sigma**2)
  9. # 将奇异值累积平方和除以总和，得到一个包含每个奇异值贡献率的数组
  10. rate = cs/cs[-1]
  11. """
  12. 找到第一个满足累积贡献率大于等于0.9的索引，并将其加1，得到奇异值的个数ind
  13. 这个值表示保留多少个奇异值，使得它们的累积贡献率达到了90%
  14. 这样做的目的是为了尽可能保留数据中的主要信息，同时实现数据降维的目标
  15. """
  16. ind = np.where(rate>=0.9)[0][0]+1
  17. # 构建一个对角矩阵，其对角线元素为前ind个奇异值，然后通过矩阵乘法得到降维后的数据b
  18. b = np.diag(sigma[:ind]) @ u.T[:ind, :] @ a
  19. # 使用np.linalg.norm()函数计算矩阵b的每列的范数，并保持维度为1
  20. c = np.linalg.norm(b, axis=0, keepdims=True)    
  21. # 计算相似度矩阵d，并归一化
  22. d = 0.5*b.T @ b/(c.T @ c)+0.5   
  23. # 其余步骤与解法一相同
  24. dd = pd.DataFrame(d)
  25. dd.to_excel('data3_6_3.xlsx', index=False)
  26. print('请输入人员编号1-18')
  27. user = int(input())
  28. n = a.shape[1]                                  
  29. no = np.where(a[user-1,:]==0)[0]                
  30. yb = set(range(n))-set(no)                      
  31. yb = list(yb)
  32. ys = a[user-1, yb]                              
  33. sc = np.zeros(len(no))                        
  34. for i in range(len(no)):
  35.     sim = d[no[i], yb]
  36.     sc[i] = ys @ sim/sum(sim)
  37. print('未评分项的编号为：', no+1)
  38. print('未评分项的分数为：', np.round(sc, 4))

2.利用SVD进行图像压缩

在数学建模中几乎用不到，了解即可。

  1. import numpy as np
  2. from numpy import linalg as LA 
  3. from PIL import Image
  4. import pylab as plt                 # 加载 Matplotlib 的 pylab 接口
  5. plt.rc('font', size=13) 
  6. plt.rc('font',family='SimHei')
  7. a = Image.open('Lena.bmp')          # 返回一个 PIL 图像对象
  8. if a.mode!='L':
  9.     a = a.convert('L')              # 转换为灰度图像
  10. b = np.array(a).astype(float)       # 把图像对象转换为数组
  11. [p, d, q] =  LA.svd(b)
  12. m, n = b.shape
  13. R = LA.matrix_rank(b)               # 图像矩阵的秩
  14. plt.figure(0)
  15. plt.plot(np.arange(1, len(d)+1), d, 'k.')
  16. plt.ylabel('奇异值'); plt.xlabel('序号')
  17. plt.title('图像矩阵的奇异值')
  18. CR = []
  19. for K in range(1, int(R/4), 10):
  20.     plt.figure(K)
  21.     plt.subplot(121)
  22.     plt.title('原图')
  23.     plt.imshow(b, cmap='gray')
  24.     I = p[:, :K+1] @ (np.diag(d[:K+1])) @ (q[:K+1, :])
  25.     plt.subplot(122)
  26.     plt.title('图像矩阵的秩='+str(K))
  27.     plt.imshow(I, cmap='gray')
  28.     src = m*n
  29.     compress = K*(m+n+1)
  30.     ratio = (1-compress/src)*100    # 计算压缩比率
  31.     CR.append(ratio)
  32.     print('Rank=%d:K=%d个:ratio=%5.2f'%(R, K ,ratio))
  33. plt.figure();
  34. plt.plot(range(1, int(R/4), 10), CR, 'bo-')
  35. plt.title('奇异值个数与压缩比率的关系')
  36. plt.xlabel('奇异值个数')
  37. plt.ylabel('压缩比率')
  38. plt.show()

声明：本文内容由网友自发贡献，转载请注明出处：【wpsshop博客】
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号