3.5.1泰勒级数与数值导数

1.泰勒级数

2.数值导数

3.5.2数值分析

1.一重积分

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

1.二分法

2.牛顿迭代法

3.用SciPy工具库求解非线性方程（方程组）

4.用fslove求非线性方程的数值解

3.5.4 函数极值点的数值解

1.一元函数的极值点

2.多元函数的极值点

3.6.1线性代数问题的符号解

1.矩阵的运算

2.解线性方程

3.特征值与特征向量

3.6.2 线性代数问题的数值解

1. 向量和矩阵的运算

2.齐次线性方程组的数值解

3.非齐次线性方程组的数值解

1）唯一解

2）多解形式

3）无解形式

4.特征值与特征向量

3.6.3求超定性方程组的最小二乘解

3.5.1泰勒级数与数值导数

1.泰勒级数

sin（x）的泰勒技术的展开式为

$sinx =\sum_{k=0}^{oo}\frac{(-1)^{k}(x)^{2k+1}}{\sum _{k=0}2k+1}$ , $x\in (-oo,+oo)$


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)  # 定义sin(x) 在k=n的函数展开def mysin(n, x): return (0 if n < 0 else mysin(n - 1, x) + item(n, x))  # 定义多项n阶求和x = np.linspace(-2 * np.pi, 2 * np.pi, 101)  # 分段从-2pi到2pi 分为101份
plt.plot(x, np.sin(x), '*-')  # 画出sin(x)的标准图像
str = ['-v', 'H--', '-.']  # 设置三种不同的图像表示  第0个，第1个，第2个
for n in [1, 2, 3]: plt.plot(x, mysin(2 * n - 1, x), str[n - 1])  # 将n=1,2,3带入，得到1，3，5阶泰勒展开式
plt.legend(['sin', 'n=1', 'n=3', 'n=5'])  # 实例
plt.savefig('figure3_16.png', dpi=500)  # 保存png形式图片，像素
plt.show()

2.数值导数

利用泰勒级数可以给出近似计算函数导数的方法.例如,若f(x)存在一阶导数,则由泰勒级数:

$f(x+\Delta x)=f(x)+f'(x)\Delta x+\circ (\Delta x)$

移项并舍弃高阶无穷小，得 $f'(x)\approx \frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}$

当函数具有更高阶导数时，如利用： $f(x+\Delta x)=f(x)+f'(x)\Delta x+\frac{f''(x)}{2!}(\Delta x)^2+\circ ((\Delta x)^2)$

以及： $f(x-\Delta x)=f(x)-f'(x)\Delta x+\frac{f''(x)}{2!}(\Delta x)^2+\circ ((\Delta x)^2)$

可得二阶导数的计算公式： $f''(x)\approx \frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{(\Delta x)^2}$

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

分析：在运动学中，速度是位移相对于时间的导数，即：

                                                                 $\nu(t)= \frac{d}{dt}r(t)$

因此，在一段很短的时间内,近似的有：

                                                         $r(t+\Delta t)\approx r(t)+v(t)\cdot \Delta t$

成立，又由于位移，速度均是矢量，因此在平面内，又有                

                                         $r_{x}(t+\Delta t)\approx r_{x}(t)+\nu (t)\cdot \Delta t\cdot cos\theta (t)$

                                 $r_{y}(t+\Delta t)\approx r_{y}+v(t)\cdot \Delta t\cdot sin\theta (t)$



联立二式，其中 $\theta (t)$ 为t时刻与x轴正向的夹角

以两个二维数组xy, xyn分别存储4个人的当前位置和下一时刻的位置,具体地说,第i个人的当前位置为xy[i]，下一时刻的位置为xyn[i],其中i取0,1,2,3时分别对应甲、乙、丙、丁。

下面语句
                        j=(i+1)%4;dxy=xy[j]-xy[i]
                        dd=dxy/ng.norm(dxy） #单位化向量

就完成了对夹角余弦值、正弦值的计算.
二维数组Txy存放4个人的所有位置信息,其中 Txyl]存放第i个人所有时刻的位置信息.具体程序代码如下:

