【学习笔记】第三章 Python在高等数学和线性代数中的应用

常用模块如下：

abc:符号变量模块;

calculus:积分相关方法;

core:基本的加、乘、指数运算等;discrete:离散数学;

functions:基本的函数和特殊函数;galgebra:几何代数;

geometry:几何实体;integrals:符号积分;

interactive:交互会话(如IPython);logic:布尔代数和定理证明;

matrices:线性代数和矩阵;ntheory:数论函数;

physics:物理学;

plotting:用 Pyglet进行二维和三维的画图;polys:多项式代数和因式分解;

printing:漂亮的打印和代码生成;series:级数;

simplify:化简符号表达式;

solvers:方程求解;

stats:统计学.

接下来会针对一些模块进行简要描述：

1.微积分模块(sympy.integrals)

微积分模块支持大量的基础与高级微积分运算功能.例如,支持导数、积分、级数展开、微分方程以及有限差分方程. SymPy还支持积分变换;在微分中,还支持数值微分、复合导数和分数阶导数。

2.离散数学模块(sympy.discrete

离散数学指变量特征是离散的数学分支,与连续变量的数学(微积分)区分开来.它主要处理整数、图形以及逻辑学中的问题.这个模块对二项式系数、乘积与求和运算有完整的支持。

3.方程求解模块(sympy.solvers)

求解器(solvers)是SymPy 中求方程解的模块.这个模块具有解复数多项式以及多项式组的能力.它可以解代数方程、常微分方程、偏微分方程和差分方程。

4.矩阵模块(sympy.matrices)

SymPy具有强大的矩阵与行列式计算的功能．矩阵属于线性代数的分支.SymPy支持矩阵创建，如全0矩阵、全1矩阵、随机矩阵以及矩阵运算.它还支持一些特殊函数，如计算黑塞矩阵(Hessian matrix)的函数、一组向量的格拉姆-施密特(Gram-Schmidt)正交化函数、朗斯基 (Wronskian)行列式计算的函数等。
另外，SymPy还支持特征值和特征向量的计算、矩阵的转置以及矩阵与行列式求解.还支持因式分解算法等.在计算中,还有零空间(null space)计算,行列式、代数余子式展开工具以及伴随矩阵。

5.物理学模块(sympy.physics)

SymPy有一个模块可以解决物理学问题.它支持力学功能,包括经典力学与量子力学以及高能物理学.它还支持一维空间与三维空间的泡利代数与量子谐振子,支持光学相关的功能.它还有一个独立的模块将物理单位集成到SymPy 里.用户可以选择相应的物理单位完成计算和单位转换.

6.统计学模块(sympy.stats)

SymPy的统计学模块支持数学计算中涉及的许多统计函数.除了常见的连续与离散随机分布函数,它还支持符号概率相关功能.SymPy库中的随机分布函数都支持随机数生成功能。

3.1.2 符号运算基础知识

使用Sympy库进行符号计算，首先要建立符号变量和符号表达式。符号变量是构成符号表达式的基本元素。可以通过库中的symbols()函数创建，如下：

>>>from sympy import *

>>>x=symbols('x') # x是符号变量的名称，'x'是符号变量的值，名称和值可以不相同

>>>y,z=symbols('y z') #构建多个符号变量时，赋值中间要用空格分隔

#构建多个符号变量还有一种方法，将m0:3传入符号函数，生成如m0,m1,m2的符号序列

在符号运算中，使用evalf()或n()方法来获得任何对象的浮点近似值，默认精度为小数点后15位。

例3.1 符号创建，类型转换以及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())  #以浮点数显示计算结果

Sympy库有很多函数可以用于处理有理数，可以对有理数进行简化，展开，合并等操作。together函数可以用于计算两个有理数的加法，apart函数可以用于处理有理数的除法。

例3.2 符号函数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)))

3.2 Scipy工具库简介

SciPy是对NumPy的功能扩展,它提供了许多高级数学函数,例如,微分、积分、微分方程、优化算法、数值分析、高级统计函数、方程求解等.SciPy是在NumPy数组框架的基础上实现的,它对 NumPy数组和基本的数组运算进行扩展，满足科学家和工程师解决问题时需要用到的大部分数学计算功能.
SciPy支持的功能包括文件处理、积分、数值分析、优化方法、统计学、信号与图像处理、聚类分析和空间分析等.下面简要介绍部分功能模块.

1. 微积分模块 (scipy.integrate)

微积分模块支持数值积分和微分方程数值解的功能。

注：以下一般不给出函数用法的详细说明,读者可以查看相关函数的帮助信息,例如,看函数romberg的帮助信息，使用命令:
>>>from scipy.integrate import romberg

>>> help(romberg)

(1) 给定函数的数值积分

quad:一重数值积分.

dblquad:二重数值积分.

tplquad:三重数值积分.

nquad:通用N重积分.

fixed-quad:使用固定阶高斯求积公式求数值积分.

quadrature:使用固定误差限的高斯求积公式求数值积分.

romberg:求函数的 Romberg数值积分.

(2) 给定离散点的数值积分

cumtrapz:用梯形法求数值积分.

simps:用辛普森法求数值积分.

romb:用 Romberg积分法求自变量均匀间隔离散点的数值积分

(3）微分方程的数值解

odeint:使用FORTRAN 库中方法求微分方程组的数值解.

ode:求一般微分方程组的数值解.

complex_ode:求复微分方程组的数值解.

2. 线性代数模块 (scipy.linalg)

与numpy.linalg 相比, scipy.linalg 函数有更高级的特征。

3. 优化模块 (scipy.optimize)

SciPy的优化模块提供了解决单变量和多变量的目标函数最小值问题的功能.它通过大量的算法解决最小化问题.优化模块支持线性回归、搜索函数的最大值与最小值、方程求根、线性规划、拟合等功能.

4. 插值模块 (scipy.interpolate)

插值模块支持一维和多维插值,例如,泰勒(Taylor)多项式插值,一维和多维样条插值.

5. 统计学模块 (scipy.stats)

统计模块提供了各种随机变量的分布、统计量的计算、分布拟合、参数检验等功能.

6. 傅里叶变换模块 (scipy.fftpack)

离散傅里叶变换和离散傅里叶逆变换可以分别用ft和ifft函数来计算.

7. 信号处理模块 (scipy.signal)

信号处理模块包含一系列滤波函数、滤波器设计函数,以及对一维和二维数据进行B-样条插值的函数.这个模块包含的函数可以进行以下操作:卷积、B-样条、滤波、滤波器设计、MATLAB式的IIR滤波器设计、连续时间的线性系统、离散时间的线性系统、线性时不变系统、信号波形、窗函数、小波分析和光谱分析等.

8. 多维图像处理模块 (scipy.cluster)

通常图像处理可以看作对二维数组的操作.这个模块提供了图像处理的各种函数,例如,图像几何变换、图像滤波等.

9. 空间分析模块 (scipy.spatial)

空间分析是一系列用于分析空间数据的算法．空间数据是指和地理空间或垂直空间相关的数据对象.这种数据包括点、线、多边形、其他几何和地理特征信息.
该模块支持 Delaunay三角剖分、Voronoi图、N维凸包等功能,支持KD树(scipy.spatial.kdtree)实现快速近邻查找算法,还可以对初始向量集合进行距离矩阵的计算.

10. 聚类模块 (scipy.cluster)

聚类是将一个大的集合分成多个组的过程. SciPy 聚类模块包括两个子模块:向量量化(vector quantization，VQ)(scipy.cluster.vq)和层次聚类(scipy.cluster.hierarchy).VQ模块支持 K均值聚类和向量量化,层次聚类模块支持分层聚类和聚合聚类.

11. 文件输入/输出模块（scipy.io）

SciPy可以使用MATLAB的“.mat”文件格式读取和写入数据，函数为load-mat 和 savemat.如果要加载数据，则可以使用如下语法:
>>>import scipy.io
>>>data=scipy.io.loadmat( 'datafile.mat ')
返回值data为一个字典，该字典包含了与“.mat”文件中保存的变量名相对应的键,对应值为 NumPy数组格式.
保存数据到“.mat”文件涉及创建一个包含要保存的所有变量的字典(变量名和值),函数为savemat,保存数组x和y 的代码如下
>>>data-0;data['x']-x; dataL'y']=y
>>>scipy.io.savemat( 'datafile.mat ' ,data)

3.3 用Sympy做符号函数画图

（与上一章中numpy利用matplotlib作图相比，sympy更方便）

1. 二维曲线画图

plot的基本使用格式为

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

多重绘制的使用格式为

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

或者

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

例3.3 在同一图形界面上画出y1=2sinx,x属于[-6,6]；y2=cos(x+pi/4),x属于[-5,5]。

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

2. 三维曲面画图

例3.4 画出三维曲面 $z=sin(\sqrt{x^2+y^2})$ 的图形

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$')

3. 隐函数画图

例3.5 绘制隐函数 $(x-1)^2+(y-2)^3-4=0$ 的图形

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=True)
pt(Eq((x-1)**2+(y-2)**3,4),(x,-6,6),(y,-2,4),xlabel='$x$',ylabel='$y$')
#转化成(x-1)^2+(y-2)^3=4的形式，Eq为隐函数画图标志

也可以使用匿名函数（lamada函数）：

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)

3.4 高等数学问题的符号解

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

1. 求极限

例3.6 验证 $\lim_{x\rightarrow 0}\tfrac{sinx}{x}=1,\lim_{x\rightarrow +\infty }(1+\tfrac{1}{x})^x=e$

from sympy import *
x=symbols('x')#定义符号变量
print(limit(sin(x)/x,x,0))#limit->取极限函数
print(limit(pow(1+1/x,x),x,oo))  #这里是两个小”o”，表示正无穷

2. 求导数

例3.7 已知： $z=sinx+x^2e^y$ ，求 $\frac{\partial^2 z}{\partial x^2},\frac{\partial z}{\partial y}$

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