有理数python_Python数学实验与建模（3）

「爱情、让人受尽委屈。」 2023-01-05 12:45 27阅读 0赞

## 第3章 Python在高等数学和线性代数中的应用 ##

1.  SymPy工具介绍
2.  SciPy工具库简介
3.  用SymPy做符号函数画图
4.  高等数学问题的符号解
5.  高等数学问题的数值解
6.  线性代数问题的符号解和数值解

### 1.SymPy工具库介绍 ###

1）sympy工具库简介

SymPy是Python版的开源计算机代数系统实现，通俗地讲，SymPy是用于符号运算的工具库，现在这个工具库包括几十个模块：

>>> help('sympy')
    Help on package sympy:
    
    NAME
        sympy
    
    DESCRIPTION
        SymPy is a Python library for symbolic mathematics. It aims to become a
        full-featured computer algebra system (CAS) while keeping the code as simple
        as possible in order to be comprehensible and easily extensible.  SymPy is
        written entirely in Python. It depends on mpmath, and other external librari
    es
        may be optionally for things like plotting support.
    
        See the webpage for more information and documentation:
    
            https://sympy.org
    
    PACKAGE CONTENTS
        abc  # 符号变量模块
        algebras (package)  # 代数
        assumptions (package) 
        benchmarks (package)
        calculus (package)
        categories (package)
        codegen (package)
        combinatorics (package)
        concrete (package)
        conftest
        core (package)  # 基本的加、乘、指数运算等
        crypto (package)  
        deprecated (package)  
        diffgeom (package)  
        discrete (package)  # 离散数学
        external (package)
        functions (package)
        galgebra  # 几何代数
        geometry (package)  # 几何实体
        holonomic (package)
        integrals (package)  # 符号积分
        interactive (package)  # 交互会话
        liealgebras (package)
        logic (package)  # 布尔代数和定理证明
        matrices (package)  # 线性代数和矩阵
        multipledispatch (package)
        ntheory (package)  # 数论函数
        parsing (package)
        physics (package)  # 物理学
        plotting (package)  # 用Pyglet进行二维和三维的画图
        polys (package)  # 多项式代数和因式分解
        printing (package)  # 漂亮的打印和代码生成
        release
        sandbox (package)
        series (package)  # 级数
        sets (package)
        simplify (package)  # 化简符号表达式
        solvers (package)  # 方程求解
        stats (package)  # 统计学
        strategies (package)
        tensor (package)
        testing (package)
        this
        unify (package)
        utilities (package)
        vector (package)

详细资料：

SymPy Modules Reference Python科学计算利器——SymPy库

几个模块的功能：

*  微积分模块（sympy.intergrals）：微积分模块支持大量的基础与高级微积分运算功能。
 *  离散数学模块（sympy.discrete）
 *  方程求解模块（sympy.solvers）
 *  矩阵模块（sympy.matrices）
 *  物理学模块（sympy.physics）
 *  统计学模块（sympy.stats）

2）符号运算基础知识

符号创建、类型转换及subs()方法代入值示例：

from sympy import *
    
    x, y, z = symbols('x  y  z')
    m0, m1, m2, m3 = symbols('m0:4')  # 创建多个符号变量
    x = sin(1)
    print("x=", x);
    print("x=", x.evalf())
    print("x=", x.n(16))  # 显示小数点后16位数字
    print("pi的两种显示格式:{},{}".format(pi, pi.evalf(3)))  # 这里不能使用n()函数
    expr1 = y * sin(y ** 2)  # 创建第一个符号表达式
    expr2 = y ** 2 + sin(y) * cos(y) + sin(z)  # 创建第二个符号表达式
    print("expr1=", expr1)
    print("y=5时，expr1=", expr1.subs(y, 5))  # 代入一个符号变量的值
    print("y=2,z=3时，expr2=", expr2.subs({
         y: 2, z: 3}))  # 代入y=2,z=3
    print("y=2,z=3时，expr2=", expr2.subs({
         y: 2, z: 3}).n())  # 以浮点数显示计算结果

运行结果：

![e1222b3ecb4f775be6ce1543da1e1b4a.png][]

符号函数together()及apart()使用示例：

from sympy import *
    
    x1, x2, x3, x4 = symbols('m1:5')
    x = symbols('x')
    print(x1 / x2 + x3 / x4)
    print(together(x1 / x2 + x3 / x4))
    print((2 * x ** 2 + 3 * x + 4) / (x + 1))
    print(simplify((2 * x ** 2 + 3 * x + 4) / (x + 1)))  # 化简没有效果
    print(apart((2 * x ** 2 + 3 * x + 4) / (x + 1)))

运行结果：

![10e1497f596e5ff83879559ab73ee616.png][]

### 2.SciPy工具库简介 ###

SciPy是对NumPy的功能扩展，它提供了许多高级数学函数，例如，微积分、微分方程、优化方法、数值分析、高级统计函数、方程求解等。SciPy是在NumPy数组框架基础上实现的，它对NumPy数组和基本的数组运算进行扩展，满足科学家和工程师解决问题时需要的大部分数学计算功能。

详细资料：SciPy v1.6.0 Reference Guide

![d5796704503d169f1fd0b2d40615481d.png][]

SciPy可以使用MATLAB的“.mat”文件格式读取和写入数据，函数为loadmat和savemat。

### 3.用SymPy做符号函数画图 ###

1）二维曲线画图

plot的基本使用格式：

plot(表达式, 变量取值范围, 属性 = 属性值)

多重绘制的使用格式为：

plot(表达式1, 表达式2, 变量取值范围, 属性 = 属性值)

或者：