# 2.数值导数import numpy as np, numpy.linalg as ng
import matplotlib.pyplot as pltN = 4
v = 1.0
d = 200.0
time = 400.0  # 的时间是随机的，可以用二分法来改变time的值来估计4人汇聚的时间
divs = 201
xy = np.array([[-d, d], [d, d], [d, -d], [-d, -d]])  # 记录4个人的当前位置
xyn = np.empty((4, 2))  # 记录4个人的下一时刻位置4*2的矩阵来记录
T = np.linspace(0, time, divs)  # 从0到指定的时间time分成200份
dt = T[1] - T[0]  # 求出时间间隔Txy = xy  # 将xy中数据保存到Txy中，Txy记录的是最终时刻的数据
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)  # 单位化向量 norm()求范数，相当于求模xyn[i] = xy[i] + v * dt * dd  # 下一跳的位置Txy = np.c_[Txy, xyn]  # 返回一个迭代器，生成数组坐标和值对。xy = xyn.copy()  # 将xyn浅拷贝给xy
for i in range(N): plt.plot(Txy[i, ::2], Txy[i, 1::2])#对于i,从0开始步长2为2，从1开始步长为2
plt.savefig("figure3_17.png", dpi=500)
plt.show()

3.5.2数值分析

1.一重积分

当一元函数在区间[a,b]不变号时，定积分的值恰好等于x=a,x=b于x轴所围成的曲边梯形的”有向面积“。

所以使用梯形法：把大的曲边梯形剖分成多个小的曲边梯形，每一个小曲边梯形的面积用一个梯形米娜及做出近似，最后累加求和得到所求定积分的数值解

下面直接给出积分 $\int_{a}^{b}f(x)dx$ 的梯形计算公式，和辛普森计算公式：

使用梯形公式时,把整个区间[a,b]进行n等分(此时步长h= (b-a)/n),分点为 $x_{i}$ (i=0,1,… ,n),对应的函数值 $f_{i}=f(x_{i})$ ,这里 $x_i=a+i\cdot h,i=0,1,\cdots ,n$ .积分 $\int_{a}^{b}f(x)dx$ 的数值解：

$I_{1}=\frac{h}{2}(f_{0}+2\sum_{i=1}^{n-1}f_{i}+f_{n})$

使用辛普森公式时，将区间[a,b]讲行2n等分，步长h= $\frac{b-a}{2n}$ ，分点 $x_i=a+i\cdot h,i=0,1,\cdots ,2n$ ，对应的函数值 $f_{i}=f(x_{i})$ .积分 $\int_{a}^{b}f(x)dx$ 的数值解:

$I_{2}=\frac{h}{3}(f_{0}+2\sum_{i=1}^{n-1}f_{2i}+4\sum_{i=1}^{n}f_{2i-1}+f_{2n})$

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

$\int_{0}^{1}\sin \left ( \sqrt{\cos +x^{2}} \right )dx$

import numpy as np
from scipy.integrate import quad  # 导入SciPy库中解决定积分的函数def trapezoid(f, n, a, b):xi = np.linspace(a, b, n)  # 将传入的上界和下界传入参数，分成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]  # 对于一维数组中的所有元素判断奇偶，然后放入xe中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.0a = 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))

关于多重疾风使用SciPy库中的函数dblquad,tplquad直接求数值解，dblquad的调用格式为

dblquad(func,a,b,gfun,hfun,args(),epsabs=1.49e-08,epsrel=1.49e-08)

其中被积函数func的格式为func(y,x),最外层x的积分区域为[a,b],内层y的积分区间为[gfun(x),hfun(x)].

tplquad的调用格式为

tplquad(func,a,b,gfun,hfun,rfun,args=(),epsabs=1.49r-08,epsrel=1.49e-8)

其中被积函数func的格式为func(z,x,y),最外层x的季风区间为[a,b]，中间层y的积分区间为[gfun(x),hfun(x)],最内层z的季风区间为[qfun(x,y),rfun(x,y)]

例3.19分别求下列积分的数值解

（1） $\int_{0}^{2}dx\int_{0}^{1}xy^{2}dy$ (2) $\iint e\tfrac{-x^{2}}{2}\sin (x^{2}+y)dxdy$

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

