python sympy库

sympy库是python的符号运算库，是电脑辅助简单数学函数计算的好工具。本文简单记录了一下有关sympy的方法。建议使用jupyter notebook，这样输出的函数很好看。

---

sympy基础

安装

pip install sympy

自变量（Symbols）

通过symbols方法将字符串声明为数学函数的变量，这些变量往往是函数的自变量。

import sympy
# 声明单个变量
x = sympy.symbols('x')

# 声明多个变量(以下方法皆可)
x,y = sympy.symbols(['x','y'])
x,y = sympy.symbols("x,y")
x,y = sympy.symbols("x y")

函数表达式（Expr）

函数表达式通过变量的运算构造具体函数，或者通过Function函数构造抽象函数。

f = sympy.sqrt(3*x*y) +x*sympy.sin(y) +y**2 +x**3  # 具体函数
u = sympy.Function('u')  # 抽象函数

变量替换与赋值

expr.subs()可以实现变量替换，替换成数字实现赋值。

g1 = f.subs(x,y) # 将f表达式中的x换成y，并将替换的结果，赋给g
g2 = f.subs({x:2*x, y:2*y}) #多次替换,字典
g3 = f.subs({x:1,y:2})

精确求值

expr.evalf(n)可以求一个表达式的保留n位有效数字的精确值

g3 = f.subs({x:1,y:2})
print(g.evalf(4))  # 8.359

微分

sympy可以实现自动求微分，方法如下

h1 = sympy.diff(f,x)   # f对x求微分
h1 = f.diff(x)   # 同上
h2 = sympy.diff(f,x,2,y,1)  # f对x求2次微分，对y求1次微分

积分

sympy可以实现自动求不定积分和定积分，区别在于是否传入积分上下限

l1 = sympy.integrate(f,x)    # 不定积分
l2 = sympy.integrate(f,(x,1,3))  #定积分

Taylor展开

一元展开

sympy可以实现泰勒展开，具体函数抽象函数都可以。但是不能对多元函数同时泰勒展开。

taylor1 = sympy.series(f,x,0,4)  # f对x在0处泰勒展开到4阶
taylor2 = sympy.series(f,x,0,4).removeO  # f对x在0处泰勒展开到4阶，去除皮亚诺余项
taylor3 = sympy.series(u(x),x,0,4)  # 抽象函数u对x在0处泰勒展开到4阶

多元展开

多元泰勒展开可以使用以下方法
参考链接：https://www.thinbug.com/q/23803320

def Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree):
    """
    Mathematical formulation reference:
    https://math.libretexts.org/Bookshelves/Calculus/Supplemental_Modules_(Calculus)/Multivariable_Calculus/3%3A_Topics_in_Partial_Derivatives/Taylor__Polynomials_of_Functions_of_Two_Variables
    :param function_expression: Sympy expression of the function
    :param variable_list: list. All variables to be approximated (to be "Taylorized")
    :param evaluation_point: list. Coordinates, where the function will be expressed
    :param degree: int. Total degree of the Taylor polynomial
    :return: Returns a Sympy expression of the Taylor series up to a given degree, of a given multivariate expression, approximated as a multivariate polynomial evaluated at the evaluation_point
    """
    from sympy import factorial, Matrix, prod
    import itertools

    n_var = len(variable_list)
    point_coordinates = [(i, j) for i, j in (zip(variable_list, evaluation_point))]  # list of tuples with variables and their evaluation_point coordinates, to later perform substitution

    deriv_orders = list(itertools.product(range(degree + 1), repeat=n_var))  # list with exponentials of the partial derivatives
    deriv_orders = [deriv_orders[i] for i in range(len(deriv_orders)) if sum(deriv_orders[i]) <= degree]  # Discarding some higher-order terms
    n_terms = len(deriv_orders)
    deriv_orders_as_input = [list(sum(list(zip(variable_list, deriv_orders[i])), ())) for i in range(n_terms)]  # Individual degree of each partial derivative, of each term

    polynomial = 0
    for i in range(n_terms):
        partial_derivatives_at_point = function_expression.diff(*deriv_orders_as_input[i]).subs(point_coordinates)  # e.g. df/(dx*dy**2)
        denominator = prod([factorial(j) for j in deriv_orders[i]])  # e.g. (1! * 2!)
        distances_powered = prod([(Matrix(variable_list) - Matrix(evaluation_point))[j] ** deriv_orders[i][j] for j in range(n_var)])  # e.g. (x-x0)*(y-y0)**2
        polynomial += partial_derivatives_at_point / denominator * distances_powered
    return polynomial

Taylor_polynomial_sympy(u(x,y),[x,y],[0,0],2)  #函数u对(x,y)在(0,0)处泰勒展开

查看展开项

taylor1.coeff(x)  # 查看taylor1中x项(x-x0项)的系数。 

极限

sympy可以实现求极限，注意极限方向

lim1 = sympy.limit(f, x, sympy.oo)  # 趋于无穷
lim2 = sympy.limit(f, x, 0)  # 趋于0，默认值 dir='+'，也就是趋于+0
lim3 = sympy.limit(f, x, 0, dir="-")  # 趋于0，默认值调整为dir='-'，也就是趋于-0

解方程

sympy可以实现解方程，方法是令Expr=0，所以在解方程时，要先构造一个等于0的左端项。返回结果是一个列表，每一项是一个解。如果是方程组，解列表每一项是一个元组，元组对应位置是对应自变量的值。

func = f-3
sympy.solve(func,x)  # 返回f=3时x的值
sympy.solve([x**2+y**2-1, x+y-1],[x,y])  #x^2+y^2=1;x+y=1

画图

sympy可以实现绘制函数图像，要先安装matplotlib库。

pip install matplotlib

sympy可以绘制一元函数、二元函数、隐函数的图像，方法如下

from sympy.plotting import plot,plot3d,plot_implicit 
from sympy import Eq  # 构造隐函数时使用

# 绘制一元函数图像
plot(x**2,x)  # 绘制以x为自变量的函数图像
plot(x**2,(x,-3,4))  # 绘制以x为自变量的函数图像，绘制区间限制在[-3,4]上
plot(x,x**2,x**3,(x,0,4))  # 绘制多个图像 

# 绘制二元函数图像，呈现3D图
plot3d(x**2*y,(x,-10,10),(y,-10,10))

# 绘制隐函数图像
plot_implicit(Eq(x**2+y**4,1),(x,-2,2),(y,-2,2))