sympy实战

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
%matplotlib inline
def logistic4(x, A, B, C, D):
    return (A-D)/(1+(x/C)**B)+D

def residuals(p, y, x): 
    A, B, C, D = p
    return y - logisctic4(x, A, B, C, D)

def peval(x, p):
    A, B, C, D = p
    return logistic4(x, A, B, C, D) 

A, B, C, D = .5, 2.5, 8, 7.3
x = np.linspace(0, 20, 20)
y_true = logistic4(x, A, B, C, D)

y_meas = y_true + 0.2 * np.random.randn(len(y_true))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

1

from sympy.parsing.sympy_parser import parse_expr
from sympy import plot_implicit
import numpy as np
1
2
3

temp3 = np.sqrt(3)
ezplot = lambda expr : plot_implicit(parse_expr(expr))
#3*sin(y)**2+cos(410*x)**2*cos(y)**2+2*1.732*cos(410*x)*sin(y)*cos(y)-4*(sin(410*x)+8)**2+sin(410*x)**2
string_math = "-3*cos(2*y)+0.5*(cos(2*x)+1)*(cos(2*y)+1)+2*{0}*cos(x)*sin(2*y)+3*cos(2*x)-1.5".format(temp3)

# string_math = "17 * x **2 - 16 *Abs(x)*y +17* y**2 - 256 "
print(string_math)
ezplot(string_math)  
1
2
3
4
5
6
7
8

-3*cos(2*y)+0.5*(cos(2*x)+1)*(cos(2*y)+1)+2*1.7320508075688772*cos(x)*sin(2*y)+3*cos(2*x)-1.5
1
2

<sympy.plotting.plot.Plot at 0x200e7cc99b0>
1
2

np.exp(1)+np.pi
1

5.8598744820488378
1
2

计算 2–√$\sqrt{2}$ 小数点后一百位。

from sympy import *
sqrt(2).evalf(100)   
1
2

1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573
1
2

len(str(sqrt(2).evalf(100)))
1

101
1
2

• 用有理数算术计算1/2 + 1/3 in rational arithmetic.

Rational(1,2)+Rational(1,3)
1

5/6
1
2

3.2.2 代数运算

• 计算(x+y)6$(x+y{)}^6$的展开。
• 化简三角表达式sin(x)/cos(x)$\sin (x)/\cos (x)$

#与其他计算机代数系统不同，在SymPy你需要显性声明符号变量:   不写会报错！！！！！！！！！！
x = Symbol('x')
y = Symbol('y')
expand((x+y)**6)
1
2
3
4

x**6 + 6*x**5*y + 15*x**4*y**2 + 20*x**3*y**3 + 15*x**2*y**4 + 6*x*y**5 + y**6
1
2

simplify(sin(x)/cos(x))
1

tan(x)
1
2

#化简是一个模糊的术语，更准确的词应该是：powsimp (指数化简)、 trigsimp (三角表达式)、logcombine、radsimp一起。
trigsimp(sin(x)/cos(x))
1
2

tan(x)
1
2

3.2.3 微积分

在SymPy中使用极限很简单，允许语法limit(function, variable, point), 因此要计算f(x)类似x→0$x\to 0$, 你应该使用limit(f, x, 0):

limit(sin(x)/x,x,0)
1

1
1
2

(求导)你可以使用diff(func, var)微分任何SymPy表达式。例如:

diff(tan(x),x)
1

tan(x)**2 + 1
1
2

你可以用下列方法检查是否正确:

limit((tan(x+y)-tan(x))/y,y,0)  #倒数的定义！！
1

tan(x)**2 + 1
1
2

用diff(func, var, n)方法来计算更高的导数:

diff(tan(x),x,2)
1

2*(tan(x)**2 + 1)*tan(x)
1
2

SymPy也知道如何计算一个表达式在一个点的Taylor序列。使用series(expr, var):

series(cos(x), x)
1

1 - x**2/2 + x**4/24 + O(x**6)
1
2

练习

• 计算limx→0sin(x)/x$\underset{x\to 0}{lim}\sin (x)/x$

• 计算log(x)对于x的导数。

limit(sin(x)/x,x,0)
1

1
1
2

diff(log(x),x)
1

1/x
1
2

SymPy支持超验基础和特殊函数的无限和有限积分，通过integrate() 功能, 使用了强大的扩展的Risch-Norman算法和启发式和模式匹配。你可以积分基本函数:

