Legendre-Galerkin方法以及Chebyshev-Legendre-Galerkin方法以及Python实现中给出了 Legendre-Galerkin 方法, 本文比较二维 Legendre-Galerkin方法一类快速实现和直接解法的效率。

考虑如下 Hemholtz 方程
−Δu+αu=f,∇u⋅n=0.\begin{aligned} -\Delta u + \alpha u = f, \\ \nabla u \cdot \mathbf{n} = 0. \end{aligned} −Δu+αu=f,∇u⋅n=0.

Dirichlet边界条件下的算法:

最近考虑Neumann 边界问题的 Legendre-Galerkin 方法求解, 由于书中只考了 Dirichlet 边界条件下的特征值方法处理，通过查阅资料，整理了Neumann 边界条件下 Legendre-Galerkin方法的快速处理方法。
对于二维的方程, Legendre-Galerkin 方法的基函数空间为两个一维的基函数空间直接张量积得到。离散后的线性系统如下：

SUM⊤+MUS⊤+αMUM⊤=F.SUM^\top + MU S^\top + \alpha M U M^\top = F. SUM⊤+MUS⊤+αMUM⊤=F.

其中SSS和MMM为一维的刚度矩阵和质量矩阵。对于 Dirichlet 边界问题, 求解上述问题的技巧如下：

Step1. 求解特征值问题:
Mx=λSx.M x = \lambda Sx. Mx=λSx.
记解得特征向量组成的矩阵为 EEE (列向量作为特征向量), 特征值组成的对角矩阵记为 Λ\LambdaΛ, 则有

ME=SEΛ.M E = SE \Lambda. ME=SEΛ.

Step2. 令 U=EUE⊤U = E \mathcal{U} E^\topU=EUE⊤, 代入到方程中:
SEUE⊤M⊤+MEUE⊤S⊤+αMEUE⊤M⊤=F.\begin{aligned} SE \mathcal{U} E^\top M^\top + M E\mathcal{U}E^\top S^\top + \alpha ME \mathcal{U} E^\top M^\top = F. \end{aligned} SEUE⊤M⊤+MEUE⊤S⊤+αMEUE⊤M⊤=F.
利用特征值和特征向量改写,注意到 Legendre-Galerkin 方法刚度矩阵是对角的,不妨记S=ΛSS = \Lambda_SS=ΛS, 从而有:
ΛSEUΛE⊤ΛS+ΛSEΛUE⊤ΛS+αΛSEΛUΛE⊤ΛS=F.\Lambda_S E \mathcal{U} \Lambda E^\top\Lambda_S + \Lambda_S E \Lambda \mathcal{U} E^\top \Lambda_S + \alpha \Lambda_S E \Lambda \mathcal{U} \Lambda E^\top \Lambda_S = F. ΛSEUΛE⊤ΛS+ΛSEΛUE⊤ΛS+αΛSEΛUΛE⊤ΛS=F.
进一步化简:
UΛ+ΛU+αΛUΛ=E⊤ΛS−1FΛS−1E.\mathcal{U} \Lambda + \Lambda \mathcal{U} + \alpha \Lambda \mathcal{U} \Lambda = E^\top \Lambda_S^{-1} F \Lambda_S^{-1} E. UΛ+ΛU+αΛUΛ=E⊤ΛS−1FΛS−1E.
上述系统已经解耦了,从而只要求解如下若干个标量方程:
(λi+λj+λiλj)Uij=∑klEik⊤FklλkSλlSElj.(\lambda_i + \lambda_j + \lambda_i \lambda_j) \mathcal{U}_{ij} = \sum\limits_{kl} E^\top_{ik} \frac{F_{kl}}{\lambda^S_{k}\lambda^S_l} E_{lj}. (λi+λj+λiλj)Uij=kl∑Eik⊤λkSλlSFklElj.
Step 3 利用
U=EUE⊤.U = E \mathcal{U}E^\top. U=EUE⊤.
去解得 UUU.

问题以及 Neumann 边界条件下的特征值算法