plot((表达式1, 变量取值范围1), (表达式2, 变量取值范围2))

例：在同一图形界面上画出

![Image 1][]

from sympy.plotting import plot
    from sympy.abc import x, pi  # 引进符号变量x及常量pi
    from sympy.functions import sin, cos
    
    plot((2 * sin(x), (x, -6, 6)), (cos(x + pi / 4), (x, -5, 5)))

![ebb8e97d9bff00de0a6f4b58a03c5a74.png][]

2）三维曲面画图

例：画出三维曲面

![Image 1][] 的图形

from pylab import rc  # pylab为matplotlib的接口
    from sympy.plotting import plot3d
    from sympy.abc import x, y  # 引进符号变量x,y
    from sympy.functions import sin, sqrt
    
    rc('font', size=16);
    rc('text', usetex=True)
    plot3d(sin(sqrt(x ** 2 + y ** 2)), (x, -10, 10), (y, -10, 10), xlabel='$x$', ylabel='$y$')

![9a9b94a6e0db7bda0ea77bd805c89453.png][]

3）隐函数画图

例：绘制隐函数

![Image 1][] 的图形。

from pylab import rc
    from sympy import plot_implicit as pt, Eq
    from sympy.abc import x, y  # 引进符号变量x,y
    
    rc('font', size=16)
    rc('text', usetex=False)
    pt(Eq((x - 1) ** 2 + (y - 2) ** 3, 4), (x, -6, 6), (y, -2, 4), xlabel='$x$', ylabel='$y$')

![848d54a08dc066089ded910f50375a32.png][]

或使用以下代码：

from sympy import plot_implicit as pt
    from sympy.abc import x, y  # 引进符号变量x,y
    
    ezplot = lambda expr: pt(expr)
    ezplot((x - 1) ** 2 + (y - 2) ** 3 - 4)

### 4.高等数学问题的符号解 ###

SymPy包括许多功能，从基本的符号算术到多项式、微积分、求解方程、离散数学和统计等。它主要处理三种类型的数据：整型数据、实数和有理数。有理数包括两个部分：分子和分母，可以用Ration类定义有理数。

1）求极限

例：验证

![Image 1][] .

from sympy import *
    
    x = symbols('x')
    print(limit(sin(x) / x, x, 0))
    print(limit(pow(1 + 1 / x, x), x, oo))  # 这里是两个小”o”，表示正无穷

运行结果：

![53c322c8cd33d3f9993d5038f7392fb4.png][]

2）求导数

例：已知

![Image 1][] ，求

![Image 1][] .

from sympy import *
    
    x, y = symbols('x y')  # 定义两个符号变量
    z = sin(x) + x ** 2 * exp(y)  # 构造符号表达式
    print("关于x的二阶偏导数为：", diff(z, x, 2))
    print("关于y的一阶偏导数为：", diff(z, y))

运行结果：

![cd1544536b01a31571a26991c9cdedf6.png][]

3）级数求和

例：验证

![Image 1][] .

from sympy import *
    
    k, n = symbols('k  n')
    print(summation(k ** 2, (k, 1, n)))
    print(factor(summation(k ** 2, (k, 1, n))))  # 把计算结果因式分解
    print(summation(1 / k ** 2, (k, 1, oo)))  # 这里是两个小o表示正无穷

运行结果：

![6d11f6358ed85a994696bc49b500b0ea.png][]

4）泰勒展开

例：写出

![Image 1][] 在0点处的3,5,7阶泰勒展开式，并在同一图形界面上画出

![Image 1][] 及它的上述各阶泰勒展开式在区间\[0, 2\]上的图形.

from pylab import rc
    from sympy import *
    
    rc('font', size=16)
    rc('text', usetex=False)
    x = symbols('x')
    y = sin(x)
    for k in range(3, 8, 2):
        print(y.series(x, 0, k))  # 等价于print(series(y,x,0,k))
    plot(y, series(y, x, 0, 3).removeO(), series(y, x, 0, 5).removeO(),
         series(y, x, 0, 7).removeO(), (x, 0, 2), xlabel='$x$', ylabel='$y$')

![026c12a4f4c7566e2aa1486307540493.png][]

5）不定积分和定积分

例：验证

![Image 1][] .

from sympy import integrate, symbols, sin, pi, oo
    
    x = symbols('x')
    print(integrate(sin(2 * x), (x, 0, pi)))
    print(integrate(sin(x) / x, (x, 0, oo)))

运行结果：

![fb327b977a3b3aac18548fab3cd6534d.png][]

6）求解代数方程（方程组）的符号解

例：求如下代数方程：（1）

![Image 1][] ; （2）

![Image 1][] .

from sympy import *
    
    x, y = symbols('x  y')
    print(solve(x ** 3 - 1, x))
    print(solve((x - 2) ** 2 * (x - 1) ** 3, x))
    print(roots((x - 2) ** 2 * (x - 1) ** 3, x))  # roots可以得到根的重数信息

运行结果：

![9eb77e708bddf859f6ba0dcff52203dc.png][]

例：求解如下代数方程组

![Image 1][]

from sympy.abc import x, y
    from sympy import solve
    
    s = solve([x ** 2 + y ** 2 - 1, x - y], [x, y])
    print("方程组的解为：", s)

运行结果：

![e21931fe83b29b685f95dcee303d7fa6.png][]

例：求函数

![Image 1][] 的驻点，并求函数在\[0, 1\]上的最大值.

驻点：函数的一阶导数为0的点(驻点也称为稳定点，临界点)。对于多元函数，驻点是所有一阶偏导数都为零的点。驻点（数学概念）\_百度百科

