Sympy 基础
Sympy是python符号计算库，有着符号求解的特点。
这篇博客内容是利用一个求解方程的例子，简单介绍Sympy符号计算的基本过程。
首先，我们要解决的问题是车流微分方程：
F为流量，
在某一时刻流量有最大值
F与
从物理意义上讲，当密度最大时拥堵最严重，车都不走了，流量为0.
什么时候流量最大呢？就是F(
那么总结起来就是：
求解方程组
在引入Sympy的时候，初始化LaTeX语法，这样就能用LaTeX语法写符号了。
import sympy sympy.init_printing()
u_max, u_star, rho_max, rho_star, A, B = sympy.symbols('u_max u_star rho_max rho_star A B')下面用这些变量表示上面三个方程：
sympy.Eq()函数将两个表达式分别放在等号两边，组成方程：
eq1 = sympy.Eq( 0, u_max*rho_max*(1 - A*rho_max-B*rho_max**2) ) eq2 = sympy.Eq( 0, u_max*(1 - 2*A*rho_star-3*B*rho_star**2) ) eq3 = sympy.Eq( u_star, u_max*(1 - A*rho_star - B*rho_star**2) )这时候输入eq1，得到：
eq2：
eq3：
假设u_star已知，那么三个未知量三个变量，方程组可解.
>>>eq2-3*eq3
sympy直接将这个字符串打印出来了，而没有转化为方程。
这时候，应该利用sympy.Eq()函数创建新的方程：
eq4 = sympy.Eq(eq2.lhs - 3*eq3.lhs, eq2.rhs - 3*eq3.rhs)
eq4：
eq2.lhs表示等号左边的代数式（left hand side）
eq2.rhs表示等号右边的代数式（right hand side）
这时候，方程eq4还没有化简。
eq4.simplify() 
注意这时候eq4本身并没有化简。
eq4依然为：
我们再看eq4.expand()函数：
eq4.expand()
expand函数将因式全部乘了起来，并且进行了化简,但是此时eq4依然没有改变。
这时候我们可以进行求解方程。
首先，从方程4求解处rho_star：
In[17]:rho_sol = sympy.solve(eq4,rho_star)
In[18]:rho_sol
Out[18]:
注意：sympy.solve()返回的是一个列表。所以我们不如采用如下代码：
rho_sol = sympy.solve(eq4,rho_star)[0] rho_sol
B_sol = sympy.solve(eq1,B)[0] B_sol
将求解的rho_sol B_sol代入方程2中：
quadA = eq2.subs([(rho_star, rho_sol), (B,B_sol)]) quadA
quadA.simplify()
A_sol = sympy.solve(quadA, A)A_sol[0]
A_sol[1]
假设
aval = A_sol[0].evalf(subs={u_star: 0.7, u_max:1.0, rho_max:10.0} ) aval
舍去负值。
bval = B_sol.evalf(subs={rho_max:10.0, A:aval} ) bval
同样，也求解得到B。
求微分：
x, nu, t = sympy.symbols('x nu t') phi = sympy.exp(-(x-4*t)**2/(4*nu*(t+1))) + sympy.exp(-(x-4*t-2*numpy.pi)**2/(4*nu*(t+1)))
phiprime = phi.diff(x)
lambdify：
fromsympy.utilities.lambdifyimport lambdify
u = -2*nu*(phiprime/phi)+4
ufunc = lambdify((t,x, nu),u)
print("The value of u at t=1, x=4, nu=3 is {}.".format(ufunc(1,4,3)))
结果为：
The value of u at t=1, x=4, nu=3 is 3.49170664206.
关闭Latex语法
sympy.init_printing(use_latex=False)
京公网安备 11010502036488号