分享我的学习心得 - Numerical Analysis
type: channel
link: http://www.pynumerical.com/category/Numerical-Analysis/
items: 计算级数
items: pygmsh的网格生成（2）：occ(OpenCASCADE)模块中完成
items: 二维Poisson方程有限元求解
items: Numerical Analysis: 改进的Euler法--显式梯形法
items: Numerical Analysis: 微分方程的欧拉法
items: Python实现传染病动力模型计算
items: Finite Differential Method: 隐式差分形式
items: Numerical Analysis: 共轭梯度法(2)--实现
items: Numerical Analysis: 共轭梯度法(1)--基本原理
items: Numerical Analysis: 什么是梯度下降

计算级数
type: item
date: 2024-03-21T20:04:56+08:00
description: 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的网格生成（2）：occ(OpenCASCADE)模块中完成
type: item
date: 2022-12-27T09:17:00+08:00
description: 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()

二维Poisson方程有限元求解
type: item
date: 2022-03-12T23:08:00+08:00
description: 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]^{\mathr…

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$意味着从不同的点出发开始移动。

Python实现传染病动力模型计算
type: item
date: 2020-02-13T17:21:00+08:00
description: $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 …

Finite Differential Method: 隐式差分形式
type: item
date: 2019-12-31T19:19:00+08:00
description: 抛物线方程的有限差分法-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)$。以抛物线方程为例，通过上述的表示方式来构造出P…

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…