from sympy import *
    
    x = symbols('x')
    y = 2 * x ** 3 - 5 * x ** 2 + x
    x0 = solve(diff(y, x), x)  # 求驻点
    print("驻点的精确解为", x0)
    print("驻点的浮点数表示为：", [x0[i].n() for i in range(len(x0))])  # 列表中的符号数无法整体转换为浮点数
    y0 = [y.subs(x, 0), y.subs(x, 1), y.subs(x, x0[0]).n()]  # 代入区间端点和一个驻点的值
    print("三个点的函数值分别为:", y0)

运行结果：

![03f7b3949fde3fd94f3d5b122c869bf6.png][]

7）求微分方程（方程组）的符号解

例：求下列微分方程的通解：

（1）齐次方程：

![Image 1][] ;

（2）非齐次方程：

![Image 1][] .

from sympy import *
    
    x = symbols('x')
    y = symbols('y', cls=Function)
    eq1 = diff(y(x), x, 2) - 5 * diff(y(x), x) + 6 * y(x)
    eq2 = diff(y(x), x, 2) - 5 * diff(y(x), x) + 6 * y(x) - x * exp(2 * x)
    print("齐次方程的解为：", dsolve(eq1, y(x)))
    print("非齐次方程的解为：", dsolve(eq2, y(x)))

运行结果：

![814ea8da32cabdc7691b73ad4f4437f9.png][]

例：求下列微分方程的解：

（1）初值问题：

![Image 1][] ;

（2）边值问题：

![Image 1][] .

from sympy import *
    
    x = symbols('x');
    y = symbols('y', cls=Function)
    eq1 = diff(y(x), x, 2) - 5 * diff(y(x), x) + 6 * y(x)
    eq2 = diff(y(x), x, 2) - 5 * diff(y(x), x) + 6 * y(x) - x * exp(2 * x)
    print("初值问题的解为：{}".format(dsolve(eq1, y(x), ics={
         y(0): 1, diff(y(x), x).subs(x, 0): 0})))
    y2 = dsolve(eq2, y(x), ics={
         y(0): 1, y(2): 0})
    print("边值问题的解为：{}".format(y2))

运行结果：

![98448912353d16b5563e7de6d72d7b22.png][]

### 5.高等数学问题的数值解 ###

大多数实际问题是无法求符号解的，只能求数值解，即近似解。

1）泰勒级数与数值导数

**泰勒级数：**

![Image 1][] 的泰勒级数展开式为

![Image 1][] .

例：画出

![Image 1][] 及它在0点处的1,3,5阶泰勒展开式在

![Image 1][] 时的图形.

编写如下函数mysin求

![Image 1][] 的近似值.调用自定义函数mysin画泰勒展开式的图形，代码如下：

import numpy as np
    import matplotlib.pyplot as plt
    
    
    def fac(n):
        return 1 if n < 1 else n * fac(n - 1)
    
    
    def item(n, x):
        return (-1) ** n * x ** (2 * n + 1) / fac(2 * n + 1)
    
    
    def mysin(n, x):
        return 0 if n < 0 else mysin(n - 1, x) + item(n, x)
    
    
    x = np.linspace(-2 * np.pi, 2 * np.pi, 101)
    plt.plot(x, np.sin(x), '*-')
    str = ['v-', 'H--', '-.']
    for n in [1, 2, 3]: plt.plot(x, mysin(2 * n - 1, x), str[n - 1])
    plt.legend(['sin', 'n=1', 'n=3', 'n=5'])
    plt.savefig('figure3_16.png', dpi=500)
    plt.show()

![b20b90205dbe4acaa4cd094b56afbd82.png][]

**数值导数**：

利用泰勒级数可以给出近似计算函数导数的方法.例如，若

![Image 1][] 存在一阶导数，则由泰勒级数

![Image 1][] ，移项并舍弃高阶无穷小，得：

![Image 1][]

这是一个常用的用于估计函数一阶导数的计算公式，它具有一阶精度.此外，还可以推出具有二阶精度的估计公式：

![Image 1][]

当函数具有更高阶的导数时，如利用

![Image 1][]

以及

![Image 1][]

可得二阶导数计算公式：

![Image 1][]

例：甲、乙、丙、丁4个人分别位于起始位置(-200,200),(200,200),(200,-200),(-200,-200)处（单位：米），并且以恒定的速率1m/s行走，在行走过程中，甲始终朝向乙的当前位置；同样，乙朝向丙、丙朝向丁、丁朝向甲.试绘制4人行走过程的近似轨迹。

import numpy as np, numpy.linalg as ng
    import matplotlib.pyplot as plt
    
    N = 4
    v = 1.0
    d = 200.0
    time = 400.0
    divs = 201
    xy = np.array([[-d, d], [d, d], [d, -d], [-d, -d]])
    T = np.linspace(0, time, divs)
    dt = T[1] - T[0]
    xyn = np.empty((4, 2))
    Txy = xy
    for n in range(1, len(T)):
        for i in [0, 1, 2, 3]:
            j = (i + 1) % 4
            dxy = xy[j] - xy[i]
            dd = dxy / ng.norm(dxy)  # 单位化向量
            xyn[i] = xy[i] + v * dt * dd  # 计算下一步的位置
        Txy = np.c_[Txy, xyn]
        xy = xyn
    for i in range(N): 
        plt.plot(Txy[i, ::2], Txy[i, 1::2])
        
    plt.savefig("figure3_17.png", dpi=500)
    plt.show()

![8b5c0b1d8f1e01fd1ba0d570733ddd9b.png][]

2）数值积分

分别使用梯形公式、辛普森公式和SciPy工具库中的函数quad求定积分

![Image 1][] 的数值解.

import numpy as np
    from scipy.integrate import quad
    
    
    def trapezoid(f, n, a, b):  # 定义梯形公式的函数
        xi = np.linspace(a, b, n)
        h = (b - a) / (n - 1)
        return h * (np.sum(f(xi)) - (f(a) + f(b)) / 2)
    
    
    def simpson(f, n, a, b):  # 定义辛普森公式的函数
        xi, h = np.linspace(a, b, 2 * n + 1), (b - a) / (2.0 * n)
        xe = [f(xi[i]) for i in range(len(xi)) if i % 2 == 0]
        xo = [f(xi[i]) for i in range(len(xi)) if i % 2 != 0]
        return h * (2 * np.sum(xe) + 4 * np.sum(xo) - f(a) - f(b)) / 3.0
    
    
    a = 0
    b = 1
    n = 1000
    f = lambda x: np.sin(np.sqrt(np.cos(x) + x ** 2))
    
    print("梯形积分I1=", trapezoid(f, n, a, b))
    print("辛普森积分I2=", simpson(f, n, a, b))
    print("SciPy积分I3=", quad(f, a, b))

运算结果：

![fd2f5c4a52d96ef84d5eed7d14066c33.png][]

例：分别求下列积分的数值解.

（1）

![Image 1][] ；（2）

![Image 1][] .

对于（2）中的二重积分，先要化成累次积分

![Image 1][]

import numpy as np
    from scipy.integrate import dblquad
    
    f1 = lambda y, x: x * y ** 2  # 第一个被积函数
    print("I1：", dblquad(f1, 0, 2, 0, 1))
    
    f2 = lambda y, x: np.exp(-x ** 2 / 2) * np.sin(x ** 2 + y)
    bd = lambda x: np.sqrt(1 - x ** 2)
    print("I2:", dblquad(f2, -1, 1, lambda x: -bd(x), bd))

运行结果：

![8b020774ae1f8e655528f4451428dfa8.png][]

例：计算

![Image 1][] ，其中

![Image 1][] 为柱面

![Image 1][] 与

![Image 1][] 两平面所围成的空间区域.

先把三重积分化成累次积分

![Image 1][] .

import numpy as np
    from scipy.integrate import tplquad
    
    f = lambda z, y, x: z * np.sqrt(x ** 2 + y ** 2 + 1)
    ybd = lambda x: np.sqrt(2 * x - x ** 2)
    print("I=", tplquad(f, 0, 2, lambda x: -ybd(x), ybd, 0, 6))

运行结果：

![fd182aab76e05cc22b7790d210172f6d.png][]

3）非线性方程（组）数值解

例：求方程

![Image 1][] 实根的近似值，使误差不超过

![Image 1][] .要求用三种方法：（1）二分法；（2）牛顿迭代法；（3）直接调用SciPy工具库求解.

import numpy as np
    from scipy.optimize import fsolve
    
    
    def binary_search(f, eps, a, b):  # 二分法函数
        c = (a + b) / 2
        while np.abs(f(c)) > eps:
            if f(a) * f(c) < 0:
                b = c
            else:
                a = c
            c = (a + b) / 2
        return c
    
    
    def newton_iter(f, eps, x0, dx=1E-8):  # 牛顿迭代法函数
        def diff(f, dx=dx):  # 求数值导数函数
            return lambda x: (f(x + dx) - f(x - dx)) / (2 * dx)
    
        df = diff(f, dx)
        x1 = x0 - f(x0) / df(x0)
        while np.abs(x1 - x0) >= eps:
            x1, x0 = x1 - f(x1) / df(x1), x1
        return x1
    
    
    f = lambda x: x ** 3 + 1.1 * x ** 2 + 0.9 * x - 1.4
    print("二分法求得的根为：", binary_search(f, 1E-6, 0, 1))
    print("牛顿迭代法求得的根为：", newton_iter(f, 1E-6, 0))
    print("直接调用SciPy求得的根为：", fsolve(f, 0))

运行结果：

![3034bd97a4952f9e6ead43e7bb4baf10.png][]

**用fsolve求非线性方程组的数值 解**

例：求下列非线性方程组的数值解.

![Image 1][]

from numpy import sin
    from scipy.optimize import fsolve
    
    f = lambda x: [5 * x[1] + 3, 4 * x[0] ** 2 - 2 * sin(x[1] * x[2]), x[1] * x[2] - 1.5]
    print("result=", fsolve(f, [1.0, 1.0, 1.0]))

运行结果：

![f1ba84bca669877793e3e7e176c1588c.png][]

或：

from numpy import sin
    from scipy.optimize import fsolve
    
    
    def Pfun(x):
        x1, x2, x3 = x.tolist()  # x转换成列表
        return 5 * x2 + 3, 4 * x1 ** 2 - 2 * sin(x2 * x3), x2 * x3 - 1.5
    
    
    print("result=", fsolve(Pfun, [1.0, 1.0, 1.0]))

4）函数极值点的数值解

求函数

![Image 1][] 在区间\[0,3\]上的极小点.

from numpy import exp, cos
    from scipy.optimize import fminbound
    
    f = lambda x: exp(x) * cos(2 * x)
    x0 = fminbound(f, 0, 3)
    print("极小点为：{}，极小值为：{}".format(x0, f(x0)))

运行结果：

![b1063701c724398e700cd0135b49bf3d.png][]

求函数

![Image 1][] 在0附近的一个极小点.

from numpy import exp, cos
    from scipy.optimize import fmin
    
    f = lambda x: exp(x) * cos(2 * x)
    x0 = fmin(f, 0)
    print("极小点为：{}，极小值为：{}".format(x0, f(x0)))

