【AI数学】极限篇与求导篇python(学习记录)

一、关键库简介

  (1)SymPy库
   SymPy库是Python的一个科学计算库，用强大的符号计算体系完成诸如多项式求值、求极限、求导、解方程、求积分、解微分方程、级数展开、矩阵运算等计算。
   (2)SciPy库
   SciPy函数库在Numpy库的基础上增加了众多的数以及工程计算中常用的库函数。可以说，Numpy是一个纯数学的计算模块，而SciPy是一个更高阶的科学计算库。

二、极限篇

  (1)$$\underset{x\to \infty }{\lim }\frac{sinx}{x}$$

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)$$\underset{x\to 1}{\lim }\frac{x^2-1}{x-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 \sqrt{sinx}$$

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$$

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+2x^4+3x^2+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+3xy+y^2$$

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编程实现梯度下降法求解下面函数的最小值，并绘图
$$minf(x,y)=x-y+2x^2+2xy+y^2$$
  给定初始值$X^0=(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

参考资料
人工智能数学基础（北京大学出版社）