第一章算法 §1二维 Laplace方程外问题 设r是平面上的一个凸多边形,它的外部区域是Ω(图1)。 我们在Ω上求解 Laplace方程的 Dirichlet问题 △=Q (x,y)∈D,(1.1 (1.2) 或者将边界条件(1.2)换成 01 (1.3) 01 (I.1),(1.3)是…个Nem 题,以后我们永远以ν表示外法 线方向。与(1.1)联系的应变能 是 dxdy (I.4 为求解问题(1.1),(1.2)或者(1.I),(1.3),我们将2剖分为 无限多个三角形单元,不妨设巫标原点O在F内部,取一个常 数5>1,以(点为相似中心,以§,52,…,5k,…为比例常数,作 F的相似形,它们分别记作,F2,…,Ik,…,每两个多边形之 间的ⅸ城称为一“层”(图2)。然后,我们将每一层进一步剖分 成单元,可以采用如下方式:在!上取若干点作为节点,其 中顶点必颁坛节点,在毎条边上还可以按需要适当地再选一些节 点,点O点出发作射线与各节点相连,这样就将每→层削分成相 似的若个四边形,再将每个四边形潮分成两个三角形,需要注意
约是,每层的剖分方式必须一致(图2)。 边界r的形状可能不是一个凸多边形,而是非常复杂的形 状。这时我们可以将厂。的外部区域分解成一个形状规则的无界 区城与一个形状复杂的有界区域。无界区域仍记作9,它可以按 照上述办法称剖分。对于有界区 域,则可以按照常规的办法作有 限单元剖分(图3)。我们还使这 两种剖分互相协调一致,即在交 界处,三角形单元的节点与边互 相吻合。在这种情况下,我们仍 然先分析区域Ω,我们将给出2 上的组合刚度矩阵Kz的箅法 Kz是一个有限阶的矩阵,在使用的时候,或者用它直接求解边 值问题,或者将它与其它单元的刚度矩阵叠加,以求解形状复杂 的区蚁的边值题。对于后者,Ω是当作一个单元来处理的。 现在考虑多边形fk。从某一点开始,按逆时针方向,将Fk
上的各节点排一个次序。设节点共有n个。(1.1)的解在节点 上的值毡排了一个次序,记作y4),k2),…,3),我们将它们排 列成…个列向量(v41,y2,…,y)“,记作Jk,这里T表示矩阵 或者向量的转置。 考虑多边形Fk1与k之间的第七层。对于每个三角形单 元,在节点上的值给定以后,作线性插值,可以按常规的方法求 单元刚度矩阼(例如参看52]),将第k层上的各单元刚度矩阵按 节点叠加,可以得到一层上的斓度矩仵。我们将这个刚度矩阵记 K (15) 即如果第k层上的应变能是Wk,则 K A (v-1, K 其山K…,K′,A都是n阶矩阵。从有限元方法的基不理论可知 (1.5)遠对称知阵、国此1与A互为转,K,K拊是对称 从(1、)和剐分的似性,不难判断、各层的飓度姸阵是札 的。将犴层的刚矩阵按节点叠加,可以得到一个无穷阶的总 树度絆阵。考虑到方程与边界条件,可以给出一个无穷阶的代数 程组。令K=K+K0、则有 Kayo -Ay,=Jos (⊥.6 A,ye-Ky1-A (1.6冫: 要鲁看BP咖D非t。。看面 Al K y (1.6) 物e专·D·卷·些脂 当考虑问题(1.1),(I.2)时,在F上不要列出方程,因此代要 (]6)。·而y是已知的。当考虑问题(1.1),(13)时,在上
要列出方程,即方程(1.6),其中∫是“等效节点力”,用常 规的有限元方法,它是可以从边界值∫求出来的。目前,我们的 目的是给出口上组合刚度矩阵的算法,因此可以暂且不管上述 边界条件的差别。 我们将在下章证明,存在一个n阶实矩阵X,使 Uk+1=Xy (1.7) X称为转移矩阵。以(1.)式代人(1.6)得 K-A+KX-AX)yo=0 但是y0可以是任意向量,所以X满足方程 A TX2-KX+A-0 (1.8) 令 K A R A 0 RI 其中是单位阵,方程(1、6)k可以写成如下形式: J R R (1,10) 设A,9分别是矩阵X的特征值和特征向量,则有 X9=Ag (1.1L) 在(1.10)中取无=1,令孙=9,并且注意到(].7).(.I1)得 g R1 nR (1.12) g 因此)与j 分别恳矩阵束R1-λR2的广义特征值与广义特 征向量1},为了求矩阵￥,我们先解上述广义特征值问题,求 出它的全部2n个特征值与特征向量,然后丛中选出ⅹ的n个特 征值与特征向量
以后,我们总是以det装示一个矩阵的行列式。我们计算 det(r-2R2)=det /x-2.AT A λ 在以上分块矩阵中,以λ乘第一列并加到第二列上得 K -A+AK-ALAT dat(R1-λR2)=d (-1)det(-A+1K-22AT) 由K的对称性,上式是一个λ的对称多项式,因此,如果九≠0是 个特征侑,1/也是一个特征值。我们将在下章中证明:X的 特祇值满足I,而且当14}=1时,必然有=【,其对应的 特征向量9;=(1,1,…,1)T,因此,我们只需将矩阵束R,-λR2的 特征值中满足!∴<I的选出,再添上λ=1、设特征佤为λ1,λ ·,九n,对巨的特征向量是只:,92…,9n,令 T=(;,9 A=diag(λ:,λ2,…,λn), 其中d表示对角矩阵或者对角矩阵。由(1.11)得 XT=TA Y: TAT (1.13 (1.13)就是计算转移矩阼X的公式。以(1.7代入(1.6)得 Axo=f 1.14) 其中Kz=A-AX,它称为组合刚度矩阵,我们将在下一章中 证明Kz是一个对称半正定矩阵 红合刚度姸阵与总刚度矩阼是有区别的。在现在的情况下, 总刚度矩阼是个无穷阶矩阼,通过它可以得到任意节点值产生 的应变能,而組合刚度矩阵仅是一个n阶矩阵,通过它只能得到 当节点值y,y1,…满足方程(I,6)4,(1.6)2,…时的应变能。从总 刚度矩阵到组合谢度矩阵的过程可以看作是一个消元过程。通过 这个过阻 都消去了 利用无限元方法求解的过程是这样的;先求出转移矩阵与组 5
合刚度矩阵,然后用与有限元方法完全同样的方式解一个代数方 程组,得到ψ,其余节点值可以由(1.⑦)确定 如果A是一个可逆矩阵,上述广义特征值问题可以化为普 通特征值问题,以R2左乘(1.12)式得 CAT)IK (AT)-IA Ag Ag 它就是一个普通特征值问题 特征值λ可能是复数,这就需要作复数运算。为避免这 点,我们可以将公式(1.13)胳作变形。特征值与特征向量必定是 成对共轭地出现的,设λ=a±,g=p土i是一对特征值与特 征向量,由(1.11)得 X(t1)=(±iB)(p±iq)。 分离实部与虛帑得 Xp=ap-Bq, Xq=Rp+ag, Y(卩,4)=(,q) 我们替换矩阵T与A中对应的两列,令 不难看出 X=711T 将各对复特征值一一作了以上替换以后,求转移矩阵X的运算就 是实的了。 §2 Fourier方法 上节中給出的方法可适用于任何的相似剖分,其运算量为解
个2n阶幃征值问题。在本节我们对一种特殊的相似剖分使用 Fourier变换,使得求矩阵X,Kz时,只需极少量的计算 我们要求图2中的剖分具有下面的性质:当图2绕O点旋转 2m/n弧度时,它在这个旋转下不变。这时0是一个正多边形,以 O点为中心,并且每个小四边形的剖分是一致的,如图4所示。 我们将在下章证坝,对于这样的特殊劑分,矩阵K0,K,A 都是循环矩阵,即它们都是如下形状 A,, ij T-罩 B 25 作单位根=e 以及两矩阵 11 1 F
则可以对B施行 Fourier变换,以F记F的共轭矩阵,通过计算 可以直接验证 FBF=die(∑b,∑“b …,∑d(t-1(-1) 令3k=Fk,k=0,1,…,以F左乘(1.6)的诸方程,得 P2。-Qz1=Pfo 2.1)a Q20+ Px1-Q22=0, (2.1) Q2k-+ P2k-gzk+I=0, 其中 FRI Q=FAF, P=FKF 它们都是对角矩阵。我们将它们记作 Po= diag(p,, po2,e, pom) P Q=diag(q,q2,…,4)。 令 zk=(z1),z2) k2) Ffa=(φ1,p2),…,p(m)T 将(2.1)写成分量的形式,得到个无穷阶代数方程组(t=1,2,… 722=0 qx},=0 下面象§1一样地作推导,存在常数x',使 k =o, 1 得到x()满足的方程 1)(x2) E (2.2) (2.2)的两个根的绝对值互为倒数。我们将在下章证明1x≤1
而且当x=1时,必有x(1=1。因此只要取 241 中绝对值较小的那一个,就是我们所需要的根x4记 Z= diag(x 则 Z2 k=0 y FZPyk, k=0,. 与§1比较,可以发现转移矩阵就是ⅹ=FzF,由此立即得到组 合刚度矩阵 Kz=Ko-atFze s3迭代法 回到§1中的任意相似剖分,虽然公式(1.13),(1.11)都是特 确的,但是求解特征值问题仍然需要通过近似计算。在这里,我 们给出另一个近似算法,它的运算量比§1中的算法的运算量要 小一些 考察由2‘(l=0,1,…)层组成的一个组合体:以 K A K 记它的组合刚度矩阵。当=0,只要取A0=A,它与§1中的符 号就是一致的。利用 Laplace方程的解使应变能(1.4)达到最小 这一事实,可以得到递推关系 K y (y,x1) K'1+ min(yT, wt)
AT (3.1) A 其中R表示n维实向量空间。计算花括号中的量关于的偏微 商,并令它等于零得 -2Aiy+2KW+2KLW-2AT2=0. 因此 (KI+K(Ay+Afz (3.2) 将(3.2)代入(3.()式并比较同类项得 K1+1EK-AT(K:+K1-AL (3.3) K4+=K-A: (K +K)-IAT (3.4) A Ai(K2+Ki (3.5) 利用(3.8)—(8.5我们可以递推地得到所有的2层组合刚度矩 阵 现在,设为已知,￠为-待定常数,并且设当k≥24, t=1,2,…,n时,ξ=C,当k=1,2,…,21-1时,方程(1.6)k 都成立。在这些限制下,由(1.4)所给出的应变能就是 K (y, cgT A[\y A yoK la (3.6 其中g1就是§I中的特征向量(1,l,….1)T, 改变常数O,使(3.6)达到最小,令n=0得 -9iAya+ciGR=0. 因此 iaI vAig 9r 91K19, 于是 W A19191A giRi 10