运行结果：

![68f0a30374845598e60b48d3aae867ed.png][]

多元函数的极值点：

求函数

![Image 1][] 的极小值.

from scipy.optimize import minimize
    
    f = lambda x: 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2
    x0 = minimize(f, [2.0, 2.0])  # [2.0, 2.0]为初始猜测值
    print("极小点为：{}，极小值为：{}".format(x0.x, x0.fun))

运行结果：

![2c44f905c0b7d6d896d919cb30fac29d.png][]

### 6.线性代数问题的符号解和数值解 ###

1）线性代数问题的符号解

**矩阵的运算**

已知

![Image 1][] ，则

![Image 1][] .

例：已知

![Image 1][] ,求

![Image 1][] （向量

![Image 1][] 的长度），

![Image 1][] .

import sympy as sp
    
    A = sp.Matrix([[1], [2], [3]])  # 列向量,即3×1矩阵
    B = sp.Matrix([[4], [5], [6]])
    print("A的模为：", A.norm())
    print("A的模的浮点数为：", A.norm().evalf())
    print("A的转置矩阵为：", A.T)
    print("A和B的点乘为：", A.dot(B))
    print("A和B的叉乘为：", A.cross(B))

运行结果：

![b40d36bac6c86e65b5fba3a413d12735.png][]

例：已知

![Image 1][] ,

![Image 1][] ,求

![Image 1][] ，A的秩R(A)，

![Image 1][] ,

![Image 1][] ,构造分块矩阵

![Image 1][] ，

![Image 1][] （提出A的左上角子块构成的矩阵），

![Image 1][] （删除A的第4行得到的矩阵）.

import sympy as sp
    import numpy as np
    
    A = sp.Matrix(np.arange(1, 17).reshape(4, 4))
    B = sp.eye(4)
    print("A的行列式为：", sp.det(A))
    print("A的秩为：", A.rank())
    print("A的转置矩阵为：", A.transpose())  # 等价于A.T
    print("所求的逆阵为：", (A + 10 * B).inv())
    print("A的平方为：", A ** 2)
    print("A,B的乘积为：", A * B)
    print("横连矩阵为：", A.row_join(B))
    print("纵连矩阵为：", A.col_join(B))
    print("A1为：", A[0:2, 0:2])
    A2 = A.copy()
    A2.row_del(3)
    print("A2为：", A2)

运行结果：

![8edf442811ffa90f2f2db13e35031e38.png][]

**解线性方程组：**

**例:求下列线性方程组的符号解.**

![Image 1][]

解 记上述线性方程组为

![Image 1][] ，可以验证系数矩阵A的秩R(A)=4，所以线性方程组有唯一解

![Image 1][] ，利用Python软件求解.

import sympy as sp
    
    A = sp.Matrix([[2, 1, -5, 1], [1, -3, 0, -6], [0, 2, -1, 2], [1, 4, -7, 6]])
    b = sp.Matrix([8, 6, -2, 2]);
    b.transpose()
    print("系数矩阵A的秩为：", A.rank())
    print("线性方程组的唯一解为：", A.inv() * b)

运行结果：

![05f5150cdd7247196bc6b421ffadbaa5.png][]

**例:求下列齐次线性方程组的基础解系.**

![Image 1][]

import sympy as sp
    
    A = sp.Matrix([[1, -5, 2, -3], [5, 3, 6, -1], [2, 4, 2, 1]])
    print("A的零空间（即基础解系）为：", A.nullspace())

运行结果：

![695cba845156af4f8016d6f96068070a.png][]

即求得的基础解系为：

![Image 1][]

**例:求下列非齐次线性方程组的通解.**

![Image 1][]

求通解要先使用rref()方法把增广矩阵化为行最简.

import sympy as sp
    
    A = sp.Matrix([[1, 1, -3, -1], [3, -1, -3, 4], [1, 5, -9, -8]])
    b = sp.Matrix([1, 4, 0])
    b.transpose()
    C = A.row_join(b)  # 构造增广矩阵
    print("增广阵的行最简形为：n", C.rref())

运行结果：

![8430c55fd00e7532f74ef81dcc5df0af.png][]

通过行最简可以写出原方程组的等价方程组为

![Image 1][]

所以方程组的通解为

![Image 1][] .

**例:求下列矩阵的特征值和特征向量.**

![Image 1][] .

import sympy as sp
    
    A = sp.Matrix([[0, -2, 2], [-2, -3, 4], [2, 4, -3]])
    print("A的特征值为：", A.eigenvals())
    print("A的特征向量为：", A.eigenvects())

运行结果：

![9ddd24c31ae39c9606f7519c75e5d4a2.png][]

求得的特征值为

![Image 1][] ,对应于-8的特征向量为

![Image 1][] .

对应于特征值1的两个线性无关 的特征向量为

![Image 1][] .

**例：把下列矩阵相似对角化，即求可逆矩阵P，使得**

**![Image 1][] 为对角阵.**

![Image 1][] .

from sympy import Matrix, diag
    
    A = Matrix([[0, -2, 2], [-2, -3, 4], [2, 4, -3]])
    if A.is_diagonalizable():
        print("A的对角化矩阵为：", A.diagonalize())
    else:
        print("A不能对角化")

运行结果：

![04755f3474af54d2d42c3247846b22f6.png][]

求得的可逆矩阵P和对角阵D分别为

![Image 1][]

2）线性代数问题的数值解

**向量和矩阵的运算：**

例：已知

![Image 1][] ,求

![Image 1][] （向量

![Image 1][] 的长度），

