这是一篇伪科学文章
文章仅仅是学习中的经验总结,可能存在大量技术错误 本文中所有代码均在jupyerlab+python3.8.5+numpy1.19.0+sympy1.6.2版本下测试
若有错误,那就这么地吧,反正没人看到。
前言 在学习lay的《线性代数及其应用》时,遇到了一道练习题是运用MATLAB等数据科学软件来计算的。由于最近正在学习Python,Python作为近几年数据处理和数据科学的爆款,很多人都证明了其地位。 于是我尝试在Python中解决这道习题,在Python中将矩阵化简为行简化矩阵(RREF)。
1. 教材中题目的简述
在原题中,提到了插值多项式即类似于p(t) = a0 +a1t + a2t2的多项式计算,根据已有的数据t和p列出线性方程组计算出几个参数a然后列出其方程式并使用其进行预测的方法,在计量统计中就是多元线性回归模型。
原题干中的数据如下:![08d16092187f259c933740edbaa1f93e.png](https://i-blog.csdnimg.cn/blog_migrate/17f4a3235cd46b37761664c4a83d832c.png)
由于题干使用的是五次多项式,所以这里我也使用五次多项式,其它次数的多项式将在文末探讨。 p(t) = a0 +a1t + a2t2+ a3t3 + a4t4 + a5t5 2. 在Python中进行操作
[Out]
3. 最终结果
成功算出当风速为750ft/s时的空气阻力约为64.8,与教材中使用MATLAB计算出的一致。
p(t) = 1.7125t – 1.1948t2 + 0.6615t3 – 0.0701t4
原题所列出的线性方程组以及最后得出的方程式如上。
4. Tips
如在使用中numpy出现显示的为科学计数值,可通过以下代码取消科学计数。
文章仅为个人学习的记录总结,若有错误,感谢指出。
文章仅仅是学习中的经验总结,可能存在大量技术错误 本文中所有代码均在jupyerlab+python3.8.5+numpy1.19.0+sympy1.6.2版本下测试
若有错误,那就这么地吧,反正没人看到。
前言 在学习lay的《线性代数及其应用》时,遇到了一道练习题是运用MATLAB等数据科学软件来计算的。由于最近正在学习Python,Python作为近几年数据处理和数据科学的爆款,很多人都证明了其地位。 于是我尝试在Python中解决这道习题,在Python中将矩阵化简为行简化矩阵(RREF)。
1. 教材中题目的简述
在原题中,提到了插值多项式即类似于p(t) = a0 +a1t + a2t2的多项式计算,根据已有的数据t和p列出线性方程组计算出几个参数a然后列出其方程式并使用其进行预测的方法,在计量统计中就是多元线性回归模型。
原题干中的数据如下:
![08d16092187f259c933740edbaa1f93e.png](https://i-blog.csdnimg.cn/blog_migrate/17f4a3235cd46b37761664c4a83d832c.png)
由于题干使用的是五次多项式,所以这里我也使用五次多项式,其它次数的多项式将在文末探讨。 p(t) = a0 +a1t + a2t2+ a3t3 + a4t4 + a5t5 2. 在Python中进行操作
首先输入数据并把数据转换成Numpy的array数组。
import numpy as np # 导入numpy库
speeds = [0, 2, 4, 6, 8, 10] # 风速的大小
resists = [0, 2.9, 14.8, 39.6, 74.3, 119] # 阻力的大小
spe_arr = np.array(speeds) # 转换成array
res_arr = np.array(resists) # 转换成array
- 由于需要的是一个矩阵,所以通过将数组spe_arr用切片的方式再转化为(6, 1)的数组。接下来为了转化为(6, 6)的数组,通过numpy的广播特性,加上一个(6, 6)的全为零的数组。
spe_arr2d = spe_arr[:, None] + np.zeros((6, 6)) # 这里的None等同于newaxis
spe_arr2d
[Out]
# 输出的是(6, 6)的数组
array([[ 0., 0., 0., 0., 0., 0.],
[ 2., 2., 2., 2., 2., 2.],
[ 4., 4., 4., 4., 4., 4.],
[ 6., 6., 6., 6., 6., 6.],
[ 8., 8., 8., 8., 8., 8.],
[10., 10., 10., 10., 10., 10.]])
- 接下来将数组中的每个元素都依次根据p(t) = a0 +a1t + a2t2+ a3t3 + a4t4 + a5t5进行处理。输出的是一个以数组表示的系数矩阵。
for i in range(6):
for j in range(6):
spe_arr2d[i, j]**=j # 依次对每个元素进行处理
spe_arr2d
- 注意:由于微信公众号的代码输入框会导致**=j出现换行,以上代码应为下图
![cc29e7d1f23d00764444e4cb901f62e7.png](https://i-blog.csdnimg.cn/blog_migrate/c57b4ededf7cdaa88308db6705567078.png)
array([
[ 1., 0., 0., 0., 0., 0.],
[ 1., 2., 4., 8., 16., 32.],
[ 1., 4., 16., 64., 256., 1024.],
[ 1., 6., 36., 216., 1296., 7776.],
[ 1., 8., 64., 512., 4096., 32768.],
[ 1., 10., 100., 1000., 10000., 100000.]])
然后就是通过numpy.concatenate()将方程的p值拼接于其后成为增广矩阵。
# 由于要在spe_arr2d数组的后面依次拼接,所以要把轴axis设置为1
# axis=1意为从外面的0轴开始数的1轴 最外层[]为0轴
mat_np = np.concatenate((spe_arr2d, res_arr[:, None]), axis=1)
mat_np
[Out]
# 这是一个6 * 7的数组(增广矩阵)
array([
[ 1. , 0. , 0. , 0. , 0. , 0. , 0.],
[ 1. , 2. , 4. , 8. , 16. , 32. , 2.9],
[ 1. , 4. , 16. , 64. , 256. , 1024. , 14.8],
[ 1. , 6. , 36. , 216. ,1296. , 7776. , 39.6],
[ 1. , 8. , 64. , 512. ,4096. , 32768. , 74.3],
[ 1. , 10. , 100. ,1000. ,10000. , 100000. ,119. ]])
- 由于numpy中似乎没有提供对矩阵进行求行简化矩阵的方法,通过网上搜索 #https://blog.csdn.net/ikoiiii/article/details/92075982 知道Python中另一个库sympy中提供了直接求行简化矩阵的方法。
- 但是其中的Matrix和numpy中的matrix似乎不兼容,所以讲刚才得到的numpy数组通过tolist()方法转化成Python列表再通过sympy.Matrix转化成矩阵再用Matrix.rref()计算行简化矩阵。
from sympy import Matrix # 从sympy库中导入Matrix
mat = Matrix(mat_np.tolist())
# Matrix中有一些参数例如类型可以设置,详见官方文档
mat
[Out]
Matrix([
[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 2.9],
[1.0, 4.0, 16.0, 64.0, 256.0, 1024.0, 14.8],
[1.0, 6.0, 36.0, 216.0, 1296.0, 7776.0, 39.6],
[1.0, 8.0, 64.0, 512.0, 4096.0, 32768.0, 74.3],
[1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0, 119.0]]
- 接下来进行简化计算:
mat_rref = mat.rref()
mat_rref
[Out]
(Matrix([[1, 0, 0, 0, 0, 0, 0], # 第一个参数a0为0
[0, 1, 0, 0, 0, 0, 1.7125], # 第二个参数a1为1.7125
[0, 0, 1, 0, 0, 0, -1.19479166666667], # 以下同上 [0, 0, 0, 1, 0, 0, 0.661458333333333],
[0, 0, 0, 0, 1, 0, -0.0700520833333333],
[0, 0, 0, 0, 0, 1, 0.00260416666666667]]),
(0, 1, 2, 3, 4, 5)) # 这一行值表示的是矩阵主元列所在的列数
# 这里可以列出方程为p(t) = 1.7125t – 1.1948t2 + 0.6615t3 – 0.0701t4
- 到此已经算出矩阵的行简化矩阵了,接下来是列出p(t) = 1.7125t – 1.1948t2 + 0.6615t3 – 0.0701t4。之后就是求出p(7.5),这里需要再次将Matrix转回numpy矩阵
- (因sympy中切片的方式与numpy不同,为更好理解,转回numpy矩阵,水平有限,暂时只能找到这种方法)。
# 这里的mat_rref必须进行切片,因为原矩阵还有一行表示主元列。
# 同时要用dtype设定类型,否则会是object类型
mat_rref_num = np.mat(mat_rref[0].tolist(), dtype='float64')
mat_rref_num
[Out]
matrix([[1, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1.71250000000000],
[0, 0, 1, 0, 0, 0, -1.19479166666667],
[0, 0, 0, 1, 0, 0, 0.661458333333333],
[0, 0, 0, 0, 1, 0, -0.0700520833333333],
[0, 0, 0, 0, 0, 1, 0.00260416666666667]], dtype=float64)
- 接下来定义函数:
# 将增广矩阵中最右列提取出来
a = []
for i in range(6):
a.append(mat_rref_num[i, 6])
# 定义p(t)函数
def p(t):
p = [a[i]*(t**i) for i in range(6)]
return np.sum(p)
p(7.5) # 计算出当t等于7.5时的p值
[Out]
64.83838274147274
3. 最终结果
成功算出当风速为750ft/s时的空气阻力约为64.8,与教材中使用MATLAB计算出的一致。
![9983e10b57a437389e181e8c844824f9.png](https://i-blog.csdnimg.cn/blog_migrate/b160985b4b60621dbfcfa092f6d9e2bb.png)
np.set_printoptions(suppress=True)
关于其它次数的多项式:
- 根据原线性方程组,不难发现总共有5个t,加上常数项总共有6个项,所以此时建立5次的多项式可以算出唯一值。
- 当项数小于5时,会发现计算出来的是一个超定方程组,方程数大于未知数个数,此时并不能得出精确解。
- 当项数大于5时,会发现是一个欠定方程组,就是方程数小于未知数个数,此时只能是无解或者无穷解。
文章仅为个人学习的记录总结,若有错误,感谢指出。