一、关键库简介
(1)SymPy库
SymPy库是Python的一个科学计算库,用强大的符号计算体系完成诸如多项式求值、求极限、求导、解方程、求积分、解微分方程、级数展开、矩阵运算等计算。
(2)SciPy库
SciPy函数库在Numpy库的基础上增加了众多的数以及工程计算中常用的库函数。可以说,Numpy是一个纯数学的计算模块,而SciPy是一个更高阶的科学计算库。
二、极限篇
(1) lim x → ∞ s i n x x { \lim_{x \to \infty} \frac{sinx}{x} } x→∞limxsinx
import sympy
from sympy import oo # 注意无穷符号表示形式为两个小写字母 o
x = sympy.Symbol('x') # 注意 Symbol 第一字母大写
f = sympy.sin(x)/x
print (sympy.limit(f,x,oo))
输出:0
(2) lim x → 1 x 2 − 1 x − 1 { \lim_{x \to 1} \frac{x^2-1}{x-1} } x→1limx−1x2−1
import sympy
x = sympy.Symbol('x')
f = (x**2-1)/(x-1)
sympy.limit(f,x,1)
print(sympy.limit(f,x,1))
输出:2
二、求导篇
(1)求导数(使用SymPy库中的diff函数) y = arcsin s i n x y={\arcsin\sqrt{sinx} } y=arcsinsinx
from sympy import *
from sympy.abc import x,y,z,f
# diff 求导函数,arcsinx 数学函数表示形式为 asin
print(diff(asin(sqrt(sin(x)))))
输出:cos(x)/(2*sqrt(1 - sin(x))*sqrt(sin(x)))
(2)求导数(使用SciPy库misc模块下的derivative函数) y = x 5 y={x^5 } y=x5
import numpy as np
from scipy.misc import derivative
def f(x): # 定义函数
return x**5
print (derivative(f, 2, dx=1e-6)) # 对函数在 x=2 处求导
输出:80.00000000230045
(3)求导数(使用numpy库里的polyld函数构造函数) y = x 5 + 2 x 4 + 3 x 2 + 5 y={x^5+2x^4+3x^2+5 } y=x5+2x4+3x2+5
import numpy as np
p = np.poly1d([1,2,0,3,0,5]) # 构造多项式
print('-------------------')
print(p) # 下面两行为多项式的显示形式,5、4、2 是下一行项数所对应的幂次
print('-------------------')
print(np.polyder(p,1)) # 求一阶导数
print(np.polyder(p,1)(1.0)) # 求一阶导数在点 x=1 的值
print('-------------------')
print(p.deriv(1)) # 求一阶导数
print(p.deriv(1)(1.0)) # 求一阶导数在点 x=1 的值
输出:
-------------------
5 4 2
1 x + 2 x + 3 x + 5
-------------------
4 3
5 x + 8 x + 6 x
19.0
-------------------
4 3
5 x + 8 x + 6 x
19.0
(4)求函数在点(1,2)处的偏导数 f ( x , y ) = x 2 + 3 x y + y 2 f(x,y)={x^2+3xy+y^2 } f(x,y)=x2+3xy+y2
from sympy import *
from sympy.abc import x,y,z,f
f = x**2 + 3*x*y + y**2
print(diff(f,x)) # 求偏导
print(diff(f,y))
fx = diff(f,x) # 求偏导并将结果赋给 fx
print(fx.evalf(subs={x:1,y:2})) # 以字典的形式传入多个变量的值,求函数值。
fy = diff(f,y)
print(fy.evalf(subs={x:1,y:2}))
输出:
2*x + 3*y
3*x + 2*y
8.00000000000000
7.00000000000000
三、提高篇
应用Python编程实现梯度下降法求解下面函数的最小值,并绘图
m
i
n
f
(
x
,
y
)
=
x
−
y
+
2
x
2
+
2
x
y
+
y
2
minf(x,y)={x-y+2x^2+2xy+y^2 }
minf(x,y)=x−y+2x2+2xy+y2
给定初始值
X
0
=
(
0
,
0
)
T
X^0={(0,0)^T }
X0=(0,0)T
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
def Fun(x, y):
return x-y + 2 * x * x + 2 * x * y + y * y
def PxFun(x, y): # 求 x 偏导
return 1 + 4 * x + 2 * y
def PyFun(x, y): # 求 y 偏导
return -1 + 2 * x + 2 * y
fig = plt.figure() # figure 对象
ax = Axes3D(fig) # Axes3D 对象
X,Y = np.mgrid[-2:2:40j, -2:2:40j] # 取样并作满射联合
Z = Fun(X, Y) # 取样点 Z 坐标打表
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap="rainbow")
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# 梯度下降
step = 0.0008 # 取步长
x = 0
y = 0
tag_x = [x]
tag_y = [y]
tag_z = [Fun(x, y)] # 三个坐标分别打入表中,该表用于绘制点
new_x = x
new_y = y
Over = False
while Over == False:
new_x -= step * PxFun(x, y)
new_y -= step * PyFun(x, y)
if Fun(x, y) - Fun(new_x, new_y) < 7e-9:
Over = True
x = new_x # 更新旧点
y = new_y
tag_x.append(x)
tag_y.append(y)
tag_z.append(Fun(x, y)) # 新点三个坐标打入表中
ax.plot(tag_x,tag_y,tag_z,'r')
plt.title('(x,y)~(' + str(x) + "," + str(y) + ')')
plt.show()
求解为:x ≈ -0.9979,y ≈ 1.4967
参考资料
人工智能数学基础(北京大学出版社)