前言

• 买了一本迪哥写的《人工智能数学基础》，在掌握了Python的基本用法后，着手学习其中的内容
• 里面很多都是考研数学的内容，让我意识到数学没有想象中的那么“没用”
• Python有很多封装好的数学工具，正文中有大量的数学公式使用了LateX的语句，参考在线LaTeX公式编辑器-编辑器 (latexlive.com)，很遗憾没能在本科期间认识到它们
• 有些LateX语句并不能正确地渲染, 先凑合着用吧

1 基础篇

2 高等数学基础

例2.6 求$\lim_{x \to 1}\frac{x ^ 2 - 1}{x - 1}$ 的极限

import sympy
from sympy import oo
import numpy as np
x = sympy.Symbol('x')
f = (x ** 2 - 1) / (x - 1)
sympy.limit(f, x, 1)  # limit(数学表达式expression, 变量variable, 要趋向的值value)

$\displaystyle 2$

例2.11 求$y = arcsin\sqrt{sinx}$的导数

$y'=\frac{1}{\sqrt{1-sinx} }\cdot\frac{1}{2\sqrt{sinx}}\cdot{cosx}=\frac{cosx}{2\sqrt{sinx-sin^2x}}$

from sympy import *
from sympy.abc import x  # 导入变量x
diff(asin(sqrt(sin(x))))  # diff求导函数

$\displaystyle \frac{\cos{\left(x \right)}}{2 \sqrt{1 - \sin{\left(x \right)}} \sqrt{\sin{\left(x \right)}}}$

例2.12 求$f(x,y)=x^2+3xy+y^2$在点$(1,2)$处的偏导数

$f_{x}(x,y)=2x+3y$

$f _{y}(x,y)=3x+2y$

$ f _{x}(1,2)=2x+3y| _{y=2}^{x=1}=8$

$f _{y}(1,2)=3x+2y| _{y=2}^{x=1}=7$

from sympy import *
from sympy.abc import x, y, f
f = x ** 2 + 3 * x * y + y ** 2
diff(f, x)  # 对x求偏导

$\displaystyle 2 x + 3 y$

diff(f, y)

$\displaystyle 3 x + 2 y$

fx = diff(f, x)
fx.evalf(subs={x:1, y:2})  # 以字典的形式传入多个变量的值, 返回计算后的数学表达式。

$\displaystyle 8.0$

fy = diff(f, y)
fy.evalf(subs={x:1, y:2})

$\displaystyle 7.0$

2.6 方向导数 2.7 梯度

函数在某点的梯度是一个向量, 它的方向与取得最大方向导数的方向一致

$gradf(x,y)=\frac{\partial f}{\partial x}\textbf{i}+\frac{\partial f}{\partial y}\textbf{j} $

$\frac{\partial f}{\partial l}=gradf\cdot \frac{\overrightarrow{l}}{\left | \overrightarrow{l} \right | }$

用梯度gradf点乘l的单位向量就得到方向导数, 这是计算方向导数最简便的方法

例2.15 用Python编程实现梯度下降法求解$f(x,y)=x-y+2x^2+2xy+y^2$的最小值

使用常规方法:

$z=x-y+2x^2+2xy+y^2$

$z'_{x}=1+4x+2y$

$z'_{y}=-1+2x+2y$

求得驻点$(-1,\frac{3}{2})$

$A=z''_{xx}| _{y=\frac{3}{2}}^{x=-1}=4$

$B=z''_{xy}| _ {y=\frac{3}{2}}^{x=-1}=2$

$C=z''_{yy}| _ {y=\frac{3}{2}}^{x=-1}=2$

$AC-B^2>0$

$A>0$

原函数在$(-1,\frac{3}{2})$处具有极小值

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


if __name__ == "__main__":
    fig = plt.figure()  # 创建自定义图象
    # MatplotlibDeprecationWarning: Axes3D(fig) adding itself to the figure is deprecated since 3.4. 
    # Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. 
    # The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6.
    ax = Axes3D(fig, auto_add_to_figure=False)  # 创建3D图象
    fig.add_axes(ax)
    # 在绘制三维图表时, 需要用到mgrid函数, 它会返回一个密集的多维网格
    # np.mgrid[开始坐标:结束坐标:步长(步长为复数表示点数，左闭右闭, 步长为实数表示间隔，左闭右开)]
    X, Y = np.mgrid[-2:2:40j, -2:2:40j]  
    Z = Fun(X, Y)
    # ax.plot_surface: https://blog.csdn.net/weixin_43584807/article/details/102331755
    # X, Y, Z 2D数组形式的数据值
    # rstride 数组行距
    # cstride 数组列距
    # cmap 曲面块颜色映射
    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.008  # 取步长
    x = 0  # 初始值
    y = 0
    tag_x = [x]  # 绘制图像的列表
    tag_y = [y]
    tag_z = [Fun(x, y)]
    new_x = x
    new_y = y
    Over = False
    while not Over:
        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+]')  # 绘制点: 'r'表示红色(好像没有作用?)
    print('(x, y)~(' + str(x) + ',' + str(y) + ')')  # 输出结果
    plt.show()  # 显示图象
    