$I=\int_{-1}^{1}dx\int_{-\sqrt{1-x^{2}}}^{\sqrt{1-x^{2}}}e\tfrac{-x^{2}}{2}sin(x^{2}+y)dy$

# 2.多重积分
import numpy as np
from scipy.integrate import dblquadf1 = 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))  # 传入参数

I1 (0.6666666666666667, 7.401486830834377e-15)#后面的返回值为误差值
I2 (0.5368603826989582, 3.696155159715886e-09)

例3.20 计算 $\int \int \int z\sqrt{x^{2}+y^{2}+1}dxdydz$ ,被积面积为柱面 $x^{2}+y^{2}-2x=0$ 于z=0,z=6两平面所谓成的空间区域

先把三重积分化为累次积分 $\int_{0}^{2}dx\int_{-\sqrt{2x-x^{2}}}^{\sqrt{2x-x^{2}}}dy\int_{0}^{6}z\sqrt{x^{2}+y^{2}+1}dz$

import numpy as np
from scipy.integrate import tplquadf = lambda z, y, x: z * np.sqrt(x ** 2 + y ** 2 + 1)  # 原积分
ybd = lambda x: np.sqrt(2 * x - x ** 2)  # z轴上的范围
print("I=", tplquad(f, 0, 2, lambda x: -ybd(x), ybd, 0, 6))  # 计算多重积分三重积分的时候需要化为累次积分进行计算

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

介绍两种常用的数值解，二分法和牛顿迭代法

1.二分法

若f(x)属于C[a,b]，([a,b]上的连续函数)且f(a)f(b)<0，则有介值定理，存在c属于[a,b]，使得f(c)=0，这时，可以使用二分法对方程进行求根

p104---p105

2.牛顿迭代法

若f(x)属于 $C^{2}$ [a,b]([a,b]上的二阶可微函数)，f(a)f(b)<0且 ${f}'$ (x)在[a,b]上不变号，则方程f(x)=0在[a,b]内必然存在某个根 $x^{*}$ .设x0是 $x^{*}$ 附近的点，则可以根据泰勒展开式

3.用SciPy工具库求解非线性方程（方程组）

例3.21 求方程 $x^{3}+1.1x^{2}+0.9x-1.4=0$ 实根的近似值，使误差不超过 $10^{-6}$ 。要求用三种方法1，二分法，2.牛顿迭代法3.直接调用SciPy工具库求解

# 用SciPy工具库求解非线性方程
import numpy as np
from scipy.optimize import fsolvedef binary_search(f, eps, a, b):  # 定义二分法c = (a + b) / 2while np.abs(f(c)) > eps:if f(a) * f(c) < 0:b = celse:a = cc = (a + b) / 2return cdef 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), x1return x1f = 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))

二分法求得的根为： 0.6706571578979492
牛顿迭代法求得的根为： 0.6706573107258097
直接调用SciPy求得的根为： [0.67065731

4.用fslove求非线性方程的数值解

例3.22 求下列非线性方程的数值解

$5x_{2}+3=0|| 4x_{1}^{2}-2\sin (x_{2}x_{3})=0||x_{2}x_{3}-1.5=0$

from numpy import sin
from scipy.optimize import fsolvef = 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]))  # 设置初始解来近似局部最优解
# 求解方程都根，对于非线性方程，通常会有多个解，因此需要设置解的大致初始值（取值范围），这样方程在初始值附近按梯度下降进行求解，可得局部最优解
# result= [-0.70622057 -0.6        -2.5       ]

上面程序中使用的是匿名函数，但python的下标从0开始，不太方便，下面利用函数定义方程，程序如下

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

3.5.4 函数极值点的数值解

下面利用SciPy库来求极值点

1.一元函数的极值点

例3.23 求函数 $f(x)=e^{x}cos(2x)$ 在区间[0,3]上的极小点

from numpy import exp, cos
from scipy.optimize import fminboundf = lambda x: exp(x) * cos(2 * x)  # 定义匿名函数
x0 = fminbound(f, 0, 3)  # 传入函数 区间
print("极小值点为：{}，极小值为：{}".format(x0, f(x0)))

极小值点为：1.8026199149262752，极小值为：-5.4251652274637722，

