返回顶部

Sympy 库的使用

参考资料:

官方文档,天下第一:https://docs.sympy.org/latest/index.html  

 

Introduction:

What is Sym?

符号计算以符号的方式处理数学对象的计算。这意味着数学对象被精确地表示,而不是近似地表示。例如,1/3 如果数值表示,那么它就是0.3333333.... 这是不精确的,如果符号表示就是1/3 ,它是精确的。

Why Sym?

首先,SymPy是完全免费的。它是开源的,并且是在自由的BSD许可下授权的,所以你可以修改源代码,如果你想的话甚至可以出售它。这与流行的商业系统如Maple或Mathematica形成鲜明对比,后者的许可证费用高达数百美元。

其次,SymPy使用Python。大多数计算机代数系统都发明了自己的语言。SymPy 则不是。SymPy完全是用Python编写的,完全是用Python执行的。这意味着,如果您已经了解Python,那么开始使用SymPy会容易得多,因为您已经知道语法(如果您不了解Python,则非常容易学习)。我们已经知道Python是一种设计良好、经过战斗考验的语言。SymPy开发人员对自己编写数学软件的能力充满信心,但编程语言设计则完全不同。通过重用现有的语言,我们能够专注于那些重要的事情:数学。 

代码

from sympy import *

x,t = symbols('x t')
init_printing(use_unicode=True)
ret = diff(sin(x) * exp(x), x)  # diff() 函数是求导函数
print(ret)  # exp(x)*sin(x) + exp(x)*cos(x)

ret2 = integrate(exp(x) * sin(x) + exp(x) * cos(x), x)  # integrate() 计算积分
print(ret2)  # sin(x)*exp(x)

ret3 = integrate(sin(x ** 2), (x, -oo, oo))  # 计算从-oo 到 oo 的积分
print(ret3)  # sqrt(2)*sqrt(pi)/2

ret4 = limit((sin(x) / x), x, 0)  # 计算x 趋近于0 的极限值
print(ret4)  # 1

ret5 = solve(x ** 2 - 2, x)  # 求解方程 x**2 -2 =0 x 的值
print(ret5) # [-sqrt(2), sqrt(2)]

# Solve the differential equation y‘’−y = e^t. 求解微分方程 y(t)
y = Function('y')
ret6 = dsolve(Eq(y(t).diff(t,t) - y(t),exp(t)),y(t))
print(ret6)  # Eq(y(t), C2*exp(-t) + (C1 + t/2)*exp(t))


# 求矩阵的特征值
ret7 = Matrix([[1,2],[2,2]]).eigenvals()
print(ret7) # {3/2 - sqrt(17)/2: 1, 3/2 + sqrt(17)/2: 1}

# 使用Latex 格式打印数学式子
ret8 = latex(sqrt(2)*sqrt(pi)/2)
print(ret8)  # \frac{\sqrt{2} \sqrt{\pi}}{2}
View Code
from sympy import *

x, y, z = symbols("x y z")

# 判断两个表达式是否相同
a = (x + 1) ** 2
b = x ** 2 + 2 * x + 1
print(simplify(a - b) == 0)

# .equals() 是通过在随机取点上的数值得到是否相等
print(a.equals(b))

# sympy 中幂用的是 ** 因为Python中^ 表示的是异或
# sympy 中分数表达要用 Rational(),因为Python中的/ 表示的是除法
print(x ** 2)
print(x + Rational(1, 2))
View Code

 

Gotchas

首先,我们应该弄清楚关于sypy的一些事情。SymPy不过是一个Python库,比如NumPy、Django,甚至是Python标准库sys或re中的模块。这意味着SymPy不向Python语言添加任何内容。Python语言固有的局限性也是SymPy固有的。这也意味着sypy会尽可能地使用Python习惯用法,这使得那些已经熟悉Python编程的人可以轻松地使用sypy编程。作为一个简单的例子,SymPy使用Python语法来构建表达式。隐式乘法(如3x或3x)在Python中是不允许的,因此在SymPy中也不允许。要与x相乘,必须输入3*。

Symbols

from sympy import * 


# 我们使用 symbols 表示符号
x = symbols("x")
print(x+1)


print(y+1) # 报错 
View Code

 

Basic Operations

Here we discuss some of the most basic operations needed for expression manipulation in SymPy. Some more advanced operations will be discussed later in the advanced expression manipulation section.

Substitution

One of the most common things you might want to do with a mathematical expression is substitution. Substitution replaces all instances of something in an expression with something else. It is done using the subs method. For example

from sympy import * 
x,y,z = symbols("x,y,z")
# 替换 
expr = cos(x)+1
print(expr.subs(x,y))  # cos(y) + 1

 

Substitution is usually done for one of two reasons:

1,Evaluating an expression at a point. For example, if our expression is cos(x) 1 and we want to evaluate it at the point 0, so that we get cos(0) 1, which is 2.

print(expr.subs(x, 0)) # 2 

2, Replacing a subexpression with another subexpression. 

from sympy import * 
x,y,z = symbols("x,y,z")
expr = x**y 
print(expr.subs(y,x**x))

 

There are two important things to note about subs.

First, it returns a new expression. SymPy objects are immutable. That means that subs does not modify it in-place. 

Second, to perform multiple substitutions at once, pass a list of (old, new) pairs to subs.

from sympy import * 
x,y,z = symbols("x,y,z")
expr = x**3 + 4*x*y -z 
ret = expr.subs([(x,2),(y,4),(z,0)])
print(ret)

'''
output:
    40
'''

 

 

 

subs 的用途:

from sympy import * 
x,y,z = symbols("x,y,z")
expr = sin(2*x)+cos(2*x)
# 我们的需求是 想要 2sin(2x)+cos(2x)
print(expand_trig(expr))  # 都会展开
print(expr.subs(sin(2*x),2*sin(x)*cos(x))) # 这时替换就显得容易一些

'''
output:
2*sin(x)*cos(x) + 2*cos(x)**2 - 1
2*sin(x)*cos(x) + cos(2*x)
'''

Converting Strings to SymPy Expressions

The sympify function (that’s sympify, not to be confused with simplify) can be used to convert strings into SymPy expressions.

from sympy import * 
x,y,z = symbols("x,y,z")

str_expr = "x*sin(x)+3*x-1/2"
expr = sympify(str_expr) # 底层用的是 eval()  
print(expr)
View Code

 

evalf

To evaluate a numerical expression into a floating point number, use evalf.

from sympy import * 
x,y,z = symbols("x,y,z")

expr = sqrt(8)
ret = expr.evalf() # 默认n=15,14位 浮点数,还有一个点
print(ret)# 2.82842712474619

ret = pi.evalf(n=100) 
print(ret)#3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
View Code

 

To numerically evaluate an expression with a Symbol at a point, we might use subs followed by evalf, but it is more efficient and numerically stable to pass the substitution to evalf using the subs flag, which takes a dictionary of Symbol: point pairs.

from sympy import * 
x,y,z = symbols("x,y,z")

expr = cos(2*x) +sin(y)
ret = expr.evalf(subs={x:1,y:2})  # 使用 evalf() 函数中的 subs 参数效率更高
print(ret) # 0.493150590278539

有时会有舍入误差小于所需精度,这些误差在计算表达式后仍然存在。用户可以通过将chop标志设置为True来删除这些数字。

from sympy import * 
x,y,z = symbols("x,y,z")

one = cos(1)**2 + sin(1)**2
print((one - 1).evalf()) #-0.e-124
print((one - 1).evalf(chop=True)) #0
View Code

 

lambdify

 

subs and evalf are good if you want to do simple evaluation, but if you intend to evaluate an expression at many points, there are more efficient ways. For example, if you wanted to evaluate an expression at a thousand points, using SymPy would be far slower than it needs to be, especially if you only care about machine precision. Instead, you should use libraries like NumPy and SciPy.

The easiest way to convert a SymPy expression to an expression that can be numerically evaluated is to use the lambdify function. lambdify acts like a lambda function, except it converts the SymPy names to the names of the given numerical library, usually NumPy. For example

import numpy
from sympy import * 

x = symbols("x")
a = numpy.arange(10)
expr = sin(x)

print(a)
f = lambdify(x,expr,"numpy")

ret = f(a)
print(ret)

b = numpy.arange(20);
ret = f(b)
print(ret)
View Code

The primary purpose of this function is to provide a bridge from SymPy expressions to numerical libraries such as NumPy, SciPy, NumExpr, mpmath, and tensorflow. In general, SymPy functions do not work with objects from other libraries, such as NumPy arrays, and functions from numeric libraries like NumPy or mpmath do not work on SymPy expressions.

 

