经典的Euler法提供了认识、求解ODE基本方法和思路，其他的求解方法均是建立在该方法基础上的一个延伸，认识、理解Euler法是非常重要的。在经典Euler法中，函数值的微分$\mathrm{d}y$是依据全微分理论构建得到的，因为ODE自身提供了求解区域中某一点处的斜率值，只要规定好步长，那么对应的函数值增量就得以确定，进而可以得到下一个点处的函数值，反复进行这一过程即可完成求解。而搜寻下一个点的关键就是斜率值，Euler中斜率值是取的当前的点，即$t_i,y_i$处的斜率值，来表示$y_i与y_{i+1}$之间的斜率，如果可以缩小步长，那么精度会上升，另一个角度来看，如果利用该区段中点处的斜率值，可能效果要更好，更具有代表性，基于这一个改进思路就得到改进的Euler法，也称为显示的梯形方法。Euler法修正的思路Euler法改进的缘由Euler法提供的思路是基于斜率场来完成的，对于某一点来说，$df = f(t_i,y_i)'$表示斜率值，那么依据全微分理论，自然可以得到对应区间上的函数值$y$的微分$\mathrm{d}y$值，它反映了$t_i, y_i$到$t_{i+1}, y_{i+1}$的函数值的变化大小，那么由此得到$y_{i+1}$的值。一般来说，$y$值表示ODE的真实值，而数值计算的值是一个估计值（近似值）,为了区别两者，设$w$表示函数值的估计值。$$ w_{i + 1} = w_{i} + h \times df(t_{i}, w_{i}) $$前述已经说明了，改进的一种方式就是对斜率值进行修正，经典法中是利用区段左端点的斜率值来作为衡量的标准，如果区段长度过大，精度下降的程度较大。斜率值的确定取区间两个端点处斜率的平均值作为该区段的斜率值，并用于计算。$$ \begin{equation} slope = \frac{1}{2}(df(t_i, w_i) + df(t_{i + 1}, w_{i+1})\\ slope = \frac{1}{2}(df(t_i, w_i) + df(t_{i} + h, w_i + h \cdot df(t_i,w_i))) \end{equation} $$改进的Euler法，依然延续了经典法的计算思路，由于经典Euler法可以计算得到下一个点处的斜率值，然后取钱一个点和后一个点处的斜率值的平均值作为斜率的修正值，再利用Euler法进行修正后的计算，这就是改进的Euler法，因为该公式与梯形面积的计算公式相似，因此也叫作显式梯形法。$$ \begin{equation} w_{0} = y_0\\ w_{i+1} = w_{i} + \frac{h}{2} \cdot (df(t_i, w_i) + df(t_i + h, w_i + hdf(t_i, w_i))) \end{equation} $$通过上述的说明，可以看出改进Euler法本质还是经典的Euler法，只不过对于一次计算，它是在Euler法基础上进行了第二次迭代计算用于修正。计算验证与对比依然以之前经典Euler法中的例子进行计算，这样方便对比两者的结果差异。$$ \begin{equation} y' = t y + t ^ 3 \quad (0\leq t \geq 1)\\ y_0 = 1 \end{equation} $$对于$t$这里共划分10段，即步长$h=0.1$进行ODE的计算。改进Euler法的计算代码如下，这里由于考虑到对比，针对同一个ODE进行计算，只不过需要利用两种方法，因此可以通过构建一个ODE类，而两个方法可以绑定为实例的方法，并且提供一个可视化方法。import numpy as np import matplotlib.pyplot as plt class Euler_Method_Analysis(object): def __init__(self, interval, y0, n, df): self.h = (interval[-1] - interval[0]) / n self.t = np.arange(interval[0], interval[-1]+self.h, self.h) self.y0 = y0 self.df = df self.classical_euler() self.modified_euler() def classical_euler(self): w = np.zeros(self.t.shape[0]) w[0] = self.y0 for i in np.arange(self.t.shape[0]-1): w[i+1] = w[i] + self.h * self.df(self.t[i], w[i]) self.w_classical = w return self.w_classical def modified_euler(self): w = np.zeros(self.t.shape[0]) w[0] = self.y0 for i in np.arange(self.t.shape[0] - 1): w[i + 1] = w[i] + 0.5 * self.h * (self.df(self.t[i], w[i]) + self.df(self.t[i] + self.h, w[i] + self.h * self.df(self.t[i], w[i]))) self.w_modified = w return self.w_modified def view_result(self): plt.style.use('ggplot') fig, ax = plt.subplots() ax.plot(self.t, self.w_classical, lw=2, marker='o') ax.plot(self.t, self.w_modified, lw=2, marker='s') ax.set_title('Result of $ df = ty + t^3$') ax.set_xlabel('t') ax.set_ylabel('y') label = ['Classical Euler Method', 'Modified Euler Method'] ax.legend(label, loc='best') fig.savefig('resfig.png', dpi=600) plt.show() # 一下为调用上述类，实现计算 def func(t, y): return t * y + t**3 res = Euler_Method_Analysis([0, 1], 1, 10, func) res.w_classical res.w_modified res.view_result()
一阶常微分方程的欧拉法一阶ODE是最简单的一种ODE，微分方程本身是一个关于函数的方程，如果说这个ODE是显式给出的，那么这实际上给出的就是这个函数的导数的表达式，从数值计算的角度来看，数值方法关心的是整个区域内各个采样点（离散点）处的函数值，而此时显式是函数的导数，那么如果说给这个式子赋值进行计算的话，得到的将是每个采样点处函数的斜率值，这个场也称为斜率场。从斜率场来理解ODE以及Euler法是一个非常直观的。假设ODE表达式为$y'=ty + t^3$，这个式子是显式表达式，其是一个ODE，同时它也是一个函数导数的表达式，由于方程右侧中不仅含有函数值$y$还包含有自变量$t$，所以这个ODE是非自相关的，同理，如果导数表达式中不含有自变量，就称为自相关的。斜率场就是对这个显式进行可视化的结果，它反映的是各个点处斜率的大小和方向。假设$t$的区间是$[0,1]$，进行斜率场的可视化，这里对于$y$的范围为了方便取$[0,1]$，即结果是反映在$[0,1]$区间上的斜率场分布。import numpy as np import matplotlib.pyplot as plt t = np.linspace(0, 1, 11) # 步长0.1 y = np.linspace(0, 1, 11) tt, yy = np.meshgrid(t, y) dy = tt * yy + tt**3 E_t, E_y = np.gradient(dy) E_norm = np.hypot(E_t, E_y) fig, ax = plt.subplots() ax.quiver(tt, yy, E_t/E_norm, E_y/E_norm, color='tab:red') ax.set_title($''$) fig.savefig('slope.png', dpi=600) plt.show()斜率场是认识、研究ODE的一个强有力的工具，也是理解Euler法的基础。试想如果ODE是显式表达的，那么就可以确定出斜率场，它意味着场中各个点处的斜率是已知的，那么对于某一点来说，它的函数值如果已知，然后利用它的斜率，可以计算出函数值变化的增量，进而可以得到下一个点，即这一思路在斜率场中表现为从一点处沿着斜率方向移动到下一个，反复的进行这一过程，即可得到问题的解答，当然这一过程完成，是建立在给定初始值的情况下，也即ODE问题实际上就是一个IVP问题。Euler法就是利用斜率与步长的乘积来确定$y$的增量，利用迭代关系完成整个过程的计算。$$ w_{0} = y_0\\ w_{i + 1} = w_i + h \times df $$其中$h$是自变量的步长，即$\mathrm{d}t$，因为ODE是显式表达的，所以$df$是ODE本身。之前已经说明了Euler法计算的思路，以及斜率场在ODE的IVP问题中的角色。以前述给出的例子进行计算来验证。若设定初始值$y(0) = 0.6$，自变量区间仍是$[0,1]$，进行Euler法的数值实现。h = 0.1 tmin = 0 tmax = 1 t = np.arange(tmin, tmax+h, h) w = np.zeros(t.shape[0]) # 初始化一个结果向量 w[0] = 0.6 for i in np.arange(w.shape[0] - 1): df = t[i] * w[i] + t[i]**3 # 计算当前点的斜率值 w[i + 1] = w[i] + h * df # 计算下一个点的函数值 print(w)为了体现出Euler法求解IVP问题的核心思想，对于上述方程给定另外一个初始值作为对比、参照，$y(0)=0.4$。结果如下。从斜率场的角度来看，当ODE是一个显式表达时，意味着场中各个离散点处的斜率时已知的，利用函数微分的概念获得函数值的增量，由此得到下一个点处的函数值。而要想完成迭代计算，那么必须提供初始的函数值$y_0$，不同$y_0$意味着从不同的点出发开始移动。
梯度下降法的缺点之前的内容已经完全学习了梯度下降算法，对于这个方法之所以在线性方程组求解领域没有更广泛地流行起来，最大的原因是每一步的迭代中近似值移动的方向是沿着负梯度方向，而相邻迭代步的残差向量(本质就是负梯度向量)$r$之间是相互正交的，这个现象会导致可能在同一个方向上搜索（移动）了多次，造成迭代效率过低，当然迭代效率好与坏还受到初始值选择的影响，但是这个现象仍然不能够被忽略。为了说明这个现象，也为了后续方法的引入，这里还是用一个可视化的方式来反应这一现象，这也有助于对于后续共轭方向的理解。这里采用的这个$2 \times 2$系数矩阵的例子在前一部分的内容中已经出现过，之所以对这个二维解向量空间的方程组进行求解与分析，主要原因在于它的解向量是一个二维空间中的向量，对于可视化较为方便与明晰。方程组如下，由于这里不在牵扯到对应的二次型，因此对于未知量设定为$x_1, x_2$:$$ \begin{bmatrix} 2 & 2\\ 2 & 5 \end{bmatrix} \begin{bmatrix} x_1\\ x_2 \end{bmatrix} = \begin{bmatrix} 6\\ 3 \end{bmatrix} \tag{1} $$梯度下降法的实现函数如下，之前的内容已经给出。import numpy as np # 每一步的残差向量r0都会被以行向量的形式存放在r矩阵中 def gd_method(A, b, x0, tol, max_time): r0 = b - np.dot(A, x0) number = 0 r = np.zeros([max_time, b.shape[0]]) x = np.zeros([max_time, b.shape[0]]) x[0:, :] = x0 r[0, :] = r0 while np.linalg.norm(r0, 2) > tol: r0 = b - np.dot(A, x0) r[number, :] = r0 alpha = np.dot(r0.T, r0) / np.dot(np.dot(r0,A), r0) x0 = x0 + alpha * r0 x[number + 1, :] = x0 number += 1 if number > max_time: print('warning:the result is not stable!') return return x, number, r求解如下，这里初始值选择为$(0, 0)$，迭代了39次达到精度要求。这里函数返回的对象是一个元组，需要进行解包，第一个元素为x矩阵，存放着每一步的解向量，以行进行存储；第二个为迭代次数，第三为每一步中的残差向量，也是按照行进行存储。迭代39次后获得了一个可以被接受的近似解$(3.99999669, -0.99999779)$。为了比较好的理解线性方程组，我们依旧将其对应的二次型函数的等值线图画出来，为了直观地感受梯度，同时提供了这个函数的在每一个采样点处的斜率方向，即斜率场。A1 = np.array([[2, 2], [2, 5]]) b1 = np.array([6, 3]) x01 = np.array([0, 0]) x_gd1 = gd_method(A1,b1, x01, 0.00001,500) print(x_gd1[1]) print(x_gd1[0][0:x_gd1[1] + 1])39 [[ 0. 0. ] [ 1.42857143 0.71428571] [ 2.04081633 -0.51020408] [ 2.74052478 -0.16034985] [ 3.04039983 -0.76009996] [ 3.38311418 -0.58874279] [ 3.52999176 -0.88249794] [ 3.69785184 -0.7985679 ] [ 3.76979188 -0.94244797] [ 3.85200907 -0.90133938] [ 3.887245 -0.97181125] [ 3.92751464 -0.95167643] [ 3.94477306 -0.98619327] [ 3.96449697 -0.97633131] [ 3.97295007 -0.99323752] [ 3.98261076 -0.98840717] [ 3.98675106 -0.99668776] [ 3.99148282 -0.99432188] [ 3.99351072 -0.99837768] [ 3.99582832 -0.99721888] [ 3.99682158 -0.99920539] [ 3.99795673 -0.99863782] [ 3.99844322 -0.99961081] [ 3.99899921 -0.99933281] [ 3.9992375 -0.99980937] [ 3.99950982 -0.99967321] [ 3.99962653 -0.99990663] [ 3.99975991 -0.99983994] [ 3.99981708 -0.99995427] [ 3.99988241 -0.9999216 ] [ 3.9999104 -0.9999776 ] [ 3.9999424 -0.9999616 ] [ 3.99995612 -0.99998903] [ 3.99997179 -0.99998119] [ 3.99997851 -0.99999463] [ 3.99998618 -0.99999079] [ 3.99998947 -0.99999737] [ 3.99999323 -0.99999549] [ 3.99999484 -0.99999871] [ 3.99999669 -0.99999779]] import numpy as np import matplotlib.pyplot as plt x = np.linspace(0, 5, 100) y = np.linspace(-2, 2, 100) xx, yy = np.meshgrid(x, y) f_xy = xx * (xx + yy) - 6 * xx + yy * (xx + 2.5 * yy) - 3 * yy g_x, g_y = np.gradient(f_xy) g_xy = np.hypot(g_x, g_y) fig, ax = plt.subplots() cb = ax.pcolormesh(xx, yy, f_xy, cmap='Spectral_r') ax.contour(xx, yy, f_xy, colors='gray') ax.quiver(xx[::10, ::10], yy[::10, ::10], g_x[::10, ::10] / g_xy[::10, ::10], g_y[::10, ::10] / g_xy[::10, ::10]) ax.set_xlabel('$x_{1}$') ax.set_ylabel('$x_{2}$') ax.set_title('the function of Ax=b') fig.colorbar(cb, ax=ax) fig.savefig('1.jpg', dpi=600) plt.show()梯度下降法中，相邻残差向量是相互正交的，从图中可以看出，每一步的向量都与前后步的向量正交，在二维平面上正交就意味着垂直，那么相隔的残差向量必定是平行的，即共线。梯度下降法的迭代或者说是搜索效率低，因为在某一个方向进行了多次的重复搜索，使得一个很简单的二元方程组的计算都需要数十次迭代，虽然这个迭代次数会受到初值的影响，但是迭代效率低是一个不争的事实。根据这个图中搜索步的显示，实际上很法线不论搜索了多少步，会发现这些向量只有两个方向，因为相隔的都是共线的，即本质上只存在两个方向，只不过没一个方向上进行了很多次的搜索。试想一下，如果对这些向量移动下，让它们首尾相接，那么似乎好像可以只进行两步就可以达到最终的点，这样岂不是效率更高。这就是共轭方向法的构想。fig, ax = plt.subplots() cs = ax.contour(xx, yy, f_xy, colors='k') ax.plot(x_gd1[0][0:x_gd1[1]][:, 0], x_gd1[0][0:x_gd1[1]][:, 1], color='tab:red', marker='o') ax.clabel(cs, inline=True, colors='k', fontsize=14) ax.set_xlabel('$x_{1}$') ax.set_ylabel('$x_{2}$') fig.colorbar(cs, ax=ax) fig.savefig('2.jpg', dpi=600) plt.show()一个改进的方向--共轭方向依据之前的阐述，我们已经明确了一个思路就是，如果是一个$n \times n$的方程组，那么其解是一个$n$维向量，之前的例子是一个二维的，这里推广到更一般的情况。一个好的想法是找到n个“方向向量”，每一次的移动搜索都能完全的搞定一个未知数$x_{i}$，这样效率应该是最高的，那么只要$n$次后，所有的$x_{i}$都将被确定！这绝对是一个好的想法，当然按照这个思路来看，点移动的方向似乎不再是负梯度方向，而是方向向量$d_{i}$了，那么自然得到，下一个解可以表示为：$$ x_{i + 1} = x_{i} + \alpha_{i} d_{i} \tag{2} $$对于这个$2 \times 2$的方程组而言，初始值$x_0$向量沿着方向向量$d_0$移动，当然还需要乘以一个步长才行，获得第一步迭代结果$x_1$向量，依据前述的原理，解向量$x$中的第一个元素也就是二维坐标轴的第一个轴的值得到确定，此时的向量$x_1$与精确解$x$之间的距离，我们称之为误差向量，$e_1 = x_1 - x$，然后从这个点按照$d_1$移动，直接得到精确解，这样一共两步就可以完成计算。按照上述的过程，可以确定一个关系，方向向量$d_i$与误差向量$e_{i + 1}$保持正交！$$ d_{i}^{\mathrm{T}} \cdot e_{i + 1} = 0\\ e_{i + 1} = x_{i + 1} - x\\ d_{i}^{\mathrm{T}} \cdot (e_{i} + \alpha_{i} d_{i}) = 0\\ \alpha_{i} = - \frac{d_{i}^{\mathrm{T}} e_{i}}{d_{i}^{\mathrm{T}} d_{i}} \tag{3} $$由于精确解本身是未知的，因此，误差向量其实也是未知的，因此，只是依据这个式子，还无法确定出步长,这里采用另一种正交方式。$$ \alpha_{i} = \frac{r_{i}^{\mathrm{T}} r_{i}}{d_{i}^{\mathrm{T}} A d_{i}} \tag{4} $$这个式子的推导利用了，矩阵与向量之间的正交！并且，依据误差向量与残差向量的相互关系$r_{i} = - A e_{i}$得出。共轭梯度法因此现在只要能找到这样的一组相互正交的方向向量$d_{i}$，那么实际上这个问题似乎就可以得到解决！而这一过程通过Gram-Schmidt过程来完成，而梯度下降法中残差向量$r_{i}$中相邻的向量互为正交，那么可以考虑使用残差作为这一组基，在此基础上进行共轭。用残差向量产生共轭后方向向量，在$n$步之后，即可得到结果。$$ r_{i + 1} = - A e_{i + 1}\\ = -A (e_{i} + \alpha_{i} d_{i})\\ = r_{i} - \alpha_{i} A d_{i} \tag{5} $$从残差向量构造方向向量，即$$ d_{i} = r_{i} + \sum_{k=0} ^{i - 1} \beta_{ik}d_{k} \quad (i > k) \tag{6} $$两侧构造矩阵向量正交,可以得到$\beta_{ij}$$$ \beta_{ij} = -\frac{r_{i}^{\mathrm{T}} A d_{j}}{d_{j}^{\mathrm{T}} A d_{j}} \tag{7} $$因为式(5)成立，可以继续进行化简，依据矩阵向量的正交性，$j$实际上只能取固定值，就可以得到关于$\beta_{ij}$的一个公式，这样的完全可以只写为$\beta_i$.$$ \beta_{i} = \frac{r_{i + 1}^{\mathrm{T}} r_{i + 1}}{r_{i}^{\mathrm{T}} r_{i}} \tag{8} $$这样，共轭梯度法的涉及到的计算过程及公式如下：$$ d_{0} = r_{0} = b - Ax_{0}\\ \alpha_{i} = \frac{r_{i}^{\mathrm{T}} r_{i}}{d_{i}^{\mathrm{T}} A d_{i}}\\ x_{i + 1} = x_{i} + \alpha_{i} d_{i}\\ r_{i + 1} = r_{i} - \alpha_{i} A d_{i}\\ \beta_{i} = \frac{r_{i + 1}^{\mathrm{T}} r_{i + 1}}{r_{i}^{\mathrm{T}} r_{i}}\\ d_{i + 1} = r_{i + 1} + \beta_{i + 1} d_{i} \tag{9} $$以上就是共轭梯度的全部过程和原理，以及实现的思想，如果你对这个过程完全了解，那么一定对这个方法很惊叹！因为，从迭代法的角度来看，它似乎有着超高的效率，有几个未知数就只需要迭代几次，极大地改进了梯度下降法的缺陷，对于一个迭代法而言，传统的迭代法都是获得的近似值，但是CG法确实通过有限步的迭代获得是精确解，最为重要的是它的迭代次数是完全确定的，这与传统的迭代法有着很大的区别，从这一个角度来看，CG法应该是一种直接法。从直接法的角度来看，它又不是完全意义上的直接法，因为该方法中需要给定的一个初始值才能求解，这一点又与迭代法有点相似。对于大型稀疏的线性方程系统来说，CG法是最好的方法！def cg_method(A, b, x0): r0 = b - np.dot(A,x0) d0 = r0 x = np.zeros([b.shape[0] + 1, b.shape[0]]) x[0,:] = x0 for i in np.arange(b.shape[0]): r0 = b - np.dot(A,x0) alpha = np.dot(r0.T, r0) / np.dot(np.dot(d0.T, A), d0) x0 = x0 + alpha * d0 x[i + 1,:] = x0 r1 = r0 - alpha * np.dot(A, d0) beta = np.dot(r1, r1) / np.dot(r0, r0) d0 = r1 + beta * d0 return xxl=cg_method(A1,b1, x01)fig, ax = plt.subplots() cs = ax.contour(xx, yy, f_xy, colors='k') ax.plot(xl[:, 0], xl[:, 1], marker='*') ax.clabel(cs, inline=True, colors='k', fontsize=14) ax.set_xlabel('$x_{1}$') ax.set_ylabel('$x_{2}$') fig.colorbar(cs, ax=ax) fig.savefig('3.jpg', dpi=600) plt.show()
共轭梯度法(1)--的基本原理之前已经搞明白了，梯度下降法的基本原理，当然解释的调度是从求函数极值的角度出发的，事实上从这个角度来理解，个人感觉是一个最为直接的理解角度，其完完全全是建立在多变量函数的微分系统中的。事实上这个方法(思想)在实际中是被用于求解线性方程的，当然单纯的梯度下降形式并没有被直接采用，被广泛用于求解对称、正系数矩阵方程组的方法是，基于梯度下降原理实现的共轭梯度法，这一方法在大型稀疏线性方程组的求解中，有极高的效率。但是，之前对于梯度下降思想的描述，好像认为这个方法是用于求函数极值的方法，似乎与线性方程组求解问题无关。因此，对于对称、正定系数矩阵的二次型的理解是很必要的，它能让我们明白矩阵系统中线性方程组与函数极值的关系，能够解释为什么一个求函数极值的思想可以被用于求解线性方程组系统。$Ax = b$对应的二次型这里就不严格去给出严格的数学表示了，因为对于非数学专业的我们来说，好像没啥必要，知晓核心的思想更为重要。注意本次针对的系数矩阵$A$都是对称、正定的，线性方程组系统$Ax = b$对应的是二次型对于自变量导数为0时的表达式，不严谨地可将其看作是函数吧，即二次型为$$ f(x) = \frac{1}{2} x^{T} A x - b^{T} x + c \tag{1} $$对于$Ax = b$来说，$c$是一个0向量。我们知道求$f(x)$极小值的方法最直接的一种就是利用其导数为0，获得极值点从而获得极小值，因此，求导得到$$ f'(x) = \frac{1}{2}A^T x + \frac{1}{2}A x - b = 0 \tag{2} $$因此，求解线性方程组$Ax = b$就等价于求解$f'(x) = 0$这个表达式，那么一阶导数为0这个表达换一种说法，不就是去求二次型$f(x)$的函数的极值嘛！这也就解释了为什么一个求解函数极值的方法却可以应用到线性方程组的求解当中来。为了更好的理解线性方程组$Ax = b$的求解与函数极值之间的关系，这里举一个简单的，可以被可视化的例子，若方程组为：$$ \begin{bmatrix} 2 & 2\\ 2 & 5 \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} 6\\ 3 \end{bmatrix} \tag{3} $$这样一个对称正定的方程组，为了待求的未知数为了更符合函数表达，这里用了$x, y$来表示，依据二次型公式，得到$$ f(x, y) = \frac{1}{2} \begin{bmatrix} x\\ y \end{bmatrix} \begin{bmatrix} 2 & 2\\ 2 & 5 \end{bmatrix} \begin{bmatrix} x & y \end{bmatrix} - \begin{bmatrix} 6 & 3 \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} $$这样一个$2 \times 2$的线性方程组对应的二次型就是一个二元函数，为了得到表达式，这里利用sympy库进行符号运算。import numpy as np from sympy import * # 先定义符号变量x, y x, y = symbols('x y') A = np.array([[2, 2], [2, 5]]) b = np.array([6, 3]) # 自变量向量这里用x1表示 x1 = np.array([x, y]) f_xy = 1/2 * x1.T @ A @ x1 - b.T @ x1 init_printing(use_unicode = True) # 更美观的显示 pprint(f_xy.expand()) 2 2 1.0⋅x + 2.0⋅x⋅y - 6⋅x + 2.5⋅y - 3⋅y 利用函数图像来直观地可视化函数极值情况。import matplotlib.pyplot as plt import numpy as np from mpl_toolkits.mplot3d import Axes3D x = np.linspace(-10, 10, 100) y = np.linspace(-10, 10, 100) xx, yy = np.meshgrid(x, y) f_xy = xx * (xx + yy) - 6 * xx + yy * (xx + 2.5 * yy) - 3 * yy # 这里用等值线云图与伪色彩图来可视化，也顺便对比下两者的区别 fig, ax = plt.subplots(1, 2) cf1 = ax[0].pcolormesh(xx, yy, f_xy, cmap='RdBu_r') cs1 = ax[0].contour(x, y, f_xy, colors='gray', levels=15) ax[0].clabel(cs1, inline=True, colors='k') fig.colorbar(cf1, orientation='vertical') ax[0].set_title('function by pcolormesh') cf2 = ax[1].contourf(x, y, f_xy, cmap='RdBu_r') cs2 = ax[1].contour(x, y, f_xy, colors='gray', levels=15) ax[1].clabel(cs2, inline=True, colors='k') ax[1].set_title('function by contourf') fig.savefig('function2.jpg', dpi=600) plt.show() # 三维的函数图像 fig1 = plt.figure() ax = Axes3D(fig1) cd = ax.plot_surface(xx, yy, f_xy, cmap='RdBu_r') ax.set_title('the function of 3D image') fig1.colorbar(cd) fig1.savefig('image3d.jpg', dpi=600) plt.show()通过函数图像的可视化，可以看到线性方程组本身是对称正定的，那么对应的二次型的函数图像则必定都位于在$z= 0 $平面以上，函数均是大于0的，且从函数的等值线及三维图像中可以看到是存在极小值的。明白了这个关系之后，对于梯度下降法或者共轭梯度法应用于方程组(系数矩阵为对称正定)的求解就不足为奇了。Gradient descent method求解方程组共轭梯度法的实现全来自于或者说是基于梯度下降法，因此，理解梯度下降法求解线性方程系统很有必要。在之前的梯度下降法的原理中，已经明白了这个方法对于函数极值的求解，是通过沿着梯度方向负向进行的，每一次的计算实质上获得是自变量的增量向量，而现在面对方程组时，由于之前的二次型的内容已经得出了线性方程组与多元函数极值之间的关系，即$f'(x_{i}) = Ax_{i} - b = 0$。那么梯度下降法可以被这样描述与构建：每一次产生的近似值$x_{i}$向量，近似解自然会有一个残差向量(residual vector)为$r_{i}$， 下表$i$表示第几步的迭代。$$ r_{i} = b - A x_{i}\\ r_{i} = - f'(x_{i}) \tag{4} $$式（4）中的$r_{i}$不就是梯度向量嘛，表示了点下一步搜索或者移动的方向，当然要想确定自变量的增量大小，还需要知道步长，这一点在梯度下降原理中已经给出了解释；在这里步长设定为$\alpha_{i}$，所以，产生的新的点$x_{i + 1}$为：$$ x_{i + 1} = x_{i} + \alpha_{i} r_{i} \tag{5} $$到这里好像似乎与之前的原理中阐述并无大的区别，实际上现在问题在于如何确定步长$\alpha_i$，而在之前的介绍中步长是一个固定的，对于迭代计算来说，固定的步长的计算效率显然是不能够令人满意的。对于函数的极小值的确定来说，在这里每一步获得的新的函数值为$f(x_{i+1})$，对于整个求解过来说，只要这个更新的函数值达到最小值那么任务就算是完成了，依据式(5)，函数值可以看做是关于步长的函数，为了获得极小值，通过导数为0来完成，即$$ \frac{\mathrm{d}f'(x_{i + 1})}{\mathrm{d} \alpha_{i}} = 0 \tag{6} $$结合式(5)来看，这是一个复合函数，利用链式求导法则，得到$$ 式(6) = f'(x_{i+1}) ^{\mathrm{T}} \cdot \frac{\mathrm{d}x_{i + 1}}{\mathrm{d}\alpha_{i}} = 0 \tag{7} $$这个式子中，$f'(x)$不就是梯度向量嘛，依据式(4)给出的定义，在这里它也称作残差向量$-r_{i + 1}$，$x_{i + 1}$求导的结果可以根据式(5)为$r_{i}$，所以，要保证函数取得极值，就要使得前后两步的梯度向量（也就是残差向量$r$）保证正交！这里写为内积形式：$$ r_{i + 1} ^ {\mathrm{T}} \cdot r_{i} = 0 \tag{8} $$那么依据这个关系，带入残差向量的表达，得到：$$ (b- A x_{i + 1}) ^ {\mathrm{T}} \cdot r_{i} = 0\\ (b - A (x + \alpha_{i} r_{i}))^{\mathrm{T}} \cdot r_{i} = 0 \tag{9} $$通过这个式子就可以得到步长$\alpha_{i}$的计算公式：$$ \alpha_{i} = \frac{r_{i} ^{\mathrm{T}} r_{i}}{r_{i} ^{\mathrm{T}} A r_{i}} \tag{10} $$到这里，所有量的表达都得到了解决，那么对于线性方程组系统来说，梯度下降法的迭代过程为：$$ repeat:\\ r_{i} = b - A x_{i}\\ \alpha_{i} = \frac{r_{i}^{\mathrm{T}} r_{i}}{r_{i} ^{\mathrm{T}} A r_{i}}\\ x_{i + 1} = x_{i} + \alpha_{i} r_{i}\\ until: x_{i}的增量幅度满足设定精度\\ end \tag{11} $$迭代计算的代码如下：import numpy as np # 每一步的残差向量r0都会被以行向量的形式存放在r矩阵中 def gd_method(A, b, x0, tol, max_time): r0 = b - np.dot(A, x0) number = 0 r = np.zeros([max_time, b.shape[0]]) x = np.zeros([max_time, b.shape[0]]) while np.linalg.norm(r0, np.inf) > tol: r0 = b - np.dot(A, x0) r[number, :] = r0 alpha = np.dot(r0.T, r0) / np.dot(np.dot(r0,A), r0) x0 = x0 + alpha * r0 x[number, :] = x0 number += 1 if number > max_time: print('warning:the result is not stable!') return return x, number, r这里对一个$3\times3$大小的对称正定矩阵进行验算，梯度下降算法迭代了67次得到结果，$[0.99983945,0.99976565, 1.99978575]$，而这个方程组的精确解为$[1,1,2]$。可以迭代结果还是比较可信与可靠的。A = np.array([[4, -2, -1], [-2, 4, -2], [-1, -2, 3]]) b = np.array([0, -2, 3]) x0 = np.array([1, 1, 1]) x_gd = gd_method(A,b, x0, 0.0001,500) print(x_gd[0][x_gd[1]-1,:]) print(x_gd[1]) print(x_gd[2]) [0.99983945 0.99976565 1.99978575] 67 [[-1. -2. 3. ] [-0.39130435 0.43478261 0.15942029] [ 0.09201888 0.02436289 0.15942029] ... [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ]]x = np.linalg.solve(A,b) print(x) [1. 1. 2.]依据上述的分析以及这个算例的计算，实际上对于线性方程组迭代的迭代计算会有一个新的认识。首先，从最为原始的迭代的角度来看，通过每一迭代步中都会产生一个残差向量，似乎在方程组的求解角度来看，残差向量仅仅也就表示一个误差吧，好像也并没有很直观的一个认识。但是从梯度下降的算法角度来看，线性方程组迭代步中的残差向量实际上就是二次型（函数）的负梯度向量，那么也就是说，这个向量$r$不断地在修正点移动的方向，使得移动方向不断的向函数的极小值点处靠近！由于gd_method函数中会记录每一步的参差向量，因为这个是一个3个未知数的方程组，因此残差向量$r$本身就是一个$R^{3}$空间中的向量，同时返回得到的$x$矩阵中存放着每一步的近似解，这样进行可视化，可以看出梯度下降中前后的参差向量确实是正交的。因为三元函数无法进行可视化，不是很容易看出梯度下降的搜寻过程。可能二元函数的化对于理解更方便吧。xs = x_gd[0][0:x_gd[1]-1,:] plt.style.use('ggplot') fig = plt.figure() ax = Axes3D(fig) ax.plot(xs[:,0], xs[:,1], xs[:,2], marker='o') ax.set_title('point in every step') ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') fig.savefig('gdimage.jpg', dpi=600) plt.show()在上述的分析中，我们会发现实际好像似乎梯度下降这个迭代算法的效率并不是很高，很简单的方程组都要迭代67次，显然这样的效率去求解大型的稀疏方程组是不能令人满意的。这个方法之所以没有被广泛地应用，主要原因是因为每一步的残差向量也就是梯度向量，都必须与相邻步的残差向量保持正交，如果为了便于理解，以二维空间为例，就是每一次的搜索(移动)方向都是相互垂直的，那么间隔的残差向量必定平行，也就说从初始值向精确解移动是呈折线的方式，一个搜索方向会被重复了好多次，这样使得迭代效率并不高。但是梯度下降法实现的思想确是很有意义的，为函数极值的确定与线性方程组的求解之间建立了联系。真正有意义的方法是在这个方法基础上实现的共轭梯度法，它对搜索方向进行了共轭处理,在理解了梯度下降的原理后，再对共轭梯度法进行学习可能更好一点。
梯度下降法的原理梯度下降法是一种数值计算方法，但是在数值分析中往往并没有直接采用它，但是在这一方法基础上发展而来的共轭梯度法却被得到应用，而梯度下降法是其核心思想与方法。目前，较为火热的深度学习方面，由于广泛采用了自动求导原理，对损失函数求极值时，梯度下降算法在该领域得到了十分广泛的应用。而我们感兴趣的是利用共轭梯度法求线性方程组，该方法的思想就是梯度下降原理，因此，很有必要对该梯度下降法有一个认识。什么是梯度下降在微积分中，很常见的一个公式就是泰勒展开，这个公式意在描述函数$f(x)$的值与$f(x+\Delta x)$之间的关系，更直白一点的说，研究的是应变量变化的幅度。对于一元函数而言，利用泰勒展开可以方便的获得函数值的变化。$$ f(x + \Delta x) = f(x) + f(x)' \cdot \Delta x + \frac{1}{2!} f(x)'' \cdot \Delta x ^2 $$当我们不考虑二阶以上的项时，便得到$f(x + \Delta x)$的一个近似值。这个式子其几何意义在很多数值方法中被应用，比如在ODE方程的数值法中，欧拉方法本质上就是采用了这个公式去完成函数值的求解，利用该点的函数值、该点的斜率以及步长，即可得到下一个点的函数值。在实际应用中往往是多元函数，尤其在深度学习中，神经网络每一层基本都是多个神经元，往往都是多元函数问题。从简单一点的情况来看，即是二元函数，此时泰勒公式的描述似乎有点繁琐，不是那么容易理解，但是从全微分的概念来看却十分清晰。二元函数的全微分公式为：$$ \mathrm{d}f(x, y) = \frac{\partial f}{\partial x} \mathrm{d}x + \frac{\partial f}{\partial y} \mathrm{d}y $$其中，$\mathrm{d}f(x, y)$表示函数的变化量，即$f(x + \Delta x, y + \Delta y) - f(x, y)$，这样很容易可以得到二元函数的泰勒展开形式。当然，这里对函数值变化量更感兴趣。对于上式的表达可以写为两个向量的内积，这样可以很直观地从几何角度来理解什么是梯度下降。函数的导数向量记为$a$$$ a = (\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}) $$而自变量的增量向量为$b$$$ b = (\mathrm{d}x, \mathrm{d}y) $$这两个向量的内积就表示函数的变化量的大小，而从向量内积的定义和几何两个方面来看，$$ a \cdot b = |a||b|\mathrm{cos} \theta $$当这两个向量平行且反向的时候，夹角为$\pi$，内积取得最小值，因为此时$\cos \pi = -1$，那么，从这个角度来看，只要这两个向量保持共线且反向，函数的变化量最小$$ b = - k a \quad (k是一个正数) $$$$ (\mathrm{d}x,\mathrm{d}y) = - k (\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}) $$这样，通过这个式子就可得到函数沿着梯度方向的一个自变量变化值。因为沿着这个方向是下降最快的一个路径，此时的$k$可以看做是步长，所以该方法也称为最速下降法。上式虽然是以二元函数为例的，但是毫无疑问的是可以推广到$n$元函数中去。更为一个通用的表达是$$ \nabla f = - k \Delta x_{i} $$Example比如函数$f(x, y) = x^2 + y ^ 2$，当然这个函数极小值很明显就是0，即当$x= 0, y=0$时取得最小值。import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 函数图像 x = np.linspace(-2, 2, 50) y = np.linspace(-2, 2, 50) xx, yy = np.meshgrid(x, y) z = xx**2 + yy**2 fig = plt.figure() ax = fig.add_subplot(111, projection='3d') cb = ax.plot_surface(xx, yy, z, cmap='RdBu_r', alpha=0.7) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') fig.colorbar(cb, ax=ax) ax.set_title('$z = x^2 + y^2$') fig.savefig('function image.jpg', dpi=600) plt.show()此时，任意选择一个初始的点,如点$(1, 1)$,即$x_0 = 1, x_1 = 1$，函数的梯度向量为$$ \nabla f = (2x, 2y) $$带入初始值后，得到$\nabla f = (2, 2)$，步长k的取值如果过大，可能会使得自变量增量过大，导致跨越过最小值点，因此，适当的较小的$k$值是获得正确解答的必要条件，这里取$k=0.1$，由此得到自变量变化量为$(-0.2, -0.2)$，即下一个新的点为$(1 - 0.2, 1 - 0.2)$，这一算法的意义为，通过梯度向量计算得到的自变量变化量，就是函数下降最快的方向，以此反复计算，知道前后计算步变化误差满足精度要求即可。对于这个函数，$x_0 = 1, y_0 = 1$，$x_1 = 0.8, y_1 = 0.8$, $x_2 = 0.64, y_2 = 0.64$，$x_3 = 0.516, y_3 = 0.516$import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 函数图像 x = np.linspace(-2, 2, 50) y = np.linspace(-2, 2, 50) xx, yy = np.meshgrid(x, y) z = xx**2 + yy**2 fig = plt.figure() ax = fig.add_subplot(111, projection='3d') cb = ax.plot_surface(xx, yy, z, cmap='RdBu_r', alpha=0.7) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') fig.colorbar(cb, ax=ax) ax.set_title('$z = x^2 + y^2$') ax.plot([1, 0.8, 0.64, 0.516], [1, 0.8, 0.64, 0.516], [2, 1.28, 0.8192, 0.532512], color='k', linewidth=2) fig.savefig('gradient.jpg', dpi=600) plt.show()import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 函数图像 x = np.linspace(-2, 2, 50) y = np.linspace(-2, 2, 50) xx, yy = np.meshgrid(x, y) z = xx**2 + yy**2 plt.style.use('ggplot') fig,ax = plt.subplots(1,1) cs = ax.contour(xx, yy, z, cmap='viridis') ax.clabel(cs, inline=True, fontsize=12) ax.set_xlabel('x') ax.set_ylabel('y') fig.colorbar(cb, ax=ax) ax.set_title('$z = x^2 + y^2$') ax.plot([1, 0.8, 0.64, 0.516], [1, 0.8, 0.64, 0.516], linewidth=2) fig.savefig('contour.jpg', dpi=600) plt.show()补充一点这里为了显示出梯度下降的过程，对之前那个函数进行编程实现，获得每一步的自变量、应变量值，并进行可视化。import numpy as np f = lambda x,y: x**2 + y**2 x0 = 4 y0 = 5 fx = lambda x: 2 * x fy = lambda y: 2 * y step = 0.1 max_time = 400 tol = 0.00001 x = np.zeros(max_time) y = np.zeros(max_time) z = np.zeros(max_time) number = 0 while f(x0, y0) > tol: dx, dy = np.array([fx(x0), fy(y0)]) x0 = x0 - dx * step y0 = y0 - dy * step z0 = f(x0, y0) x[number] = x0 y[number] = y0 z[number] = z0 number += 1 if number > max_time: break print(number) print(x[0:number]) print(y[0:number]) print(z[0:number])结果为35 [3.20000000e+00 2.56000000e+00 2.04800000e+00 1.63840000e+00 1.31072000e+00 1.04857600e+00 8.38860800e-01 6.71088640e-01 5.36870912e-01 4.29496730e-01 3.43597384e-01 2.74877907e-01 2.19902326e-01 1.75921860e-01 1.40737488e-01 1.12589991e-01 9.00719925e-02 7.20575940e-02 5.76460752e-02 4.61168602e-02 3.68934881e-02 2.95147905e-02 2.36118324e-02 1.88894659e-02 1.51115727e-02 1.20892582e-02 9.67140656e-03 7.73712525e-03 6.18970020e-03 4.95176016e-03 3.96140813e-03 3.16912650e-03 2.53530120e-03 2.02824096e-03 1.62259277e-03] [4.00000000e+00 3.20000000e+00 2.56000000e+00 2.04800000e+00 1.63840000e+00 1.31072000e+00 1.04857600e+00 8.38860800e-01 6.71088640e-01 5.36870912e-01 4.29496730e-01 3.43597384e-01 2.74877907e-01 2.19902326e-01 1.75921860e-01 1.40737488e-01 1.12589991e-01 9.00719925e-02 7.20575940e-02 5.76460752e-02 4.61168602e-02 3.68934881e-02 2.95147905e-02 2.36118324e-02 1.88894659e-02 1.51115727e-02 1.20892582e-02 9.67140656e-03 7.73712525e-03 6.18970020e-03 4.95176016e-03 3.96140813e-03 3.16912650e-03 2.53530120e-03 2.02824096e-03] [2.62400000e+01 1.67936000e+01 1.07479040e+01 6.87865856e+00 4.40234148e+00 2.81749855e+00 1.80319907e+00 1.15404740e+00 7.38590339e-01 4.72697817e-01 3.02526603e-01 1.93617026e-01 1.23914897e-01 7.93055338e-02 5.07555416e-02 3.24835466e-02 2.07894698e-02 1.33052607e-02 8.51536685e-03 5.44983478e-03 3.48789426e-03 2.23225233e-03 1.42864149e-03 9.14330553e-04 5.85171554e-04 3.74509795e-04 2.39686269e-04 1.53399212e-04 9.81754956e-05 6.28323172e-05 4.02126830e-05 2.57361171e-05 1.64711150e-05 1.05415136e-05 6.74656869e-06]可视化编程： from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt xi = np.linspace(-6, 6, 50) yi = np.linspace(-6, 6, 50) xx, yy = np.meshgrid(xi, yi) zi = xx**2 + yy**2 fig = plt.figure() ax = Axes3D(fig) ax.plot_surface(xx, yy, zi, cmap='RdBu_r', alpha=0.5) ax.plot(x[0:number], y[0:number], z[0:number], color='tab:blue') fig.savefig('demo.jpg',dpi=600) plt.show()
Numerical Analysis: Iteration Method of Equation Systems
对角占优线性方程组的迭代法线性方程组的求解是多数数值问题最终的一个环节，很多问题最终都变为了一个求解线性方程组的问题。求解方程组最广为人知的方法是高斯消元法，这是一种直接方法，通过构造某种形式，可以直接得到最终的结果，其矩阵形式就对应着LU分解法。现代的计算机语言对于科学计算的解决，其已经包含了直接法，并且实现的效率极高，而迭代法计算效率低下，一般计算中往往不会采用它。但是，迭代法求解问题的思路却是十分有价值的，迭代思想也不仅仅局限于方程组的求解。迭代的直观感受以一个二元一次方程组为例，$$ a_{11}x_{1} + a_{12}x_{2} = b_{1}\\ a_{21}x_{1} + a_{22}x_{2} = b_{2} $$从方程组的第一式得到：$$ x_{1} = \frac{b_{1} - a_{12}x_{2}}{a_{11}} $$第二式得到：$$ x_{2} = \frac{b_{2} - a_{21}x_{1}}{a_{21}} $$以上构造的关系式，再附加上一组初始值，通过反复的计算就可以得到方程组的近似解。上述的迭代公式的构造称为不动点迭代，该院里在方程求解中有过应用，这里该方法也称为雅克比法，Jacobian法。Jacobian法前述的实例可以看出，jacobian法本质上就是对每一个方程进行变换，依次得到一个$x_i$的表达式，然后依据初始值就可以玩车呢个迭代计算，且这一方法中，每一个式均是对应$b$减去其他项的和，然后除以$x_i$的系数。这样可以方便地得到其矩阵形式。方程组为$Ax = b$，对系数矩阵$A$进行分解，将其分解为对角线一下元素组成的矩阵$L$，只有对角线元素的矩阵$D$，以及对角线元素以上的元素组成的矩阵$U$，这样矩阵$A = L + D + U$，由此得到：$$ (L + D + U) x = b\\ Dx = b - (L + U) x\\ x = D^{-1} b - D^{-1}(L + U)x $$得到一个迭代公式：$$ X_{k + 1} = D^{-1} b - D^{-1}(L + U) X_{k} \quad (k = 0,1,2,3,....) $$这样第一次的给定的初始值记做$X_{0}$，利用它可以得到一组近似值为$X_{1}$，反复迭代直到精度达到事先设定的程度。这一方法可以很方便地利用代码实现。import numpy as np # 定义一个函数来完成jacobian法 # 输入参数包括方程组本身、初始值、精度 # 这里为了避免陷入死循环，还需要一个控制最大迭代次数 def Jacobian_Method(A, b, x0, eps=0.0001, max_times=400): L = np.tril(A, -1) U = np.triu(A, 1) D = np.diag(np.diag(A, 0)) number= 0 # 记录迭代次数 while np.linalg.norm(np.dot(A, x0) - b, np.inf) > eps: F = np.linalg.solve(D, b) G = np.linalg.solve(D, L + U) x0 = F - np.dot(G, x0) number += 1 #此时，迭代完成一次，则次数加一 if number > max_times: print('Warning! The solution is not stable!') return return x0, number这是雅克比法的数值实现代码。该方法的收敛速度并不快，并且并非对任意的线性方程组都能够进行求解，迭代法的适用范围主要是针对对角占优矩阵的。Gauss-Seidel迭代从Jcobian方法中可以看出的是，每一次的计算中，第一个未知数必定会被先得到确定，那么接下来第二个未知数的计算中人然会用第一个未知数，只不过对于Jcobian来说，它的计算依据是采用的上一步的值，而Gauss-Seidel法采用了更新的值用于后面的未知数的计算，这样可以大大提高收敛的速度。其矩阵形式为：$$ (L + D + U)x = b\\ （L + D）x = b - Ux\\ x = (L + D)^{-1}B - (L + D)^{-1} U x $$迭代公式为：$$ X_{k + 1} = (L + D)^{-1} b - (L + D)^ {-1} U X_{k} $$def Gauss_Seidel_Method(A, b, x0, eps=0.0001, max_times=400): L = np.tril(A, -1) U = np.triu(A, 1) D = np.diag(np.diag(A, 0)) number= 0 # 记录迭代次数 while np.linalg.norm(np.dot(A, x0) - b, np.inf) > eps: F = np.linalg.solve(L + D, b) G = np.linalg.solve(L + D, U) x0 = F - np.dot(G, x0) number += 1 if number > max_times: print('Warning: The result is not stable!') return return x0, numberGauss-Seidel迭代法的收敛速度也明显快于Jacobian法，所以这一方法使用的更为广泛。但是，有没有收敛速度更快的方法呢？这就是连续松弛迭代法，它本质是对Gauss-Seidel法的一个改进，通过松弛参数$\omega$。连续松弛迭代法SOR在gauss-seidel法的基础上，通过对$Ax=b$两边乘以一个松弛参数$\omega$，再进行迭代公式的构建。$$ \omega A x = \omega b\\ \omega (L + D + U) X = \omega b\\ \omega Lx + \omega Dx + \omega Ux + Dx = \omega b + Dx\\ (\omega L + D)x = \omega b + Dx - \omega D x - \omega U x\\ x = (\omega L + D)^{-1}[(1 - \omega)D - \omega U]x + \omega (\omega L + D)^{-1} b $$得到迭代公式$$ X_{k+1} = (\omega L + D)^{-1}[(1 - \omega)D - \omega U]X_{k} + \omega (\omega L + D)^{-1} b $$这一方法的收敛速度与$\omega$有关，当大于1时为过松弛迭代，小于1时为欠松弛迭代，等于1是就完全退化为gauss-seidel法。而且对于松弛参数的确定并没有一个确定的方法，需要通过不断的试算采纳呢个估计到一个最优的值。def SOR_Method(A, b, x0, w, eps=0.0001, max_times=400): L = np.tril(A, -1) U = np.triu(A, 1) D = np.diag(np.diag(A, 0)) number= 0 # 记录迭代次数 while np.linalg.norm(np.dot(A, x0) - b, np.inf) > eps: F = w * np.linalg.solve((w * L + D), b) G = np.linalg.solve(w * L + D, (1 - w) * D - w * U) x0 = F + np.dot(G, x0) number += 1 if number > max_times: print('Warning: The result is not stable!') return return x0, number收敛条件与效率对比上述的迭代法是否一定收敛？这个会收到系数矩阵$A$的影响，如果系数矩阵一个严格的对角占优矩阵，那么上述的三个迭代法一定收敛，但是收敛并不意味着这个矩阵一定是严格的对角占优矩阵，严格对角占优只是收敛的充分条件。对角占优表示每一行的主元比该行其他所有元素之和的绝对值还要大，这就是对角占优。在线性方程组的求解中，迭代法现在一般不会被使用，主要原因是迭代法的计算效率差，没有直接法快速，计算机硬件的普及，使得较大规模的直接发求解很容易实现，并且数值计算语言中也已经设置好了算法，matlab中的"\"就是直接法，这是一个高级的LU分解法的实现，并且利用底层的c++进行了预编译，速度十分快，而自己定义的迭代法在没有编译的情况下，速度慢大约近100倍，在python中，numpy库实现了线性代数的基本运算，这里的solve方法相当于matlab的"\"运算。A = np.array([[3,1,-1],[2,4,1],[-1,2,5]]) b= np.array([4,1,1]) x0 = [0,0,0] w = 1.1 a1 = Jacobian_Method(A, b, x0) a2 = Gauss_Method(A, b, x0) a3 = SOR_Method(A, b, x0, w) print(a1) print(a2) print(a3)结果如下：(array([ 1.99995 , -0.99994562, 0.99995301]), 26) (array([ 1.99996414, -0.99997146, 0.99998141]), 13) (array([ 1.9999659 , -0.99997676, 0.9999868 ]), 9)雅克比法迭代了26次才达到精度，高斯赛德尔法用了13次，SOR中，当松弛参数为1.1时，连续过松弛方法的收敛用了9次。
Numerical Analysis: Numerical Integration
Newton-Cotes积分对于定积分的计算，解析解的存在特别受制于被积函数的形式，往往原函数难以找到，使得积分计算无法进行。积分的意义是对所在区域进行无穷小划分并求和，因此，从数学上来说，即使是找不到原函数，但是积分表示的物理意义依然存在。对于积分的数值方法，Newton-Cotes法是一种最容易理解的方法，这一方法实际上是一系列的该类型数值方法的总称，因为这些方法的数值积分实现的思路是一样的。数值方法对于经典微积分的实现，一个解决的思路在于寻求近似解以替换精确解，其实在定积分的定义过程中，黎曼的定义就已经给出了这种寻求近似的思路，只不过他通过无穷小求和的形式给出了定积分的精确定义。而Newton-Cotes法就是依据这一种切割成小块然后求和的思想来实现数值解答的，其中最为关键的就是利用多项式插值来近似表示小区段内原始的函数。我们知道，这里采用的拉格朗日插值多项式，可以依据所给的插值点个数，来获得n-1次的多项式；对于2个点的区段来说，插值时只能得到一个一次多项式，这是线性插值，不管原函数在该区间上的形状是什么，都用一个一次多项式来替代；而当区间内取3个点时，拉格朗日多项式为2次，，当然也可以利用4点进行插值近似。通常，当采用两点插值的一次多项式来替换原函数时，被称为梯形法(trapz)；采用三点拉格朗日插值时，被称为Simpson法；4点时称为Simpson 3/8法。Newton-Cotes法的实现思路示意图。所以，对于Trapz法来说，在$[x_0, x1]$这个小区间上，采用两点拉格朗日插值公式得到线性函数$f_1(x)$，那么近似解就可以表示为一次函数在这个区间上的积分。$$ \int_{x_0} ^ {x_1} (\frac{x - x_1}{x_0 - x_1} \cdot y_0 + \frac{x - x_0}{x_1 - x_0} \cdot y_1 )\mathrm{d}x $$因为这个积分是可以直接计算的，设区间长度（步长）为$h$，结果为：$$ S = \frac{h}{2} \cdot y_0 + \frac{h}{2} y_1=\frac{y_0 + y_1}{2}h $$对于这个结果，其正好时梯形的面积公式，因此当为线性插值时，该方法亦称之为梯形法。同理，采用二阶多项式时，依据拉格朗日插值多项式的要求，此时需要3个点，设为$x_0, x_1, x_2$，且取整个区段的长度为$2h$，且连个子区间长度相等，则有：$$ \int_{x_0} ^{x_2} \frac{(x - x_1)(x - x_2)}{(x_0 - x_1)(x_0 - x_2)} \cdot y_0 + \frac{(x - x_0)(x - x_2)}{(x1 - x_0)(x_1 - x_2)}\ cdot y_1 + \frac{(x - x_0)(x - x_1)}{(x_2 - x_0)(x_2 - x_1)} \cdot y_2\\ $$结果为：$$ \frac{h}{3}(y_0 + 4y_1 + y_2) $$这只是对一小段积分区间的计算，当每一小段的积分近似结果都得到确定后，然后进行求和，由此完成积分的数值计算。对于梯形法来说，分析的对象是由两个点组成的一个区间，当对每一段求和后，可以知道的是，整个积分区域中的第一个子区间上的$y_0$，以及最后一个区间上$y_1$均只出现一次，对于其他的点对应的$y$值来说，它们即是上一个梯形的下底，也是下一个梯形的上底，所以求和后会出现两次，这样真正在使用时的梯形法表示为：$$ \int f(x)\mathrm{d}x = \frac{h}{2} \cdot (y_{0} + y_{n} + 2 \cdot \sum_{i = 1} ^ {n - 1} y_{i}) $$这里，离散函数值向量元素个数为$n+1$个，index从0开始。对于Simpson法来说，特别的一点是在研究“单个”对象时，由于其实二阶的插值，需要三个点，使得区间被分为2段，因此，就整体求解来说，该方法的子区间数目必须为偶数才能获得计算。设离散点数为$2n + 1$，则区间总的数量为$2n$。$$ \int f(x) \mathrm{d}x = \frac{h}{3} \cdot (y_{0} + y_{2_n} + 4 \sum_{i=1} ^{n} y_{2 i -1} + 2 \cdot \sum_{i= 1} ^ {m - 1} y_{2i}) $$这里的求和，注意区分下标为奇数以及偶数的区别。每一段中，只有中间点的函数值系数是4，并且它与前后的区段不接触，因此，单纯地保持求和就可以；除去首尾区间，区段边界都是连续接触的，所以求和会出现两次。编程实现时，函数向量$y$的获得在外部单独计算，并且这里的计算都是等步长的。import numpy as np # 梯形法 def trapz_method(y, x): h = x[1] - x[0] s = h / 2 * (y[0] + y[-1] + 2 * y[1:-1].sum()) return s # 以x**2在[0, 1]上的积分为例 x = np.linspace(0, 1, 201) y = x**2 res = trapz_method(y, x)'''Simpson法中必须保证离散的点个数为奇数，才能确保 区间段数为偶数。 ''' def simpson_method(y, x): h = x[1] - x[0] s = h / 3 * (y[0] + y[-1] + 4 * y[1:-1:2].sum() + 2 * y[2:-1:2].sum()) return s
不动点迭代不动点迭代法是利用函数的不动点原理来构造迭代关系式的一种方法。函数的不动点的定义为：对于函数$f(x)$来说，如果存在一个值$x_0$，并且这个值带入函数后，满足这个关系：$$ f(x_0) = x_0 $$那么就称$x_0$ 为函数$f(x)$的一个不动点。比如，对于余弦函数$\cos(x)$来说，0.739085是它的不动点。对函数来说，不动点原理形成了一个迭代关系式，因此，借助于这个思路，在方程的求解时，努力地构造出这样的形式，把方程变换为一侧为$x$，另一侧为一个函数即可。那么就可以得到迭代关系式：$$ x_{k+1} = f(x_{k})\quad k=0,1,2,3,....\infty $$但是，FPI中非常重要的一点就是，这样的迭代关系能否保证收敛？在这里，只要保证迭代式中的函数在不动点附近的斜率绝对值小于1即可，即导数值小于1，这样就能保证收敛。例如，求方程$x^3 + x - 1 = 0$的根。如果我们构造这样的形式：$$ x = 1 - x ^ 3 $$通过计算，我们会发现这一迭代计算是发散的，不会得到一个有用的估计值。因为函数$1 - x ^ 3$不能完全的保证它的导数的绝对值小于1。所以，要想得到有效的形式，必须要是一侧函数的斜率绝对值小于1。得到：$$ x = \sqrt[3]{1 - x} $$import numpy as np import matplotlib.pyplot as plt # 定义一个fpi的函数用于计算方程的根 def myfpi(f, x0, iter_max, eps): # x0为初始值 # f为构造后的函数 # iter_max最大的迭代次数，避免陷入死循环 # eps为阀值，当误差小于这个这个值停止计算 number = 0 # 用于记录迭代的次数，先初始化为0 while np.abs(x0 - f(x0)) > eps: print('the error is: ',np.abs(x0 - f(x0))) x1 = f(x0) x0 = x1 number += 1 if number > iter_max: print('Warning: The result is not stable!') return return x0 # 定义迭代式中的函数 def func1(x): # 这是另外一种，收敛速度更快 return (1 + 2 * x**3) / (1 + 3 * x**2) def func2(x): # (1-x)**(1/3) return (1 - x)**(1/3) x0 = 0.4# 初始值 max_time = 400 eps = 0.0001 res1 = myfpi(func1, x0, max_time, eps) print(res1) res2 = myfpi(func2, x0, max_time, eps) print(res2)结果的输出为：the error is: 0.36216216216216224 the error is: 0.07470645955965782 the error is: 0.005105524662041616 0.6823501779404628 the error is: 0.4434326653017492 the error is: 0.30445961169099023 the error is: 0.2335452355366372 the error is: 0.1620700735059284 the error is: 0.11988614320526159 the error is: 0.08427085993241235 the error is: 0.06129859706965046 the error is: 0.04345059424243991 the error is: 0.03135485804934435 the error is: 0.022330818149177367 the error is: 0.01605160099338865 the error is: 0.011460910899801835 the error is: 0.008222159662768846 the error is: 0.00587845620934957 the error is: 0.004213100804023595 the error is: 0.0030142409755113952 the error is: 0.0021592285329744554 the error is: 0.0015453565272561809 the error is: 0.0011067207210083696 the error is: 0.0007922226584874403 the error is: 0.0005672828725480361 the error is: 0.00040611545134028315 the error is: 0.00029078551365002703 the error is: 0.00020818212995765695 the error is: 0.00014905678931631527 the error is: 0.00010671688227037457 0.6822832757973981误差不论是之前的二分法，还是上述的FPI法，它们都是一个迭代的方法，只不过形成迭代关系时用到的构造原理不同而已，但是仍然都是迭代法的一种。从计算的角度来看，迭代必须采用循环来让其重复计算，知道获得满意的的结果。我们通过设置一个误差标准，当每一次计算的结果的偏离值满足设定的值时，就停止计算。对于while循环来说，必须有入口条件，而数值结果的误差就是描述、控制循环停止强有力的一个标准。后向、前向误差后向误差，指的是依找计算的输入的过程来看的，因为计算最终的步骤就是获得一个解$x_a$，因此，当我们使用$f(x_a)$这一项时，称这个为后向误差。那么只要使用了输入值，就是前向误差，前向误差一般都是指，当前的结果与真实值$r$之间的偏离度，因为真实值往往不知道，所以使用前向误差很有逻辑上的困难。迭代的收敛当我们使用后向误差来定义误差时，即$e_{i} = |r - x_{i}|$，那么有$$ \lim_{i \rightarrow \infty} \frac{e_{i + 1}}{e_{i}} = S $$当这一结果存在，并且收敛到S，则称为线性收敛，收敛速度为$S$。因此，可以知道的是，二分法的收敛速度是$\frac{1}{2}$，全程保持这个速度；而FPI由于收敛速度与函数的导数值有关，收敛速度并不确定。当这个速度小于$\frac{1}{2}$时，其速度比二分法要高效。而牛顿拉夫逊法，它的收敛速度$S$为0，计算效率极高，并且使用非常广泛。牛顿拉夫逊法牛顿法计算的原理，在于利用当前估计值处的函数导数（斜率）来构造函数的切线，然后利用这个切线方程与$x$轴的交点，来确定一个估计值，一直迭代循环计算，直到满足误差要求。若我们所要求的方程为$f(x) = 0$，设定一个初始的估计值$x_0$，对应的函数值为$f(x_0)$，由此确定了一个函数曲线上的一个点$(x_0, f(x_0))$，则该点处的斜率（导数值）为$f'(x_0)$，由此可以确定处经过该点的直线方程，然后利用这个直线方程去确定新的估计值。$$ f'(x_0)(x - x_0) = 0 - f(x_0)\\ x - x_0 = - \frac{f(x_0)}{f'(x_0)}\\ x = x_0 - \frac{f(x_0)}{f'(x_0)} $$由此得到迭代关系为：$$ x_{k+1} = x_{k} - \frac{f(x_k)}{f'(x_k)} \quad k=0,1,2,3,\cdots $$牛顿法的求解方式思路清晰，很容易理解。编程实现时，比较麻烦的是必须要确定出函数的导数形式。def newton_method(f, fd, x0, iter_max, eps): number = 0 while np.abs(f*(x0)) > eps: x1 = x0 - f(x0) / fd(x0) x0 = x1 number += 1 if number > iter_max: print('Warning: The result is not stable!') return return x0
数据的拟合、回归分析最为常用的一个原理就是最小二乘原理，这一原理不仅在数据的处理中得到广泛的应用，而且在深度学习中，线性回归中OLS损失函数是最为常用的模型，这里的OLS就是最小二乘原理。最小二乘原理应用的目的是为了控制最优化问题，即以什么样的标准才能确定问题的最优解。对于这一原理（方法）的认识，更为有价值的是为什么要研究这个原理，它能够解决什么样的问题，从这一角度去理解学习，可以避免枯燥的数学公式推导以及证明。现实中遇到了什么问题？现实中，对于一个理工科学生来说，在实验室进行实验研究是很重要的科学研究手段，而实验研究过程数据的记录以及对于数据的处理分析是十分重要的环节。假设，我做了一个实验，得到了一系列的数据（二维的数据，x和y值），可能几百个（也可能上千个），假设然后我发现这些数据的分布接近线性，也就是说x与y之间的变化规律比较符合线性关系，那么我自然会依据数学函数的相关知识，将数据的模型建立为$y = a x + c$，其中$a$和$b$这个一次函数的系数，它们都是未知的；既然我已经观测到了那么多的数据点，将每一个数据点都带入事先设定的这个线性模型中，一定会得到很多个方程（有多少点就有多少个方程），对于我而言，只需要确定两个未知的系数即可，即$a$和$c$；但问题在于，现在方程的数量远远大于未知数的数量，可以取任意两个方程组成方程组进行求解得到一组结果，而这种结果可以有很多个，那么究竟选择哪一个才最合适？我想这个结果，一定是要最优的，最能反映出每一个方程的要求的，这就需要利用最小二乘原理来进行解答。总之，现实中我遇到的问题，现在转为了一个方程组求解的问题，并且具体情况是这样的：得到的方程组中方程的个数远远大于未知数的个数。从向量空间来看待最小二乘原理前述已经阐明了面临的问题，以及这个问题具体的情形是什么。不从线性代数中线性方程组解的特性来理解，因为这样往往过于抽象和公理化，很难透彻地理解最小二乘原理的意义。这里从向量空间的角度去理解，因为从根本上来说，线性代数研究的本质就是研究$N$维向量空间中的变换问题。前面已经说了，这个问题依然是一个线性方程组求解的问题，那么自然可以用以下的式子来表示：$$ A x = b \tag{1} $$比较特别的是，这里的矩阵$A$不是一个方阵，它的行数实际上就是我们实测的数据点的个数，由于我们之前的模型为线性模型，只有两个待确定的系数，因此这个矩阵的列数只有2列。$x$就是我们要求的未知系数，此时它是一个长度为2列向量，具有两个元素即$a$和$c$；向量$b$就是实测数据点的$y$值。这里不用纠结用什么变量表示，因为问题的本质特性不会因为变量名字的改变而发生变化。为了更好地从向量空间的角度理解，我可以将式(1)，依据矩阵乘法的运算法则进行展开，得到：$$ \begin{bmatrix} a_{11}\\ a_{21}\\ a_{31}\\ \vdots\\ a_{n1} \end{bmatrix} x_{1} + \begin{bmatrix} a_{12}\\ a_{22}\\ a_{32}\\ \vdots\\ a_{n2} \end{bmatrix} x_{2} = \begin{bmatrix} b_{1}\\ b_{2}\\ b_{3}\\ \vdots\\ b_{n} \end{bmatrix} \tag{2} $$从这个形式来看，实际上就是n维向量之间的线性关系的研究。上述形式中，第一个列向量和第二个列向量组成的n维向量空间进行线性运算后得到$b$向量。有一点需要说明的是，对于空间来说，人只能认识到三维空间，再高维度的空间无法被可视化，列向量中元素的个数表明了这个向量的空间维度，为了可视化表示这里假定图中的向量是$R^2$中的向量，如下图 。这两个向量各自都是一条线，可以确定一个平面，那么这个平面中任何其他的向量可以用它们的线性组合来表示，也就是$Ax$，也必定是与他们处于同一平面上的（或者说是同以向量子空间中），而真实的列向量$b$肯定不在这个平面上，那么这个问题所要寻求的最佳结果，就是在这个平面上找到一个$b'$向量，它由$x'$来确定，即我们寻找的一个近似的解。使得其与向量$b$最接近！这一要求，自然是向量$b$在这个平面上的投影就是，这里向平面做一个垂线，由此得到投影$b'$向量；由于$Ax$与垂线$b - Ax'$垂直，利用向量的正交性质，内积必定为0，所以可以写为$$ (Ax)^{T}(b - Ax') = 0 \tag{3} $$转置运算后，得到$$ x^{T}A^{T} \cdot (b - Ax') = 0 \tag{4} $$由此得到，满足$$ A^{T}Ax' = A^{T}b \tag{5} $$即可，求解该方程组，得到最小二乘原理的解答。式（3）也被称为法线方程。对于数据拟合的流程从上述对于最小二乘原理的解释以及法线方程的意义，可以得出做数据拟合的基本步骤，此时不借助于数据分析软件，也能完成这一项工作。首先，具有一系列的数据点$(t_1, y_1), (t_2, y_2),(t_3, y_3),\cdots,(t_n, y_n)$，共有$n$个点。选择一个合适的数据模型。比如$y =a t + b$ 等。强制数据点满足模型，即直接带入到模型中，则获得n个方程，形成$Ax=b$的形式。利用法线方程进行求解。$A^{T}Ax = A^{T}b$例子利用随机数构造出拟合数据设定我们理论上精确的模型关系式为$y = 1.25 x^2 + 0.7x + 2.1$。利用该式子构造出一些数据点，然后添加noise进行扰动，得到“实测”的数据点供拟合用import numpy as np import matplotlib.pyplot as plt # x的范围设定在-2,2之间 x = np.linspace(-2, 2, 50) # 50个x的值 # 定义出1.25 x^2 + 0.7x + 2.1这个原始的函数 def fun(x): return 1.25 * x**2 + 0.7 * x + 2.1 # 产生对应的y y = fun(x) # 产生干扰波动，利用随机数 wave = np.random.randn(y.shape[0]) * 0.7 # 构成拟合所需的ydata ydata = y + wave # 数据点的可视化 plt.style.use('ggplot') # ggplot风格 fig,ax = plt.subplots() ax.scatter(x, ydata, marker='o', label='data points') ax.plot(x, y, label='real point') ax.set_xlabel('x') ax.set_ylabel('y') ax.set_title('The points') ax.legend(loc='best') plt.show()# 对数据进行拟合 # 1.依据数据图像，确定一个比较合适的模型 # 这里选用y = a x^2 + b x + c这个模型 # 2. 将所有数据点x,ydata带入，强制构造方程组， # 主要任务是确定系数矩阵A A1 = x**2 A2 = x A3 = np.ones(x.shape[0]) A = np.c_[A1, A2, A3] # b向量就是ydata b = ydata # 3. 获得法线方程，并求解 result = np.linalg.solve(np.dot(A.T, A), np.dot(A.T, b)) print(result)结果为：$$ \begin{bmatrix} 1.23244467 & 0.7380719 & 2.10100332 \end{bmatrix} $$对拟合结果于原数据进行对比。# 参数解包 a, b, c = result # fitting result y_fit = a * x**2 + b * x + c plt.style.use('ggplot') # ggplot风格 fig,ax = plt.subplots() ax.scatter(x, ydata, marker='o', label='data points') ax.plot(x, y, label='real point', color='tab:blue') ax.plot(x, y_fit, linestyle='-.', color='tab:green') ax.set_xlabel('x') ax.set_ylabel('y') ax.set_title('The points') ax.legend(loc='best') plt.show()在scipy库中，optimize模块的curve_fit函数可以进行数据的拟合，只不过，需要先构造一个含有未知数的函数，具体例子如下。import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit def func(x, a, b, c): return a * np.exp(-b * x) + c xdata = np.linspace(0, 4, 50) y = func(xdata, 2.5, 1.3, 0.5) np.random.seed(1729) y_noise = 0.2 * np.random.normal(size=xdata.size) ydata = y + y_noise plt.style.use('ggplot') plt.plot(xdata, ydata, color='tab:blue', label='data') popt, pcov = curve_fit(func, xdata, ydata) plt.plot(xdata, func(xdata, *popt), 'r-', label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt)) plt.xlabel('x') plt.ylabel('y') plt.legend() plt.show()如果你觉得这些对你还有些收获！当然可以打赏我。请点击右下方Tip进行打赏！
Numerical analysis: Solving Equations with Method of Bisection
二分法方程的求解是科学计算中不可或缺的一个部分，这里主要对单一的方程进行求解，针对超越方程这种非多项式构成的方程，由于其复杂度较高，获得解析解基本不太可能。一般常用的数值方法有二分法、不动点原理迭代法和最为著名的牛顿-拉夫逊法。同时还会指出每一种方法的代价是什么，谈到迭代技术，则收敛速度必须要得到关注；此外，对应的代码会利用Python给出。二分法是一种十分重要的计算或查找逻辑，也许二分法这一种方法不是最优秀的求解方法，但是二分法的这一种求解或搜寻的逻辑思路是一种应用广泛的思想，不仅仅在数值计算领域，在计算机编程算法中，二分思想是一个很有效的算法。在数值分析中，二分法的产生的数学基础是依赖于函数的中值定理以及函数的连续性。更通俗一点来看，二分法的思路就是将正确的解“括住”，寻找到一个足够小的区间，满足事先设定好的精度即可。函数的中值定理：若$f(x)$是闭区间$[a, b]$上的连续函数，如果说该函数$f(x)$在区间边界处满足这个关系式:$f(a) \cdot f(b) < 0$，那么该函数在该区间上必定存在一个根，即使得函数为0的自变量取值必定在区间$(a, b)$内。思路与过程这一定理的本质意义可以这样理解。首先，函数的连续性可以让函数曲线的变化是平滑、不间断的，就排除了函数值跳跃的可能性；再者，这一定理中给出了一个先决条件$f(a) \cdot f(b) < 0$，即左端点，右端点处的函数值是异号的，前面由于连续性的要求，对于这一情况的描述之可能存在唯一的一种形式，就是：如果从一个函数值小于0的端点处开始，函数的变化必须平滑、连续，那么函数值要向从负值，变化到正值，一定会穿过$f(x) = 0$这一水平线，那么此时的$x$值就是对应方程的根。上述只是阐述了“二分法”的计算思路，还没有体现出“二分”的概念；让想得到足够精确的解，还需对这个区间的搜索做一点控制工作，此时“二分”的概念得以体现。在寻找“括住根”的区间过程中，很难一步得到一个很精确的区间，因此，迭代思想必定会在这个过程中得到应用，第一步区间可以大一点，但是必须满足中值定理，只要中值定理满足，则这个区间就是具有意义的，为了减小区间范围，对区间长度进行对半分，此时一个关键问题不可避免，这两个区间到底使用哪一个？这一定是一个选择性问题。当我们对二分法求解的数学思路有了认识后，应当可以感觉到，要想把根括住，必须要让函数在区间端点的值是异号，也是说如何选择两个区间中的哪一个，唯一的判断依据就是“继续保持中值定理的满足”，只要这个区间还满足中值定理，那么就选择它。在进行多次的循环后，由于每一次区间长度都减半，一定会得到一个较小的区间。每一次完整的计算后，解的不确定性都会减小一半，得到的最后一个区间，取这个区间的中点作为估计值即可。精度与收敛速度对于二分法的精度的确定，完全可以从整个求解过程中得到，每一次计算都会使得区间长度减半。如果初始的区间是$[a, b]$，则完成第一次完整的计算后，区间长度变为$(b - a) / 2$；第二次计算完成后区间长度变为$(b - a) / 4$，以此类推可以得到进行了n次迭代后的区间长度为$(b - a) / 2^{n}$；注意在最后得到的区间长度是$(b - a) / 2^{n}$，而方程的根是取这个区间的中点作为估计值，由此可以得到二分法的精度为:$$ eps = \frac{b - a}{2 ^{n + 1}} $$对于收敛速度或者效率来说，可以自然得到为$(b - a) / 2$。二分法的实现代码前述二分法过程中，包含了两个重要信息；第一需要多次的迭代计算，使得估计值越来越精确，这意味着需要使用循环结构来实现；此外，在任何的一次计算中，需要对当前区间的二分后，要取选择下一次的区间时，要利用中值定理的条件去做判断，所以分支结构也一定会被利用。伪代码对求解的描述初始区间为$[a, b]$，其已经满足端点处函数值为异号的要求。最外层：利用循环结构控制:循环控制条件为(b - a) / 2 与预先设定的tolerance的关系 获得初始区间的中点：c = (b - a) / 2（中点作为估计值） 进行二分后区间的选择：利用分支语句 if f(c) == 0: 这个就是正确解，直接可以结束循环 if f(a) * f(c) < 0: 左侧端点和中点c满足中值定理的要求，则选择它 b = c 直接利用c作为右侧端点即可 else: 此时意味着左侧的a处函数值与c处的值不满足中值定理，选择另一半 a = c c作为左侧端点使用即可Input的确定:初始区间（这个区间必须满足中值定理的描述），以端点值最为输入即可tolerance值得给定，用于迭代计算得终止所要求解的方程对应的函数$f(x)$，以def定义的函数形式引入即可Output：最终的估计值是必须要返回的。其他的可以根据自己的需求而设定CODE# 库的调用 import numpy as np # 导入numpy库，且重命名为np import matplotlib.pyplot as plt # 定义了一个二分法求解方程的函数 def bisection_method(f, a, b, tol=0.0001): # f:待求的方程对应的函数 # a:初始区间的左侧端点值 # b:初始区间的右侧端点值 # tol: 当区间长度小于这个 # 值时，结果可以接受，可结 # 束计算，默认值为0.0001 # 确定端点对应的函数值 fa = f(a) fb = f(b) # 利用while循环控制 while (b - a) / 2 > tol: # 入口条件满足就要执行 mid = (b + a) / 2 # 这一步实际上就是进了区间二分 # 中点处的函数值f(mid) fmid = f(mid) if fmid == 0: return mid # 此时c就是真正的解 # 利用if语句进行区间的选择 # 注意，区间的左、右始终用变量a，b来表示 if fa * fmid < 0: # 选择中点左侧的区间 b = mid fb = fmid else: a = mid fa = fmid # 上述循环结束后，区间a，b的值均已经发生了 # 变化，取此时的中点作为估计值即可。 return (a + b) / 2eg:求方程$\sin(x) - 6x -5=0$的根。首先确定出方程根大致的区间# 先定义出函数 def fun(x): return np.sin(x) - 6 * x - 5 # 利用函数图像来确定初始区间范围 # 产生x的离散点，范围可自行尝试。 # 此为-10，10之间的，产生100个点 x = np.linspace(-10, 10, 100) # 带入定义好的函数fun y = fun(x) fig, ax = plt.subplots() ax.plot(x, y,linestyle='--', linewidth=2, color='tab:blue',label='fun') ax.set_title('$\sin(x) - 6x -5$') ax.set_xlabel('X') ax.set_ylabel('Y') ax.grid('True') ax.legend(loc='best') plt.show() 由此，先可以大致确定出初始区间可以选择为$[-2.5, 0]$。# 调用bisection函数进行求解 # 前述的代码已经完成了函数的定义以及方程对应函数的确定 # 此时，要对方程进行计算即可 sol = bisection_method(fun, -2.5, 0, 0.0001) print(sol)结果为$-0.9708404541015625$。OOP方式的编程实现为什么要采用面向对象的编程？上述采用面向过程的编程范式，通过函数定义的实现计算过程，虽然也很好地解决了问题，但是更进一步地想，我们真正面对的问题是“求方程的根”，至于用什么数值方法可能不是我们的目的，并且方法的种类又不唯一。如果按照之前定义函数的方式来实现，求解对象与数值方法函数之间时存在割裂的。而采用OOP时，可以通过类的定义，确定出方程这个对象，然后将方程的具体形式（也就是这个方程对应的函数）作为对象的属性，初始区间、精度都作为对象的属性，而方程的解（或者说是求解）定义为方程对象的方法，而方法的数量又不是唯一的，因此，可以定义多个方法，如二分法、不动点迭代法，以及牛顿-拉夫逊法，这样的封装形式更加方便。import numpy as np # 定义一个方程对象 class myequation(object): def __init__(self, fun, a, b, tol=0.0001): self.fun = fun self.a = a self.b = b self.tol = tol def bisection(self): a_ = self.a b_ = self.b fa = self.fun(a_) fb = self.fun(b_) while (b_ - a_) / 2 > 0: mid_ = (b_ + a_) / 2 fmid = self.fun(mid_) if self.fun(mid_) == 0: return mid_ if fa * fmid < 0: b_ = mid_ fb = fmid else: a_ = mid_ fa = fmid return (a_ + b_) / 2 # 应用 def fun(x): return np.sin(x) - 6 * x - 5 eq = myequation(fun, -2.5, 0) sol1 = eq.bisection() print(sol1)结果为$-0.9708989235042559$。对于二分法的理解核心是二分这种逼近、缩小范围的思路，对于方程的求解，二分的实现方式被限定在满足函数中值定理的条件下，通过多次的迭代计算获取一个能够接受的小区间，真实解是存在于这个区间上。如果你觉得这些对你还有些收获！当然可以打赏我。请点击右下方Tip进行打赏！