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