分享我的学习心得 - Numerical Analysis
type: channel
link: http://www.pynumerical.com/tag/Numerical-Analysis/
items: Numerical Analysis: 改进的Euler法--显式梯形法
items: Numerical Analysis: 微分方程的欧拉法
items: Numerical Analysis: 共轭梯度法(2)--实现
items: Numerical Analysis: 共轭梯度法(1)--基本原理
items: Numerical Analysis: 什么是梯度下降
items: Numerical Analysis: Iteration Method of Equation Systems
items: Numerical Analysis: Numerical Integration
items: Numerical Analysis: 方程求解的不动点迭代与牛顿拉夫逊法
items: Numerical Analysis: 最小二乘原理
items: Numerical analysis: Solving Equations with Method of Bisection

Numerical Analysis: 改进的Euler法--显式梯形法
type: item
date: 2020-02-18T21:43:00+08:00
description: 经典的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…

Numerical Analysis: 微分方程的欧拉法
type: item
date: 2020-02-15T11:32:00+08:00
description: 一阶常微分方程的欧拉法一阶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$意味着从不同的点出发开始移动。

Numerical Analysis: 共轭梯度法(2)--实现
type: item
date: 2019-11-01T22:07:00+08:00
description: 梯度下降法的缺点之前的内容已经完全学习了梯度下降算法，对于这个方法之所以在线性方程组求解领域没有更广泛地流行起来，最大的原因是每一步的迭代中近似值移动的方向是沿着负梯度方向，而相邻迭代步的残差向量(本质就是负梯度向量)$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.997956…

Numerical Analysis: 共轭梯度法(1)--基本原理
type: item
date: 2019-10-29T15:53:00+08:00
description: 共轭梯度法(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, …

Numerical Analysis: 什么是梯度下降
type: item
date: 2019-10-27T21:41:00+08:00
description: 梯度下降法的原理梯度下降法是一种数值计算方法，但是在数值分析中往往并没有直接采用它，但是在这一方法基础上发展而来的共轭梯度法却被得到应用，而梯度下降法是其核心思想与方法。目前，较为火热的深度学习方面，由于广泛采用了自动求导原理，对损失函数求极值时，梯度下降算法在该领域得到了十分广泛的应用。而我们感兴趣的是利用共轭梯度法求线性方程组，该方法的思想就是梯度下降原理，因此，很有必要对该梯度下降法有一个认识。什么是梯度下降在微积分中，很常见的一个公式就是泰勒展开，这个公式意在描述函数$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_t…

Numerical Analysis: Iteration Method of Equation Systems
type: item
date: 2019-10-19T10:13:00+08:00
description: 对角占优线性方程组的迭代法线性方程组的求解是多数数值问题最终的一个环节，很多问题最终都变为了一个求解线性方程组的问题。求解方程组最广为人知的方法是高斯消元法，这是一种直接方法，通过构造某种形式，可以直接得到最终的结果，其矩阵形式就对应着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…

Numerical Analysis: Numerical Integration
type: item
date: 2019-10-07T21:32:00+08:00
description: 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$的获得在外部单独计算，并且…

Numerical Analysis: 方程求解的不动点迭代与牛顿拉夫逊法
type: item
date: 2019-09-21T15:37:00+08:00
description: 不动点迭代不动点迭代法是利用函数的不动点原理来构造迭代关系式的一种方法。函数的不动点的定义为：对于函数$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 err…

Numerical Analysis: 最小二乘原理
type: item
date: 2019-09-20T13:24:00+08:00
description: 数据的拟合、回归分析最为常用的一个原理就是最小二乘原理，这一原理不仅在数据的处理中得到广泛的应用，而且在深度学习中，线性回归中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）也被…

Numerical analysis: Solving Equations with Method of Bisection
type: item
date: 2019-09-20T13:21:00+08:00
description: 二分法方程的求解是科学计算中不可或缺的一个部分，这里主要对单一的方程进行求解，针对超越方程这种非多项式构成的方程，由于其复杂度较高，获得解析解基本不太可能。一般常用的数值方法有二分法、不动点原理迭代法和最为著名的牛顿-拉夫逊法。同时还会指出每一种方法的代价是什么，谈到迭代技术，则收敛速度必须要得到关注；此外，对应的代码会利用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：最终的…