(x, y)~(-0.9993608094022046,1.498965767887478)

2.9.5 高手点拨: 求导的三种方式

已知$f(x)=x^5+2x^4+3x^2+5$, 求$f'(1)$

使用Sympy的diff函数

import sympy
from sympy.abc import x, f

f = x ** 5 + 2 * x ** 4 + 3 * x ** 2 + 5
fx = sympy.diff(f)
fx.evalf(subs={x:1})

$\displaystyle 19.0$

使用scipy.misc模块下的derivative函数

import numpy as np
from scipy.misc import derivative


def f(x):
    return x ** 5 + 2 * x ** 4 + 3 * x ** 2 + 5

derivative(func=f, x0=1, dx=1e-6, n=1)  # 函数为f, 取值为1, 间距为1e-6, 一阶导数

18.999999999991246

使用NumPy模块里的poly1d构造$f(x)$

import numpy as np

p = np.poly1d([1, 2, 0, 3, 0, 5])  # 构造多项式
np.polyder(p, 1)(1.0)  # 一阶导数在x=1.0时的取值

19.0

p.deriv(1)(1.0)  # 一阶导数在x=1.0时的取值

19.0

2.10 习题

$lim_{x\to1}sin(lnx)$

import sympy
from sympy.abc import x, f

f = sympy.sin(sympy.log(x))
sympy.limit(f, x, 1)

$\displaystyle 0$

$lim_{x\to8}\frac{\sqrt[3]{x}-2}{x-8}$

import sympy
from sympy.abc import x, f

f = (x ** (1/3) - 2) / (x - 8)
sympy.limit(f, x, 8)

$\displaystyle \frac{1}{12}$

求$y=x^4-2x^3+5sinx+ln3$的导数

import sympy
from sympy.abc import x, y

y = x ** 4 - 2 * x ** 3 + 5 * sympy.sin(x) + sympy.log(3)
diff(y)

$\displaystyle 4 x^{3} - 6 x^{2} + 5 \cos{\left(x \right)}$

求$z=(3x^2+y^2)^{4x+2y}$在点$(1,2)$的偏导数

$lnz=(4x+2y)\cdot ln(3x^2+y^2)$

$\frac1zdz=\left[4ln(3x^2+y^2)+\frac{24x^2+12xy}{3x^2+y^2}\right]dx+\left[2ln(3x^2+y^2)+\frac{8xy+2y^2}{3x^2+y^2}\right]dy$

import sympy
from sympy.abc import x, y, z

z = (3 * x ** 2 + y ** 2) ** (4 * x + 2 * y)
zx = sympy.diff(z, x)
zx

$\displaystyle \left(3 x^{2} + y^{2}\right)^{4 x + 2 y} \left(\frac{6 x \left(4 x + 2 y\right)}{3 x^{2} + y^{2}} + 4 \log{\left(3 x^{2} + y^{2} \right)}\right)$

zx.evalf(subs={x:1, y:2})

$\displaystyle 84401203.0927369$

zy = sympy.diff(z, y)
zy

$\displaystyle \left(3 x^{2} + y^{2}\right)^{4 x + 2 y} \left(\frac{2 y \left(4 x + 2 y\right)}{3 x^{2} + y^{2}} + 2 \log{\left(3 x^{2} + y^{2} \right)}\right)$

zy.evalf(subs={x:1, y:2})

$\displaystyle 48788945.5463684$

求方向导数和梯度

求函数$z=x^2+y^2$在点$(1,2)$处沿点$(1,2)$到点$(2,2+\sqrt{3})$方向的方向导数, 以及在点$(1,2)$的梯度

import sympy as sp
import numpy as np
from sympy.abc import x, y, z

z = x ** 2 + y ** 2
zx, zy = z.diff(x), z.diff(y)
gradz = np.array([zx.evalf(subs={x:1, y:2}), zy.evalf(subs={x:1, y:2})], dtype=float)  # 求梯度
gradz

