python做微积分 | 第10期

      我们最开始使用python基本都是做数值计算的，而在处理数学问题时常常也需要进行一些符号运算，python能否胜任这种需求呢？当然是可以的，我们就以微积分为例，来见识一下python处理符号运算和数值运算的双重能力。需要安装一款python的第三方库sympy，可以直接使用pip3 install sympy来进行安装，也可以通过其他方式安装。本文无意于让大学生完成高等数学作业时投机取巧，不过作为验证自己计算是否正确的工具倒是不错。sympy更有意义的用途是用于科学计算模拟系统运行状况。

1. 变量定义

     数学问题中，我们常常需要定义一些自变量，然后使用一组映射作用在这些自变量上构成数学方程。sympy提供了数学符号相关的数据结构Symbol，在进行符号运算之前需要先定义一些自变量符号。单个符号直接使用sympy.Symbol()，多个符号就可以使用sympy.symbols()来进行定义。

import sympy as sp
# 单个变量定义
x = sp.Symbol('x')
print('x={}, type(x)={}'.format(x, type(x)))
# 多个变量定义
x, y, z = sp.symbols('x y z')
print('x={}, type(x)={}'.format(x, type(x)))

Output:

x=x, type(x)=<class 'sympy.core.symbol.Symbol'>
x=x, type(x)=<class 'sympy.core.symbol.Symbol'>

      从以上演示可以看出，sympy中定义了一个sympy.core.symbol.Symbol类，这个类就是sympy进行符号运算最基础的数据结构。

2. 微积分基础

    事实上，sympy可以完成高等数学中几乎所有的操作，当然也包括一些基础的运算。而计算微积分的时候我们常常需要具备一些数学运算的基础，比如极限、表达式展开和合并等等，这里顺便演示使用sympy一下。

x = sp.Symbol('x')
# 极限
print(sp.limit(x**5, x, 2))
print(sp.limit(x/sp.sin(x), x, 0))
# 展开
print(sp.expand((x+1)**2))
# 化简
print(sp.simplify(x**2 + 2*x + 1 + x+(1/2)*x**2 + 0.5))

Output:

32
1
x**2 + 2*x + 1
1.5*x**2 + 3*x + 1.5

      从以上演示看出，sympy在处理极限、表达式展开和化简等数学基础运算时非常直观，基本与人工手动计算的形式一致。如果你有一个很复杂的带有同类项的表达式，可以使用sympy.simplify函数试试。

3. 一元微积分

      我们先来看看一元微积分，即只有一个自变量的微积分。由于普通的python库基本都是进行数值计算的，因此在构建符号函数时不能直接使用，比如numpy；而应该使用sympy中提供的基础符号函数来构建复杂的符号函数，比如sympy.sin(), sympy.exp()等等。计算符号函数导数的方法也很简单，直接使用sympy.diff()，而计算积分可以直接调用sympy.integrate()函数。

x = sp.Symbol('x')
f = sp.diff(x**5)
print("f(x)={}, f'(x)={}".format("x^5", f))
print("f(x)={}, f'(x)={}".format("xe^x", sp.diff(x*sp.exp(x))))
print("f(x)={}, f'(x)={}".format("xln(x)", sp.diff(x*sp.ln(x))))
print("f(x)={}, f'(x)={}".format("xsin(x)cos(x)", sp.diff(x*sp.sin(x)*sp.cos(x))))
print("f(x)=x^5, f'(2)={}".format(f.evalf(subs={'x':2})))
print("f(x)=5x^4, 不定积分integrate(f)={}".format(sp.integrate(f, x)))
print("f(x)=5x^4, 定积分integrate(f, 0, 2)={}".format(sp.integrate(f, (x, 0, 2))))

Output:

