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}

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

 

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) # 報錯 

 

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

 

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

 

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

 

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)

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)

 

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)

 

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

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)

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

 

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

 

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)

 

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