array([2., 4.])

A = np.array([1, 2])
B = np.array([2, 2 + 3 ** 0.5], dtype=float)
gradz * (B - A) / np.linalg.norm(B - A)  # 计算方向导数

array([1.        , 3.46410162])

3 微积分

例3.6 定积分

应用SciPy科学计算库求$\int_{0}^{3} cos^2(e^x)dx$

import numpy as np
from scipy.integrate import quad

func = lambda x:np.cos(np.exp(x)) ** 2  # 定义被积分函数
quad(func, 0, 3)  # 调用quad积分函数

(1.296467785724373, 1.397797133112089e-09)

输出结果(积分值, 误差)

例3.7 二重积分

求$\iint_{D}e^{-x^x-y^2}dxdy$, 其中$D=\left \{ (x,y) | 0 \le x \le 10,0 \le y \le 10\right \}$

import numpy as np
from scipy.integrate import dblquad  # 二重积分


def integrand(x,y):
    return np.exp(-x ** 2 - y ** 2)

x_a = 0
x_b = 10
y_a = 0
y_b = 10
"""
scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)
参数：
func： 可调用的
至少有两个变量的 Python 函数或方法：y 必须是第一个参数，x 必须是第二个参数。

a, b： 浮点数
x:a < b 的积分极限

gfun： 可调用或浮点数
y 中的下边界曲线，它是一个函数，采用单个浮点参数 (x) 并返回浮点结果或表示恒定边界曲线的浮点数。

hfun： 可调用或浮点数
y 中的上边界曲线(与 gfun 的要求相同)。

args： 顺序，可选
传递给 func 的额外参数。

epsabs： 浮点数，可选
绝对容差直接传递到内部一维正交积分。默认值为 1.49e-8。 dblquad 尝试获得 abs(i-result) <= max(epsabs, epsrel*abs(i)) 的精度，其中 i = 从 gfun(x) 到 hfun(x) 的 func(y, x) 的内部积分，而 result 是数值近似值。请参阅下面的 epsrel。

epsrel： 浮点数，可选
内部一维积分的相对容差。默认值为 1.49e-8。如果epsabs <= 0, epsrel 必须大于 5e-29 和50 * (machine epsilon).看易胜宝更多。

返回：
y： 浮点数
结果积分。

abserr： 浮点数
误差的估计。
"""
dblquad(integrand, x_a, x_b, lambda x:y_a, lambda x:y_b)

(0.7853981633974476, 1.375309851021853e-08)

例3.8 定积分近似求解

用定义法求$\int_{0}^{3} cos^2(e^x)dx$的近似解

from numpy import *

a, b = 0, 3


def f(x):
    return cos(exp(x)) ** 2


def trape(n):
    h = (b - a) / n  # 将区间分成n等分后的长度
    x = a
    sum = 0
    for i in range(1, n):
        x2 = a + i * h
        sum += (f(x) + f(x2)) * h / 2  # 小梯形的值(上底加下底的和乘高除以二)
        x = x2
    return sum

trape(10)

0.944822326405313

trape(100)

1.2843391540917448

trape(1000)

1.2960750567338157

不定积分

$ \int lnx $

from sympy import *
from sympy.abc import x

expr = log(x)
integrate(expr,x)

$\displaystyle x \log{\left(x \right)} - x$

3.8习题

$\int _{1}^{2} x^2 + \frac {1}{x^4}dx$

$ =(\frac {1}{3}x^3 - \frac {1}{3}x^{-3}) | _{1}^{2}$

$ = \frac {21}{8}$

import numpy as np
from scipy.integrate import quad

func = lambda x: x ** 2 + x ** (-4)
quad(func, 1, 2)

(2.625, 2.914335439641036e-14)

$ \int _{-1}^{0}\frac {3x^4+3x^2+1}{x^2+1}dx$

$ =\int _{-1}^{0}3x^2+\frac {1}{1+x^2}dx$

$ = (x^3 + arctanx)| _{-1}^{0} $

$ = 1 + \frac{\pi}{4} $

import numpy as np
from scipy.integrate import quad

func = lambda x: (3 * x ** 4 + 3 * x ** 2 + 1) / (1 + x ** 2)
quad(func, -1, 0)

(1.7853981633974483, 1.9821901491273144e-14)

利用定积分的定义计算极限:

$ lim _{n \to \infty}\frac {1^{p}+2^{p}+...+n^{p}}{n^{p+1}}$

$ = lim _{n \to \infty}\frac {\sum_{i=1}^{n}i^p}{n^{p+1}} $