integrate(6*x**5, x)
1

x**6
1
2

integrate(sin(x), x)
1

-cos(x)
1
2

也可以计算一下有限积分:

integrate(sin(x),(x,0,2*pi))
1

0
1
2

integrate(sin(x),(x,0,pi))
1

2
1
2

不标准积分也支持:

integrate(exp(-x), (x, 0, oo))
1

1
1
2

SymPy可以求解线性代数方程，一个或多个变量:

solve(x**4 - 1, x)
1

[-1, 1, -I, I]
1
2

复合函数求解

solve([x + 5*x*y - 2, -3*x + 6*y - 15], [x, y])
1

[(-27/10 + sqrt(809)/10, 23/20 + sqrt(809)/20),
 (-sqrt(809)/10 - 27/10, -sqrt(809)/20 + 23/20)]
1
2
3

求解超越方程（有限的）:

solve(exp(x) + 1, x)
1

[I*pi]
1
2

多项式方程的另一个应用是factor。factor将多项式因式分解为可化简的项，并且可以计算不同域的因式:

f = x**4 - 3*x**2 + 1
factor(f)
1
2

(x**2 - x - 1)*(x**2 + x - 1)
1
2

#SymPy也可以解布尔方程，即，判断一个布尔表达式是否满足。对于这个情况，我们可以使用satisfiable函数:


satisfiable(x & y)
1
2
3
4

{x: True, y: True}
1
2

#这告诉我们(x & y)是真，当x和y都是True的时候。如果一个表达式不是True，即它的任何参数值都无法使表达式为真，那么它将返回False:


satisfiable(x & ~x)
1
2
3
4

False
1
2

3.2.4.1 练习

• 求解系统方程x+y=2$x+y=2$, 2⋅x+y=0$2\cdot x+y=0$ （求解方程组）
• 是否存在布尔值，使(~x | y) & (~y | x)$\text{(~x | y) \& (~y | x)}$为真?

方程组有一组数值解

原方程组如下：

2x+3y+z=4
4x+2y+3z =17
7x+y-z=1
1
2
3
4

求解方程组的代码Python代码如下：

from sympy import *

a = Matrix([[2,3,1],[4,2,3], [7,1,-1]])
b = Matrix([[4],[17],[1]])
x = symarray(‘x’, 3)
solve(a*x-b)

执行以上代码得到

{x_1: -1, x_2: 5, x_0: 1}
1
2

求解存在自由变量的方程组

代码如下：

from sympy import *

a = Matrix([[3,0,-1,0],[8,0,0,-2], [0,2,-2,-1]])
x = symarray(‘x’, 4)
solve(a*x)

得到：

{x_1: 5x_3/4, x_2: 3x_3/4, x_0: x_3/4}
1
2

3.2.5.2 微分方程

SymPy可以解 (一些) 常规微分。要求解一个微分方程，使用dsolve。首先，通过传递cls=Function来创建一个未定义的符号函数:

In [38]:

f, g = symbols(‘f g’, cls=Function)

f 和 g是未定义函数。我们可以调用f(x), 并且它可以代表未知的函数:

In [39]:

f(x)

Out[39]:

f(x)

In [40]:

f(x).diff(x, x) + f(x)

Out[40]:

f(x) + Derivative(f(x), x, x)

In [41]:

dsolve(f(x).diff(x, x) + f(x), f(x))

Out[41]:

f(x) == C1*sin(x) + C2*cos(x)

关键词参数可以向这个函数传递，以便帮助确认是否找到最适合的解决系统。例如，你知道它是独立的方程，你可以使用关键词hint=’separable’来强制dsolve来将它作为独立方程来求解:

In [42]:

dsolve(sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x), f(x), hint=’separable’)

Out[42]:

[f(x) == -asin(sqrt(C1/cos(x)**2 + 1)) + pi,
f(x) == asin(sqrt(C1/cos(x)**2 + 1)) + pi,
f(x) == -asin(sqrt(C1/cos(x)**2 + 1)),
f(x) == asin(sqrt(C1/cos(x)**2 + 1))]

练习

求解Bernoulli微分方程

$\frac{d f(x)}{x} + f(x) - f(x)^2=0$

使用hint=’Bernoulli’求解相同的公式。可以观察到什么?
1
2
3
4
5
6

1