（计算数学专业论文）边值问题的微分求积法及其应用.pdf

同济大学硕士学位论文 摘要 微分求积( d i f f e r e n t i a lq u a d r a t u r e ，简称d q ) 方法是一种非常准确有效的数 值计算方法。但是传统的d q 方法的使用仅限于规则区域上。要处理不规则区 域问题，坐标变换不可缺少。钟 3 7 - 3 9 1 提出了三角微分求积( t r i a n g u l a rd i f f e r e n - t i a lq u a d r a t u r e ，简称t d q ) 方法，避免了坐标变换。为求解五边形上的椭圆型 微分方程边值问题，我们借助了区域分裂法( d o m a i nd e c o m p o s i t i o nm e t h o d ， 简称d d m ) ，将大区域分成几个小区域，在每个小区域上，用t d q 方法这 种方法我们称之为三角微分求积区域分裂法( t r i a n g u l a rd i f f e r e n t i a lq u a d r a t u r e d o m a i nd e c o m p o s i t i o nm e t h o d ，简称t d q d d m 。本文还引用了边界归化技 术( b o u n d a r yr e d u c t i o nt e c h n i q u e ，简称b r t ) ，消除内点，那么只需要求解拟边 界点所满足的方程即可，大大减少了计算量。数值实验表明，t d q d d m 对不规 则区域问题简单易行。 本文还研究了另外一个专题当求解微分方程问题时，数值微分是常用的 方法，将微分方程转化成代数方程组求解。但是数值微分对误差非常敏感，因 此我们想到用对误差不敏感的数值积分代替数值微分。在d q 方法的基础上， 本文提出了一种新的方法一一基于最高阶导函数插值的微分求积( d i f f e r e n t i a l q u a d r a t u r em e t h o db a s e do nt h eh i g h e s td e r i v a t i v e ，简称d q l h d ) 方法。这种方 法的出发点是通过对未知函数最高阶导数进行插值，然后积分得到对低阶导数 和原未知函数的逼近。并用这种方法研究了各向同性弹性薄板的弯曲问题，数 值实验表明这种方法准确性高，收敛性好，计算量少，边界条件的处理简单。 关键词：微分求积法，三角微分求积法，区域分裂法，边界归化技术，边值问 题 v i a b s t r a c t a b s t r a c t t h ed i f f e r e n t i a lq u a d r a t u r em e t h o d ( d q m ) i sa na t t r a c t i v en u m e r i c a lm e t h o d w i t hh i g he f f i c i e n c ya n da c c u r a cy b u tt h ec o n v e n t i o n a ld q mi sl i m i t e di ni t s a p p l i c a t i o nt 0r e g u l a rr e g i o n s t od e a lw i t hp r o b l e m so ni r r e g u l a rg e o m e t r i cd 0 _ m a i u s c o o r d i n a t et r a n s f o r m a t i o nh a st ob ec o n d u c t e d t h et r i a n g u l a rd i f f e r e n t i a l q u a d r a t u r em e t h o d ( t d q m ) p r o p o s e db yz h o n g 4 ，6 】 a v o i dt h ec o o r d i n a t et r a n s - f o r m a t i o n i nt h i sp a p e r ，t h ed o m a i nd e c o m p o s i t i o nm e t h o d ( d d m ) i su s e df o r t h ee l l i p t i e a lb o u n d a r yp r o b l e m s0 na p e n t a g o n a lr e g i o n i ne v e r ys u b - d o m a i n ，w e s o l v et h ep a r t i a ld i f f e r e n t i a le q u a t i o n sw i t ht d qm e t h o d ，n a m e l y , t r i a n g u l a rd i f - f e r e n t i a lq u a d r a t u r ed o m a i nd e c o m p o s i t i o n ( t d q d d ) m e t h o d w i t hb o u n d a r y r e d u c t i o nt e c h n i q u e t h ef u n c t i o n a lv a l u e so ni n t e r n a lp o i n t sc a nb ee l i m i n a t e d t h es y s t e mo fe q u a t i o n sw h i c hs a t i s f i e db yt h eb o u n d a r yp o i n t sc a l lb eo b t a i n e d n u m e r i e a lr e s u l t ss h o wt h a ti se a s ya n de f f e c t i v ef o rt r e a t i n gt h ep r o b l e m so n i r r e g u l a rr e g i o n a n o t h e rs u b j e c ti sa l s or e s e a r c h e di nt h i sp a p e r t h en u m e r i c a | d i f f e r e n t i a t i o n i so f t e nu s e dw h e nd e a l i i l gw i t ht h ed i f i e r e n t i a le q u a t i o n s u s i n gt h en u m e r i c a l d i f f e r e n t i a t i o n ，t h ed i f f e r e n t i a le q u a t i o l l sc a nb et r a n s f o r m e di n t oa l g e b r a i ce q u 争 t i o n s t h e nw ec a ng e tt h en u m e r i c a l8 0 l u t i o nf r o mt h ea l g e b r a i ce q u a t i o n s b u t t h eu n m e r i c a ld i f i e r e n t i a t i o np r o c e s si sv e r ys e n s i t i v et oe v e nas m a i i1 e v e lo fe r - r o l s i nc o n t r a s ti ti se x p e c t e dt h a to na v e r a g et h en u m e r i c a li n t e g r a t i o np r o c e s s i sm u c h1 e s ss e n s i t i v et oe r r o r s i nt h i sp a p e r w ep r o v i d ean e wm e t h o du s i n g t h ed qm e t h o db a s e do nt h ei n t e r p o l a t i o no ft h eh i g h e s td e r i v a t i v e ( d q i h d ) f o r t h ed i f f e r e n t i a ie q u a t i o n s t 】舱o r i g i n a lf u n c t i o ni st h e no b t a i n e db yi n t e g r a t i o n i nt h i sp a p e r ，t h ed q i h dm e t h o dw a sa p p l i e dt ot h eb u c k l i n ga n a l y s i so ft h i n i s o t r o p i cp l a t e s ，t h en u m e r i c a lr e s u l t ss h o wt h a to u fm e t h o di so fh i g ha c c u r a c y , o fg o o dc o n v e r g e n c ew i t hl i t t l ec o m p u t a t i o n a le f f o r t s a n di ti se a s yt od e a lw i t h t h eb o u n d a r yc o n d i t i o 璐 k e y w o r d s ：d i f f e r e n t i a lq u a d r a t u r em e t h o d ；t r i a n g u l a rd i f f e r e n t i a lq u a d r a t u r e m e t h o d ；d o m a i nd e c o m p o s i t i o nm e t h o d ；b o u n d a r yr e d u c t i o nt e c h n i q u e ；b o u n d a r y p r o b l e m i 同济大学学位论文原创性声明 本人郑重声明；所呈交的学位论文，是本人在导师指导下，进行 研究工作所取得的成果。除文中已经注明引用的内容外，本学位论文 的研究成果不包含任何他人创作的、已公开发表或者没有公开发表的 作品的内容。对本论文所涉及的研究工作做出贡献的其他个人和集 体，均己在文中以明确方式标明。本学位论文原创性声明的法律责任 由本人承担。 学位论文作者签名： 年月日 学位论文版权使用授权书 本人完全了解同济大学关于收集、保存、使用学位论文的规定， 同意如下各项内容；按照学校要求提交学位论文的印刷本和电子版 本；学校有权保存学位论文的印刷本和电子版，并采用影印、缩印、 扫描、数字化或其它手段保存论文；学校有权提供目录检索以及提供 本学位论文全文或者部分的阅览服务；学校有权按有关规定向国家有 关部门或者机构送交论文的复印件和电子版；在不以赢利为目的的前 提下，学校可以适当复制论文的部分或全部内容用于学术活动。 学位论文作者签名： 年月日 第一章前言 第一章前言 众所周知，对工程系统的分析包括两个重要步骤，即合理反映物理现象的 数学模型的建立和数学方程的求解。一般来说，工程问题是由线性或非线性 的偏微分方程描述的。因此，很多工程学科中用数值逼近法求解偏微分方程。 最常用的有：有限元法、有限差分法、边界元法等，只要选择合适的网格点 数就可求得所要求精度的数值解。有限元法是迄今为止应用最广泛，最有效的 数值计算方法，但是所用的网格点比较多，计算起来不方便。有必要寻找更为 有效的数值计算工具f 1 3 1 。d q 方法自提出以来，已被成功运用到许多工程物理 中去。这个方法数学原理简单，计算精度高。计算量少，使用方便，不依赖泛 函和变分原理，边界条件不用另外考虑。当问题具有全局性光滑解时，d q 方 法因为其精度高，且只需较少的网格点而成为优于有限元和有限差分的更佳选 择 6 - 1 2 ，1 4 ，2 7 - 3 4 ，4 0 - 4 4 ，4 6 - 5 9 i ，d q 方法的原理将在第二章中阐述。 在了解国内外最新研究状况的基础上，我们将d q 方法应用到一般情况， 把d q 方法、d d m 和b r t 等结合起来，提出t t d q d d m 方法，用来求解五边 彤上的椭圆型微分方程边值问题，这在第三章中作了详细介绍。 数值微分法对工程物理问题求解数值解非常重要。例如，我们需要对导数 的估计以便将控制方程转换成离散方程组形式。对一个确定的微分方程，出发 点是通过插值得到未知函数的估计。然后求导得到导数估计。但是众所周知， 数值微分对误差非常敏感，因此我们希望用对误差不太敏感的数值积分来代替 数值微分【2 3 ，2 4 】。 在f 2 3 ，2 4 1 中，作者给出了m u l t i - q u a d r i c 径向基函数网格法。径向基函数逼近 有砖种类型：直接径向基函数法和问接径向基函数法。直接径向基函数法的出 发点是对未知函数插值：然而，间接径向基函数法的出发点是对未知函数的最 高阶导函数的插值，然后通过积分得到低阶导函数和原函数的逼近。在间接径 向基函数法中，最后的数值解通过广义最小二乘法得到。间接径向基函数法比 直接的径向基函数法准确性更高。但是，如上所述，在间接径向基函数法中要 用到广义最小二乘法，这会给原问题的求解带来不便。因此我们提出了基于最 高阶导函数插值的微分求积法，克服了上述困难。 虽然d q 方法有很多优点，但是传统的d q 方法只适用于不超过两阶的微分 方程，因为每个离散点只有一个边界条件。但是结构力学中的控制方程一般都 是四阶微分方程，每条边上给出了两个边界条件【4 ，5 ，9 ，1 0 】。如果仍用原来的方 法离散，非常麻烦f 3 1 1 。 鉴于此，我们将d q 方法和数值积分相结合，提出了一种新的求解方法一 - - d q i h d 方法。这种方法的出发点是对未知函数最高阶导数插值，然后积分得 到对低阶导数和未知函数的逼近。通过这种方法导数的精度与未知函数的精度 类似，这是现存的其他方法无法实现的，这是第四章的主要内容。 在第五章中，我们介绍了用d q i h d 方法求解两维四阶微分方程问题，并分 析了各向同性的弹性薄板的小挠度弯曲问题，数值试验表明，这种方法有很好 的收敛性，对边界条件的处理非常容易。 l 同济大学硕士学位论文 第二章微分求积法和区域分裂法 d q 方法是最近几年引起广泛注意的一种数值方法，已在工程和数学的许 多领域得到成功应用。这种方法具有数学原理简单、使用方便、计算精度高、 计算量和内存需求量少，适合微型机等优点，是瑞利里兹法、迦辽金法、配 点法和拟谱法的一种有竞争力的替代方法有很好的发展前景，对规则区域问 题，d q 方法比有限元，有限差分和边界元法有高得多的效率 7 1 4 】。 2 1d q 方法的提出及发展 d q 方法是r i c h a r db e l l m a n 和他的同事们【l ，3 1 在2 0 世纪7 0 年代初期提出的求 解非线性偏微分方程的一种新方法。自提出以来d q 方法已被成功运用到许多 工程物理中去。这个方法数学原理简单，计算精度高，计算量少，使用方便， 不依赖泛函和变分原理，边界条件不用另外考虑。当问题具有全局性光滑解 时，d q 方法因为其精度高，且只需较少的网格点而成为优于有限元和有限差分 的更佳选择睁1 2 ，1 4 ，2 7 - 3 4 ，4 0 - 4 4 ，4 6 - 5 9 1 。近年来，对d q 方法的理论和应用的研 究，国内外都取得了很大进展。 美国曼哈顿k a n s a ss t a t e 大学的m i n g l e 应用d q 方法求解非线性扩散方 程 5 0 1 ；美国o k l a h o m a 大学的c i v a n 和s l i e p c e v i c h 应用广义的d q 方法处理输运 过程、t h o m a s - f e r m i 方程、泊松方程、以及多维问题等f 5 l 一5 4 1 ：近几年，美国 的b e r t 和他的同事们大力推崇d q 方法，他们第一次把d q 方法作为一种数学工 具用到结构力学中进行结构分析 5 - 7 ，2 9 ，3 0 1 。l i e w 和他的合作者对d q 方法的理 论和实际应用都作了深刻的研究，他们最近的研究主要是对板的三维振动，弯 曲和稳定性进行分析2 8 ，4 9 ，5 5 1 。d o k u z 大学的5 m e rc i v a l e k 把d q 方法用到弹性 杆的稳定性、振动、弯曲分析f 1 3 ，6 0 ，6 1 1 。 近年来对d q 方法的研究还体现在适当选取基函数，使其应用更为广泛且 得到更为准确的数值解。在清华大学的钟洪志的文章中提到用广义拉格朗日插 值函数作为d q 逼近的基函数。精度高，所需网格点少，运算比较简单 3 7 - 3 9 。 我们知道d q 方法是一种导数逼近的数值离散方法，它起源于传统的积分求积思 想。传统的d q 方法实际上是基于一维函数的逼近，而在文献| 4 7 1 中，s h u 和他 的同事又研究使用径向基函数作为插值基函数，将d q 逼近的思想推广到一般情 形，打破了传统d q 方法中沿网格线使用函数值的限制。根据这个方法，任意空 间导数都可以用整个物理空间中所有函数值的加权线性和来逼近。 在d q 方法的基础上，还研究出了很多新的d q 方法。l i e w 和他的同事提出 了移动最小二乘微分求积( m o v i n gl e a s ts q u a r e sd i f f e r e n t i a lq u a d r a t u r em e t h o d ， 简称m l s d q m ) 方法5 6 1 。l i e w 提出了微分求积元法( d i f f e r e n t i a lq u a d r a t u r ee l - e m e n tm e t h o d ，简称d q e m ) 用这种方法对r e i s s n e r - m i n d l i n 板进行了静态分 析 4 s l 。l i e w 还提出了d q u ，d q n ，d q z 等方法，参考文献f 1 5 2 2 ，2 s 展示了这 些方法的有效性。吴雄华、丁志宏、沈烨、孔文彬等人用d q 方法来研究了美 式期权定价问题，体t 觅- j d q 方法在金融领域的应用和发展4 0 ，5 7 1 。根据d q 方 2 第二章微分求积法和区域分裂法 法的思想，在传统的d q 方法的基础上并结合d d m ，吴雄华老师和他的学生李 晨，吴芸、沈烨、刘书亭等人提出了微分求积区域分裂法( d i f f e r e n t i a lq u a d r a - t t t r ed o m a i nd e c o m p o s i t i o n ，简称d q d d ) 方法f 3 3 ，4 9 ，5 8 ，5 9 1 ，来求解三角形区域 上的偏微分方程边值问题、非线性奇异摄动问题及抛物型方程自由边界问题。 他们的工作表明这种方法十分高效，很好地发挥- j d q 方法和d d m 各自的优 点除此之外，还有很多的国内外专家学者，对d q 方法的研究也有着重要的显 著的成果。 从上述专家学者发表的专著论文来看，d q 方法被认为是一种非常有效的数 值方法，用很少的点就可求得准确的数值解。但是这种方法也有其局限性，仅 局限于规则区域上维问题和多维问题的求解，对复杂的不规则区域问题的求 解存在一定困难f 弱l ，需要用坐标变换将求解区域化成规则区域，这导致微分方 程变得非常复杂。 2 2d q 方法的基本原理 d q 方法的基本思想是函数对某变量的偏导数可以由在此变量方向上所有 离散点处的函数值的加权线性求和来逼近。而其中的权系数不依赖于任何具体 问题，只与网格剖分有关。因此偏微分方程可以简化成用这些权系数表示的代 数方程组。考虑，( z ) 用n 个离散点的插值函数逼近，有 ，( 司= 乏：a x ) i ( * a ， ( 2 2 1 ) j = l 其中巧是离散节点，( ) 是函数在离散点上的值，乃 ) 是l a g r a n g e 插值多项 式( 2 6 】对式( 2 2 1 ) 求一阶、二阶导数得： 掣i 。= o , i j f ( x j ) ，江1 n ( 2 删 。 j = 1 c 0 2 a f 。( 2 x ) i 。：n 地) ，l ：1 仃 ( 2 删 一 j = l 其中啦j 与6 j 是z 方向上对应的一阶与二阶偏导数的权系数， = o f 如j ( x ) k 岵= 掣k 方程( 2 。2 2 ) 与( 2 2 3 ) 对不超过n 一1 次的多项式必须是正合的。为确定上述的权 系数。常用的试函数有很多，例如多项式、正弦、余弦函数、及各种正交多 项式。其中拉格朗日插值函数应用广泛，不受样点选取的影响，在计算权系数 时，不会导致范德蒙矩阵的病态即， ( z ) = 肇k ：兰= 警 ( 2 ，2 4 ) 将此式代入方程( 2 2 2 ) ，可得到直接计算一阶导数权系数的公式： = f _ l ，忙l ，2 ，靠 = 磊一。i 一。 3 同济大学硕士学位论文 2 击0 毛嚣烈2 ，川“，2 ，n 令a = ( ) ，b = ( ) ，则易得b = a 2 得到权系数矩阵a 和b 后，可以建立d q 公式。对于只有一个变量的函 数，) ，令x = ( t 1 , z 2 一，蕾。) t 是n l 的列向量，则有 ，( x ) = a f ( x ) ，”( x ) = b ( x ) ( 2 2 5 ) 其中 ，( x ) = ( ，0 - ) ，( 如) ) 7 ，( x ) = ( ，7 ( z 1 ) ，7 ( z 。) ) r ， ，( x ) = ( ，( z t ) ，”( z 。) ) r d q 方法的精确性取决于权系数的精确性，因此d q 方法的精确性受离散节点的 选择的影响，故权系数的确定多用正交多项式的零点。 同样对两个变量的函数，( z ，y ) ，令x = ( 。l ，z 。) ，y = ( y a ，) 。 f = ( ，( ，) ) 。则有 筹= 鼠e 雾= f 垛 筹吐f 筹一，t ，畿= 允脚t 其中也，a ，e ，玩分n x ，方向上的权系数矩阵。d q 公式的引进，使原来包含 函数各阶导数的微分方程( 组) 转化为以离散点函数值为未知数的代数方程组， 数学原理简单，操作起来非常方便。 2 3 区域分裂法 区域分裂法( d d m ) 是八十年代崛起的一种新方法。随着计算机技术的发 展，大型工程问题的求解求助于通过计算机来实现，而计算机并行技术的发 展，使得效率大大提高。区域分裂法作为并行理论的一种，受到人们的广泛重 视。由于该方法的许多优点，如区域分裂方式的任意性，区域分裂后物理问题 的数学描述的多样性以及问题求解方法的灵活性，因而它已成为并行计算机上 进行大型科学计算问题的算法设计的重要方法1 5 ，2 6 1 。1 9 8 5 年后研究渐趋活 跃。1 9 8 7 年之后，每年召开一次国际会议，美、俄、法、意、中国等数值分析 学家竟相参加此项研究。进入九十年代，区域分裂法已成为当今计算数学的热 门领域，其趋势方兴未艾。 d d m 的发展史，最原始的思想可以追溯到1 8 7 0 年由德国数学 家h a s c h w a r z 提出的著名的s c h w a r z 交替法f 4 5 l ，其本意是借用交替法的思 想来论证两个相互重叠区域的和集上l a p l a c e 方程d i r i c h l e t 问题解的存在 性。1 8 9 0 年p i c a r d 进一步发展了s c h w a r z 的思想，并用之于论证非线性椭圆型 偏微分方程解的存在性和唯一性。由于s c h w a r z 算法可以把复杂区域分解为若干 相互重叠的子区域，在子区域上可以采用诸如f o u r i e r 快速变换( f f t ) 、交替方 向法等快速算法求解，因此引起了广泛的关注。近十年来，随着计算机技术和 4 第二章微分求积法和区域分裂法 并行算法的迅猛发展，s c h w a r z 交替法由于其并行的可实现性，极大地刺激了研 究者的兴趣。d d m 最初是为了在并行机上解决p d e s ，特别是在m i m d 上，现 在它不仅是一种纯粹的数值方法，而且用于觎决各种各样的问题。 与以往的方法不同，区域分裂法并不直接对问题本身进行处理，而是将计 算区域划分为多个规则的小区域分别单独处理。简而言之，区域分裂法是把计 算区域分解为若干个子区域，子区域的形状应尽量规则，于是原问题的求解转 化为在子区域上求解。因此，区域分裂法可以极大地缩小计算规模。同时由于 计算规模的缩小和子区域的规则形状为在子区域上使用快速算法提供了可能 d d m 备受关注是因为它具有其他方法无法比拟的优点： ( 1 ) 它把大问题化为若干小问题，缩小计算规模。 f 2 ) 子区域形状如果规则( 如长方形) ，其上或者允许使用熟知的快速算法，如 快速f o u r i e r 变换( f f t ) ，谱方法r 方法等：或者已经有解这类问题的高效率 软件备用。 ( 3 ) 允许使用局部拟一致网格，无需用整体拟一致网格。甚至各子区域可以用不 f - j 离散方法进行计算。这对于形状极不规则的闯题，如建筑结构中的叛，梁组 合结构，具有很大的灵活性。 ( 4 ) 允许在不同的子区域使用不同的数学模型，以便整体模型更适合于工程物理 实际情况。再如对数学中颇为棘手的混合型方程，如果把区域的椭圆型部分和 双曲型部分作为两个子区域分别考虑，就简单多了。 ( 5 ) 算法是高度并行的，即计算的主要步骤是在各子区域内独立进行的。 ( 6 ) 对对称区域有更简单的区域分裂算法。 d d m 主要有三种类型：( 1 ) 重叠型区域分裂法( o v e r l a p p i n gd o m a i nd e - c o m p o s i t i o n m e t h o d ，d d m 1 ) ；( 2 ) 非重叠型区域分裂法( n o n o v e r l a p p i n g d o - m a i nd e c o m p o s i t i o nm e t h o d ，d d m 一2 ) ；( 3 ) 虚拟区域分裂法( f i c t i t i o u sd o - m a i nm e t h o d d d m 一3 ) 。 虽然d d m 一1 比较流行，但是为了消除奇性，不跨间断计算，我们借鉴吴 雄华等人的做法f l 一5 1 ，采用d d m - 2 ，并引进边界归化的思想来消除内部虚拟边 界。d d m 。2 求解方法之一是子结构消元法。先消去子区域内网点上的未知量， 再解虚拟边界网点上未知量组成的低阶方程组，然后回代求出各子区域内网点 上的解。如果虚拟边界网点上的方程采用直接方法求解而不采用迭代方法，贝 j 无需考虑算法的收敛性，故可用来处理很广泛的一类线性问题。但是要想消去 内点得到虚拟边界点所满足的方程组是比较困难的。边界元中边界归化思想， 利用分块矩阵导出边界点和内点的关系，得到内点用边界点表示的表达式， 消除内点。从而得到边界点上的函数值所满足的方程，求解即得到边界点函数 值，然后回代得到子区域内点上的函数值。 d d m 2 的实质是引进两个变量：内部变量啦和虚拟边界变量“b ，先在子区 域上并行求解子问题，将内部变量啦消去( a o 用u 口表示t ) 得到u b 所满足的方 程；再求解最后回代得到原问题的解。为简单起见，我们以两个子区域情形为 倒来说明其主要思想。考虑如下形式d i r i c h l e t 问题： l 社u - - - f , 裟 整个大区域q 被分成互不重叠的两个子区域q l ，q 2 ，虚拟边界为b = a q ln a q 2 5 同济大学硕士学位论文 ( 如图2 1 所示) 。 a l b 图2 1 ：区域分裂示意图 令u 1 ，“2 ，u b 分别表示整个区域q 上的解在子区域q 1 。q 2 ，及虚拟边 界b 上的值，那么u = ( t 1 ，缸2 ，t b ) r ，并且在虚拟边界b 处满足连续性条件： u l2 t 2 o u l o v a+ o m 概 这里m 表示区域q 在边界处b 的单位外法向量。 d d m 2 的关键是如何用u b 表示子区域内函数值t “，然后在整个区域上的问 题就分裂为两个子区域上的予问题，可以完全独立并行求解。 6 第三章 三角微分求积区域分裂法求解多边形问题 第三章三角微分求积区域分裂法求解多边 形闯题 在上一章中，已经对微分求积法作了简单的介绍。我们知道微分求积法对 规则区域上的一维和多维问题能够用很少的计算量求得准确的数值解，但求解 复杂的不规则区域问题存在一定困难另外，微分求积法是基于插值函数进行 函数逼近所以如果在整个区域上使用微分求积法，则应通过坐标变换将非矩 形的求解区域转换到规范化的计算区域上，结果简单的控制方程转变成复杂的 微分方程，尤其是对高阶的情况，且在转换过程中可能产生奇性。 对处理三角形区域上的问题，清华大学的钟宏志提出了三角微分求积法， 克服了上述困难，我们把这_ 种方法与区域分裂法相结合，形成三角微分求积区 域分裂法，对处理多维问题非常有效。这一章我们结合边界归化技术，来求解五 边形上的椭圆型微分方程边值问题。 3 ，1 三角微分求积法概述 t d q 方法的定义 如z h o n g ( 2 0 0 0 ) 所述，首先将三角形区域剖分网格( 见图3 1 ) ，即一个 顶点确定条对边，该边法向相应地确定一个方向。用平行线将顶点1 和对 边1 按方向1 等分，每条线段用从0 到m 的数字来表示，线段0 代表边1 ，线段m 代表通过顶点1 的线段在l 方向上的任何一条线段可记怖，另外两个方向采 取同样的措施。很明显，用“q ，r 分别表示三个方向的任何一条线段，则网格中 任意一点可用q ，r ) 来表示。任意一点的面积坐标可表示为扣m ，q m ，r m ) 且有： p + 口+ r = m ，0sp ，g ，rs m ( 3 1 1 ) 整个三角形共有 m = ( m + 1 ) ( n + 2 ) 2( 3 1 2 ) 个网格点。任意一点( z ，! ，) 的三个面积坐标可通过 l t = ( 啦+ 以z + q y ) ( 2 ) ，i = 1 ，2 ，3 ( 3 1 3 ) 计算。建直角坐标系中三角形的面积，我们知道 i 1z 1 l l 2 a = l1x 2 抛1 ( 3 1 4 ) 1 1x 3 驰l 其中( 3 1 3 ) 中的系数即上述行列式中相应代数余子式的值，例如： o l = z 2 珈一x 3 抛，b l = y 2 一船， c 1 = 一( z 2 一黝) ，( 3 1 5 ) 其余的系数可以通过交换下标1 ，2 ，3 相应得到。三个顶点应按逆时针排列，以保 证( 3 1 4 ) 中行列式值为正 7 同济大学硕士学位论文 图3 1 ：任意三角形的网格剖分 雄2 y 2 ( 0 1 1 1o 根据d q 方法，三角形内任一点空间变量的偏导数可通过三角形内所有点的 函数值的加权线性和来表示。即， mm - j 巩 m ，y ) ) 蚋= c 端，o ( 3 1 6 ) j = o t ：o 其中巩是n 阶微分算子；带下标的风 ，( z ，) 蚋表示点( q ，屈们处的导数 值：c ：臻，是与点g ，r ) 处的函数值有关的系数。求和指标，g r ) 满足 ( p ，q ，r ) = ( m 一 一j ，j )( 3 1 7 ) 与d q 方法类似，对微分方程问题引进t d q 方法，则建立了以三角形内各点 的函数值为未知量的代数方程组，写成矩阵乘积形式为： 警如k 雾“k 差。如k 等“b w ( 3 1 8 ) 此处的y 代表所有网格点的函数值组成的列向量，如，a 。，岛，z 是权系数的矩 阵形式。 权系数的确定及其性质 与d q 方法类似，只要确定了权系数，则问题就可以解决了。检查网格，发 现可用有限元中对三角形的区域常用的广义拉格朗日形状函数为插值基函数。 因此，我们取，为下面的m 个基函数： ，p 旷= 瓦( l 1 ) 了g ( l 2 ) 7 7 ( 岛) ，0sp ，q ，rs 仇 ( 3 1 9 ) 其中辅助函数为： - p 溉) = ￥一厶“+ 1 v t 譬孑 m ； ( 3 1 1 0 ) 8 第三章三角微分求积区域分裂法求解多边形问题 类似的有7 口( l 2 ) 和- r ( 岛) 我t r i e r 容易看出这些基函数的性质， 旆f h l2 乏管? 竺鬈：鬻m 幻柳z 3 ) i 如咖2 知 ( ) ,( a ，卢，y ) = 扫，g ，r ) ； ( 3 1 1 1 ) 一、o ，o t h e r w i s e 这里的j 是k r o n e c k o r 积。对基函数求一阶- 9 数有， 警l 嘶= 唱。o 扣略。轧= 嘲， = 0 卸 = 0 = o ( 3 1 1 2 ) 啦秤键等等， = b ，( 2 a ) k ( 2 , ) b 3 c 2 a ) 1 袈i r - p 笼圣f 矗7 。畿j 蚋 ( 3 1 1 3 ) ：二三兰：? ( 臻巍# 盟- - - - ) 蚋 t ，卟川训圳训叫懒孰p 1 。1 4 鼽加_ 够豢， l 0 ， 同理，无( l 1 ) 和- r ( l 1 ) 分别对l 2 和l 3 的一阶导数为 跏。：陋m l l q 瓢 f ，叫， 2 p o r ； o ? s p 一1 ； ( 3 1 1 5 ) p = l ： p = 0 2 p ； 0 卢q l ； ( 3 1 1 6 ) 口= 1 ； q = 0 一f ，n ：= l 厶3 - 幽k - t 1 ， 2 r 竹 瓦a l i b 咖= 磐蛐卅m “l 警g 。；( 3 “7 ) l0 ， r = 0 数盥盥 同济大学硕士学位论文 。， 0 l o ，1 喻 1 1 1 ，二l b j 一 e ( - a , b ) ，n 5 n o 一1 ( 4 ) 的” ， 夕 弋? ( 3 )v ， n o ( - c ，d l n 6 图3 2 ：例2 的求解区域( 口= c o s ( 口l o ) ；b = s i n o r l o ) ；c = s ( 3 丌1 0 ) ；d = 一s i n ( 3 7 r l o ) ) 根据t d q 方法的定义，高阶导数的权系数矩阵可通过一阶导数的权系数矩 阵的自乘得到，例如 仇, n - j c 端，= c 黯，。酬阿， j - - - o = o mm - j c 拱，= c 纷。础妒 3 2 三角微分求积区域分裂法 要在五边形上用t d q d d m ，我们首先将五边形划分成五个三角形。从 五边形的中心连接五个顶点，分成五个小三角形区域( 见图3 2 ) 。设第一个区 域( d a b ) 上的控制方程为 。磐+ 6 b 百i ，2 v 1 + c 娑+ d 娑+ e v i ：f 4 - ( 3 ，2 1 )o 石- 了+ 6 彳f + c i 一+ d i =- 【j ，z l j 上式里的口l ，o ，6 ，c ，d ，e 和，都是茁，v 的函数。根据t d q f f 式，把上式写成矩阵乘 积的形式 a o a “y 1 + b o b w y l - i - c o a 。v 1 - - i - do 岛v 1 + e o v l = f ( 3 2 2 ) 这里的“o ”代表矩阵的h a d a m a r d 乘积，a 。，如，马是mxm 矩阵( 已 在( 3 ，1 ，8 ) 中给出) ，v 1 = ( ? 3 1 ( 以，u a ) m ，a = ( o ( 甄，v j ) ) m ，b = ( b ( 瓤，y j ) ) m ，c ? ( c ( ，y a ) m ，d = ( d ( x 。，协) ) ，e = ( e ( ，船) ) m ，f = ( ，( 霉，协) ) j i f 将j 4 ，a z ，v 1 分 第三章 三角微分求积区域分裂法求解多边形问题 a 一( 委 4 产( 奏萎萎萎萎萎萎 k ( 善 3 3 边界归化技术 根据分块矩阵的乘法，矿满足下式 ( d i a g ( a ) a 。+ d i a g ( b ) b + d i a g ( - c ) - a , ：+ d i a g ( d ) 岛+ d i a g ( e ) ) v 1 = f 一( - a o ( a x x l 2 曙+ a x x l 3 晗+ a 4 呀+ a x z l 5 w + a x x l 6 v g + 屯w ) + 百o ( b w l 2 呀 + b w l 3 曙+ b y y l 4 呀+ b 1 5 w + b 1 6 曙+ b u y l 7 w ) + 虿o ( a m 嵋+ a m 曙 + a 略+ a m 增+ a 圳增+ a m w ) + 五。( 岛1 2 晗+ 易1 3 日+ 岛1 4 曰 + b u l 5 w + 岛1 6 曙+ 黾1 7 w ) ) ( 3 3 1 ) 从而找到了边界点和内点的关系，所有内点的函数值可用边界点的函数值来表 示： 驴1 = p p 0 1 + p p 2 1 曙+ p p 3 1 曙+ p p 4 1 曙+ p p 5 1 曙+ 删1 咯+ p p 7 1 曰( 3 3 2 ) 消除内点，再解以拟边界点上函数值为未知量组成的低阶方程组，然后回代求 出各子区域内网点上的解。 3 4 拟边界上满足的方程 接下来我们推导五条拟边界上的外法向导数公式以及在顶点处给定方向的 导数公式，通过连续性条件，建立拟边界上点所满足的方程。 对任意所给的方向n ，有 娑：芸胡+ 娑5 触 ( 3 4 1 ) 这里的0 是z 轴正方向与所给方向n 间的夹角。根据t d q 公式，上式可写成 = 畦。o ( 屯曙) + 畦。o ( 岛曙) ( 3 4 2 ) 这里的明“= ( c o s o ( x j 【，协) ) 。一1 ，w 1 _ n = ( s i n o ( x l ，协) ) 。一1 设n l ，m ，n 3 ，n 4 ，n 5 分别是边o a ，o b ，o c ，o d ，o e 的外法向，n o 是顶点o 处给定的方向( 见图3 2 ) 。 设是五条拟边界的公共点的函数值，磁，恐，五，恐，五分别是区域( 1 ) 和区 同济大学硕士学位论文 域( 2 ) ，区域( 2 ) 和区域( 3 ) ，区域( 3 ) 和区域( 4 ) ，区域( 4 ) 和区域( 5 ) ，区域( 5 ) 和区域( 1 ) 的 公共边上点的函数值，则有 咯= 曙= x 2 ，砰= 馏= 恐，婶= 嵋= 置 昭= 憎= 恐，孵= w = 墨曙= 谤= x o 在拟边界上，由法向导数的连续性建立如下方程： 鬻o n = 鬻，型o n 3 = 丝o n 3 ，塑o n 4 = 丝o n 4 ， ( 3 a 3 ) 2 舰 鬻= 筹，筹= 筹。筹= 筹 c s 删 肌5 锄5 加l 挑l 矾 根据分块矩阵的乘法，以筹为例，( 3 4 2 ) 可写成： 型1 = w _ 7 l ( l 矿1 + a 。7 2 曙+ a 挪w + a = 7 4 讶+ 也7 5 馆+ a 。t 6 咯+ a 川w ) + h 盘7 1 ( 岛7 l 矿1 + 玩7 2 曙+ b 7 3 订+ b z 7 4 呀+ 见7 5 w + 耽7 6 曙+ 毋竹w ) 由( 3 ，3 ，2 ) 式矿7 可由边界点的函数值曙，曙，曙。嵋，曙，w 表示，因此筹也可由 它们表示： ；导= p o + p 2 曙+ b 曙+ 只w 十r 嵋+ r 曙+ p 7 w 这表明( 3 4 3 ) 和( 3 4 4 ) 中的导数都可由边界点的函数值的加权线性和来表示。 对于d i r i c h l e t 边界问题，最后只要求解包含x l ，恐，甄，恐和j 岛为未知 量的代数方程组即可，如后面的例1 。对混合边界问题( 见图3 ，2 ) ，不妨设己 知c d 边上有边界条件焦+ 口= g ，那么j 岛，墨，恐，妫，蜀，x 5 和x c d 都是未知 量。c d 边上的导数的妊理方法同上，如例2 。通过求解代数方程组，可以求 得拟边界点和公共的顶点的函数值，回代入( 3 3 2 ) 式，整个区域上的问题即可 解决。 3 5 数值实验 例1 1 - - a v ( + 嘲2 v = 2 列裟 这是一个具有d i r i c h l e t 边界的椭圆型方程，求解域见图3 2 ，它的解析解为 0 ，可) 窖矿却+ x y 从解析解中，我们可求出 ( z ，) ，采用前面讲述的t d q d d 方法，将求得的 解与解析解比较，表3 1 给出了各子区域上的绝对误差的最大值( 用e r r o r ( i ) 表 示