2.多元函数的极值点

使用的是scipy.optimize模块下的minimize函数，该函数求多元化函数的极小值和极小点，其调用函数形式为

minimize(fun,x0,args=(),method=None)

例3.25 求函数 $f(x)=100(x_{2}-x_{1}^{2})^{2}+(1-x_{1})^{2}$ 的极小值

from scipy.optimize import minimizef = lambda x: 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2;
x0 = minimize(f, [2.0, 2.0]) # 取值不一样，近似值也不一样
print("极小点为：{},极小值为:{}".format(x0.x, x0.fun))

极小点为：[0.99999565 0.99999129],极小值为:1.8932893809017893e-11

3.6.1线性代数问题的符号解

SymPy线代模块的函数和矩阵操作包括：求行列式的得值，特殊矩阵的构建，求矩阵的特征值，特征向量的转置和逆矩阵利用eye,zero和ones等函数，可以快速构建特殊矩阵，还可以删除矩阵中某些选中的行或列，基本运算符，入+,-,*,**也可以用于矩阵

1.矩阵的运算

基础知识：

已知 $\alpha =[a_1,a_2,a_3]^{T},\beta =[b_1,b_2,b_3]^T$ ，则 $\alpha \cdot \beta =a_1b_1+a_2b_2+a_3b_3$ ， $\alpha \times \beta =\begin{vmatrix} i& j &k \\ a_1& a_2 & a_3\\ b_1&b_2 & b_3 \end{vmatrix}$ (T为矩阵的转置)

例3.26 已知 $\alpha =[1,2,3]^{T},\beta =[4,5,6]^{T}$ ，求 $||\alpha ||_{2}$ ,（向量α）， $\alpha ^{T}$ ， $\alpha \cdot \beta$ ， $\alpha \times \beta$ 。

import sympy as spA = sp.Matrix([[1], [2], [3]])
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))

例3.27 已知 $A=\begin{bmatrix} 1 & 2 & 3 &4 \\ 5& 6 &7 &8 \\ 9& 10 & 11 & 12\\ 13& 14 & 15 & 16 \end{bmatrix}$ ， $B=\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 &0 \\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix}$ ，求 $\left | A \right |$ ，A的秩R（A）， $A^T$ , $(A+10E)^{-1}$ , $A^2,AB$ ，构造分块矩阵[A,B], $\begin{bmatrix} A\\ B \end{bmatrix}$ , $A_1=\begin{bmatrix} 1 & 2\\ 5 & 6 \end{bmatrix}$ (提出A的左上角子块构成的矩阵)， $A_2=\begin{bmatrix} 1 & 2 & 3 & 4\\ 5& 6 & 7 &8 \\ 9 & 10 & 11 & 12 \end{bmatrix}$ （删除A的第四行得到的矩阵）

import sympy as sp
import numpy as npA = sp.Matrix(np.arange(1, 17).reshape(4, 4))  # 建立从1 到 16的数组，转换成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)