![Image 1][] .

from numpy import arange, cross, inner
    from numpy.linalg import norm
    
    a = arange(1, 4);
    b = arange(4, 7)  # 创建数组
    print("a的二范数为：", norm(a))
    print("a点乘b=", a.dot(b))  # 行向量a乘以列向量b
    print("a,b的内积=", inner(a, b))  # a,b的内积,这里与dot(a,b)等价
    print("a叉乘b=", cross(a, b))

运行结果：

![57c655a23529ba212e3e7e659203ec6c.png][]

例：已知

![Image 1][] ,

![Image 1][] ,求

![Image 1][] ，A的秩R(A)，

![Image 1][] ,

![Image 1][] ,构造分块矩阵

![Image 1][] ，

![Image 1][] （提出A的左上角子块构成的矩阵），

![Image 1][] （删除A的第4行得到的矩阵）.

import numpy as np
    import numpy.linalg as LA
    
    A = np.arange(1, 17).reshape(4, 4)
    B = np.eye(4)
    print("A的行列式为：", LA.det(A))
    print("A的秩为：", LA.matrix_rank(A))
    print("A的转置矩阵为：n", A.transpose())  # 等价于A.T
    print("所求的逆阵为：n", LA.inv(A + 10 * B))
    print("A的平方为：n", A.dot(A))
    print("A,B的乘积为：n", A.dot(B))
    print("横连矩阵为：", np.c_[A, B])
    print("纵连矩阵为：", np.r_[A, B])
    print("A1为：", A[0:2, 0:2])
    A2 = A.copy();
    A2 = np.delete(A2, 3, axis=0)
    print("A2为：", A2)

运行结果：

A的行列式为： 4.7331654313261276e-30
    A的秩为： 2
    A的转置矩阵为：
     [[ 1  5  9 13]
     [ 2  6 10 14]
     [ 3  7 11 15]
     [ 4  8 12 16]]
    所求的逆阵为：
     [[ 0.11277778  0.00333333 -0.00611111 -0.01555556]
     [-0.005       0.09       -0.015      -0.02      ]
     [-0.02277778 -0.02333333  0.07611111 -0.02444444]
     [-0.04055556 -0.03666667 -0.03277778  0.07111111]]
    A的平方为：
     [[ 90 100 110 120]
     [202 228 254 280]
     [314 356 398 440]
     [426 484 542 600]]
    A,B的乘积为：
     [[ 1.  2.  3.  4.]
     [ 5.  6.  7.  8.]
     [ 9. 10. 11. 12.]
     [13. 14. 15. 16.]]
    横连矩阵为： [[ 1.  2.  3.  4.  1.  0.  0.  0.]
     [ 5.  6.  7.  8.  0.  1.  0.  0.]
     [ 9. 10. 11. 12.  0.  0.  1.  0.]
     [13. 14. 15. 16.  0.  0.  0.  1.]]
    纵连矩阵为： [[ 1.  2.  3.  4.]
     [ 5.  6.  7.  8.]
     [ 9. 10. 11. 12.]
     [13. 14. 15. 16.]
     [ 1.  0.  0.  0.]
     [ 0.  1.  0.  0.]
     [ 0.  0.  1.  0.]
     [ 0.  0.  0.  1.]]
    A1为： [[1 2]
     [5 6]]
    A2为： [[ 1  2  3  4]
     [ 5  6  7  8]
     [ 9 10 11 12]]

**齐次方程的数值解：**

例:求下列齐次线性方程组的基础解系.

![Image 1][]

import numpy as np
    from scipy.linalg import null_space
    
    A = np.array([[1, -5, 2, -3], [5, 3, 6, -1], [2, 4, 2, 1]])
    print("A的零空间（即基础解系）为：", null_space(A))

运行结果：

![2f9c970bcdab360c281db28160508664.png][]

**非齐次线性方程组的数值解：**

例:求下列线性方程组的数值解.

![Image 1][]

解 记上述线性方程组为

![Image 1][] ，可以验证系数矩阵A的秩R(A)=4，所以线性方程组有唯一解

![Image 1][] ，利用Python软件求解.

import numpy as np
    import numpy.linalg as LA
    
    A = np.array([[2, 1, -5, 1], [1, -3, 0, -6], [0, 2, -1, 2], [1, 4, -7, 6]])
    b = np.array([[8, 6, -2, 2]]);
    b = b.reshape(4, 1)
    print("系数矩阵A的秩为：", LA.matrix_rank(A))
    print("线性方程组的唯一解为：", LA.inv(A).dot(b))  # 使用逆矩阵
    print("线性方程组的唯一解为：", LA.pinv(A).dot(b))  # 使用伪逆
    print("线性方程组的唯一解为：", LA.solve(A, b))  # 利用solve求解

运行结果：

![d2120905594cd94f751db040805c73e0.png][]

例:下列非齐次线性方程组有无穷多解，求它的最小范数解.

![Image 1][]

from numpy import array
    from numpy.linalg import pinv
    
    A = array([[1, 1, -3, -1], [3, -1, -3, 4], [1, 5, -9, -8]])
    b = array([1, 4, 0])
    b.resize(3, 1)
    x = pinv(A).dot(b)  # 求最小范数解
    print("最小范数解为：", x)

运行结果：

![9a27c059061fd945103814f1f873f3ed.png][]

例：求下列矛盾方程组的最小二乘解.

![Image 1][]

from numpy import array
    from numpy.linalg import pinv
    
    A = array([[1, 1], [2, 2], [1, 2]])
    b = array([1, 3, 2])
    b.resize(3, 1)
    x = pinv(A).dot(b)  # 求最小二乘解
    print("最小二乘解为：", x)

