如何在固有值问题中导入Dirichlet条件

考虑一般固有值问题

我们经常遇到要求某一固有矢量的分量为零。即 νi=0 ,这样的Dirichlet约束条件。

首先，有些要求用户控制计算流程的固有值求解器，如Arpack允许在计算过程中调整约束条件，但是这个在实现上要求比较高。本文讨论的是先调整好K, M阵，然后导入如scipy的固有值求解器

evals, evecs = scipy.linalg.eigh(K, M)

这样的做法。此时最直接的方法是

• 删除K阵和M阵的第i行和第i列。

但是这种方法一般实装不便，尤其在大型并行计算中伴随的大量内存，通信操作，效率不高。另一种方法则是

• 给 Kii 乘上一个大数P，M阵不变。

也就是所谓的罚函数法。这是实装最简单的方法。但是专业的计算软件开发者一般尽量避免采用这种方法。其导入的数值误差还在其次，它的主要问题是使矩阵的条件数变差，进而导致迭代法求解时间增加，甚至不收敛。

另一种方法是

• 把K阵的第i行和第i列清零，并设 Kii=1 . 把M阵的第i行和第i列清零。

该方法的有效性取决于固有值求解器的算法。比如说采用shift-invert计算 (K−σM)−1 的话没有问题。但是如采用Regular Inverse Mode计算 M−1K 的则会遇到矩阵奇异的错误。

下面的这种方法出自https://scicomp.stackexchange.com/questions/14429/applying-dirichlet-b-c-to-the-eigenvalue-problem

其做法是

• 把K阵的第i行和第i列清零，并设 Kii=c . 把M阵的第i行和第i列清零，并设 Mii=1 .

在这里c为用户设定的值。该方法导入了一个非现实(unphysical)的固有值c，其对应的固有矢量为 (0,0,0,......,νi=1,....0,0) .而在其他所有的固有值下对应的固有矢量中 νi=0 . 非常聪明的方法！虽然它导入了一个多余的固有值c，但是只要我们把c值取到感兴趣的固有值范围之外，该结果就不会在计算结果中出现。

在多物理场有限元软件

四元gitee.com/hillyuan_be3b/siyuan

中导入了这种方法。在用户定义的文件中

Eigen:      Method: block krylov schur      Which: SR      Sigma: 0.0      Num Eigenvalues: 10      Block Size: 10      Convergece Tolerance: 1.e-10      Maximum Iterations: 100000      K pivot: 1000.0

K pivot被用于定义上述c值。