摘要：在学习与科研中，经常会遇到一些数学运算问题，使用计算机完成运算具有速度快和准确性高的优势。python的numpy包具有强大的科学运算功能，且具有其他许多主流科学计算语言不具备的免费、开源、轻量级和灵活的特点。本文使用python语言的numpy库，解决数学运算问题中的线性方程组问题、积分问题、微分问题及矩阵化简问题，结果准确快捷，具有一定的借鉴意义。
1.sympy库简介
sympy一个用于符号型数学计算（symbolic mathematics）的python库。它旨在成为一个功能齐全的计算机代数系统（computer algebra system，cas），同时保持代码简洁、易于理解和扩展。sympy完全是用python写的，并不需要外部的库。
本文选择sympy库的原因在于：
- 免费：该库基于bsd开源许可，免费且开源；
- 基于python：该库完全是用python写就，并以python作为该库操作语言；
- 轻量级：为了使sympy简单易用，该库仅基于mpmath库（一个纯python库，用于浮点运算）；
- 灵活性：除了用作交互工具，还可插入其他应用或软件拓展功能中。
具体说来，如果x与y未曾赋值，那么下列语句就会报错
|
1
2
|
#测试语句
print(x+y)
|
而符号运算则不同，符号运算多用于公式推导，不需要赋值，此时使用sympy进行符号运算便具有方便快捷的优势，如下述语句便不会报错。
|
1
2
3
4
|
#测试语句
x=symbol('x')
y=symbol('y')
print(x+y)
|
2 sympy库解决数学运算问题实现
2.1 求解线性方程组
解方程的功能主要是使用sympy中solve函数实现。以式(1)为例，求解过程如下：
（1） 符号表示
sympy库中使用symbol函数定义符号变量，
|
1
2
3
4
5
6
|
from sympy import *
x=symbol('x')
y=symbol('y')
#或者用如下语句
x,y=symbol('x y')#第二个用空格隔开
|
（2）方程表示
使用代码表示数学符号与手写体的数学运算符号存在一定的差异，下面列举常用的运算符：
- 加号
- 加号 +
- 减号 -
- 除号 /
- 乘号 *
- 指数 **
- 对数 log()
- e的指数次幂 exp()
对于长的表达式，如果不确定运算符的优先级，可以加入小括号提升其优先级。由于需要将表达式都转化成右端等于0,这里把常数3和7移到等式左边。题目中表达式可表示为：
|
1
2
|
2*x-y-3=0
3*x+y-7=0
|
(3)使用solve函数解方程
在使用solve函数解方程之前，我们先来看一下solve函数的定义。solve函数的第一个参数是要解的方程，要求右端等于0，第二个参数是未知数。
对于式(1)的求解，代码如下：
|
1
2
3
4
|
from sympy import *
x = symbol('x')
y = symbol('y')
print(solve([2*x-y-3,3*x+y-7],[x,y]))
|
2.2 求解微积分问题
2.2.1 求解极限问题
在2.1中通过一个简单的二元一次方程组求解熟悉了该库求解数学问题的基本过程，下面本文通过示例，讲解使用sympy库求解微积分的过程。
求解式（2）所示的极限问题，需要用到limit函数求极限。
(1)符号及方程表示
引入sympy库并定义n为符号变量与2.1中一致。
|
1
2
3
|
from sympy import *
n = symbol('n')
s = ((n+3)/(n+2))**n
|
(2)利用limit函数求极限
首先我们介绍limit函数的调用格式：limit(e, z, z0, dir='+'),e为任意表达式，表示求取e(z)在点z0处的极限，dir='+'表示取右极限，die='-'则表示取左极限。则上式的求解代码可表示如下：
|
1
2
3
|
from sympy import *
n = symbol('n')
s = ((n+3)/(n+2))**n
|
print(limit(s,n,oo)) #无穷的表示方法是两个小写的字母o
2.2.2 求解定积分
(1)符号表示
|
1
2
3
|
from sympy import *
t = symbol(t)
x = symbol(x)
|
(2)方程表示
|
1
2
|
m = integrate(sin(t)/(pi-t),(t,0,x))
n = integrate(m,(x,0,pi))
|
完整代码如下：
|
1
2
3
4
5
6
|
from sympy import *
t = symbol(t)
x = symbol(x)
m = integrate(sin(t)/(pi-t),(t,0,x))
n = integrate(m,(x,0,pi))
print(n)
|
2.2.3 求解微分问题
如求取的通解
（1）符号表示
这里与之前不同的是增加了函数的表示(用f(x)表示y)，即例题中的y还有微分表示
|
1
2
3
|
from sympy import *
f = function('f')
x = symbol('x')
|
y'的表示方法由以下代码组成
|
1
|
diff(f(x),x)
|
这里对diff函数稍作介绍：
上面是求一阶导的方法，求解高阶导的方法如下所示：
|
1
2
3
4
5
6
7
8
9
10
|
>>> diff(x**3,x)
3*x**2
>>> diff(x**3,x,1)
3*x**2
>>> diff(x**3,x,2)
6*x
>>> diff(x**3,x,3)
6
>>> diff(x**3,x,4)
0
|
即改变第三个参数即可。
下面继续我们的解题过程。
|
1
2
3
4
5
6
7
8
9
10
11
12
|
#左端
diff(f(x),x)
#看一下
print(diff(f(x),x))
#result
#d
#--(f(x))
#dx
#右端
2*f(x)*x
|
用dsolve函数解微分方程
dsolve函数是用来解决微分方程（differential equation）的函数。
函数的一个用法为：
|
1
|
dsolve(eq, f(x))
|
第一个参数为微分方程（要先将等式移项为右端为0的形式)。第二个参数为要解的函数(在微分方程中)
举个例子:
|
1
2
3
4
5
6
7
8
9
10
|
>>> from sympy import *
>>> f = function('f')
>>> x = symbol('x')
>>> pprint(2*x-diff(f(x),x))
d
2*x - --(f(x))
dx
>>> dsolve(2*x - diff(f(x),x), f(x))
#result
#eq(f(x), c1 + x**2)
|
这样，我们可以将我们要解的题目，用以下代码表示。
|
1
|
dsolve(diff(f(x),x) - 2*f(x)*x, f(x))
|
结果为：
eq(f(x), c1*exp(x**2))
#即f(x) = c1*exp(x**2)
对比答案可以发现正确。
完整代码:
|
1
2
3
4
|
from sympy import *
f = function('f')
x = symbol('x')
print(dsolve(diff(f(x),x)-2*f(x)*x,f(x))
|
2.2.4 矩阵化简
平时线性代数问题中我们会遇到化简问题，虽然不算难，但着实麻烦。而且，出一点错就会导致
结果出错。不过好运的是sympy提供了相关的支持。
例题：
符号表示与矩阵表示
|
1
2
3
4
5
6
|
from sympy import *
x1,x2,x3 = symbols('x1 x2 x3')
a11,a12,a13,a22,a23,a33 = symbols('a11 a12 a13 a22 a23 a33')
m = matrix([[x1,x2,x3]])
n = matrix([[a11,a12,a13],[a12,a22,a23],[a13,a23,a33]])
v = matrix([[x1],[x2],[x3]])
|
注意m的表示，需要有两个中括号
化简实现
|
1
|
print(m*n*v)
|
得到的是：
matrix([[x1*(a11*x1 + a12*x2 + a13*x3) + x2*(a12*x1 + a22*x2 + a23*x3) + x3*(a13*x1 + a23*x2 + a33*x3)]])
使用
|
1
2
|
f = m * n * v
print f[0]
|
可以进一步得到化简后的式子
也许你要问我要化简后在计算怎么办？下面我就举个例子。
如果上式中x1，x2，x3均等于1，则可这样代入。
|
1
2
3
4
5
6
7
8
|
from sympy import *
x1,x2,x3 = symbols('x1 x2 x3')
a11,a12,a13,a22,a23,a33 = symbols('a11 a12 a13 a22 a23 a33')
m = matrix([[x1, x2, x3]])
n = matrix([[a11, a12, a13], [a12, a22, a23], [a13, a23, a33]])
v = matrix([[x1], [x2], [x3]])
f = m * n * v
print f[0].subs({x1:1, x2:1, x3:1})
|
可得
a11 + 2*a12 + 2*a13 + a22 + 2*a23 + a33
以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持服务器之家。
原文链接：https://blog.csdn.net/whcxytj/article/details/75270873