运行结果：

![f40cb7aba488201afa289ebf38398ee4.png][]

**特征值与特征向量：**

例:求下列矩阵的特征值和特征向量.

![Image 1][] .

import numpy as np
    from numpy.linalg import eig
    
    A = np.array([[0, -2, 2], [-2, -3, 4], [2, 4, -3]])
    values, vectors = eig(A)
    print("A的特征值为：", values)
    print("A的特征向量为：", vectors)

运行结果：

![f76666b7537b4f89bc25d0c2b382d624.png][]

例：对于矩阵

![Image 1][] ,验证

![Image 1][] .

from numpy import array, dot
    from numpy.linalg import eig, inv
    
    A = array([[0, -2, 2], [-2, -3, 4], [2, 4, -3]])
    values, vectors = eig(A)
    check = dot(inv(vectors), A).dot(vectors)
    print("check=n", check)

运行结果：

![e04f73d13d0cd562771cfad056ce0294.png][]

3）求超定线性方程组的最小二乘解 参考：pyhton scipy最小二乘法（scipy.linalg.lstsq模块） numpy的ones\_like函数\_christne1225i的专栏-CSDN博客 np.c\_和np.r\_的用法解析\_demonszjin-CSDN博客

求超定线性方程组

![Image 1][] 的最小二乘解，相当于给定x的观测值

![Image 1][] ,y的观测值

![Image 1][] ,拟合经验函数

![Image 1][] .

import numpy as np
    import numpy.linalg as LA
    from matplotlib.pyplot import plot, rc, legend, show, savefig
    
    x = np.array([0, 1, 2, 3])
    y = np.array([-1, 0.2, 0.9, 2.1])
    A = np.c_[x, np.ones_like(x)]
    m, c = LA.lstsq(A, y, rcond=None)[0]
    print(m, c)
    rc('font', size=16)
    plot(x, y, 'o', label='原始数据', markersize=5)
    plot(x, m * x + c, 'r', label='拟合直线')
    rc('font', family='SimHei')  # 用来正常显示中文标签
    rc('axes', unicode_minus=False)  # 用来正常显示负号
    legend()
    savefig("figure3_41.png", dpi=500)
    show()

运行结果：

![c961cdf1d694b365aca8f27e7f60f8cc.png][]

![3b3677705effa4f6e304b8f3ed68035f.png][]

例：为了测量刀具的磨损速度，做这样的实验：经过一段时间（如每隔1小时），测量一次刀具的厚度，得到一组实验数据

![Image 1][] 如下表所示，试用实验数据拟合经验公式

![Image 1][] .

import numpy as np
    import numpy.linalg as LA
    import matplotlib.pyplot as plt
    
    t = np.arange(8)
    y = np.array([27.0, 26.8, 26.5, 26.3, 26.1, 25.7, 25.3, 24.8])
    A = np.c_[t, np.ones_like(t)]
    ab = LA.lstsq(A, y, rcond=None)[0]  # 返回值为向量
    print(ab)
    plt.rc('font', size=16)
    plt.plot(t, y, 'o', label='Original data', markersize=5)
    plt.plot(t, A.dot(ab), 'r', label='Fitted line')
    plt.legend()
    plt.savefig("figure3_42.png", dpi=500)
    plt.show()

运行结果：

![b9ff20dbad7ccb58240bc4ca70e38aa7.png][]

![f1d5f9bff12874d607ee13518fcee466.png][]