$ = \frac{1}{n} lim _{n \to \infty}\sum_{i=1}^{n}(\frac{i}{n})^p $

$ = \int _{0}^{1}x^pdx $

$ = \frac {1}{p+1} $

4 泰勒公式与拉格朗日乘子法

指对连, 三角断,
三角对数隔一换, 三角指数有感叹

例4.7 根据$e^x$的$n$次泰勒多项式展开式, 计算$e$的近似值

$ e\approx 1+1+\frac{1}{2!}+\frac{1}{3!}+...+\frac{1}{n!}$

import numpy as np
import pandas as pd


def f(n):
    sum1 = 1
    if n == 0:
        sum1 = 1
    else:
        m = n + 1
        for i in range(1, m):
            sum2 = 1.0
            k = i + 1
            for j in range(1, k):
                sum2 = sum2 * j
            sum1 = sum1 + 1.0 / sum2
    return sum1


num = 10
pd.DataFrame(np.array([[i,f(i)] for i in range(1, num + 1)]), columns=['n', 'e'])

	n	e
0	1.0	2.000000
1	2.0	2.500000
2	3.0	2.666667
3	4.0	2.708333
4	5.0	2.716667
5	6.0	2.718056
6	7.0	2.718254
7	8.0	2.718279
8	9.0	2.718282
9	10.0	2.718282

例 4.13 $sinx$的$n$阶泰勒多项式

$sinx = x - \frac {x^3}{3!} + \frac {x^5}{5!} - \frac {x^7}{7!} + ...+ (-1)^{m-1}\frac {x^{2m-1}}{(2m-1)!}+R_{2m}(x)$

import numpy as np
import pandas as pd


def fsin(x):
    m = 20  # 项数
    sum = 0.0
    for i in range(1, m+1):
        n = 2 * i - 1  # 2m-1
        temp1, temp2, temp3 = 1, 1, 1
        for j in range(1, i): # (-1)^(m-1)
            temp1 = -temp1
        for j in range(1, n + 1):
            temp2 = temp2 * x  # x^(2m-1)
            temp3 = temp3 * j  # (2m-1)!
        sum += temp1 * temp2 / temp3
    return sum


pd.DataFrame({'np.sin(x)': np.array([np.sin(x) for x in range(-20, 1)]),
              'fsin(x)': np.array([fsin(x) for x in range(-20, 1)]), 
              'error': np.array([fsin(x) - np.sin(x) for x in range(-20, 1)])},
              index=np.array([x for x in range(-20, 1)]))

	np.sin(x)	fsin(x)	error
-20	-0.912945	5364.411846	5.365325e+03
-19	-0.149877	666.994385	6.671443e+02
-18	0.750987	74.739042	7.398806e+01
-17	0.961397	8.185042	7.223645e+00
-16	0.287903	0.899283	6.113793e-01
-15	-0.650288	-0.606249	4.403901e-02
-14	-0.990607	-0.987967	2.640574e-03
-13	-0.420167	-0.420039	1.282664e-04
-12	0.536573	0.536578	4.880595e-06
-11	0.999990	0.999990	1.394300e-07
-10	0.544021	0.544021	2.831679e-09
-9	-0.412118	-0.412118	3.790196e-11
-8	-0.989358	-0.989358	2.858824e-13
-7	-0.656987	-0.656987	3.441691e-15
-6	0.279415	0.279415	3.497203e-15
-5	0.958924	0.958924	-2.775558e-15
-4	0.756802	0.756802	-6.661338e-16
-3	-0.141120	-0.141120	1.110223e-16
-2	-0.909297	-0.909297	0.000000e+00
-1	-0.841471	-0.841471	0.000000e+00
0	0.000000	0.000000	0.000000e+00

4.9 习题

$lim _{x \to 0}(\frac {sinx-xcosx}{sin^3x})$

$=lim _{x \to 0}\frac {x-\frac{1}{6}x^3-x(1-\frac{1}{2}x^2)}{x^3}$

$=lim _{x \to 0}\frac {-\frac{1}{6}x^3+\frac{1}{2}x^3}{x^3}$

$=\frac{1}{3}$

import sympy
from sympy import sin, cos, limit

x = sympy.Symbol('x')
f = (sin(x) - x * cos(x)) / (sin(x) ** 3)
limit(f, x, 0)

$\displaystyle \frac{1}{3}$

上一篇：Python-人工智能数学基础(5-6)

下一篇：Python-莫烦python学习笔记（multithreading）