import numpy as np
from sympy import *
sympy.__version__
'1.4'
E
I
pi
1. 欧拉公式
E**(I*pi)
x = symbols("x")
expand( E**(I*x), complex=True)
x = symbols("x", real=True)
expand( E**(I*x), complex=True)
1.1 泰勒展开
证明欧拉公式展开的正确性
tmp = series(exp(I*x), x, 0, 10)
tmp
re(tmp)
im(tmp)
series(sin(x), x, 0, 10)
series(cos(x), x, 0, 10)
series(I*sin(x)+cos(x), x, 0, 10)
2. 球体积分
integrate() 计算符号积分
integrate()计算不定积分:
integrate(x*sin(x), x)
integrate()计算定积分:
integrate(x*sin(x), (x,0,2*pi))
计算圆面积:
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随着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
integrate(circle_area, (x,-r,r))
- 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
f_diff = f(x).diff(x)
f_diff
f_diff = f(x).diff(x,1)#对x求1阶导数
f_diff
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)
expr_diff = Derivative.as_finite_difference(f_diff, [x, x-h, x-2*h, x-3*h])
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))
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
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 比较点数与误差之间的关系