A的行列式为： 0
A的秩为： 2
A的转置矩阵： Matrix([[1, 5, 9, 13], [2, 6, 10, 14], [3, 7, 11, 15], [4, 8, 12, 16]])
所求逆矩阵： Matrix([[203/1800, 1/300, -11/1800, -7/450], [-1/200, 9/100, -3/200, -1/50], [-41/1800, -7/300, 137/1800, -11/450], [-73/1800, -11/300, -59/1800, 16/225]])
A的平方为： Matrix([[90, 100, 110, 120], [202, 228, 254, 280], [314, 356, 398, 440], [426, 484, 542, 600]])
A,B的乘积为： Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
横连矩阵为 Matrix([[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]])
纵连矩阵为 Matrix([[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为： Matrix([[1, 2], [5, 6]])
A2为： Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])

2.解线性方程

3.28

$\left\{\begin{matrix} 2x_1+x_2-5x_3+x_4=8\\ x_1-3x_2-6x_4=6 \\ 2x_2-x_3+2x_4=-2\\ x_1+4x_2-7x_3+6x_4=2 \end{matrix}\right.$

解：记上述线性方程组的解为 Ax=b 可以验证系数矩阵 A的秩R（A）=4,所以以线性方程组有唯一解 $x=A^{-1}b$ ,利用Python软件可得

import sympy as spA = 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)

系数矩阵A的秩为： 4
线性方程组的唯一解为： Matrix([[4], [-14/9], [-2/9], [4/9]])

3.29 求下列齐次方程组的基础解系

$\left\{\begin{matrix} x_1-5x_2+2x_3-3x_4=0\\ 5x_1+3x_2+6x_3-x_4=0\\ 2x_1+4x_2+2x_3+x_4=0 \end{matrix}\right.$

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

A的零空间（即基础解系）为： [Matrix([
[-9/7],
[ 1/7],
[ 1],
[ 0]]), Matrix([
[ 1/2],
[-1/2],
[ 0],
[ 1]])]

3.30

$\left\{\begin{matrix} x_1+x_2-3x_3-x_4=1\\ 3x_1-x_2-3x_3+4x_4=4\\ x_1+5x_2-9x_3-8x_4=0 \end{matrix}\right.$


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())

增广矩阵的行最简形为：
(Matrix([
[1, 0, -3/2, 3/4, 5/4],
[0, 1, -3/2, -7/4, -1/4],
[0, 0, 0, 0, 0]]), (0, 1))
求通解首先要用rref()方法把增广阵化成行最简形：

通过上述返回值可以看出，增广阵的列向量组的最大无关组由第1列，第2列组成，对应的行最简形的列是单位坐标向量。求得的行最简形矩阵为： $\begin{bmatrix} 1 & 0 &-\frac{3}{2} & \frac{3}{4} &\frac{5}{4} \\ 0 & 1& -\frac{3}{2} & -\frac{7}{4} &-\frac{1}{4} \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix}$ ，通过行最简形可以写出原方程组的等价方程组为：

$\begin{Bmatrix} x_1=\frac{3}{2}x_3-\frac{3}{4}x_4+\frac{5}{4}\\ x_2=\frac{3}{2}x_3+\frac{7}{4}x_4-\frac{1}{4} \end{Bmatrix}$ 所以，方程组的通解为：

$\begin{bmatrix} x_1\\ x_2\\ x_3\\ x_4 \end{bmatrix}=c_1\begin{bmatrix} \frac{3}{2}\\ \frac{3}{2}\\ 1\\ 0 \end{bmatrix}+c_2\begin{bmatrix} -\frac{3}{4}\\ \frac{7}{4}\\ 0\\ 1 \end{bmatrix}+\begin{bmatrix} \frac{5}{4}\\ -\frac{1}{4}\\ 0\\ 0 \end{bmatrix}$

3.特征值与特征向量

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

$A=\begin{bmatrix} 0 & -2 & 2\\ -2 & -3 & 4\\ 2 & 4 & -3 \end{bmatrix}$


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

A的特征值为： {-8: 1, 1: 2}
A的特征向量为： [(-8, 1, [Matrix([
[-1/2],
[ -1],
[ 1]])]), (1, 2, [Matrix([
[-2],
[ 1],
[ 0]]), Matrix([
[2],
[0],
[1]])])]

求得的特征值为 $\lambda _{1}=-8$ ， $\lambda _{1}=\lambda _{2}=1$ ，对应于特征值-8的特征向量为 $\xi =[-1/2,-1,1]^{T}$

对应于特征值1的两个线性无关的特征向量为 $\xi =[-2,1,0]^{T}$ , $\xi =[2,0,1]^{T}$

例 3.32 把下列矩阵相似对角化，即求得可逆矩阵P，使得 $P^{-1}AP=D$ 为对角阵

$A=\begin{bmatrix} 0 & -2 & 2\\ -2 & -3 &4 \\ 2& 4 & -3 \end{bmatrix}$

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

A的对角化矩阵为：
(Matrix([
[-1, -2, 2],#可逆矩阵P
[-2, 1, 0],
[ 2, 0, 1]]), Matrix([#对角矩阵D
[-8, 0, 0],
[ 0, 1, 0],
[ 0, 0, 1]]))

3.6.2 线性代数问题的数值解

使用Numpy的子模块linalg（线性代数的缩写）

注：在NumPy库的array数组中，*和multiply运算等价，都表示矩阵的逐个元素相乘；@和dot()函数都表示矩阵乘法

1. 向量和矩阵的运算

例 3.33（续3.26）已知 $\alpha = [1,2,3]^T,\beta =[4,5,6]^T$ ，求 $\left \| \alpha \right \|_2$ （向量 $\alpha$ 的长度）， $\alpha ^T,\alpha \cdot \beta ,\alpha \times \beta$ 程序如下

from numpy import arange, cross, inner
from numpy.linalg import norma = arange(1, 4)
b = arange(4, 7)
print("a的二范数：", norm(a))
print("a点乘b=", a.dot(b))
print("a,b的内积为=", inner(a, b))  # a,b的内积，这里和dot(a,b)等价
print("a叉乘b=", cross(a, b))

a的二范数： 3.7416573867739413
a点乘b= 32
a,b的内积为= 32
a叉乘b= [-3 6 -3]

例3.34（续3.27）已知 $A=\begin{bmatrix} 1 & 2 & 3 &4 \\ 5& 6 &7 &8 \\ 9& 10 & 11 & 12\\ 13& 14 & 15 & 16 \end{bmatrix}$ ， $B=\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 &0 \\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix}$ ，求 $\left | A \right |$ ，A的秩R（A）， $A^T$ , $(A+10E)^{-1}$ , $A^2,AB$ ，构造分块矩阵[A,B], $\begin{bmatrix} A\\ B \end{bmatrix}$ , $A_1=\begin{bmatrix} 1 & 2\\ 5 & 6 \end{bmatrix}$ (提出A的左上角子块构成的矩阵)， $A_2=\begin{bmatrix} 1 & 2 & 3 & 4\\ 5& 6 & 7 &8 \\ 9 & 10 & 11 & 12 \end{bmatrix}$ （删除A的第四行得到的矩阵）