上述算法不能直接应用到 Neumann 边界问题, 因为在 Neumann 边界条件下,刚度矩阵的特征向量有一个为0,无法直接像上面那样求解特征值问题得到。需要改进上面这个算法。下面描述文章中给出的方法。该方法的思路是首先通过修正原来的基函数,使得参考矩形下的Legendre谱方法的刚度矩阵变成如下形式：
S=(000I⋆).S = \begin{pmatrix} 0 & 0 \\ 0 & I_\star \end{pmatrix}. S=(000I⋆).
为了达到这个目的,回忆Neumann边值问题的 Legendre-Galerkin 谱方法的基函数为
ϕk=Lk(x)−bkLk+2(x),bk=k(k+1)(k+2)(k+3).\phi_k = L_k(x) - b_k L_{k+2}(x), \ b_k = \frac{k(k+1)}{(k+2)(k+3)}. ϕk=Lk(x)−bkLk+2(x), bk=(k+2)(k+3)k(k+1).
而由文章中计算出的该基函数下的刚度矩阵为对角矩阵,并且对角线元素为:
Skk=(4k+6)bk.S_{kk} = (4k+6)b_k. Skk=(4k+6)bk.
因此定义如下修正的基函数:
ϕk⋆(x)=1(4k+6)bkϕk(x)=ak⋆Lk(x)−bk⋆Lk+2(x).\begin{aligned} \phi^\star_k(x) &= \frac{1}{\sqrt{(4k+6)b_k}} \phi_k(x) \\ &= a^\star_k L_k(x) - b_k^\star L_{k+2}(x). \end{aligned} ϕk⋆(x)=(4k+6)bk1ϕk(x)=ak⋆Lk(x)−bk⋆Lk+2(x).
其中 ak⋆a_k^\starak⋆ 以及 bk⋆b_k^\starbk⋆ 为修正后的基函数组合系数,为:
ak⋆={1,k=0,(k+2)(k+3)2k(k+1)(2k+3),k>0,bk⋆={0,k=0,k(k+1)2(k+2)(k+3)(2k+3),k>0.a_k^\star = \left\lbrace \begin{aligned} &1, & k = 0, \\ &\sqrt{\frac{(k+2)(k+3)}{2k(k+1)(2k+3)}}, & k >0, \end{aligned} \right. \quad b_k^\star = \left\lbrace \begin{aligned} &0, \ &k = 0, \\ &\sqrt{ \frac{k(k+1)}{2(k+2)(k+3)(2k+3)} }, & k> 0. \end{aligned} \right. ak⋆=⎩⎪⎪⎨⎪⎪⎧1,2k(k+1)(2k+3)(k+2)(k+3),k=0,k>0,bk⋆=⎩⎪⎪⎨⎪⎪⎧0, 2(k+2)(k+3)(2k+3)k(k+1),k=0,k>0.
显然, 在这组修正的基函数下, 刚度矩阵变为 diag{0,1,1,⋯,1}diag\{0, 1, 1, \cdots ,1 \}diag{0,1,1,⋯,1}.下面重新计算新基函数下的质量矩阵:
mjk=∫ϕk⋆(x)ϕj⋆(y)dx=∫(ak⋆Lk(x)−bk⋆Lk+2(x))(aj⋆Lj(x)−bj⋆Lj+2(x))dx=ak⋆aj⋆∫Lk(x)Lj(x)dx−ak⋆bj⋆∫Lk(x)Lj+2(x)dx−bk⋆aj⋆∫Lk+2(x)Lj(x)dx+bk⋆bj⋆∫Lk+2(x)Lj+2(x)dx.\begin{aligned} m_{jk} &= \int \phi^\star_k(x) \phi^\star_j(y) dx \\ &= \int (a^\star_k L_k(x) - b^\star_kL_{k+2}(x)) (a^\star_j L_j(x) - b^\star_j L_{j+2}(x)) dx \\ &= a_k^\star a_j^\star \int L_k(x) L_j(x) dx - a_k^\star b_j^{\star} \int L_k(x) L_{j+2}(x) dx \\ & - b_k^\star a_j^\star \int L_{k+2}(x) L_j(x) dx + b_k^\star b_j^\star \int L_{k+2}(x) L_{j+2}(x)dx. \end{aligned} mjk=∫ϕk⋆(x)ϕj⋆(y)dx=∫(ak⋆Lk(x)−bk⋆Lk+2(x))(aj⋆Lj(x)−bj⋆Lj+2(x))dx=ak⋆aj⋆∫Lk(x)Lj(x)dx−ak⋆bj⋆∫Lk(x)Lj+2(x)dx−bk⋆aj⋆∫Lk+2(x)Lj(x)dx+bk⋆bj⋆∫Lk+2(x)Lj+2(x)dx.
从而
mkk=(ak⋆)21k+12+(bk⋆)21k+52={2,k=0,(k+2)(k+3)k(k+1)(2k+1)(2k+3)+k(k+1)(k+2)(k+3)(2k+3)(2k+5),k>0,\begin{aligned} m_{kk} &= (a_k^\star)^2 \frac{1}{k+\frac{1}{2}} + (b_k^\star)^2 \frac{1}{k+\frac{5}{2}} \\ &= \left\lbrace \begin{aligned} &2, & k=0, \\ &\frac{(k+2)(k+3)}{k(k+1)(2k+1)(2k+3)} + \frac{k(k+1)}{(k+2)(k+3)(2k+3)(2k+5)}, & k > 0, \end{aligned} \right. \end{aligned} mkk=(ak⋆)2k+211+(bk⋆)2k+251=⎩⎪⎨⎪⎧2,k(k+1)(2k+1)(2k+3)(k+2)(k+3)+(k+2)(k+3)(2k+3)(2k+5)k(k+1),k=0,k>0,

mk,k+2=mk+2,k=−k(k+1)(k+4)(k+5)(2k+3)(2k+7)1(k+2)(k+3)(2k+5).m_{k,k+2} = m_{k+2, k} = - \sqrt{ \frac{k(k+1)(k+4)(k+5)}{(2k+ 3)(2k + 7)}} \frac{1}{(k+2)(k+3)(2k + 5)}. mk,k+2=mk+2,k=−(2k+3)(2k+7)k(k+1)(k+4)(k+5)(k+2)(k+3)(2k+5)1.
有了上述的过程,我们重新考虑离散系统:
SUM+MUS+αMUM=F.SUM + MU S + \alpha M U M = F. SUM+MUS+αMUM=F.
此时
S=(000I⋆),M=(m000M⋆).S = \begin{pmatrix} 0 & 0 \\ 0 & I_\star \end{pmatrix}, \ M = \begin{pmatrix} m_0 & 0 \\ 0 & M_\star \end{pmatrix}. S=(000I⋆), M=(m000M⋆).
注意,此时和之前最大的区别就是刚度矩阵和质量矩阵可被同一正交矩阵对角化. 考虑特征值问题:
M⋆x=λxM_\star x = \lambda x M⋆x=λx
则有
M⋆E⋆=E⋆Λ⋆M_\star E_\star = E_\star \Lambda_\star M⋆E⋆=E⋆Λ⋆
即
M⋆=E⋆Λ⋆E⋆⊤M_\star = E_\star \Lambda_\star E^\top_\star M⋆=E⋆Λ⋆E⋆⊤
从而有
M=(100E⋆)(m000Λ⋆)(100E⋆⊤):=EΛME⊤.S=(100E⋆)(000I⋆)(100E⋆⊤):=EΛSE⊤.\begin{aligned} &M = \begin{pmatrix} 1 & 0 \\ 0 & E_\star \end{pmatrix} \begin{pmatrix} m_0 & 0 \\ 0 & \Lambda_\star \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & E^\top_\star \end{pmatrix} := E \Lambda_M E^\top \quad. \\ &S = \begin{pmatrix} 1 & 0 \\ 0 & E_\star \end{pmatrix} \begin{pmatrix} 0 & 0 \\ 0 & I_\star \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & E^\top_\star \end{pmatrix} := E \Lambda_S E^\top. \end{aligned} M=(100E⋆)(m000Λ⋆)(100E⋆⊤):=EΛME⊤.S=(100E⋆)(000I⋆)(100E⋆⊤):=EΛSE⊤.
从而有:
EΛSE⊤UEΛME⊤+EΛME⊤UEΛSE⊤+αEΛME⊤UEΛME⊤=F.E \Lambda_S E^\top U E \Lambda_M E^\top + E \Lambda_M E^\top U E \Lambda_S E^\top + \alpha E\Lambda_ME^\top U E\Lambda_ME^\top = F. EΛSE⊤UEΛME⊤+EΛME⊤UEΛSE⊤+αEΛME⊤UEΛME⊤=F.
令 U‾=E⊤UE\overline{U} = E^\top U EU=E⊤UE 以及 F‾=E⊤FE\overline{F} = E^\top F EF=E⊤FE, 则上述系统化简为
ΛSU‾ΛM+ΛMU‾ΛS+αΛMU‾ΛM=F‾.\Lambda_S \overline{U} \Lambda_M + \Lambda_M \overline{U} \Lambda_S +\alpha \Lambda_M \overline{U} \Lambda_M = \overline{F}. ΛSUΛM+ΛMUΛS+αΛMUΛM=F.
从而求解标量方程
[(λS)j(λM)k+(λM)j(λS)k+α(λM)j(λM)k]U‾jk=F‾jk[(\lambda_S)_j(\lambda_M)_k + (\lambda_M)_j (\lambda_S)_k + \alpha (\lambda_M)_j (\lambda_M)_k] \overline{U}_{jk} = \overline{F}_{jk} [(λS)j(λM)k+(λM)j(λS)k+α(λM)j(λM)k]Ujk=Fjk
解得 U‾\overline{U}U, 再利用变换 U=EU‾E⊤U = E \overline{U} E^\topU=EUE⊤ 来获得 UUU.

Neumann边界条件下二维Legendre-Galerkin方法快速计算相关推荐

php 二维数组根据键值合并二维数组_php数组实现根据某个键值将相同键值合并生成新二维数组的方法详解...
这篇文章主要介绍了php数组实现根据某个键值将相同键值合并生成新二维数组的方法,涉及php数组的遍历.赋值相关运算技巧,需要的朋友可以参考下本文实例讲述了php数组实现根据某个键值将相同键值合并生成 ...
二维码推广方法20种
二维码推广方法20种二维码,又称二维条码,最早发明于日本,它是用某种特定的几何图形按一定规律在平面(二维方向上)分布的黑白相间的图形记录数据符号信息的.在移动互联网发展中顺势成长,已成为当前手机端 ...
二维码制作方法有哪些？教你简单的二维码制作方法
二维码是怎么制作的呢?二维码是用某种特定的几何图形按照一定规律在平面(二维方向)分布的黑白相间的图形记录数据符号信息的.现如今,随着智能手机的广泛普及和技术的不断改进,二维码已经被广泛应用于商业领域中 ...
java二维数组的创建,java二维数组创建方法
java动态创建二维数组,从零学java笔录-第31篇图解二位数组在内存中存储,java二维数组动态赋值,java二维数组创建方法二维数组的定义 type arrayName[ ][ ]; typ ...
怎样生成二维码？分享几种轻易生成二维码的方法
怎样能够生成二维码呢?在日常中,使用二维码拥有很多的便利之处,比如,二维码可以被轻松地扫描和分享,使得信息的传递更加便捷.二维码还可以被用于实现物品追踪.防伪和溯源等功能,保证商品的质量和安全性.总之 ...
python append函数二维_python创建与遍历List二维列表的方法
python创建与遍历List二维列表的方法 python 创建List二维列表 lists = [[] for i in range(3)] # 创建的是多行三列的二维列表 for i in ran ...
编程示例：表格程序开发的EXCEL方法，以二维码的数据容量计算为例
编程示例:表格程序开发的EXCEL方法,以二维码的数据容量计算为例在二维码的计算中,它的第一个表格是以版本号为参数,计算该版本下的数据容量. 表1如下: 在EXCEL中以公式的形式生成与上图一致的表 ...
二维数组更改vue_使用vue中的v-for遍历二维数组的方法
如下所示: {{itemss}} 其中,data数据为: this.data = [ [ { type: '', name: '资产', start: '期末余额', end: '期初余额' }, { ...
在OpenCV环境下写的灰度图像二维傅里叶换,幅值计算,频谱平移和将数值归一化到0到255区间的四个函数
图像处理开发需求.图像处理接私活挣零花钱,请加微信/QQ 2487872782 图像处理开发资料.图像处理技术交流请加QQ群,群号 271891601 灰度图像的二维傅里叶变换(cv_gray_fft ...

Neumann边界条件下二维Legendre-Galerkin方法快速计算

Dirichlet边界条件下的算法:

问题以及 Neumann 边界条件下的特征值算法

Neumann边界条件下二维Legendre-Galerkin方法快速计算相关推荐

最新文章

热门文章