f(x)=x^5, f'(x)=5*x**4
f(x)=xe^x, f'(x)=x*exp(x) + exp(x)
f(x)=xln(x), f'(x)=log(x) + 1
f(x)=xsin(x)cos(x), f'(x)=-x*sin(x)**2 + x*cos(x)**2 + sin(x)*cos(x)
f(x)=x^5, f'(2)=80.0000000000000
f(x)=5x^4, 不定积分integrate(f)=x**5
f(x)=5x^4, 定积分integrate(f, 0, 2)=32

      一元微积分的符号运算以及对应的数值运算如上图演示。大家可以发现，在计算不定积分时，省略了一个常数项。通常在未给定初值的情况下，我们也不关心这个常数项，与日常计算微积分的情形基本一致。

4. 多元微积分

     了解了一元微积分的计算，自然就更想见识一下多元微积分的计算，sympy在这方面也是非常地强大。所使用的接口仍然是diff和integrate这两个函数。通过变量与基础函数的拼接所形成的复杂函数在sumpy中是一个sympy.core.add.Add对象，使用这个对象直接调用diff函数并指定自变量就可以求解对应的偏导数。多重积分的计算与一元积分是相似的，只是需要传入想要进行积分运算的符号自变量。

x, y, z = sp.symbols('x y z')
# 多元微分
f = (x**2) + x*sp.ln(y) + sp.exp(z)
print("type(f)={}".format(type(f)))
print("f对x的偏导数:{}".format(f.diff(x)))
print("f对y的偏导数:{}".format(f.diff(y)))
print("f对z的偏导数:{}".format(f.diff(z)))
# 多元积分
f = x**2 + 2*y*z
print("三重不定积分:{}".format(sp.integrate(f, x, y, z)))  # 三重不定积分
# 三重定积分
print("三重定积分:{}".format(sp.integrate(f, (x,0,2), (y,0,2), (z,0,2))))

Output:

type(f)=<class 'sympy.core.add.Add'>
f对x的偏导数:2*x + log(y)
f对y的偏导数:x/y
f对z的偏导数:exp(z)
三重不定积分:x**3*y*z/3 + x*y**2*z**2/2
三重定积分:80/3

5. 微分方程

      基本的微积分运算都了解了，但是微积分中一类非常常见的问题还是有必要演示一下，即微分方程的求解。sympy可以求解普通的方程或者方程组是理所当然的，那么对于微分方程的求解是否也那么顺利呢？sympy对方程组的求解函数与matlab是一致的，这对于从matlab迁移到python的朋友来说是个很不错的消息。手动解算过微分方程的朋友都应该知道，微分方程的基础是建立在一个未知函数的导数之上的，那么这个未知的函数该如何表示这是个问题。在sympy中，使用sympy.Function()接口创建这个未知的函数，然后利用这个函数构建微分方程，使用dsolve求解。

x, y, z = sp.symbols('x y z')
# 求解一元方程组
print(sp.solve([x**2 + 2*x + 1]))
# 求解多元方程组
print(sp.solve([x+y+z-1, x+y-2, x+z-3]))
# 求解微分方程
f = sp.Function('f')
i_f = sp.dsolve(f(x).diff(x) + f(x)**2 + f(x), f(x), ics={f(0):1})
print("f'(x)+f(x)^2+f(x)=0, f(0)=1求解f(x):{},\r\ntype(f(x))={}".format(i_f, type(i_f)))
# 方程求值
print('f(x)={}'.format(i_f.subs([(x, 1)])))

Output:

[{x: -1}]
{x: 4, y: -2, z: -1}
f'(x)+f(x)^2+f(x)=0, f(0)=1求解f(x):Eq(f(x), -1/(2*(1/2 - exp(x)))),
type(f(x))=<class 'sympy.core.relational.Equality'>
f(x)=Eq(f(1), -1/(2*(1/2 - E)))

      到此，使用sympy模块进行微积分符号和数值计算的演示完毕。最开始我们处理这类问题的首选是matplab，等到掌握sympy之后，就没有再使用matlab的动力了。sympy不仅是免费的，而且轻量级，使用过程与matlab一样方便。

       欢迎对python和人工智能感兴趣的朋友加入微信群讨论和交流，本群也可提供pythn编程咨询：