其他例子:

import numpy
from sympy import * 

x = symbols("x")
expr = sin(x)+cos(x)
f = lambdify(x,expr,"math")
ret = f(0.1)
print(ret)

ret2 = expr.subs(x,0.1)
print(ret2)
View Code

 

Printing

As we have already seen, SymPy can pretty print its output using Unicode characters. This is a short introduction to the most common printing options available in SymPy.

Printers

There are several printers available in SymPy. The most common ones are

  • str

  • srepr

  • ASCII pretty printer

  • Unicode pretty printer

  • LaTeX

  • MathML

  • Dot

Setting up Pretty Printing

If all you want is the best pretty printing, use the init_printing() function. This will automatically enable the best printer available in your environment.

 

省略...

 

simplification

simplify

Now let’s jump in and do some interesting mathematics. One of the most useful features of a symbolic manipulation system is the ability to simplify mathematical expressions. SymPy has dozens of functions to perform various kinds of simplification. There is also one general function called simplify() that attempts to apply all of these functions in an intelligent way to arrive at the simplest form of an expression. Here are some examples

from sympy import * 
init_printing(use_unicode=True)  
x = symbols("x")

ret = simplify(sin(x)**2+cos(x)**2)
print(ret) # 1

ret = simplify(sin(2*x) + cos(2*x))
print(ret)  # sqrt(2)*sin(2*x + pi/4)
View Code

 

But simplify() has a pitfall. It just applies all the major simplification operations in SymPy, and uses heuristics to determine the simplest result. But “simplest” is not a well-defined term. For example

from sympy import * 
init_printing(use_unicode=True)  
x = symbols("x")

ret = simplify(x**2+2*x+1)
print(ret) # x**2+2*x+1
View Code

We did not get what we want. There is a function to perform this simplification, called factor(), which will be discussed below.

Another pitfall to simplify() is that it can be unnecessarily slow, since it tries many kinds of simplifications before picking the best one. If you already know exactly what kind of simplification you are after, it is better to apply the specific simplification function(s) that apply those simplifications.

Applying specific simplification functions instead of simplify() also has the advantage that specific functions have certain guarantees about the form of their output. These will be discussed with each function below. For example, factor(), when called on a polynomial with rational coefficients, is guaranteed to factor the polynomial into irreducible factors. simplify() has no guarantees. It is entirely heuristical, and, as we saw above, it may even miss a possible type of simplification that SymPy is capable of doing.

factor 因式分解

factor()接受一个多项式,并将其分解为有理数上的不可约因子。

from sympy import * 
init_printing(use_unicode=True)  
x = symbols("x")

ret = simplify(x**2+2*x+1) # 不能简化
print(ret)

ret = factor(x**2+2*x+1) # 可以简化
print(ret)
View Code

For polynomials, factor() is the opposite of expand()

factor()在有理数上使用一个完整的多元因子分解算法,这意味着factor()返回的每个因子都保证是不可约的。

 

expand 展开

from sympy import * 
init_printing(use_unicode=True)  
x = symbols("x")

ret = expand((x+1)**2) 
print(ret)  # x**2 + 2*x + 1
View Code

 

collect 合并同类项

from sympy import * 
init_printing(use_unicode=True)  
x,y,z = symbols("x,y,z")

expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3
ret = collect(expr,x)
print(ret)  # x**3 + x**2*(2 - z) + x*(y + 1) - 3
View Code

 

cancel 通分

from sympy import * 
init_printing(use_unicode=True)  
x,y,z = symbols("x,y,z")

expr = 1/x + (3*x/2 - 2)/(x - 4)
ret = cancel(expr)
print(ret) #(3*x**2 - 2*x - 8)/(2*x**2 - 8*x)
View Code

 

apart   部分   分数分解

from sympy import * 
init_printing(use_unicode=True)  
x,y,z = symbols("x,y,z")

expr = (4*x**3 + 21*x**2 + 10*x + 12)/(x**4 + 5*x**3 + 5*x**2 + 4*x)
ret = apart(expr)
print(ret)  # (2*x - 1)/(x**2 + x + 1) - 1/(x + 4) + 3/x
View Code

 

posted @ 2020-07-16 22:52  Zcb0812  阅读(1178)  评论(0编辑  收藏  举报