import numpy as np
import numpy.linalg as LAA = np.arange(1, 17).reshape(4,4)
B = np.eye(4)  # 单位矩阵
print("A的行列式为：", LA.det(A))
print("A的秩为：", LA.matrix_rank(A))
print("A的转置矩阵：", A.transpose())  # 等价于  A.T
print("所求逆矩阵：", LA.inv(A + 10 * B))
print("A的平方为：", A.dot(A))
print("A,B的乘积为：", 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]]

2.齐次线性方程组的数值解

使用scipy.linalg模块的null_space函数可求齐次线性方程组的基础解系

$\left\{\begin{matrix} x_1-5x_2+2x_3-3x_4=0\\ 5x_1+3x_2+6x_3-x_4=0\\ 2x_1+4x_2+2x_3+x_4=0 \end{matrix}\right.$

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

A的零空间（基础解系）为： [[-0.35455452 0.71506457]
[ 0.41182294 0.03091606]
[-0.05010987 -0.65273305]
[-0.83796298 -0.24832727]]

3.非齐次线性方程组的数值解

1）唯一解

例3.36 求下列非齐次线性方程组的数值解

$\left\{\begin{matrix} 2x_1+x_2-5x_3+x_4=8\\ x_1-3x_2-6x_4=6 \\ 2x_2-x_3+2x_4=-2\\ x_1+4x_2-7x_3+6x_4=2 \end{matrix}\right.$

import numpy as np
import numpy.linalg as LAA = 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))  # 使用逆矩阵  inv()
print("线性方程组的唯一解为", LA.pinv(A).dot(b))  # 使用伪逆   pinv()
print("线性方程组的唯一解为", LA.solve(A,b))  # 使用solve求解

系数矩阵A的秩为： 4
线性方程组的唯一解为 [[ 4. ]
[-1.55555556]
[-0.22222222]
[ 0.44444444]]
线性方程组的唯一解为 [[ 4. ]
[-1.55555556]
[-0.22222222]
[ 0.44444444]]
线性方程组的唯一解为 [[ 4. ]
[-1.55555556]
[-0.22222222]
[ 0.44444444]]

2）多解形式

例3.37（3.30）

