python解zuobiaoxi方程_在学习线性代数时所探索的Python运用

这是一篇伪科学文章
文章仅仅是学习中的经验总结，可能存在大量技术错误 本文中所有代码均在jupyerlab+python3.8.5+numpy1.19.0+sympy1.6.2版本下测试
若有错误，那就这么地吧，反正没人看到。

---

前言 在学习lay的《线性代数及其应用》时，遇到了一道练习题是运用MATLAB等数据科学软件来计算的。由于最近正在学习Python，Python作为近几年数据处理和数据科学的爆款，很多人都证明了其地位。 于是我尝试在Python中解决这道习题，在Python中将矩阵化简为行简化矩阵(RREF)。
1. 教材中题目的简述
在原题中，提到了插值多项式即类似于p(t) = a₀ +a₁t + a₂t²的多项式计算，根据已有的数据t和p列出线性方程组计算出几个参数a然后列出其方程式并使用其进行预测的方法，在计量统计中就是多元线性回归模型。
原题干中的数据如下：
由于题干使用的是五次多项式，所以这里我也使用五次多项式，其它次数的多项式将在文末探讨。 p(t) = a₀ +a₁t + a₂t²+ a₃t³ + a₄t⁴ + a₅t⁵ 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) # 转换成arrayres_arr = np.array(resists) # 转换成array

• 由于需要的是一个矩阵，所以通过将数组spe_arr用切片的方式再转化为(6, 1)的数组。接下来为了转化为(6, 6)的数组，通过numpy的广播特性，加上一个(6, 6)的全为零的数组。

[In]

spe_arr2d = spe_arr[:, None] + np.zeros((6, 6)) # 这里的None等同于newaxisspe_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) = a₀ +a₁t + a₂t²+ a₃t³ + a₄t⁴ + a₅t⁵进行处理。输出的是一个以数组表示的系数矩阵。

[In]

for i in range(6):    for j in range(6):        spe_arr2d[i, j]**=j # 依次对每个元素进行处理spe_arr2d

• 注意：由于微信公众号的代码输入框会导致**=j出现换行，以上代码应为下图

[Out]

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值拼接于其后成为增广矩阵。

[In]

# 由于要在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()计算行简化矩阵。

[In]

from sympy import Matrix # 从sympy库中导入Matrixmat = 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]]

• 接下来进行简化计算：

[In]

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矩阵，水平有限，暂时只能找到这种方法)。

[In]

# 这里的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)

• 接下来定义函数:

[In]

# 将增广矩阵中最右列提取出来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计算出的一致。     p(t) = 1.7125t – 1.1948t2 + 0.6615t3 – 0.0701t4 原题所列出的线性方程组以及最后得出的方程式如上。 4. Tips 如在使用中numpy出现显示的为科学计数值，可通过以下代码取消科学计数。

np.set_printoptions(suppress=True)

     关于其它次数的多项式：

• 根据原线性方程组，不难发现总共有5个t，加上常数项总共有6个项，所以此时建立5次的多项式可以算出唯一值。
• 当项数小于5时，会发现计算出来的是一个超定方程组，方程数大于未知数个数，此时并不能得出精确解。
• 当项数大于5时，会发现是一个欠定方程组，就是方程数小于未知数个数，此时只能是无解或者无穷解。

---

文章仅为个人学习的记录总结，若有错误，感谢指出。