#可以进行多重校正(扩展子空间)方法的软件包 ------ ##总体介绍 这里主要用来介绍可以实现多重校正算法的软件包,为了更容易使人理解算法的主要思想,我们这里采用Matlab语言来进行实现,具体的软件包名称为 FEM_MATLAB, 可以在这里下载这个软件包: [FEM_MATLAB下载地址](http://lsec.cc.ac.cn/~hhxie/Files/FEM_MATLAB.rar "FEM_MATLAB下载地址")。当然为了避免Matlab的速度瓶颈,这里尽量避免去使用网格单元循环的方式,所以里面形成刚度矩阵的程序不是那么直观。但是如果有必要需要去修改组装刚度矩阵的程序的话,花上一定的时间应该也能看明白的。为了能大致了解一下这个软件包,这里大概做一个介绍。 与大多数的有限元软件包的结构一样,也都是包括网格、有限单元(有限元空间)、刚度矩阵组装、线性方程求解、一些必要的后处理计算。 如果想尽快接触这个软件包进行具体的学习,可以先进入软件包,然后执行一下下面命令(每次要运行这个软件包的时候都需要执行一下): InitFEM_MATLAB 然后执行一下下面这个例子: [Err_1,Err_0,N] = Elliptic_Problem(4) 之后会得到相应的结果(二次元的结果): Err_1 = 0.12956744189696 0.03340054568577 0.00842058913824 0.00210982743517 Err_0 = 0.00378404501965 0.00047617309076 0.00005960498432 0.00000745458994 N = 32 128 512 2048 同时会有如下一幅展示收敛速度的图像输出: ![误差收敛图](http://data.xyzgate.com/9f5201bc16b25c6d2de37b40767e592b.jpg "误差收敛图") 具体的过程可以查看相应的程序: edit elliptic_problem **下面来具体结合有限元软件来讲述多重校正算法及其实现。** ##有限元网格 这里主要支持的是三角形网格,区域采用的是Matlab自带的偏微分方程计算包pdetool中的生成网格软件。为了生成网格,需要首先定义描写区域的文件,关于pdetool中区域的定义规则可以查看相关的说明书,当然也可以直接查看一下我们给出的环形区域的描述文件(在1_mesh文件夹中)。有了区域的描述文件之后,可以调用如下的程序产生相应的网格: mesh=Initial_mesh_From_Tool(Domain,hsize) 下面来介绍一下网格的数据结构:elem保存的是在每个单元上的节点编号,逆时针排列;edge是保存每个单元上的边的编号,也是逆时针排,并且满足单元局部的第一个节点的对边是单元局部的第一条边。node4edge保存了每条边的节点编号。dirichlet node和dirichlet edge保存的是的dirichlet边条件的节点信息和边的信息。其他边界条件的边用neumann来表示,因为ditichlet边条件的处理与neumann边条件的处理不相同,所以将二者区分开来。在软件包中,我们提供了对网格进行一致加密和自适应加密的功能,一致加密只需要运行如下的命令即可: mesh = refinement(mesh); 自适应加密将会结合具体的自适应算法来介绍。 另外为了之后算法设计的方便,我们同时也产生了网格边的信息,每个单元包含的三条边的信息,同时也产生了边界点的信息和边界边的信息。为了突出处理Dirichlet本质变条件,我们这里特意把这一部分的边界信息储存下来。当然同时为了后面形成刚度矩阵的计算,在网格部分也提供了计算每个单元与标准单元之间仿射变换的Jacobian矩阵和逆矩阵,同时也提供了单元的面积,边的长度、法向量和切向量的计算。 另外为了实现多重网格算法和自适应算法,我们也提供了对网格进行一致加密和自适应加密的程序,同时为了应用多重网格算法、多重校正算法,在加密的程序中同时也提供了对一般的Lagrange单元产生插值矩阵的功能。 ##单元的构造 有限单元介绍 在建立网格之后,下面介绍如何建立相应的有限元空间。简言之,有限元空间就是在有限元网格上赋予一定的单元定义和边界信息。单元的定义是通过标准单元上基函数的定义以及基函数在一般单元和标准单元之间的关系,这里采用的Lagtange插值基函数只需要仿射变换,比较适合于Matlab语言的实现。 目前我们主要实现了基于插值的Lagrange型的基函数。基函数的定义结合网格就形成有了我们通常认为的有限元空间。下面是一个最简单的协调线性有限元的定义: %triangle Base.b_3 = @base_Lg_3_1; Base.gb_3 = @grad_base_Lg_3_1; Base.Hessian_3 = @Hessian_3_1; Base.gb_multinode_3 = @grad_base_Lg_3_1_multinode; Base.DOF_3 = [1 0 0]; Base.scale_3 = [0 0 0]; Base.nodal_3 = nodal_Lg_3_1; Base.nodal_funct_3 = @nodal_funct_Lg_3_1; %-------------Standard Lagrange P_1 element------------------- function v=base_Lg_3_1(x) v=[1-x(:,1)-x(:,2),x(:,1),x(:,2)]; function v=grad_base_Lg_3_1(x) %give the gradient of the basis of function v=[-1 1 0; -1 0 1]; function [vx,vy] = grad_base_Lg_3_1_multinode(x) %give the gradient of the basis of function for multi nodes dim =size(x,1); vx=ones(dim,1)*[-1 1 0]; vy=ones(dim,1)*[-1 0 1]; function v= Hessian_3_1(x) %give the gradient of the basis of function v=[0 0 0; 0 0 0; 0 0 0; 0 0 0]; % give the nodal functional of this type of base function v = nodal_Lg_3_1 %nrNodes = size(mesh.node,1); v = [ 0 0; 1 0; 0 1]; % give the nodal functional of this type of base function v = nodal_funct_Lg_3_1(f,mesh) v = f(mesh.node); 其中的基函数是定义在参考单元上的基函数,再接合仿射变换就可以定义在一般的三角形单元上的基函数了,进而可以定义出有限元空间。 为了使得处理有限元单元的统一,这里在给单元的定义的时候除了给出通常的单元基函数的多项式描述,自由度在单元点线面上分布的同时也给出了与单元基函数相对应的泛函定义,这样比较方便地进行有限元插值和在多重网格算法中给出插值矩阵。同时也有便于实现多重网格方法的时候对一般有限元方法形产生插值矩阵和限制矩阵。当然为了以后升级的可能,我们这里也对单元的类型进行了一定的藐视,如何单元基函数是按照积分的方式来定义的,就取负数,如果是按照微分的方式来定义的就取相应的正数,如果是按照函数值的方式定义的就设定为零。 ##形成刚度矩阵: 之前形成的网格和有限元空间都是为了可以进行有限元刚度矩阵的组装操作。为了实现多重校正算法,我们需要在取值函数空间和试探函数空间不同的情况下组装矩阵,这个不同不仅指的是所采用的基函数不同,同时所依赖的网格有可能不同。在多重校正算法中往往是最细网格和最粗网格联合形成刚度矩阵。为了达到这个目的,我们在加密网格的时候,在每个加密的单元上保存父亲单元上的祖先单元的信息,同时也产生加密后的单元节点在祖先单元中的重心坐标。在本软件包中,我们实现了这种矩阵的组装功能,具体可以查看组装刚度矩阵程序: A = Matrix_Assemble(meshs,Matrix) ##代数求解: 形成刚度矩阵和右端项之后接下来就是进行代数求解了。当然最直接的求解方法是利用线性方程组的直接求解方法进行求解。当然我们这里也提供了多重网格方法来求解线性方程组。 ##具体的算例介绍 上面只是介绍了软件包的大概,如果有兴趣学习使用和研究多重校正算法的研究者可以结合下面的实例来进行学习。发布本软件包的主要原因是可以设计求解非线性问题和特征值问题的多重校正算法的。对于具体的算法介绍可以参看后面的参考文献。 ---------- ## 线性特征值问题 ## 为了简单起见,我们求解如下的线性特征值问题 ```math \begin{cases}{} -\Delta u=\lambda u,\ \ \ {\rm in}\ \Omega,\\ u=0,\ \ \ \ \ {\rm on}\ \partial\Omega,\\ \int_\Omega u^2d\Omega =1. \end{cases} ``` 为了进行有限元离散,我们先定义上面特征值问题相应的变分形式: 求$$(\lambda,u)\in\mathbb R\times H_0^1(\Omega)$$使得$$\\| u\\|_0=1$$且满足如下的方程 ```math (\nabla u,\nabla v) = \lambda(u,v),\ \ \ \forall v\in H_0^1(\Omega). ``` 采用最简单的线性有限元进行离散,用$$V_h$$来表示线性有限元空间。则具体的离散格式如下: 求$$(\bar\lambda_h,\bar u_h)\in\mathbb R\times V_h$$使得$$\\|\bar u_h\\|_0=1$$且满足如下的方程 ```math (\nabla \bar u_h,\nabla v_h) = \bar\lambda_h(\bar u_h,v_h),\ \ \ \forall v_h\in V_h. ``` 这里我们采用$$(\bar \lambda_h,\bar u_h)$$来表示精确的有限元离散逼近(与后面的利用迭代方法得到的逼近值$$(\lambda_h,u_h)$$相区别), 整个网页其实都是在讨论如何能够快速地求解出这个精确的有限元离散特征对,或者说目的是构造出相应的代数解法器。 为了更好的理解后面的算法,这里有必要介绍一下最基本的特征值问题有限元离散的误差估计: ```math \|u-\bar u_h\|_1 \leq (1+C\eta_a(V_h))\inf\limits_{v_h\in V_h}\|u-v_h\|_1, ``` ```math \|u-\bar u_h\|_0 \leq C\eta_a(V_h)\|u-u_h\|_1, ``` ```math 0<\bar\lambda_h-\lambda <\|u-u_h\|_1^2. ``` 其中$$\eta_a(V_h)$$的定义如下(与通常边值问题有限元方法中Aubin-Nitsche技巧一样) ```math \eta_a(V_h) = \sup\limits_{f\in L^2(\Omega),\|f\|_0=1}\inf\limits_{v_h\in V_h}\|Tf-v_h\|_1. ``` 这里的算子$$T$$定义如下 ```math (\nabla Tf,v) = (f,v),\ \ \ \forall v\in H_0^1(\Omega). ``` 当然上面的结论对于特征值问题一般的有限维逼近空间逼近往往也成立,在后面的论述中我们就会用到在$$V_h$$的一个子空间中逼近特征对$$(\bar\lambda_h, \bar u_h)$$的误差估计,其结果跟这里的表述完全类似。 为了进行多重校正算法,我们需要采用与求解边值问题的多重网格方法一样的方式来产生嵌套的网格。假设现有一个粗网格$$\mathcal T_H$$, 然后对其进行一致加密,得到一系列的嵌套网格$$\mathcal T_{h_1}, \cdots, \mathcal T_{h_n}$$,加密到足够细的网格停止。 由于是一致加密,所以网格尺寸满足如下的关系: ```math \frac{h_k}{h_{k+1}}=\beta,\ \ \ \ \ k=1, \cdots, n-1. ``` 为了能够设计扩展子空间方法,我们需要保存每个细网格单元在粗网格单元中的重心坐标,同时也在最粗网格和细网格上定义有限元空间$$V_H$$和$$V_{h_1}, \cdots, V_{h_n}$$。对于每个有限元空间定义如下对特征值问题的逼近精度: ```math \delta_{h_k}(\lambda) = \sup\limits_{u\in M(\lambda)}\inf\limits_{v_{h_k}\in V_{h_k}}\|u-v_{h_k}\|_1, ``` 其中$$M(\lambda)$$表示相对于特征值$$\lambda$$的特征空间. 为了简单起见,我们这里考虑在凸区域上的特征值问题,所以有限元的逼近精度满足如下的性质: ```math \frac{\delta_{h_k}(\lambda)}{\delta_{h_{k+1}}(\lambda)}=\beta,\ \ \ \ \ k=1, \cdots, n-1. ``` 依据定义的有限元空间系列,我们可以定义求解特征值问题的扩展子空间迭代步。假设现有一个定义在细有限元空间内的特征函数逼近$$(\lambda_{h_k},u_{h_k})\in \mathbb R\times V_{h_k}$$。扩展子空间迭代步定义如下: **算法1: 扩展子空间迭代算法** 1. 求解如下边值问题:求$$\widehat u_{h_{k+1}}\in V_{h_{k+1}}$$满足 ```math (\nabla\widehat u _{h_{k+1}},\nabla v_{h_{k+1}})= \lambda_{h_k}b(u_{h_k},v_{h_{k+1}}), \ \ \ \forall v_{h_{k+1}}\in V_{h_{k+1}}. ``` 2. 定义子空间 **$$V_{H,h_{k+1}}=V_H+{\rm span}\\{\widehat u_{h_{k+1}}\\}$$**, 然后在这个低维子空间上求解如下的特征值问题:求$$(\lambda_{h_{k+1}},u_{h_{k+1}})\in \mathbb R\times V_{H,h_{k+1}}$$使得$$\|u_{h_{k+1}}\|_0=1$$且满足如下方程 ```math (\nabla u_{h_{k+1}},\nabla v_{H,h_{k+1}}) = \lambda_{h_{k+1}}(u_{h_{k+1}},v_{H,h_{k+1}}),\ \ \ \forall v_{H,h_{k+1}}\in V_{H,h_{k+1}}. ``` 得到新的特征对$$(\lambda_{h_{k+1}},u_{h_{k+1}})$$,为了后面行文的方便,我们把上面算法记为: ```math (\lambda_{h_{k+1}},u_{h_{k+1}}) = {\tt Correction}(\lambda_{h_k},u_{h_k},V_{h_k}, V_{h_{k+1}},V_H). ``` 从算法1中我们可以看到,由于扩展子空间$$V_{H,h_{k+1}}=V_H+{\rm span}\\{\widehat u_{h_{k+1}}\\}$$的维数相对细有限元空间$$V_{h_{k+1}}$$的维数很小,所以求解第二步的特征值问题的计算量将会很小,即迭代算法的计算量主要是在第一步中的求解定义在细有限元空间的线性边值问题。由于这个问题是一个标准的Laplace方程,可以应用几乎所有的线性解法器来求解(最容易想到的自然是多重网格算法)。 **需要特别指出的是,这里的扩展子空间算法与一般的扩展方法不一样的。一般的扩展方式可以说是Krylov式的,即通过矩阵乘向量产生的空间, i.e., $$P_k(A)W$$, 这里$$P_k$$表示$$k$$次多项式, $$W$$表示某些向量组成的矩阵。这里的扩展方式没有采取Krylov的方式,而是采用了多重网格方法中的粗空间来扩展。一方面其具有稀疏性,另一方面,粗空间又能对特征空间进行一定的逼近。** 详见文章:Yunhui He, Qichen Hong, Hehu Xie, Meiling Yue and Chunguang You, [Energy Error Estimates of Subspace Method and Multigrid Algorithm for Eigenvalue Problems](https://arxiv.org/abs/1705.03038 "Energy Error Estimates of Subspace Method and Multigrid Algorithm for Eigenvalue Problems"), arXiv:1705.03038. 下面可以通过简单地分析来揭示这个迭代步对近似特征函数逼近的精度改进。 **定理 1:** 假设给定的特征对$$(\lambda_{h_k},u_{h_k})\in\mathbb{R}\times V_{h_k}$$满足如下的条件 ```math \|u-u_{h_k}\|_1 \lesssim \varepsilon_{h_k}(\lambda), ``` ```math \|u-u_{h_k}\|_{-1} \lesssim \eta_a(H)\|u-u_{h_k}\|_1, ``` ```math |\lambda-\lambda_{h_k}| \lesssim \varepsilon_{h_k}^2(\lambda). ``` 则经过扩展子空间迭代之后得到的特征对, $$(\lambda_{h_{k+1}},u_{h_{k+1}})\in\mathcal{R}\times V_{h_2}$$ 具有如下的误差估计 ```math \|u-u_{h_{k+1}}\|_1 \lesssim \varepsilon_{h_2}(\lambda), ``` ```math \|u-u_{h_{k+1}}\|_{-1} \lesssim \eta_a(H)\|u-u_{h_{k+1}}\|_1,\ \ \ ({\rm important}) ``` ```math |\lambda-\lambda_{h_{k+1}}|\lesssim \varepsilon_{h_{k+1}}^2(\lambda), ``` 其中 $$\varepsilon_{h_{k+1}}(\lambda):=\eta_a(H)\varepsilon_{h_k}(\lambda)+\delta_{h_{k+1}}(\lambda)$$. 关于上面定理的证明请参加后面的参考文献。 基于上面定义的扩展子空间迭代算法,我们可以在嵌套有限元空间系列上定义求解特征值问题的多重网格算法。具体的算法定义如下: **算法2: 求解特征值问题的多重网格算法** 1. 在初始有限元空间$$V_{h_1}$$上求解特征值问题: ```math (\nabla u_{h_1},\nabla v_{h_1}) = \lambda_{h_1}(u_{h_1},v_{h_1}),\ \ \ \forall v_{h_1}\in V_{h_1}. ``` 2. 对 $$k=1,\cdots, n-1$$进行如下的迭代 ```math (\lambda_{h_{k+1}},u_{h_{k+1}}) = {\tt Correction}(\lambda_{h_k},u_{h_k},V_{h_k}, V_{h_{k+1}},V_H). ``` 最后得到在最细有限元空间$$V_{h_n}$$上的特征函数逼近$u_{h_n}$和相应的特征值逼近$$\lambda_{h_n}$$. **定理 2** 假设条件 $$C\beta \bar\lambda \eta_a(V_H)<1$$成立, 运行算法2之后得到的特征对逼近$$(\lambda_{h_n}, u_{h_n})$$满足如下的误差估计 ```math \|u-u_{h_n}\|_1 \lesssim \delta_{h_n}(\lambda), ``` ```math \|u-u_{h_n}\|_0 \lesssim \eta_a(V_H)\|u-u_{h_n}\|_1, ``` ```math |\lambda-\lambda_{h_n}|\lesssim \|u-u_{h_n}\|_1^2. ``` 在软件包可以运行如下的程序看看上面定义的多重网格算法的效果: ```` [LAM,LAM_Dir,Dis_1,Dis_0,N] = Eigen_Multigrid(4) ```` 下面是多重网格得到特征值的逼近: ``` LAM = 20.50549180291402 19.92984722056495 19.78679654626770 19.75110057778336 ``` 为了比较,我们同时计算了直接求解方法得到的特征值逼近: ``` LAM_Dir = 20.50549180291402 19.92977699513403 19.78678910433666 19.75110004218717 ``` 输出直接方法和多重网格方法得到的特征函数之间的误差:$$H_1$$范数 ``` Dis_1 = 0 0.00425240074569 0.00139879989711 0.00037663782628 ``` $$L^2$$范数 ``` Dis_0 = 1.0e-003 *[ 0 0.19353715290493 0.07177703203547 0.02019913961686]``` 每层网格上的单元个数 ``` N = 128 512 2048 8192 ``` 同时也会输出如下的误差收敛图: ![多重网格误差图](http://data.xyzgate.com/51f5832c9a668196f4a5ebcacd2181c2.jpg "多重网格误差图") 由此可以看到,多重网格方法得到的特征值和特征函数与直接算法得到误差差不多的逼近值,并且同样具有最优的收敛速度。 ### 做为一个特征值代数解法器的形式 上面设计的扩展子空间算法及由之够造的多重网格迭代算法的目的是为了得到与直接求解有限元离散特征值问题相同的收敛速度,或者说是为了得到具有相同收敛速度的特征对的逼近。我们知道最后求解出来的特征对的逼近与精确解之间的误差其实是包含了离散误差$$\\|u-\bar u_{h_n}\\|_1$$和代数误差$$\\|\bar u_{h_n}-u_{h_n}\\|_1$$这两部分的。上面的算法其实是为了去求解一个逼近有限元离散解的解法,那么自然的问题是上面的算法是否可以用来设计出一个代数解法器使得代数精度可以任意小以至于达到机器精度。观察**算法1**可以发现我们最重要的一步就是在扩展子空间$$V_{H,h_{k+1}}$$上求解一个特征值问题,这样使得特征函数逼近在$$L^2(\Omega)$$范数下的误差比在$$H_0^1(\Omega)$$范数下的误差具有更高的收敛苏。基于这样的理解我们可以在同一个网格层上执行多次的扩展子空间迭代算法。为了表示多次迭代,我们用$$u_{h_k}^{(\ell)}$$表示在空间$$V_{h_k}$$上第$$\ell$$次迭代的特征函数逼近。 假如已经有一个特征对逼近$$(\lambda_{h_k}^{(\ell)}, u_{h_k}^{(\ell)}$$, 则具体的迭代算法构造如下: **算法3: 扩展子空间迭代算法** 1. 求解如下边值问题:求$$\widehat u_{h_{k}}^{(\ell+1)}\in V_{h_{k+1}}$$满足 ```math (\nabla\widehat u _{h_k}^{(\ell+1)},\nabla v_{h_k})= \lambda_{h_k}^{(\ell)}b(u_{h_k}^{(\ell)},v_{h_k}), \ \ \ \forall v_{h_k}\in V_{h_k}. ``` 我们这里采用多重网格迭代(以$$u_{h_k}^{(\ell)}$$为初值)求解上面的线性方程得到一个逼近解$$\widetilde u_{h_k}^{(\ell+1)}$$. 2. 定义子空间 **$$V_{H,h_k}=V_H+{\rm span}\\{\widetilde u_{h_k}^{(\ell+1)}\\}$$**, 然后在这个地位子空间上求解如下的特征值问题:求$$(\lambda_{h_k}^{(\ell+1)},u_{h_k}^{(\ell+1)})\in \mathbb R\times V_{H,h_k}$$使得$$\|u_{h_k}^{(\ell+1)}\|_0=1$$且满足如下方程 ```math (\nabla u_{h_k}^{(\ell+1)},\nabla v_{H,h_k}) = \lambda_{h_k}^{(\ell+1)}(u_{h_k}^{(\ell+1)},v_{H,h_k}),\ \ \ \forall v_{H,h_k}\in V_{H,h_k}. ``` 得到新的特征对$$(\lambda_{h_{k+1}},u_{h_{k+1}})$$,为了后面行文的方便,我们把上面算法记为: ```math (\lambda_{h_k}^{(\ell+1)},u_{h_k}^{(\ell+1)}) = {\tt CorrectionStep}(\lambda_{h_k}^{(\ell)},u_{h_k}^{(\ell)},V_{h_k}, V_H). ``` **定理 3:** 假设给定的特征对$$(\lambda_{h_k}^{(\ell)},u_{h_k}^{(\ell)})\in\mathbb{R}\times V_{h_k}$$满足如下的条件 ```math |\bar\lambda-\lambda_{h_k}^{(\ell)}|+\|\bar u-u_{h_k}^{(\ell)}\|_{-1} \lesssim \eta_a(H)\|\bar u_{h_k}-u_{h_k}^{(\ell)}\|_1, ```` 则经过扩展子空间迭代之后得到的特征对, $$(\lambda_{h_k}^{(\ell+1)},u_{h_k}^{(\ell+1)})\in\mathcal{R}\times V_{h_k}$$ 具有如下的误差估计 ```math \|\bar u_{h_k}-u_{h_k}^{(\ell+1)}\|_1 \leq \gamma\|\bar u_{h_k}-u_{h_k}^{(\ell)}\|_1, ``` ```math \|\bar u_{h_k}-u_{h_k}^{(\ell+1)}\|_{-1} \lesssim \eta_a(H)\|\bar u_{h_k}-u_{h_k}^{(\ell+1)}\|_1,\ \ \ ({\rm important}) ``` ```math |\lambda-\lambda_{h_{k+1}}|\lesssim \|\bar u_{h_k}-u_{h_k}^{(\ell+1)}\|_1^2 \leq \eta_a(H)\|\bar u_{h_k}-u_{h_k}^{(\ell+1)}\|_1, ``` 其中$$ \gamma:= \big(\theta+(1+\theta)C\eta_a(V_H)\big)\big(1+C\eta_a(V_H)\big)$$. 关于上面定理的证明请参加后面的参考文献。 基于上面定义的扩展子空间迭代算法,我们可以在嵌套有限元空间系列上定义求解特征值问题的多重网格算法。具体的算法定义如下: **算法4: 求解特征值问题的多重网格算法** 1. 在初始有限元空间$$V_{h_1}$$上求解特征值问题: ```math (\nabla u_{h_1},\nabla v_{h_1}) = \lambda_{h_1}(u_{h_1},v_{h_1}),\ \ \ \forall v_{h_1}\in V_{h_1}. ``` 2. 对 $$k=2,\cdots, n$$进行如下的迭代 2.1. 定义 $$(\lambda_{h_k}^{(0)},u_{h_k}^{(0)}):=(\lambda_{h_{k-1}},u_{h_{k-1}})$$ 2.2. 对$$\ell=0,\cdots, p_k$$进行如下的循环 ```math (\lambda_{h_k}^{(\ell+1)},u_{h_k}^{(\ell+1)}) = {\tt CorrectionStep}(\lambda_{h_k}^{(\ell)},u_{h_k}^{(\ell)},V_{h_k}, V_H). ``` 2.3. 把$$(\lambda_{h_k}^{(p)},u_{h_k}^{(p)})$$作为本层的特征对输出. 最后得到在最细有限元空间$$V_{h_n}$$上的特征对逼近$$(\lambda_{h_n},u_{h_n})$$. **定理 4** 假设条件 $$\gamma^p\beta<1$$成立, 运行算法2之后得到的特征对逼近$$(\lambda_{h_n}, u_{h_n})$$满足如下的误差估计 ```math \|\bar{u}_{h_n}-u_{h_n}\|_1 \leq \frac{2\gamma^p\beta}{1-\gamma^p\beta}\delta_{h_n}(\lambda), ``` ```math \|\bar u_{h_n}-u_{h_n}\|_{-1} \lesssim \eta_a(V_H)\|\bar u_{h_n}-u_{h_n}\|_1, ``` ```math |\bar\lambda_{h_n}-\lambda_{h_n}|\lesssim \|\bar u_{h_n}-u_{h_n}\|_1^2. ``` 从上面的定理可以看出算法4收敛的条件$$\gamma^p\beta<1$$中是关于$$\gamma$$的$$p$$次方,所以我们可以控制$$p$$使得这个$$\gamma^p$$到我们所需要的收敛速度。其实在实际计算中$$p$$往往取$$1$$或$$2$$就已足够,因为在算法3中求解线性方程组的时候用的是多重网格迭代,我们知道其具有一致的收敛速度,所以$$\theta$$是会比较小的(如有必要可以多做几次多重网格迭代,但其实往往做一次就可以了)。由此可以知道条件$$\gamma^p\beta<1$$是很容易满足的。另外,代数精度$$\|\bar{u}_{h_n}-u_{h_n}\|_1$$有如下控制 ```math \frac{2\gamma^p\beta}{1-\gamma^p\beta}\delta_{h_n}(\lambda), ``` 则可以知道我们提高迭代次数$$p$$就可以使得代数精度任意小。 在软件包中可以运行如下的程序看看上面定义的多重网格算法的效果: ### 特征值代数解法器更直接的形式 从上面算法的定义中可以看出来,在网格层$$V_{h_1},\cdots,V_{h_{n-1}}$$上即使我们把代数误差降低到机器精度,然后将所得的特征函数逼近插值到最细空间$$V_{h_n}$$上作为初值,这个时候初值的代数精度也是由$$h_n$$(在能量范数意义下)决定,并不会因为上一层的代数精度很小而使得本层初值的精度提高。所以,更有效率的方式就是在粗网格层上$$\mathcal T_{h_1},\cdots,\mathcal T_{h_{n-1}}$$只进行一次的扩展子空间迭代,只需要能得到在最细网格层上具有精度$$h_n$$的初值即可。然后直接在最细网格层$$\mathcal T_{h_n}$$上进行多次的扩展子空间迭代来提高代数精度。如此我们可以定义一种更有效率的多重网格算法如下: **算法5: 作为特征值代数解法器的多重网格算法** 1. 在初始有限元空间$$V_{h_1}$$上求解特征值问题: ```math (\nabla u_{h_1},\nabla v_{h_1}) = \lambda_{h_1}(u_{h_1},v_{h_1}),\ \ \ \forall v_{h_1}\in V_{h_1}. ``` 2. 对 $$k=2,\cdots, n-1$$进行如下的迭代 2.1. 定义 $$(\lambda_{h_k}^{(0)},u_{h_k}^{(0)}):=(\lambda_{h_{k-1}},u_{h_{k-1}})$$ 2.2. 执行如下的扩展子空间迭代步 ```math (\lambda_{h_k}^{(1)},u_{h_k}^{(1)}) = {\tt CorrectionStep}(\lambda_{h_k}^{(0)},u_{h_k}^{(0)},V_{h_k}, V_H). ``` 2.3. 把$$(\lambda_{h_k}^{(1)},u_{h_k}^{(1)})$$作为本层的特征对输出. 3. 在最细的网格层上进行如下的运算 3.1. 定义 $$(\lambda_{h_n}^{(0)},u_{h_n}^{(0)}):=(\lambda_{h_{n-1}},u_{h_{n-1}})$$. 3.2. 对$$\ell=0,\cdots, p_k$$进行如下的循环 ```math (\lambda_{h_n}^{(\ell+1)},u_{h_n}^{(\ell+1)}) = {\tt CorrectionStep}(\lambda_{h_n}^{(\ell)},u_{h_n}^{(\ell)},V_{h_n}, V_H). ``` 3.3. 把$$(\lambda_{h_n}^{(p)},u_{h_n}^{(p)})$$作为本层的特征对输出. 最后得到在最细有限元空间$$V_{h_n}$$上的特征对逼近$$(\lambda_{h_n},u_{h_n})$$. 在上面的算法定义中,为了提高最后特征对逼近$$(\lambda_{h_n},u_{h_n})$$的代数精度,我们只需要在最细网格层上多做扩展子空间迭代即可。 在软件包可以运行如下的程序看看上面定义的多重网格算法的效果: ```` [LAM,LAM_Dir,Dis_1,Dis_0,Err,Err_dir,N]=Full_Eigen_Multigrid(5) ```` 这里取得$$V_H$$为一直三角形剖分且$$H=1/8$$, 到最细网格层之前,在每层网格上制作一次的扩展子空间迭代,然后在最细网格层上做7次的扩展子空间迭代。 下面是得到特征值的逼近: ``` LAM = 20.50554489770677 19.93269612856454 19.78882369326495 19.75169551743347 19.74218157151798 ``` 为了比较,我们同时计算了直接求解方法得到的特征值逼近: ``` LAM_Dir = 20.50554489770677 19.92978984221579 19.78679229019024 19.75110083703540 19.74218157147764 ``` 输出直接方法和在最细层网格上迭代得到的特征函数之间的误差:$$H_1$$范数 ``` Dis_1 = 0.00597295493595 0.00102759047565 0.00027726384660 0.00008746233116 0.00002858902923 0.00000950914893 0.00000320004740 ``` 相应的 $$L^2$$范数 ``` Dis_0 = 1.0e-004 *[0.52153135178134 0.09590753872280 0.02164326307480 0.00580079567563 0.00171249352253 0.00053881759789 0.00017475733101] ``` 每层网格上的多重网格算法的$$H^1_0(\Omega)$$误差: ``` Err = 0.43298803620732 0.21952099806761 0.11134133416743 0.05587227968197 0.02726042207640``` 每层网格上精确特征函数逼近的$$H^1_0(\Omega)$$误差: ```` Err_dir = 0.43298803620732 0.21769522738474 0.10899564658228 0.05451624413096 0.02726042191283 ``` 每层网格上的单元个数 ``` N = 128 512 2048 8192 32768 ``` 从上面的数值结果可以发现,在最细网格层上多进行扩展子空间迭代就可以提高代数精度。即使在最细网格层上代数精度很高,但是离散误差并没有提高。 同时也会输出如下代数精度与迭代次数之间的关系图: ![代数解法器的误差图](http://data.xyzgate.com/fb6e49ef987faa4c02b522c7de02c90b.jpg "代数解法器的误差图") 由此可以看到,在最细网格层上进行扩展子空间迭代可以提高特征值和特征函数的代数精度。 当然我们也可以通过数值实验发现,一方面如果我们在每次扩展子空间迭代的时候,对线性方程组做更多次的迭代求解(或更多次的前后光滑),比如在每次迭代的时候进行1次的MG迭代,而在每次迭代的时候前后光滑都用2次的CG迭代,这样就可以得到如下的代数误差收敛图: ![更多次CG迭代的代数解法器误差图](http://data.xyzgate.com/e24ea66bb509d7aba0c8b2f294c8440c.jpg "更多次CG迭代的代数解法器误差图") 另一方面,也可以把网格$$\mathcal T_H$$的尺寸变小,收敛速度也会更快。比如我们把$$\mathcal T_H$$的网格尺寸变成$$\frac{1}{16}$$,在每一次迭代中,执行MG迭代求解线性方程组的时候,前后光滑均做两步CG迭代,则会得到如下的代数精度与迭代次数图像: ![更快的代数解法器误差图](http://data.xyzgate.com/50bf91a1c600fc985c10a000aba0804e.jpg "更快的代数解法器误差图") 如果更进一步,每次求解线性方程组的时候进行2次MG迭代且前后光滑的时候做2次CG迭代,收敛速度将会更进一步提高。数值结果如下图: ![更多MG迭代更快的的代数解法器误差图](http://data.xyzgate.com/4ae000aeb1e204636d1d7230a240a257.jpg "更多MG迭代更快的的代数解法器误差图") 虽然上面的数值算例展示的是只计算一个特征值,但是上面的算法都是很自然的可以推广到计算多个特征值的(对于一般的算法这其实是很自然的)。比如可以使用如下的命令计算Laplace特征值问题的前$$13$$个特征值: ``` [LAM,LAM_Dir,Dis_1,Dis_0,Dis_LAM,N]=Full_Eigen_Multigrid_Multi(5) ``` 这里我们取$$\mathcal T_H$$的网格尺寸为$$\frac{1}{16}$$, 下图展示的是前$$13$$个特征值的代数精度的收敛速度: ![前13个特征值的收敛速度](http://data.xyzgate.com/172529589e89e57959b50ec38b448c95.jpg "前13个特征值的收敛速度") 当然把网格$$\mathcal T_H$$变小同时也会增加每一次迭代步中在低维空间中求解特征值问题的计算量,但是这一部分的时间如果相对于最细网格上的计算量是小量的话,这样的代价是值得的。尤其是在并行计算中,根据我们的计算经验,粗网格尺寸应该尽量取得小一些,这样才能保证并行的效率。否则,粗网格很小,然后这么小的维数分布在很多的进程中,这样的维数分布下,数据传输的时间将占大比率而使得并行效率下降。所以并行计算的时候$$\mathcal T_H$$的网格尺寸尽量取得小一些,这样一方面可以提高并行的效率,同时也可以加快扩展子空间迭代的收敛速度。 ##参考文献: 18. Xiaole Han, Hehu Xie and Fei Xu, A cascadic multigrid method for eigenvalue problem, Journal of Computational Mathematics, 35(1) (2017), 7-90. 17. Shanghui Jia, Hehu Xie, Manting Xie and Fei Xu, A full multigrid method for nonlinear eigenvalue problems, Sci China Math, 59 (2016), 2037–2048. 16. Hehu Xie and Manting Xie, A multigrid method for ground state solution of Bose-Einetein condensates, Commun. Comput. Phys., 19(3) (2016), 648-662. 15. Hongtao Chen, Hehu Xie and Fei Xu, A full multigrid method for eigenvalue problems, Journal of Computational Physics,322 (2016),747-759. 14. Xiaole Han, Yu Li, Hehu Xie and Chunguang You, Local and parallel finite element algorithm based on multilevel discretization for eigenvalue problems, International Journal of Numerical Analysis & Modeling, 13(1) (2016),73-89. 13. Hehu Xie and Tao Zhou, A multilevel finite element method for Fredholm integral eigenvalue problems,Journal of Computational Physics,303 (2015),173-184 12. Hongtao Chen,Yunhui He, Yu Li and Hehu Xie, A multigrid method for eigenvalue problems based on shifted-inverse power technique. European Journal of Mathematics, 1(1) (2015), 207-228. 11. Xiaole Han, Yu Li and Hehu Xie, A multilevel correction method for Steklov eigenvalue problem by nonconforming finite element methods, Numerical Mathematics: Theory, Methods and Applications, 8(03) (2015), 383-405. 10. 谢和虎, 非线性特征值问题的多重网格算法,中国科学:数学,45(8) (2015),1193-1204. Hehu Xie, A multigrid method for nonlinear eigenvalue problems (in Chinese), Sci Sin Math, 45 (2015), 1193-1204 9. Qun Lin, Hehu Xie and Fei Xu, Multilevel correction adaptive finite element method for semilinear elliptic equation, Application of Mathematics, 60(5) (2015), 527-550. 8. Wei Gong, Hehu Xie, Ningning Yan, A multilevel correction method for optimal controls of elliptic equation, SIAM J. Sci. Comput., 37(5) (2015), A2198-A2221, 2015. 7. Qun Lin, Fusheng Luo and Hehu Xie, A multilevel correction method for Stokes eigenvalue problems and its application, Mathematical Methods in the Applied Sciences, 38(18) (2015), 4540-4552. 6. Qun Lin and Hehu Xie, A multi-level correction scheme for eigenvalue problems, Mathematics of Computation, 84(291) (2015), 71-88. 5. Hehu Xie, A type of multilevel method for the Steklov eigenvalue problem, IMA J. Numer. Anal., 34(2) (2014), 592-608. 4. Hehu Xie, A multigrid method for eigenvalue problem, Journal of Computational Physics, 274(1), 2014, pages 550-561. 3. Xia Ji, Jiguang Sun and Hehu Xie, A multigrid method for Helmholtz transmission eigenvalue problems, Journal of Scientific Computing, 60(2), pages: 276-294, 2014. 2. Qun Lin and Hehu Xie, A multilevel correction type of adaptive finite element method for Steklov eigenvalue problems, Proceedings of the International Conference Applications of Mathematics 2012, pp. 134-143. 1. Qun Lin and Hehu Xie, An observation on Aubin-Nitsche Lemma and its applications, Mathematics in Practice and Theory, 41(17) (2011), 247-258.(In Chinese, English title and abstract)
/
Go to slide: