import numpy as np
from sympy import *
sympy.__version__
'1.4'
E
e \displaystyle e e
I
i \displaystyle i i
pi
π \displaystyle \pi π
1. 欧拉公式
E**(I*pi)
− 1 \displaystyle -1 −1
e i x \displaystyle e^{i x} eix
x = symbols("x")
expand( E**(I*x), complex=True)
i e − im ( x ) sin ( re ( x ) ) + e − im ( x ) cos ( re ( x ) ) \displaystyle i e^{- \operatorname{im}{\left(x\right)}} \sin{\left(\operatorname{re}{\left(x\right)} \right)} + e^{- \operatorname{im}{\left(x\right)}} \cos{\left(\operatorname{re}{\left(x\right)} \right)} ie−im(x)sin(re(x))+e−im(x)cos(re(x))
x = symbols("x", real=True)
expand( E**(I*x), complex=True)
i sin ( x ) + cos ( x ) \displaystyle i \sin{\left(x \right)} + \cos{\left(x \right)} isin(x)+cos(x)
1.1 泰勒展开
证明欧拉公式展开的正确性
tmp = series(exp(I*x), x, 0, 10)
tmp
1 + i x − x 2 2 − i x 3 6 + x 4 24 + i x 5 120 − x 6 720 − i x 7 5040 + x 8 40320 + i x 9 362880 + O ( x 10 ) \displaystyle 1 + i x - \frac{x^{2}}{2} - \frac{i x^{3}}{6} + \frac{x^{4}}{24} + \frac{i x^{5}}{120} - \frac{x^{6}}{720} - \frac{i x^{7}}{5040} + \frac{x^{8}}{40320} + \frac{i x^{9}}{362880} + O\left(x^{10}\right) 1+ix−2x2−6ix3+24x4+120ix5−720x6−5040ix7+40320x8+362880ix9+O(x10)
re(tmp)
x 8 40320 − x 6 720 + x 4 24 − x 2 2 + re ( O ( x 10 ) ) + 1 \displaystyle \frac{x^{8}}{40320} - \frac{x^{6}}{720} + \frac{x^{4}}{24} - \frac{x^{2}}{2} + \operatorname{re}{\left(O\left(x^{10}\right)\right)} + 1 40320x8−720x6+24x4−2x2+re(O(x10))+1
im(tmp)
x 9 362880 − x 7 5040 + x 5 120 − x 3 6 + x + im ( O ( x 10 ) ) \displaystyle \frac{x^{9}}{362880} - \frac{x^{7}}{5040} + \frac{x^{5}}{120} - \frac{x^{3}}{6} + x + \operatorname{im}{\left(O\left(x^{10}\right)\right)} 362880x9−5040x7+120x5−6x3+x+im(O(x10))
series(sin(x), x, 0, 10)
x − x 3 6 + x 5 120 − x 7 5040 + x 9 362880 + O ( x 10 ) \displaystyle x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \frac{x^{9}}{362880} + O\left(x^{10}\right) x−6x3+120x5−5040x7+362880x9+O(x10)
series(cos(x), x, 0, 10)
1 − x 2 2 + x 4 24 − x 6 720 + x 8 40320 + O ( x 10 ) \displaystyle 1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} - \frac{x^{6}}{720} + \frac{x^{8}}{40320} + O\left(x^{10}\right) 1−2x2+24x4−720x6+40320x8+O(x10)
series(I*sin(x)+cos(x), x, 0, 10)
1 + i x − x 2 2 − i x 3 6 + x 4 24 + i x 5 120 − x 6 720 − i x 7 5040 + x 8 40320 + i x 9 362880 + O ( x 10 ) \displaystyle 1 + i x - \frac{x^{2}}{2} - \frac{i x^{3}}{6} + \frac{x^{4}}{24} + \frac{i x^{5}}{120} - \frac{x^{6}}{720} - \frac{i x^{7}}{5040} + \frac{x^{8}}{40320} + \frac{i x^{9}}{362880} + O\left(x^{10}\right) 1+ix−2x2−6ix3+24x4+120ix5−720x6−5040ix7+40320x8+362880ix9+O(x10)
2. 球体积分
integrate() 计算符号积分
integrate()计算不定积分:
integrate(x*sin(x), x)
− x cos ( x ) + sin ( x ) \displaystyle - x \cos{\left(x \right)} + \sin{\left(x \right)} −xcos(x)+sin(x)
integrate()计算定积分:
integrate(x*sin(x), (x,0,2*pi))
− 2 π \displaystyle - 2 \pi −2π
计算圆面积:
x, y = symbols("x, y")# symbols() 可以一次创建多个符号
r = symbols("r", positive=True)
circle_area = 2 * integrate(sqrt(r**2-x**2), (x,-r,r) )
circle_area
π r 2 \displaystyle \pi r^{2} πr2
计算球形体积:
r随着x的变化而变化,r = sqrt(R2-x2)
x, y = symbols("x, y")# symbols() 可以一次创建多个符号
r = symbols("r", positive=True)
circle_area = 2 * integrate(sqrt(r**2-x**2), (x,-r,r) )
circle_area = circle_area.subs(r, sqrt(r**2-x**2))
circle_area
π ( r 2 − x 2 ) \displaystyle \pi \left(r^{2} - x^{2}\right) π(r2−x2)
integrate(circle_area, (x,-r,r))
4 π r 3 3 \displaystyle \frac{4 \pi r^{3}}{3} 34πr3
- expression.subs(x,y): 将算式中的x替换为y
- expression.subs({x:y, u:v}): 字典多次替换,依次将算式中的x,u 替换为y,v
- expression.subs([(x,y), (u,v)]): 列表多次替换
3. 数值微分
as_finite_diff() : 自动计算N点公式
3.1 符号微分
x = symbols("x", real=True)
h = symbols("h", positive=True)
f = symbols("f", cls=Function)
f_diff = f(x).diff()
f_diff
d d x f ( x ) \displaystyle \frac{d}{d x} f{\left(x \right)} dxdf(x)
f_diff = f(x).diff(x)
f_diff
d d x f ( x ) \displaystyle \frac{d}{d x} f{\left(x \right)} dxdf(x)
f_diff = f(x).diff(x,1)#对x求1阶导数
f_diff
d d x f ( x ) \displaystyle \frac{d}{d x} f{\left(x \right)} dxdf(x)
3.2 符号微分: 转化为多点共同运算
expr_diff = as_finite_diff(f_diff, [x, x-h, x-2*h, x-3*h])
expr_diff
F:\Anaconda\Anaconda3_201910\lib\site-packages\sympy\core\decorators.py:38: SymPyDeprecationWarning:
_as_finite_diff has been deprecated since SymPy 1.1. Use
Derivative.as_finite_difference instead. See
https://github.com/sympy/sympy/issues/11410 for more info.
_warn_deprecation(wrapped, 3)
11 f ( x ) 6 h − f ( − 3 h + x ) 3 h + 3 f ( − 2 h + x ) 2 h − 3 f ( − h + x ) h \displaystyle \frac{11 f{\left(x \right)}}{6 h} - \frac{f{\left(- 3 h + x \right)}}{3 h} + \frac{3 f{\left(- 2 h + x \right)}}{2 h} - \frac{3 f{\left(- h + x \right)}}{h} 6h11f(x)−3hf(−3h+x)+2h3f(−2h+x)−h3f(−h+x)
expr_diff = Derivative.as_finite_difference(f_diff, [x, x-h, x-2*h, x-3*h])
expr_diff
11 f ( x ) 6 h − f ( − 3 h + x ) 3 h + 3 f ( − 2 h + x ) 2 h − 3 f ( − h + x ) h \displaystyle \frac{11 f{\left(x \right)}}{6 h} - \frac{f{\left(- 3 h + x \right)}}{3 h} + \frac{3 f{\left(- 2 h + x \right)}}{2 h} - \frac{3 f{\left(- h + x \right)}}{h} 6h11f(x)−3hf(−3h+x)+2h3f(−2h+x)−h3f(−h+x)
expr_diff.args
(-3*f(-h + x)/h, -f(-3*h + x)/(3*h), 3*f(-2*h + x)/(2*h), 11*f(x)/(6*h))
3.3 比较 数值求导 和 符号求导的误差
f(x) = x*exp(-x^2)
f_diff.subs(f(x), xexp(-x^2)) : 将f(x)替换为 xexp(-x^2)
doit(): 计算导数公式(symbols,符号,非数值)
sym_dexpr = f_diff.subs(f(x), x*exp(-x**2)).doit()
sym_dexpr
− 2 x 2 e − x 2 + e − x 2 \displaystyle - 2 x^{2} e^{- x^{2}} + e^{- x^{2}} −2x2e−x2+e−x2
1. 符号微分公式转化为数值计算的函数
为后面计算 符号求导的结果做准备
sym_dfunc = lambdify([x], sym_dexpr, modules="numpy")
sym_dfunc(np.array([-1,-.5, 0, 0.5, 1]) )
array([-0.36787944, 0.38940039, 1. , 0.38940039, -0.36787944])
- 通配符与模板:
求得 coefficient,为数值求导 做准备
expr_diff
expr_diff.args
(-3*f(-h + x)/h, -f(-3*h + x)/(3*h), 3*f(-2*h + x)/(2*h), 11*f(x)/(6*h))
w = Wild("w")
c = Wild("c")
pattern = [arg.match(c*f(w)) for arg in expr_diff.args]
pattern
[{w_: -h + x, c_: -3/h},
{w_: -3*h + x, c_: -1/(3*h)},
{w_: -2*h + x, c_: 3/(2*h)},
{w_: x, c_: 11/(6*h)}]
coefficients = [t[c] for t in sorted(pattern, key=lambda t:t[w])]
coefficients
[-1/(3*h), 3/(2*h), -3/h, 11/(6*h)]
将系数表达式列表中的h替换为0.001,
得到系数数组时Sympy.Float对象
需要使用float将其转化为python的float格式
coeff_arr = np.array([np.float(coeff.subs(h, .001)) for coeff in coefficients])
coeff_arr
array([ -333.33333333, 1500. , -3000. , 1833.33333333])
3. 计算数值求导的结果
def moving_window(x, size):
from numpy.lib.stride_tricks import as_strided
x = np.ascontiguousarray(x)
return as_strided(x, shape=(x.shape[0]-size+1, size),
strides=(x.itemsize, x.itemsize))
x_arr = np.arange(-2,2, 1e-3)
y_arr = x_arr * np.exp(-x_arr*x_arr)
num_res = (moving_window(y_arr, 4)*coeff_arr).sum(axis=1)
num_res
array([-0.12931154, -0.12968029, -0.13004973, ..., -0.12931154,
-0.12894349, -0.12857613])
4. 求符号求导计算得到的结果
sym_res = sym_dfunc(x_arr[3:])
sym_res
array([-0.12931154, -0.12968029, -0.13004973, ..., -0.12931154,
-0.12894349, -0.12857613])
5. 比对 符号求导 和 数值求导 的误差
print( np.max(np.abs( num_res - sym_res)) )
4.089441674182126e-09
3.4 比较点数与误差之间的关系