$\left\{\begin{matrix} x_1+x_2-3x_3-x_4=1\\ 3x_1-x_2-3x_3+4x_4=4\\ x_1+5x_2-9x_3-8x_4=0 \end{matrix}\right.$

A = np.array([[1, 1, -3, -1], [3, -1, -3, 4], [1, 5, -9, -8]])
b = np.array([1, 4, 0])
b.resize(3, 1)
x = pinv(A).dot(b)
print("最小范数解为", x)

最小范数解为 [[ 0.35040431]
[-0.0916442 ]
[-0.38814016]
[ 0.42318059]]

3）无解形式

例3.38

$\left\{\begin{matrix} x_1+x_2=1\\ 2x_1+2x_2=3\\ x_1+2x_2=2 \end{matrix}\right.$

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

最小二乘解为 [[0.8]
[0.6]]

4.特征值与特征向量

例。3.39求特征值和特征向量

$A=\begin{bmatrix} 0 & -2 & 2\\ -2 & -3 & 4\\ 2 & 4 & -3 \end{bmatrix}$

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

A的特征值为： [ 1. -8. 1.]
A的特征向量为 [[ 0.94280904 -0.33333333 -0.2981424 ]
[-0.23570226 -0.66666667 0.74535599]
[ 0.23570226 0.66666667 0.59628479]]

例3.40（续3.39）对于矩阵

$A=\begin{bmatrix} 0 & -2 & 2\\ -2 & -3 & 4\\ 2& 4 & -3 \end{bmatrix},P=\begin{bmatrix} 0.9428 &-0.3333 &-0.4254 \\ -0.2357 & -0.6667 &0.7454 \\ 0.2357 & 0.6667 &0.5247 \end{bmatrix},\Lambda =\begin{bmatrix} 1 & 0 &0 \\ 0 & -8 &0 \\ 0 & 0 & 1 \end{bmatrix}$

验证 $P^{-1}AP=\Lambda$

rom numpy import dotA = np.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)

check=
[[ 1.00000000e+00 -2.77555756e-16 -1.11022302e-16]
[ 0.00000000e+00 -8.00000000e+00 0.00000000e+00]
[-5.55111512e-17 3.88578059e-16 1.00000000e+00]]

##P矩阵不唯一

3.6.3求超定性方程组的最小二乘解

建模中常用线性最小二乘法，就是求超定线性方程组（未知数个数少，方程个数多）的最小二乘解

例3.41 求超定线性方程组 $\begin{bmatrix} 0 & 1\\ 1& 1\\ 2 &1 \\ 3 & 1 \end{bmatrix}\begin{bmatrix} m\\ c \end{bmatrix}=\begin{bmatrix} -1\\ 0.2\\ 0.9\\ 2.1 \end{bmatrix}$ 的最小二乘解，相当于给定x的观测值 $x=[0,1,2,3]^T$ ，y的观测值 $y=[-1,0.2,0.9,2.1]^T$ ，拟合经验函数y=mx+c

>>> np.ones_like(x)

array([[1, 1, 1],

[1, 1, 1]]

np.c_:将传入的两个参数将切片对象沿第二个轴（按列）连接。

import numpy as mp
import numpy.linalg as LA
from matplotlib.pyplot import plot, rc, legend, show, savefigx = 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的值为：{}".format(m,c))
rc('font', size=16)
plot(x, y, 'o', label='原始数据', markersize=15)
plot(x, m * x + c, 'r', label='拟合直线')
rc('font', family='SimHei')
rc('axes', unicode_minus=False)
legend()
savefig("figure3_41.png", dpi=500)
show()

0.9999999999999999 -0.9499999999999997#分别是m和c的值

例3.42 为了测量刀具的磨损速度,做这样的实验:经过一定时间(如每隔一小时),测量一次刀具的厚度,得到一组实验数据 $(t_i,y_i)$ (i = 1,2,… ,8)如表3.3所示.试用实验数据拟合经验公式y = at＋b.

import numpy as np
import numpy.linalg as LA
import matplotlib.pyplot as pltt = np.arange(8)
y = 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()

[-0.30357143 27.125 ] ##y=-0.30357143t+27.125

高登数学，线性代数问题的数值解（SciPy第三方库，近似解）相关推荐

Python常用第三方库
Python常用第三方库一. 文件读写二.网络抓取和解析三.数据库连接四.数据清洗转换五.数据计算和统计分析六.自然语言处理和文本挖掘七.图像和视频处理八.音频处理九.数据挖掘/机器 ...
python的标识库和第三方库_Python 标准库、第三方库
Python数据工具箱涵盖从数据源到数据可视化的完整流程中涉及到的常用库.函数和外部工具.其中既有Python内置函数和标准库,又有第三方库和工具.这些库可用于文件读写.网络抓取和解析.数据连接.数清 ...
python concat函数多张表_最全Python数据工具箱：标准库、第三方库和外部工具都在这里了 - Mr_YJY...
导读:Python数据工具箱涵盖从数据源到数据可视化的完整流程中涉及到的常用库.函数和外部工具.其中既有Python内置函数和标准库,又有第三方库和工具.这些库可用于文件读写.网络抓取和解析.数据连接 ...
线性代数知识点总结_考研数学线性代数部分怎样复习
前言: 对于线性代数而言,数一.数二.数三的差别并不是很大,所以在这里我就不区分了.在线性代数中,线性方程组和矩阵的相似是考察的重点,并且大家还要注意线性方程组和向量之间的相结合,矩阵的相似和二次型的 ...
数学/线性代数 {子式,余子式,代数余子式,拉普拉斯展开}
数学/线性代数 {子式,余子式,代数余子式,拉普拉斯展开/行列式按{行/列}展开} @LOC_COUNTER=2 子式定义给定行列式A, 从中任意的选择k行和k列, 这些行列的交点, 组成新的 k ...
三阶行列式的题目_考研数学 | 线性代数中的行列式重难点分析
巧儿姐姐叨一叨线性代数是考研高数试卷中的一个重要组成部分,大约占整个试卷中22%的比例.据历年的考察情形来看,线代的题型变化不大,比较容易拿分.这样也就要求大家在线性代数的地方一定要把该拿的分拿到! ...
数学-线性代数3(相关性、基、维数、四个基本子空间)
目录: 九.线性相关性.基.维数 1.线性无关与线性相关 1)背景知识 2)线性无关与线性相关 3)零空间的作用 4)生成空间 2.基 3.维数 4.总结 ---------------- ...
数学/线性代数 {矩阵初等变换,[阶梯形/最简形]矩阵,初等矩阵}
数学/线性代数 {矩阵初等变换,[阶梯形/最简形]矩阵,初等矩阵}; @LOC_COUNTER: 3; 矩阵的初等变换定义矩阵的初等变换和行列式的变换是完全一样的; . LINK: (http ...
数学/线性代数 {行列式, 行列式变换,行列式操作,行列式计算}
数学/线性代数 {行列式, 行列式变换,行列式操作,行列式计算} @LOC_COUNTER: 5 行列式定义给定方形矩阵S [ a b c d ] \begin{bmatrix} a & ...
深度学习中的数学-线性代数
深度学习中的数学-线性代数 1 矩阵和向量相乘 1.1 标准乘积 1.2 元素对应乘积 2 线性相关和生成子空间 3 特征分解 4 奇异值分解推荐书目参考 1 矩阵和向量相乘 1.1 标准乘积如 ...

高登数学，线性代数问题的数值解（SciPy第三方库，近似解）

3.5.1泰勒级数与数值导数

1.泰勒级数

2.数值导数

3.5.2数值分析

1.一重积分

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

1.二分法

2.牛顿迭代法

3.用SciPy工具库求解非线性方程（方程组）

4.用fslove求非线性方程的数值解

3.5.4 函数极值点的数值解

1.一元函数的极值点

2.多元函数的极值点

3.6.1线性代数问题的符号解

1.矩阵的运算

2.解线性方程

3.特征值与特征向量

3.6.2 线性代数问题的数值解

1. 向量和矩阵的运算

2.齐次线性方程组的数值解

3.非齐次线性方程组的数值解

1）唯一解

2）多解形式

3）无解形式

4.特征值与特征向量

3.6.3求超定性方程组的最小二乘解

高登数学，线性代数问题的数值解（SciPy第三方库，近似解）相关推荐

最新文章

热门文章