PDE采用解析方法求解后，通常解的形式都是傅里叶级数形式表示的，计算各个采样点的值需要完成级数的求和计算。$$ u(z,t) = \frac{4 u_0}{\pi} \sum_{n=1}^{\infty} \frac{1}{2n-1} \sin(\frac{2n-1}{2H} \pi z) \exp(- \frac{((2n-1)\pi}{2H})^2 C_v t) $$一维空间中的非稳态问题，场函数必定是一个二元函数，自变量有$x$和$t$，级数的计算还需要考虑$n$的取值，需要计算每一个$n$对应下的结果，因此，需要不断地更新$n$生成新的函数，然后求和就可以得到$z,t$对你应下的结果。利用sympy实现计算，其中summation函数能够实现级数的计算，因此，只需要每一次更新$z,t$的值即可。import numpy as np import sympy as sp import matplotlib.pyplot as plt n = sp.symbols('n') H = 1.0 zv = np.linspace(0, H, 11) t = 1e5 tv = np.linspace(0, t, 50) cv = 1e-5 u0 = 1000 # 初始化结果矩阵，时间为竖向，空间值按照列进行放置 sol = np.zeros([tv.shape[0], tv.shape[0]]) for row in range(tv.shape[0]): for col in range(zv.shape[0]): f = 1 / (2 * n - 1) * sp.sin((2 * n - 1)/ (2 * H) * sp.pi * zv[col]) * sp.exp(-1 * ((2*n-1)*sp.pi/ (2*H))**2 * 1e-5 * tv[row]) sol[row, col] = 4 * u0 / np.pi * sp.sumation(f, (n, 1, 50)).evalf()evalf()方法，能够将计算结果转换为浮点数进行表示，sympy中类似的还有sympify和lambdify，都能实现将结果以数值的形式表示出来。
pygmsh中除了最基本的geo模块，为了适应更加复杂以及通用的CAD那样的建模方式，还提供了occ模块。从源代码看，geo以及occ中均有Geometry()类，但是他们的类方法具有不同，当然，需要注意的是，这个类不管是是在geo，还是在occ中，这个类都是继承于一个父对象的，父对象本身包含了一些最基本的方法，比如add_polygon 等方法，所在在使用的时候既有差异，也存在相同的操作。occ中Geometry对象存在的意义，在大程度上体现在能够进行几何对象的布尔操作，当然一些方法也需要注意，比如disk、rectangle等等，这些是基本图形操作。生成以后，可以用boolean_difference(c0, c1) 来进行差集运算。boolean_union([c0, c1])来进行交运算。通过布尔运算能够生成更加复杂的几何对象。当然，这个过程中一定要注意mesh_size的施加方式，多种布尔操作下，最好通过geom.characteristic_length_max来设定。其余的，像拉伸、旋转等操作与geo中的操作都是一样的，只不过布尔运算必须是occ中的几何对象才可以。import pygmsh with pygmsh.occ.Geometry() as geom: geom.add_disk([0, 0], 2, mesh_size=0.1) mesh = geom.generate_mesh() with pygmsh.occ.Geometry() as geom: geom.characteristic_length_max = 0.05 rec = geom.add_polygon([[0, 0], [1, 0], [1, 1], [0, 1]]) cr = geom.add_disk([0.5, 0.5], 0.15) geom.boolean_difference(rec, cr) mesh = geom.generate_mesh() with pygmsh.occ.Geometry() as geom: geom.characteristic_length_max = 0.02 rec = geom.add_polygon([[0, 0], [1, 0], [1, 1], [0, 1]]) r1 = geom.add_disk([0.5, 0.5], 0.15) rr = geom.boolean_difference(rec, r1) d3 = geom.extrude(rr, [0, 0, 1], num_layers=10) mesh = geom.generate_mesh() with pygmsh.occ.Geometry() as geom: geom.characteristic_length_max = 0.2 disk1 = geom.add_disk([0, 0], 5) disk2 = geom.add_disk([0, 0], 2) flat = geom.boolean_difference(disk1, disk2) geom.extrude(flat, [0, 0, 0.3]) mesh = geom.generate_mesh()
1. 方程二维Poisson方程的形式为：$$ \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} = f(x, y) $$这里不给出边界条件的原因是，采用数值方法（尤其是有限元方法）对于模型几何域的不敏感，任何形状都能够完成求解。边界条件的施加可以在控制方程得到离散后进行。2. Galerkin法进行PDE方程的离散Poisson方程属于边值问题，即物理场函数仅是空间坐标的函数，与时间无关。为了获得其有限元形式，需要首先构建出PDE对应的等效泛函形式。这里采用加权残余方法，当权函数取形函数时，该方法被称为Galerkin法。首先，有限元法应当先确定出单元的类型，通常为了计算方便选择低阶三角形单元，一方面在单元积分计算时能够直接进行显式表达，另一方面，三角形单元的保形能力较好，对于几何域的逼近效果要好。2.1 单元类型单元类型为三节点三角形单元，选定单元类型意味着形函数就会确定，这里本质要做的工作就是，在一个单元内，用这个单元三个节点处的场函数值，利用插值的思想去确定出单元内任意位置处的值。形函数的形式通常选择线性插值，利用单元三个节点的值，能够反向求解出插值系数，从而得到形函数的形式。即$$ N_i = \frac{1}{2A} (\alpha_i + \beta_i x + \gamma_i y)\\ N_j = \frac{1}{2A} (\alpha_j + \beta_j x + \gamma_j y)\\ N_m = \frac{1}{2A} (\alpha_m + \beta_m x + \gamma_m y) $$此时，单元内的任意点的场函数值可以表示为$$ u(x, y) = N_i \cdot u_i + N_j \cdot u_j + N_m \cdot u_m $$表示为矩阵形式为：$$ u(x, y) = \begin{bmatrix} N_i & N_j & N_m \end{bmatrix} \cdot \begin{bmatrix} u_i \\ u_j\\ u_m \end{bmatrix}=[N] \cdot \{u_e\} $$2.2 PDE等效积分形式利用加权残余方法构建方程的积分形式。这里我们寻求的是数值解，是一个近似解答，如果记数值近似解为$\omega(x,y)$，既然为近似，那么代入原始的控制方程中必定是不满足的，因此，会存在一个残差$R$，为$$ R = \frac{\partial ^2 \omega}{\partial x ^ 2} + \frac{\partial ^2 \omega}{\partial y^2} - f(x,y) $$寻求近似解，自然要找到误差或者偏差最小的那一个，因此，为了找到误差最小的，强制使得的残差为0，应当是最好的。注意残差应当是相对于区域而言的，此时针对的是一个单元，所以需要函数内积的概念，加权函数选择形函数相同的，得到$$ \iint_\Omega [N]^\mathrm{T} (\frac{\partial ^2 \omega}{\partial x^ 2} + \frac{\partial ^2 \omega}{\partial y^2} - f(x,y))\mathrm{d}x \mathrm{d}y $$这个形式其实就是PDE所对应的等效积分强形式，其中函数仍然要求保证二阶导的连续，对于场函数的连续并没有降低，因此，称之为强形式。此时进行求解仍然具有较大的困难。并且，数学上二阶连续性的要求对于现实中的情况而言过高，降低连续性通常会得到更加真实的解答。因此，寻求弱形式似乎很有必要。此时利用高斯定理，对于二维问题高斯定律则退化为格林公式，本质上就是积分的分部积分。$$ \iint_{\Omega}[N]^{\mathrm{T}} \frac{\partial ^ 2 \omega}{\partial x^2}\mathrm{d}x\mathrm{d}y= \int[\int[N]^\mathrm{T} \frac{\partial ^ 2 \omega}{\partial x ^ 2}\mathrm{d}x]\mathrm{d}y $$进行分部积分，有$$ \int[\int[N]^\mathrm{T} \frac{\partial ^ 2 \omega}{\partial x ^ 2}\mathrm{d}x]\mathrm{d}y= \int [{[N]^\mathrm{T}\frac{\partial \omega}{\partial x}|_{\partial \Omega}} - \int\frac{\partial [N]^{\mathrm{T}}}{\partial x}\frac{\partial \omega}{\partial x}\mathrm{d}x]\mathrm{d}y $$注意到，在离散方程的时候，此时没有考虑边界条件，因此边界均是无约束的情况，此时分部积分的第一项的值应当为0，由此得到$$ \iint_{\Omega}[N]^{\mathrm{T}} \frac{\partial ^ 2 \omega}{\partial x^2}\mathrm{d}x\mathrm{d}y= -\iint_{\Omega}\frac{\partial [N]^{\mathrm{T}}}{\partial x} \frac{\partial \omega}{\partial x}\mathrm{d}x\mathrm{d}y $$关于$y$的偏导项进行同样的处理后，得到等效积分弱形式。$$ \iint_{\Omega}[\frac{\partial [N]^{\mathrm{T}}}{\partial x}\frac{\partial \omega}{\partial x}+\frac{\partial [N]^{\mathrm{T}}}{\partial y}\frac{\partial \omega}{\partial y} - [N]^{\mathrm{T}}f(x,y)]\mathrm{d}x\mathrm{d}y=0 $$上式即为等效积分弱形式，其中场函数要求一阶连续，连续性得到降低。为了得到有限元形式，在弱形式的基础上，利用形函数、节点处函数值表示的场函数代入后，得到$$ \iint_{\Omega}\frac{\partial [N]^{\mathrm{T}}}{\partial x}\frac{\partial [N]}{\partial x} + \frac{\partial [N]^{\mathrm{T}}}{\partial y}\frac{\partial [N]}{\partial y}\mathrm{d}x\mathrm{d}y \cdot \{u_e\} = \iint_{\Omega}[N]^{\mathrm{T}}f(x,y)\mathrm{d}x\mathrm{d}y $$此时，得到了离散后的有限元形式，左端将形函数矩阵以及形函数表达式代入后，可以发现，形函数的一阶导是域空间坐标无关的，因此，这个积分是直接可以进行计算的，积分就是这个单元的面积$A$。而形函数矩阵相乘后得到左端为$$ [B] = \frac{1}{4A} \begin{bmatrix} \beta_i^2 + \gamma_i^2 & \beta_i \beta_j + \gamma_i \gamma_j & \beta_i \beta_m + \gamma_i \gamma_m\\ \beta_j \beta_i + \gamma_j \gamma_i & \beta_j^2 + \gamma_j^2 & \beta_j \beta_m + \gamma_j \gamma_m\\ \beta_m \beta_i + \gamma_j \gamma_i & \beta_m \beta_j + \gamma_m \gamma_j & \beta_m^2 + \gamma_m^2 \end{bmatrix} $$现在只剩下右端的源项的积分了，因为形函数与$x,y$有关，并且源项$f(x,y)$也是$x,y$的函数，因此，这个积分不能直接进行计算，它是在三角形区域上的一个积分，采用数值积分应当是最好的形式。即需要用高斯积分进行计算。由于整个推导采用的都是基于全局坐标进行的，并没有引入等参描述，采用高斯积分时还需要利用等参描述，太过于麻烦。这样的情况下一般会采用近似积分来代替，即利用单元中心处的积分值乘以单元面积进行近似计算$$ x_c = \frac{x_i + x_j + x_m}{3}\\ y_c = \frac{y_i + y_j + y_m}{3} $$有：$$ \iint_{\Omega}[N]^{\mathrm{T}}f(x,y)\mathrm{d}x\mathrm{d}y=[N(x_c, y_c)]^{\mathrm{T}} \cdot f(x_c, y_c) \cdot A $$从而得到右侧的单元荷载向量，长度为3，记作$F$，则单元的线性方程组为$$ [B]\{\omega_e\} = \{F\} $$在获得单元刚度矩阵以及荷载向量后，通过直接集成的方法得到整个模型的刚度矩阵和向量，然后利用乘大数方法或者置1法进行边界条件的施加，然后通过线性方程组的求解，获得最终的解答。
经典的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$意味着从不同的点出发开始移动。
$SI$模型$N$表示一个城市的总人口，在不考虑出生人口和死亡人口的情况下，$I$代表已经感染人的个数，称为患者数（或者僵尸数）；这一模型将人群分为两类，一类是感染的人，另外一类是未被感染的人，也称为易感染人，数量用$S$表示。自然可以有$N = I + S$。每一天一个患者能够传染的速率为$\lambda$，那么已有的感染人的数量乘以这个感染传播速率就得到新增的感染人数，$I \times \lambda$。这些都是对于每一天来说的，需要注意的是 前一天的情况会得到下一天的情况，因此，这一个过程是一种完全的迭代实现。令感染者在人群中的比例为$i = \frac{I}{N}$，易感染者在人群中的比例为$s = \frac{S}{N}$，因为对于人群的划分只分为两类，所以得到$1 = i + s$。依据上述的表达，可以计算出每一天新增的感染人数为$\lambda \times iN \times s$，这里乘以$s$是因为，只要感染者增加，那么易感人数必定减少，而新增的感染者一定是从未感染者中来的。这样新增感染者比例的增加量为$\lambda i s$。同时借助于$1 = s + i$，得到$s = 1- i$，所以后一天感染者所占的比例为$i + \lambda i (1 - i)$。由此通过该迭代关系可以完成预测。以一个例子实现，假设总人口为三百七万人(3.7e6)，初始时刻感染的人为10个，感染速率$\lambda$设定为1，即每一天一个感染者会感染一个人，模型总时间为100天。这里不讨论参数如何确定，只是单纯地为了计算，结果并不能代表现实中的情况。import numpy as np import matplotlib.pyplot as plt class SI(object): def __init__(self, N, T, m, alpha): self.N = N self.T = T self.m = m self.alpha = alpha self.get_infective_ratio() def get_infective_ratio(self): ii = np.zeros(self.T) ii[0] = self.m / self.N for t in np.arange(self.T - 1): ii[t+1] = ii[t] + self.alpha * ii[t] * (1 - ii[t]) self.s = 1 - ii self.i = ii return self.i, self.s def visulization(self): plt.style.use('ggplot') fig, ax = plt.subplots() ax.plot(self.i, color='tab:red', linewidth=2, label='I') ax.plot(self.s, color='tab:blue', linewidth=2, label='S') #ax.grid(True) ax.set_xlabel('Days/ day') ax.set_ylabel('Ratio') ax.set_title('SI Model') ax.legend(shadow=True,loc='best') fig.savefig('si.png', dpi=600) plt.show() demo = SI(3.7e6, 100, 10, 1) demo.visulization()$SI$模型是一个最简单的模型，这一模型中人群分为感染者和未感染者两种，因为总人口不变，所以这两者之间的关系一定是“此消彼长”的，感染者数量的增加一定会导致未感染者数量的下降，所以，直观地可以意识到，感染者所占的比例已开始会缓慢增加，随后会随着已有感染者数量的速度加快，但是这不是无限的，因为未感染者数量在持续的下降，直到两者数量相等时，出现分水岭，感染比例速度开始下降。这一模型过于理想化，考虑的因素也过于简单。$SIR$模型$SI$模型中没有考虑死亡和治愈的人的影响，实际情况中，死亡的人不会再进行疾病的传播，如果患者被治愈其体内会有抗体存在，应当来说是不会再被感染的，这两种情况造成的影响都是相同的，都是“去除”的效果，这也是$R$的意义(remove)。$SIR$模型是在$SI$模型上的一个改进，基本的原理不变，只是加入了$R$的影响。移出者$R$一定是来自于已有的感染人群，不论是死亡还是治愈，设定移出率为$\gamma$。这一模型中人群被分为了3类。$S$易感染者：由于疾病一直在传播，所以易感染的人数量一直在下降，它的变化数量就是增加的感染者数量，$\lambda \times N \times i \times s$，则对应的变化比例为$\lambda i s$。$I$感染者：不考虑移除为$\lambda N i s$，现在加入了移除$R$，那么需要减去$R$的数量，得到$\lambda N i s - \gamma N i$，得到对应的比例为$\lambda is - \gamma i$。$R$移出者：移出者是来自于已经感染的人群，即$\gamma N i$，所占的比例为$\gamma i$。这样基于上述的描述，利用迭代计算可以计算出设定的天数内各个人群比例的变化值。为了演示计算，设定总人口依然为三百七十万人，初始时刻的感染者为6人，$\lambda=1$，Remove速率为0.3，分析的总时间是100天。class SIR(object): def __init__(self, N, T, m, alpha, gamma): self.N = N self.T = T self.m = m self.alpha = alpha self.gamma = gamma self.get_result() def get_result(self): i0 = np.zeros(self.T) s0 = np.zeros(self.T) r0 = np.zeros(self.T) i0[0] = self.m / self.N s0[0] = 1 - i0[0] for t in np.arange(self.T-1): i0[t+1] = i0[t] + self.alpha * i0[t] * s0[t] - i0[t] * self.gamma s0[t+1] = s0[t] - self.alpha * i0[t] * s0[t] r0[t+1] = r0[t] + self.gamma * i0[t] self.i = i0 self.s = s0 self.r = r0 return self.i, self.s, self.r def visulization(self): plt.style.use('ggplot') fig, ax = plt.subplots() ax.plot(self.i, color='tab:red', linewidth=2, label='I') ax.plot(self.s, color='tab:blue', linewidth=2, label='S') ax.plot(self.r, color='tab:orange', linewidth=2, label='R') #ax.grid(True) ax.set_xlabel('Days/ day') ax.set_ylabel('Ratio') ax.set_title('SIR Model') ax.legend(shadow=True,loc='best') fig.savefig('SIR.png', dpi=600) plt.show() demo1 = SIR(3.7e6, 100, 6, 1, 0.3) demo1.visulization()这里只是理想地进行一个计算，具体情况影响因素很多，上述的模型都是过于理想化的，但是能反映出感染率等因素对于疾病传播规模的影响，其中$\lambda$值在实际情况中就体现为防控措施。
抛物线方程的有限差分法-Dirichlet条件下偏微分方程求解问题中，是否有解析解很大程度上受到边界条件、初始条件以及方程形式的影响，在复杂边界条件下，寻求解析解是不可能完成的任务，因此，偏微分方程的数值解法具有很强的实用价值，对于物理场问题的描述，离不开数值解法。PDE的数值求解方法中，有限差分法是最基本的、也是求解思路最为清晰的一种数值方法，这一方法构造简洁，是在对场问题描述、理解的基础上得来的；在一维问题中，FDM要比FEM有着更高的精度和计算效率。对于PDE的数值求解，以FDM作为学习的开始，是一个不错的选择。有限差分形式PDE是是一种特殊的方程，在方程中往往包含着场函数的导数，而最终的求解是需要得到场函数本身，那么从结果出发来看，对于一个一维的非稳态问题来说，场函数的是一个二元函数，对于它的数值描述，实际可以这样理解，即在这个域上，离散后会得到许多个离散的点，每一个点的坐标实际上是已知的，这样，对于这个二元函数来说，它在每一个离散点处函数值，就是我们需要的数值解的描述，一般我们更愿意依据离散点的分布形式，将它的函数值放在一个同形状的矩阵中。数值方法是一种尽可能近似的、有限度的描述场函数分布的方法，理论上的真解是连续函数，而数值解描述的是区域上某点处的场函数值，因此，数值方法的实现往往第一步就是需要对问题的离散，有限差分法也不例外。这里针对的是二阶线性偏微分方程，并且是各向同性的，即是一个均匀方程，那么不管怎样，PDE一定会有导数项，如果能将导数表示为场函数，那么方程就完全转化为了场函数本身，如果能够进行求解，那么所得结果正好就是我们所需要的数值解。对于导数差分形式的理解，从泰勒展开的角度来看，个人感觉这是一个比较清晰和简洁的方式。一个函数的泰勒展开形式为：$$ f(x+\Delta x) = f(x)+\frac{1}{1!}f'(x)(\Delta x)^1+\frac{1}{2!}f''(x)(\Delta x)^2+o(x^3) $$泰勒展开式提供了一种思路，即函数在$x + \Delta x$处的函数值可以利用函数在$x$处的函数值以及这个函数的导数去估计，当略去高阶项时，可以得到函数一阶导的表示方式，$$ f'(x)=\frac{f(x+\Delta x) - f(x)}{\Delta x} $$这就是一阶导数的差分形式，此时为两点前向差分公式，它将函数的一阶导数表示为两点处的函数值。当采用点$x$和点$x - \Delta x$这两个点来构造时，得到的一阶导数的两点后向差分公式，后向差分形式是一种很重要的形式，虽然好像由此得到的方程的差分形式没有前向的简洁，但是它所得的差分迭代关系式在数值计算中通常是无条件稳定的，这一点复杂度的代价换来的求解优势还是很值得的。$$ f'(x)=\frac{f(x)-f(x-\Delta x)}{\Delta x} $$上述是一阶导数的差分公式，而PDE中二阶导数是最为常见的，通常需要三个点的函数值来描述二阶导数的值，当略去二阶以上的项时，$x+\Delta x$ 、$x-\Delta x$分别对应的泰勒展开为$$ f(x + \Delta x) = f(x) + f'(x)\Delta x + \frac{1}{2}f''(x)(\Delta x)^2\\ f(x - \Delta x) = f(x) - f'(x)\Delta x + \frac{1}{2}f''(x)(\Delta x)^2 $$上述两式相加，消去一阶导项，可以得到二阶导的三点中心差分公式：$$ f''(x) = \frac{f(x-\Delta x) -2f(x)+f(x+\Delta x)}{(\Delta x)^2} $$以上是最基本的差分公式，多数问题依据它们便可以得到解决，当然，可以通过泰勒展开来构造出更为复杂的形式。抛物线方程的有限差分形式当PDE中的导数项利用差分形式进行表示后，此时方程中是不含有导数项的，完全是由场函数值来组成的。而不论是几维问题，偏导数永远都是针对某一方向而言的，因此，前述的差分形式都是通用的，为了便于计算，可以这样定义：$x、x + \Delta、、x - \Delta$可以理解为中间点，后一个点、前一个点，一个PDE问题，其空间域往往是得到确定的，一维问题中x的长度是已知的，当采用均匀的划分后，$\Delta x$实际上已经得到了确定，如果说$x$方向上，区间左端点为$x_0$，那么下一个点记为$x_1$，其坐标必定为$x_0 + 1 \cdot \Delta x$，以此类推得到$x_i =x_0 + i \cdot \Delta x$，为了表达式的简洁，我们可以直接用点号来表示点，即$u(x_i)$记作$u(i)$。以抛物线方程为例，通过上述的表示方式来构造出PDE对应的有限差分形式。PDE形式如下：$$ D\frac{\partial ^ 2 u(x, t)}{\partial x ^ 2} = \frac{\partial u(x,t)}{\partial t} $$这个方程是一维抛物线型方程的最简形式，它通常的物理意义可以表示一维的非稳态温度场、一维的固结问题等。对于其中的二阶导采用三点中心差分公式来替换，一阶导采用前向差分公式，且对于$x$方向来说点编号记作$i$，时间$t$方向上点编号记作$j$，得到$$ D \cdot \frac{u(i-1, j) - 2u(i, j) + u(i+1, j)}{(\Delta x)^2} = \frac{u(i,j+1) - u(i,j)}{\Delta t} $$该方程是一个非稳态问题控制方程，可以理解为场函数在前一时刻的分布情况下得到下一时刻场的分布，当位于初始时刻时，由于初始条件的存在，便可以得到下一个时刻的状态分布，而每一个时刻中，边界条件都是需要强制满足的，我们只需要将$j$时刻与$j+1$时刻的项分别放在等号左右两侧即可反映出这一种计算思路。令网格比为$\sigma = D \frac{\Delta t}{(\Delta x)^2}$$$ u(i, j+1) = \sigma u(i-1,j)+(1-2\sigma)u(i,j)+\sigma u(i+1,j) $$这是前向差分得到的迭代式，这里先不写出矩阵形式，因为矩阵形式需要考虑到边界条件的影响，确切地来说，边界条件的形式会决定自由点的个数，导致矩阵的大小不同。对于时间的导数，采用后向差分法时，得到后向差分形式为：$$ D \cdot \frac{u(i-1, j) - 2u(i, j) + u(i+1, j)}{(\Delta x)^2}=\frac{u(i,j) - u(i,j-1)}{\Delta t}\\ -\sigma u(i-1,j) +(2\sigma +1)u(i,j)-\sigma u(i+1,j) = u(i,j-1) $$对比前向差分形式与后向差分形式可以发现，前向差分形式实际上是一个显式形式，因为前一时刻的状态都是已知的，所以下一时刻的$i$处的值由前一时刻的3个点完全可以得到确定；而对于后向形式而言，却恰好相反，已知的一个点与下一时刻的3个节点值有关，这是一种隐式。前向形式简洁明了，但是存在一个较大的隐患，即对网格比有要求，为了迭代的收敛，$\sigma$必须小于等于0.5才行，这样使得时间步长$\Delta t$与空间步长$\Delta x$的取值都提出了限制，这会极大地限制数值解的描述范围。而后向差分形式得到的迭代关系式是一个隐式，它是一种无条件稳定的，所以对于问题的解的描述自由度很大；造成这一种结果的根本原因是迭代过程中矩阵的谱半径必须要小于1才能够保证误差不会被放大。从求解精度来看，不论是前向还是后向，它们的精度都是一样的，为$O(\Delta t + (\Delta x)^2$。除了上述的两种形式外还有很多中差分形式，但是最为著名和常用的应该是Crank-Nicolson法，它是显式与隐式的一种结合体，对于空间项而言，该方法采用了平均的思想，即取$j-1$时刻与$j$时刻的二阶导的平均值来表示二阶导项，而时间项依然采用后向差分形式。$$ D(\frac{u(i-1,j-1) - 2 u(i, j - 1) + u(i+1, j-1)}{(\Delta x)^2} + \frac{u(i-1,j) - 2u(i, j)+u(i+1, j)}{(\Delta x)^2})\\ =\frac{u(i, j) - u(i, j-1)}{\Delta t} $$将$j-1$与$j$项分开，得到Crank-Nicolson法的迭代关系式为：$$ -\sigma u(i-1,j)+(2 \sigma +2) u(i, j) - \sigma u(i+1, j) = \sigma u(i-1,j-1)+ (2-2 \sigma)u(i, j-1)+ \sigma u(i+1, j-1) $$C-N法也是一种隐式形式，所以需要通过方程组的求解才能确定下一个时刻的状态值，该形式下精度为$O((\Delta t)^2+ (\Delta x)^2)$，精度比后向、前向差分形式要高，并且也是无条件稳定的。边界条件的引入PDE的求解问题一般也称为边值问题BVP（稳态）或者IBVP初边值问题（非稳态），实际上这反映了PDE的求解的核心是边界条件、初始条件。PDE反映了一类具有共性的物理场问题，而边界条件、初始条件的引入使得问题转为某一特定的问题或者研究对象。常见的边界条件有三种。第一类边界条件，也称为本质边界条件，亦称为Dirichlet条件，这一种边界条件反映了物理场状态量函数本身的取值情况，这一种边界条件的描述使得边界处的场函数值是已知的；第二类边界条件，即自然边界条件、Nuemann条件，它描述了场函数在边界处的导数值，但是边界处场函数的值来说是未知的。而有限差分在得到差分形式后，所要求解的就是各个差分点处的场函数值，而Nuemann条件使得在边界处的差分点场函数值与内部的点一样，都是未知的，所以都是自由点，这样从每一个时刻要求解的自由点的个数要比第一类边界条件下的要多1个。后向差分在Dirichlet条件下的矩阵形式对于一维空间的抛物线型PDE来说，当两个端点的边界条件均为Dirichlet条件时，如果空间步数为$n$那么，那么实际$x$方向的离散点的个数为$n+1$个，因为第一类边界条件的含义就是指定边界点的场函数取值，所以就整个过程而言，这些点的值都是已知的。由于后向差分、C-N形式都是隐式差分形式，无法显式地表示出下一时刻与上一时刻的关系，将他们的迭代关系写为矩阵形式，那么最终的求解就是反复地线性方程组求解。本质上，都是利用上一时刻的各点处的函数值去计算得到下一时刻的点的函数值。后向差分形式时，迭代矩阵可以表示为：$$ A \cdot u_{j} = u_{j - 1} + \sigma \cdot s_{j} $$其中：$$ A = \begin{bmatrix} 1 + 2 \sigma & -\sigma & 0 & \dots & 0\\ -\sigma & 1 + 2 \sigma & -\sigma & \ddots & \vdots\\ 0 & -\sigma & 1 + 2 \sigma & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & -\sigma\\ 0 & \dots & 0 & -\sigma & 1 + 2 \sigma \end{bmatrix} $$$$ s_{j}= \begin{bmatrix} u_{0j}\\ 0\\ \vdots\\ 0\\ u_{nj} \end{bmatrix} $$注意到，右侧的一维向量就是第一类边界条件的体现，采用$n$段划分时，第一个节点为0号节点，那么一共具有$n + 1$个节点，最后一个端点（右端点）为$n$号节点，那么自由点的点号范围时$1$~$n-1$，共$n - 1$个自由点。Crank-Nicolson法的矩阵形式:C-N法的迭代关系要比后向差分形式复杂，前一时刻的状态量本身也会成为一个矩阵，但是最终的计算都是归结于线性方程组的求解。$$ A \cdot u_{j} = B \cdot u_{j-1} + \sigma \cdot s $$其中：$$ A = \begin{bmatrix} 2 \sigma + 2 & -\sigma & 0 & \dots & 0\\ -\sigma & 2 + 2\sigma & -\sigma & \ddots & \vdots\\ 0 & -\sigma & 2 + 2\sigma & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & -\sigma\\ 0 & \dots & 0 & -\sigma & 2 + 2\sigma \end{bmatrix} $$$$ B = \begin{bmatrix} 2 - 2\sigma & \sigma & 0 & \dots & 0\\ \sigma & 2 - 2\sigma & \sigma & \ddots & \vdots\\ 0 & \sigma & 2 - 2\sigma &\ddots & 0\\ \vdots & \ddots & \ddots & \ddots & \sigma\\ 0 & \dots & 0 & \sigma & 2 - 2\sigma \end{bmatrix} $$$s$向量用于表示第一类边界条件的影响。$$ s = \begin{bmatrix} u_{0, j-1}\\ 0\\ \vdots\\ 0\\ u_{n, j-1} \end{bmatrix} + \begin{bmatrix} u_{0, j}\\ 0\\ \vdots\\ 0\\ u_{n, j} \end{bmatrix} $$第一类边界条件也成为本质条件，它完完全全地针对场状态变量进行限定，因此，从求解问题的角度看，第一类边界条件使得边界处的点是已知的，它们会作为系数的形式参与到计算中。正式因为这一种特性，使得在不同类型的边界条件下，隐式差分形式的计算会些许有点不同，主要是影响到自由点的个数，Nuemann条件下的计算会在以后进行介绍。实现import numpy as np # 边界条件为完全的Dirichlet条件 class Parabolic_CnMethod(object): def __init__(self, D, x_range, tmax, xn, tn): self.D = D self.x_range = x_range self.tmax = tmax self.xn = np.int32(xn) self.tn = np.int64(tn) self.dx = (x_range[-1] - x_range[0]) / self.xn self.dt = tmax / self.tn self.sigma = (self.D * self.dt) / self.dx**2 self.x = np.arange(x_range[0], x_range[-1] + self.dx, self.dx) self.t = np.arange(0, tmax + self.dt, self.dt) self.left_matrix() self.right_matrix() def left_matrix(self): # 因为边界全部为第一类边界，所以两端边界处的点是已知的 x_number = self.x.shape[0] - 2 # 自由点个数 # 自由点的个数决定了矩阵的形状大小 l_m0 = np.diag(np.ones(x_number) * (self.sigma*2 + 2), 0) l_ml = np.diag(np.ones(x_number - 1) * -self.sigma, -1) l_mu = np.diag(np.ones(x_number - 1) * -self.sigma, 1) self.l_m = l_m0 + l_ml + l_mu return self.l_m def right_matrix(self): r_m0 = np.diag(np.ones(self.x.shape[0] - 2) * (2 - 2*self.sigma), 0) r_ml = np.diag(np.ones(self.x.shape[0] - 2 - 1) * self.sigma, -1) r_mu = np.diag(np.ones(self.x.shape[0] - 2 - 1) * self.sigma, 1) self.r_m = r_m0 + r_ml + r_mu return self.r_m def solve(self, lbc, rbc, ic): # 初始化一个自由点的结果矩阵w w = np.zeros([self.t.shape[0],self.x.shape[0]-2]) w[0, :] = ic s = np.zeros(self.x.shape[0] - 2) for iters in np.arange(self.tn): s[0] = lbc[iters] + rbc[iters + 1] s[-1] = lbc[iters] + rbc[iters + 1] b = np.dot(self.r_m, w[iters, :]) + self.sigma * s w[iters + 1, :] = np.linalg.solve(self.l_m, b) ww = np.c_[lbc, w, rbc] return ww针对的方程为：$$ \frac{\partial ^2 u}{\partial x ^ 2}=\frac{\partial u}{\partial t} $$其中$D=1$，边界条件为$$ u(0, t) = 0\\ u(1, t)=0 $$初始条件：$u(x, 0)=100, (0<x<1)$。上述创建的这个类，是Crank-Nicolson法，针对的是两个边界均为第一类条件，且边界值为一个常数，如果边界条件为一个函数，那么需要稍微改造下即可。这里的计算总时间取0.3s。$x$方向划分的段数为100，总时间步取1000。demo = Parabolic_CnMethod(1, [0, 1], 0.3, 100, 1000) u = demo.solve(np.zeros(demo.t.shape[0]), np.zeros(demo.t.shape[0]), np.ones(demo.x.shape[0] - 2) * 100) import matplotlib.pyplot as plt cb = plt.pcolormesh(res.x, res.t, ww, cmap='Spectral_r', shading='gouraud') plt.contour(res.x, res.t, ww, colors='gray', levels = 10) plt.colorbar(cb) plt.show()2020！！新年快乐！
梯度下降法的缺点之前的内容已经完全学习了梯度下降算法，对于这个方法之所以在线性方程组求解领域没有更广泛地流行起来，最大的原因是每一步的迭代中近似值移动的方向是沿着负梯度方向，而相邻迭代步的残差向量(本质就是负梯度向量)$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()