[e1222b3ecb4f775be6ce1543da1e1b4a.png]: /images/20221119/f2778d0e1be94173876314759890051c.png
[10e1497f596e5ff83879559ab73ee616.png]: https://img-blog.csdnimg.cn/img_convert/10e1497f596e5ff83879559ab73ee616.png
[d5796704503d169f1fd0b2d40615481d.png]: https://img-blog.csdnimg.cn/img_convert/d5796704503d169f1fd0b2d40615481d.png
[Image 1]: 
[ebb8e97d9bff00de0a6f4b58a03c5a74.png]: https://img-blog.csdnimg.cn/img_convert/ebb8e97d9bff00de0a6f4b58a03c5a74.png
[9a9b94a6e0db7bda0ea77bd805c89453.png]: https://img-blog.csdnimg.cn/img_convert/9a9b94a6e0db7bda0ea77bd805c89453.png
[848d54a08dc066089ded910f50375a32.png]: https://img-blog.csdnimg.cn/img_convert/848d54a08dc066089ded910f50375a32.png
[53c322c8cd33d3f9993d5038f7392fb4.png]: https://img-blog.csdnimg.cn/img_convert/53c322c8cd33d3f9993d5038f7392fb4.png
[cd1544536b01a31571a26991c9cdedf6.png]: https://img-blog.csdnimg.cn/img_convert/cd1544536b01a31571a26991c9cdedf6.png
[6d11f6358ed85a994696bc49b500b0ea.png]: https://img-blog.csdnimg.cn/img_convert/6d11f6358ed85a994696bc49b500b0ea.png
[026c12a4f4c7566e2aa1486307540493.png]: https://img-blog.csdnimg.cn/img_convert/026c12a4f4c7566e2aa1486307540493.png
[fb327b977a3b3aac18548fab3cd6534d.png]: https://img-blog.csdnimg.cn/img_convert/fb327b977a3b3aac18548fab3cd6534d.png
[9eb77e708bddf859f6ba0dcff52203dc.png]: https://img-blog.csdnimg.cn/img_convert/9eb77e708bddf859f6ba0dcff52203dc.png
[e21931fe83b29b685f95dcee303d7fa6.png]: https://img-blog.csdnimg.cn/img_convert/e21931fe83b29b685f95dcee303d7fa6.png
[03f7b3949fde3fd94f3d5b122c869bf6.png]: https://img-blog.csdnimg.cn/img_convert/03f7b3949fde3fd94f3d5b122c869bf6.png
[814ea8da32cabdc7691b73ad4f4437f9.png]: https://img-blog.csdnimg.cn/img_convert/814ea8da32cabdc7691b73ad4f4437f9.png
[98448912353d16b5563e7de6d72d7b22.png]: https://img-blog.csdnimg.cn/img_convert/98448912353d16b5563e7de6d72d7b22.png
[b20b90205dbe4acaa4cd094b56afbd82.png]: https://img-blog.csdnimg.cn/img_convert/b20b90205dbe4acaa4cd094b56afbd82.png
[8b5c0b1d8f1e01fd1ba0d570733ddd9b.png]: https://img-blog.csdnimg.cn/img_convert/8b5c0b1d8f1e01fd1ba0d570733ddd9b.png
[fd2f5c4a52d96ef84d5eed7d14066c33.png]: https://img-blog.csdnimg.cn/img_convert/fd2f5c4a52d96ef84d5eed7d14066c33.png
[8b020774ae1f8e655528f4451428dfa8.png]: https://img-blog.csdnimg.cn/img_convert/8b020774ae1f8e655528f4451428dfa8.png
[fd182aab76e05cc22b7790d210172f6d.png]: https://img-blog.csdnimg.cn/img_convert/fd182aab76e05cc22b7790d210172f6d.png
[3034bd97a4952f9e6ead43e7bb4baf10.png]: https://img-blog.csdnimg.cn/img_convert/3034bd97a4952f9e6ead43e7bb4baf10.png
[f1ba84bca669877793e3e7e176c1588c.png]: https://img-blog.csdnimg.cn/img_convert/f1ba84bca669877793e3e7e176c1588c.png
[b1063701c724398e700cd0135b49bf3d.png]: https://img-blog.csdnimg.cn/img_convert/b1063701c724398e700cd0135b49bf3d.png
[68f0a30374845598e60b48d3aae867ed.png]: https://img-blog.csdnimg.cn/img_convert/68f0a30374845598e60b48d3aae867ed.png
[2c44f905c0b7d6d896d919cb30fac29d.png]: https://img-blog.csdnimg.cn/img_convert/2c44f905c0b7d6d896d919cb30fac29d.png
[b40d36bac6c86e65b5fba3a413d12735.png]: https://img-blog.csdnimg.cn/img_convert/b40d36bac6c86e65b5fba3a413d12735.png
[8edf442811ffa90f2f2db13e35031e38.png]: https://img-blog.csdnimg.cn/img_convert/8edf442811ffa90f2f2db13e35031e38.png
[05f5150cdd7247196bc6b421ffadbaa5.png]: https://img-blog.csdnimg.cn/img_convert/05f5150cdd7247196bc6b421ffadbaa5.png
[695cba845156af4f8016d6f96068070a.png]: https://img-blog.csdnimg.cn/img_convert/695cba845156af4f8016d6f96068070a.png
[8430c55fd00e7532f74ef81dcc5df0af.png]: https://img-blog.csdnimg.cn/img_convert/8430c55fd00e7532f74ef81dcc5df0af.png
[9ddd24c31ae39c9606f7519c75e5d4a2.png]: https://img-blog.csdnimg.cn/img_convert/9ddd24c31ae39c9606f7519c75e5d4a2.png
[04755f3474af54d2d42c3247846b22f6.png]: https://img-blog.csdnimg.cn/img_convert/04755f3474af54d2d42c3247846b22f6.png
[57c655a23529ba212e3e7e659203ec6c.png]: https://img-blog.csdnimg.cn/img_convert/57c655a23529ba212e3e7e659203ec6c.png
[2f9c970bcdab360c281db28160508664.png]: https://img-blog.csdnimg.cn/img_convert/2f9c970bcdab360c281db28160508664.png
[d2120905594cd94f751db040805c73e0.png]: https://img-blog.csdnimg.cn/img_convert/d2120905594cd94f751db040805c73e0.png
[9a27c059061fd945103814f1f873f3ed.png]: https://img-blog.csdnimg.cn/img_convert/9a27c059061fd945103814f1f873f3ed.png
[f40cb7aba488201afa289ebf38398ee4.png]: https://img-blog.csdnimg.cn/img_convert/f40cb7aba488201afa289ebf38398ee4.png
[f76666b7537b4f89bc25d0c2b382d624.png]: https://img-blog.csdnimg.cn/img_convert/f76666b7537b4f89bc25d0c2b382d624.png
[e04f73d13d0cd562771cfad056ce0294.png]: https://img-blog.csdnimg.cn/img_convert/e04f73d13d0cd562771cfad056ce0294.png
[c961cdf1d694b365aca8f27e7f60f8cc.png]: https://img-blog.csdnimg.cn/img_convert/c961cdf1d694b365aca8f27e7f60f8cc.png
[3b3677705effa4f6e304b8f3ed68035f.png]: https://img-blog.csdnimg.cn/img_convert/3b3677705effa4f6e304b8f3ed68035f.png
[b9ff20dbad7ccb58240bc4ca70e38aa7.png]: https://img-blog.csdnimg.cn/img_convert/b9ff20dbad7ccb58240bc4ca70e38aa7.png
[f1d5f9bff12874d607ee13518fcee466.png]: https://img-blog.csdnimg.cn/img_convert/f1d5f9bff12874d607ee13518fcee466.png