（计算数学专业论文）连续有限元求解非线性动力系统的高效算法.pdf

摘要 用有限元方法求解各种微分方程，在科学研究、工程技术等方 面有广泛的应用。本文研究用二次有限元方法求解非线性动力系统 我们利用四级显示r - k 公式，提供单元右节点上的初值，并寻找了 不用重新计算新的函数值( 致) ，只需要用在计算右节点时已经得到 的( i = 1 ，2 ，3 ，4 ) 值的一个新的线性组合就可得到单元中点具有三 阶精度的值并证明了在利用四级显示r - k 公式时这是能达到的最 高阶精度。而且得出对任意次元，我们都可用类似的方法得到单元 内部具有同样精度的初值。这样利用这些初值和单元左端点的已知 值对用有限元方法( 或插值系数有限元法) 离散后的非线性方程组进 行几次简单迭代修正即可从而在求解离散后的非线性方程组时可 以不用牛顿法这样能大大地减少了计算量 最后我们还提出了一个可以提供初值的方法，即利用前面已经 计算的两个单元上的五个节点值( 含中点值) 构造四次插值。利用这 个插值函数为后面的单元提供中点和右节点上的初值而最前面两 个单元可以利用任意方法启动 数值实验结果表明上述计算方法有效而且此方法还可以用于 半线性抛物型方程的计算 关键词：动力系统，有限元方法，r o u n g e - k u t t a 方法，数值解 a b s t r a c t i n t h i sp a p e r ，w eg i v eo u rc o n s i d e r a t i o n si ns o l v i n gt h en o n l i n e a rd y - n a m i c a ls y s t e m sb yq u a d r a t i cf i n i t ee l e m e n tm e t h o d f i r s t l y , w eu s e4 - s t a g e e x p l i c i tr u n g e - k u t t am e t h o dt og e tt h er i g h tn o dv a l u eo ft h ee l e m e n t ，a n d t h e nw ep r o p o s ean e wf o r m u l at og e tt h em i d n o dv a l u eo ft h ee l e m e n tb yj u s t g i v i n gan e wl i n e a rc o m b i n a t i o no ft h ek sw h i c hw e r ec a l c u l a t e db e f o r e s o w en e e d n tc a l c u l a t en e wf u n c t i o nv a l u e l a s t l y , w ea l s op r o v et h a tt h en e w f o r m u l ah a st h r e eo r d e ro fa c c u r a c ya n dt h i si st h eb e s tw ec a ng e tt h r o u g h 4 - s t a g ee x p l i c i tr u n g e - k u t t am e t h o d i nt h el i g h to ft h i sm e t h o dw ec a r la l s o g e ts i m i l a rf o r m u l af o ra n yn o di nt h ee l e m e n ta n dt h e yh a v et h es o m ep r e - c i s i o na st h em i d n o df o r m u l aw eg i v e t h i sm e a n sw ec a nu s ef i n i t ee l e m e n t m e t h o do fa r b i t r a r yp o l y n o m i a lo r d e r f o rs o l v i n gt h en o n - l i n e a re q u a t i o n s w h i c ha x ed i s c r e t ef o r m u l a t i o no ft h en o n - l i n e a rd i f f e r e n t i a le q u a t i o nb yf i n i t e e l e m e n tm e t h o do ri n t e r p o l a t e dc o e f f i c i e n tf i n i t ee l e m e n tm e t h o d ，u s u a l l yw e u s en e w t o nm e t h o d b u tn o ww i t ht h em i d n o da n dt h er i g h tn o dv a l u e si n h a n d ，t o g e t h e rw i t ht h el e f tn o dv a l u eo ft h ee l e m e n tw ea l r e a d yk n o w ，w e u s et h i sv a l u e sa si n i t i a lv a l u e sf o ri t e r a t i v ea l g o r i t h m ，t h e nw eo n l yn e e dt o u s es i m p l ei t e r a t i v em e t h o df o rt h ed i s c r e t en o n l i n e a re q u a t i o n ss e v e r a lt i m e s t h i sm e t h o ds a v e sal o to fc o m p u t a t i o nt i m ea n dq u a n t i t yo fc a l c u l a t i o n w ec a na l s ou s ee x t r a p o l a t i o nm e t h o dt op r o v i d ei n i t i a lv a l u e s w ec a l - c u l a t et h ef i r s tt w oe l e m e n t sw i t ha n yk i n dm e t h o d ，f r o mt h et h i r de l e m e n t o n ，w eu s et h et w oe l e m e n t s ( i na l lf i v ep o i n t s ) j u s tc a l c u l a t e di m m e d i a t e l y b e f o r et h i se l e m e n t w ec o n s t r u c taq u a r t i ci n t e r p o l a t i o np o l y n o m i a la n dl e t i tp r o v i d e st h em i d n o d a n dt h er i g h tn o dv a l u e so ft h en e x te l e m e n t n u m e r i c a le x p e r i m e n t si n d i c a t et h a tt h ea b o v em e t h o dw ep r o p o s e di s e f f e c t i v e a n dt h i sm e t h o dc a na l s ob eu s e df o rs o l v i n gs e m i - l i n e a rp a r a b o l i c e q u a t i o n k e yw o r d s ：d y n a m i c a ls y s t e m ，f e m ，r u n g e - k u t t am e t h o d ，n u m e r i c a ls o - l u t i o n l i 湖南师范大学学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下，独 立进行研究工作所取得的成果。除文中已经注明引用的内容外，本 论文不含任何其他个人或集体已经发表或撰写过的作品成果对本 文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明 本人完全意识到本声明的法律结果由本人承担。 学位论文作者签名：奄卅扩铭如l ，年莎月厂日 湖南师范大学学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定， 研究生在校攻读学位期间论文工作的知识产权单位属湖南师范大学 同意学校保留并向国家有关部门或机构送交论文的复印件和电子版， 允许论文被查阅和借阅。本人授权湖南师范大学可以将本学位论文 的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印 或扫描等复制手段保存和汇编本学位论文 本学位论文属于 1 、保密口，在z年解密后适用本授权书 2 、不保密f ( 请在以上相应方框内打“ ”) 作者签名：洲扩够日期：扣尸年占月 日 导师签名2 纩白仍蜡日期：o 产月易日 4 1 连续有限元求解非线性动力系统的高效算法 1 引言 现代科学技术中出现了大量的非线性微分方程，构成了非线性 科学和非线性数学研究的核心内容之一现实生活中许多现象的物 理背景就是非线性动力系统非线性动力系统的研究以其理论的深 刻性和应用的广泛性成为国际数学、力学、物理以及金融经济领域 中最活跃的前沿 与线性系统理论不同，非线性系统目前尚无统一的研究方法由 于绝大多数确定性非线性系统的微分方程没有求得精确解的有效方 法，自2 0 世纪2 0 年代以来，人们发展了多种近似方法，如小参数 法、坐标变换法，多尺度法、慢变参数法，k a m 法、等效线性化、谐 波平衡法，r i t z - g a l e r k i n 等【7 , 8 ，2 3 】 本文主要是用二次有限元方法求解非线性动力系统的初值问题 在具体处理非线性项时，用到了插值系数有限元方法。它对常微和 偏微都有效，此法可以使求解工作量显著减少。而且还保持着最佳阶 收敛性和超收敛性。对这部分工作研究中，国内外有很多学者做了大 量工作，国内的有陈传淼、张乃莹等国外的有z l a m a l ，l a r s s o n ，t h o m e e 等。见f 1 , 3 ，4 ，6 ，1 7 对于非线性常微分方程初值问题的有限元求解方法【3 ，1 3 】。现考 虑如下方程 仳t = f ( t ，u ) ，乱( o ) = u o ，0 t 。( 1 - 1 ) 在单元( t k ，t k + - ) 上用m 次有限元逼近厶u = 凳oc j x j ，可得如下有 关墨，j = 1 ，2 ，m 的非线性方程组 ( d , ( i h u ) ，已) = ( f ( t ，h u ) ，已) ，已p m 一1 ，i = 1 ，2 ，m 其中凰是单元的已知左端点值。一般取6 = ( t t k ) 卜1 ，i = 1 ，2 ，m ， 用n e w t o n 法求解，需要多次计算导数矩阵 眠川= 警k 卸= 1 2 ，m 工作量很大。而如果采用插值系数有限元法m z l a m a l ( 1 9 8 0 在讨论 1 硕士学位论文 半线性热传问题时提出，即对f ( t ，) 采用同样的分片插值函数) h f ( t ，u ) = c j f ( t j ，码) j = o 则得到关于未知数x j ，j = 1 ，2 ，仇的仇阶非线性方程组 mm ( 蟛( ) ，引= ( 咖，鳓f ( 如，x j ) ，i ，歹= 1 ，2 ，m j = oj = o 写成矩阵形式为 k x = m f ( x ) 】+ b ，【f ( x ) 】= ( f ( t l ，x 1 ) ，f ( t 2 ，x 2 ) ，f ( t m ，叉- m ) ) t ( 1 - 2 ) 这里7 n m 阶方阵k = ( 以( ) ，】与m = ( 勿，鳓】可以一次性计 算好，而n e w t o n 迭代只是在玛与f ( t ，恐) 之间进行，计算切矩阵 m d f ( x ) 】= m d i a g ( d j f ( t j ，玛) ) 的运算比较简单 求解非线性方程组的主要方法是n e w t o n 法，它的计算量很大从 上面可知在用插值系数有限元法离散时已经能很好的减少计算量 而在本文我们研究了一种更直接的计算方法这里以常用的2 次连 续有限元法为例。如果能以某种方式提供单元右节点和单元中点的 精度比较好的初值，而且计算代价比较小，然后直接代人方程( 1 - 2 ) 作几次简单迭代修正即可，从而可以不用n e w t o n 法求解( 因而不 必求切矩阵) 这样就能更加减少计算量这是本文的出发点 在本文中我们首先考虑了用4 级4 阶显式r u n g e - k u t t a 法得到单 元右节点值巩+ 1 它具有4 阶精度，比如用如下的经典公式 + h j ( k 1 + 2 k 2 + 2 k 3 + k 4 ) ，心) + h j 2 ，u j + h j k l f 2 、 + h | 2 ，钍l 七h j k 2 | + ，+ h j k 3 ) 我们研究了用k ( i = 1 ，2 ，3 ，4 ) 的某种线性组合来提供中点初值+ 。z ， 使其精度尽量高我们验证了常用的几种4 级r - k 公式，并重新考 察了4 级显示r - k 公式，得出了在此情况下能得到的最好精度是三 2 嘶似似泐 =， r = = = = 啪虬段 ，i_-i-，、_ii【 连续有限元求解非线性动力系统的高效算法 阶要使中点也具有4 阶精度，就算在改变尬的情况下也不成立 但如果改变，甄，则能达到4 阶精度，不过其计算量差不多是对中 点重新用了一次r - k 法计算，计算开销较大 4 级显示r - k 公式的一般形式为 iu ( t + h ) 一钍( t ) = h ( e l k l + e 2 k 2 + e 3 k a + e 4 k , ) ik 1 = f ( t ，u ) k 2 = m + a 2 h ，a 2 k l 危) l 蚝= f ( t + a 3 h ，仳+ ( n 3 6 3 2 ) 风+ b 3 2 k ： h ) li 0 ，q = q ( x ) 0 记 小= 6 帅 + q u v d x 加川= z 6 加d x v = 乱l t | ，钍7 l a ，6 】，u ( a ) = t | ( 6 ) = o ) 这样两点边值问题可化为：求u v 且满足下式 a ( u ，u ) = ( ，u ) ，帕v( 2 - 1 ) 问题( 2 - 1 ) 的有限元逼近是构造y 的有限维子空间s ncv ，考虑如 下有限维问题( 或称离散问题) ： 记j = 陋，6 】，给定j 的一个任意剖分：口= x o x l x n - 1 x n = b 。e i = 【x i - 1 鼢】称为一个单元，h i = 甄一孔- 1 h = m a x l 。 丸) 其中 踟是由剖分区域上的简单函数( 例如分片多项式) 组成( 小支集基函 数构成的有限维子空间) 我们用这些简单函数去逼近原问题的解 对应于上面的例子，其中 8 n = 乱 c a ，b l u , , i 岛p m ( e i ) ，1 i n ) 有限元方法中关健的一步是：把边值问题化为变分问题之后，对 求解区域q 作剖分，使q 成为有限个“单元”的和，然后在每一个 单元上作未知函数的某种多项式插值，并且使它们在相邻单元的公 共边界上满足某种连续性条件，以保证用这种分片插函数组成的有 限维函数空间踟是未知函数解空间v 的子空间。这种子空间满足 微分方程边值问题的本质边界条件，从而克服了古典的r i t z - g a l e r k i n 方法在构造子空间鼬时所遇到的困难此外，由于采用了分片插值 的办法来构造子空间鼬的基函数，可以使得形成的近似变分方程 6 得蛳 使踟卜 批 ，_，1-一， 连续有限元求解非线性动力系统的高效算法 ( 即线性代数方程组) 的系数矩阵时其计算量大为减少特别，系数 矩阵还是一个稀疏矩阵。如果对结点编号采取某些措施，还可以使 它成为一个带状的稀疏矩阵，这对实际计算是大有好处的，因为可 以节省计算机的存贮量【4 ，1 0 ，1 4 ，1 9 有限元的计算是逐单元进行的。当所有单元都计算完以后，再进 行组装成总矩阵，并处理本质边界条件。最后求解所得的线性方程 组。它很容易形成程序在计算机上实现。基于上述原因，有限元方 法现在已经成为一种求解微分方程的重要方法之一详见【1 ，2 】 2 2 插值系数有限元 有限元方法完全可以用于解非线性问题，且具有解线性问题相 同的好精度。但数值求解非线性方程组时，由于n e w t o n 法的迭代过 程中要多次对不同的值u j 计算切矩阵，其工作量很大为此人们提 出了一种叫做插值系数有限元方法的简化算法它用在离散未知函 数时的基把非线性项也进行插值可以显著简化计算切矩阵从而 节约大量计算时间而且它还保持着最佳阶收敛性与超收敛性此 法对椭圆，抛物与双曲问题都有效【1 ，4 ，1 5 ，2 0 】 为了定义插值系数有限元方法，先考察如下形式的半线性二阶 椭圆问题 a u + ( u ) = g ( z ) i nq ，u = 0mf( 2 - 2 ) 其弱形式为寻求u 日1 满足 a ( u ，u ) + ( ，( u ) ，u ) = ( 9 ( z ) ，u ) ，v v 月吾( q ) 其中qcr d 是边界为r 的有界域，且双线性型a ( 乱，u ) = f a ( a i j ( x ) d u d j v + a ( x ) u v ) d x 是岛一强制的设s 是维有限元子空间且，【) 搀。 是它的一组基因此( 2 - 2 ) 的有限元解u s 可以表示为v ( x ) = 羔1 伤s ，再取钉= i = 1 ，2 ，n ，则有 n n a ( ，仇) v j + ( ，( 仍( z ) ) ，协) = ( 夕，忱) ，i = 1 州2 ，n ( 2 - 3 ) j = lj = l 当用n e w t o n 迭代法求解非线性方程组( 2 - 3 ) 时，每一次迭代，都必需 7 硕士学位论文 计算它的切矩阵 a ( ，协) n n + 妒( ( z ) u j ) v j ，忱) n n ( 2 - 4 ) j = l 但这个切矩阵的计算却依赖于所选定的u 值。因此用n e w t o n 法求解 ( 2 - 3 ) 时，人们不得不多次计算这些切矩阵的值，其工作量是相当大 的 为了简化算法，可以采用一种朴实但非常优美的思想，即对函 数f ( u ) 本身作插值i h f ( u ) ，用i h f ( u ( x ) ) = 竺。妒j f ( u j ) s 台代替 ，( u ( z ) ) 在q 上定义次插值系数有限元u 满足 nn a ( ，似) + f ( u j ) ( q o j ，忱) = ( 9 ，忱) ，i = 1 ，2 ，n ( 2 - 5 ) j = lj = l 用n e w t o n 迭代法求解非线性方程组( 2 5 ) 时，切矩阵为 陋( 仍，仍) n x n + ( 仍，仇) i ，1 。x l r t d i a g ( f 7 ( 巩) ，7 ( ) ，7 ( ) ) ( 2 6 ) 其中刚度矩阵 a ( v j ，仇) 与质量矩阵【( 仍，协) 】可以一次性计算好，而 计算第二项只是简单的乘法因此迭代计算只是简单地在节点值 和，( ) 之间完成，因此整个计算工作量将大大减少以上方法即为 插值系数有限元法f 2 0 ，2 l 】 2 3 有限元求解常微分方程 下面我们介绍用连续有限元方法求解常微分方程的初值问题考 虑如下初值问题 乱7 = f ( t ，u ) ，0 t 。，u ( o ) = u o( 2 - 7 ) 其中f ( t ，仳) 可以是线性( f ( t ，u ) = a ( t ) u + b ( t ) 一种重要的特殊情况) 或 非线性 将区间j = ( 0 ，t ) 作剖分磊：t o = 0 t l t 2 t l v = t 记单 元厶= ( t j ，t j + 1 ) ，步长为b = t j + l t j ，最大步长为h = m a x h j 记连续 且分片m 次有限元空间 s = u ： c ( j ) ， i 尸i m ，j = 0 ，1 ，2 ，一1 ) 8 连续有限元求解非线性动力系统的高效算法 解初值问题的特点( 与边值问题不同) 是可以逐个单元求解。连续有 限元空间s 的元素u 在乃虽然是m 次多项式，但是左端点的值 u ( t j ) = u j 已知，因而在易上只有m 个自由度现定义m 次连续有 限元解u s 在每个单元易上满足 ， ( u 7 一f ( t ，u ) ) r l d t = 0 ，j = 0 ，1 ，2 ，一，一1 ，u ( o ) = 咖 ( 2 - 8 ) j1 3 这里，7 是任何m 一1 次多项式，特别地可取 町= ( t t j ) 4 ，i = 0 ，1 ，2 ，仇一1 于是得到一个m 阶非线性方程组，它们唯一地确定了u 的值，特别 地确定了节点值u ( o ) = u j + l 的值见【1 ，1 3 】 用连续有限元求解常微分方程的初值问题它是单步法，可以 改变步长，自开始如果是解线性问题，那么它每步只需要求解一 个m 个未知数的方程组( 非线性情况是求解非线性方程组，可以通 过预估效正格式来减轻计算量) 。而本文就是针对非线性情况，并结 合插值系数有限元方法。在用二次元求解时，通过较少的计算代价 提供较好的单元中点和右节点初值，然后对用有限元离散的方程做 几次简单迭代修正，以此来减少计算量 9 连续有限元求解非线性动力系统的高效算法 3 关于r _ k 法提供中点值的研究 由上面的介绍我们可知：在用有限元方法或插值系数有限元法 求解非线性动力系统时最后离散得到的方程组是非线性的我们 对非线性方程组的求解，一般是用n e w t o n 法求解，而n e w t o n 法需 要相当大的计算量因为每次迭代都需要计算切矩阵不过用插值 系数有限元法离散具有相对简单的切矩阵，从而能很好地减少计算 量。在本文，我们的目的是在用二次元求解常微分方程时：以小代 价提供单元右节点和单元中点上的预估初值，结合单元左节点的已 知值，然后把这些值代人离散后的非线性方程组作几次简单迭代修 正即可。从而不用n e w t o n 法求解非线性方程组。这样可以大大减 少计算量。为此我们考虑了以下两种提供预估初值的方式： u pu 川li j | t j lt j 刈1t 赣 图3 - 1 ：二次元单元示意图 一、我们利用4 级显示r - k 法得到单元右节点值巩，再利用在 计算时的4 个k ( 江1 ，2 ，3 ，4 ) 的值的某种线性组合，或者允许甄 改变，以此希望能得到v j 一。：的一个尽量精确的值，以它们作为单 元右节点和中点的预估初值此法可自开始，可取变步长 二、我们还考虑了用前面两个单元已经计算出来的值作一个四 次插值函数。利用此函数的外插值给出下一个单元的右节点和中点 值作为预估初值。此方法为多步法计算量也很少但是它不能自 开始必需利用其它的方法提供前面两个单元上的计算 下面我们就来考虑怎样得到单元中点的初值，而且使其具有尽 量高的精度。这样在求解非线性方程组时，我们希望：1 、如果用牛 顿法，以此为初值可以减少迭代次数2 、也可以直接以此初值对方 硕士学位论文 程做简单迭代。这样都能减少计算量 3 14 阶r - k 法的重新考察 对于微分方程 u t = f ( t ，让)( 3 - 1 ) 我们可以用4 级4 阶显式r u n g e - k u t t a 法给出右节点值，然后考虑给 硷，鲍，蚝，致一个新的线性组合，希望其能使单元中点也具有4 阶 精度即希望下式也有4 阶精度，或使其尽量精确 u ( t + 百h ) 一u ( t ) ：鲁( 9 1 k 1 + 夕2 鲍+ 夕3 k 3 + 9 4 k 4 ) 常用的4 级4 阶公式有： 1 、4 级4 阶古典显式r u n g e - k u t t a 公式： iu j + i = + h ( k 1 + 2 1 , ：2 + 2 1 ( 3 + k 4 ) ik 1 = f ( t j ，吩) k 2 = ，( 巧+ 九2 ，嘶+ h k l 2 ) ik 3 = ，( 岛+ 2 ，+ h k 2 2 ) ik 4 = f ( t j + h ，u j + 娲) 2 、4 级4 阶显式k u t t a 公式： i 呦+ 1 = u j + h ( k 1 + 3 k 2 + 3 k 3 + k 4 ) ik 1 = ，( 岛，吻) k 2 = ，( 如+ 危3 ，吻+ h k l 3 ) i 9 3 = f ( t j + 2 h 3 ，一h k l 3 + h k 2 ) 【甄= f ( t j + h ，嘶+ h k l 一h 9 2 + h k 3 ) 3 、4 级4 阶显式g i l l 公式： f + l = 嘶+ l h ( k 1 + ( 2 一以) 鲍+ ( 2 + 钜) 硒+ 甄) lk 1 = f ( t j ，) t ( 2 = f ( t j + 2 ，+ h k l 2 ) i 蚝= ，( 岛+ 2 ，u j + k l + ( 1 一孚) ) 【甄= ，( 岛+ ，一乎 鲍+ ( 1 + 乎) 蚝) 1 2 连续有限元求解非线性动力系统的高效算法 我们用这些公式来得到单元右节点值的同时，下面我们考虑对已经 给定的k 。，凰，尬一个新的线性组合使单元中点也具有4 阶精度， 或使其尽量精确 h 吩+ 昙= u j + 等( 夕1 硒+ 9 2 k 2 + 夕3 配+ 9 4 ) ( 3 - 2 ) 1111 + 一吩2 壶u 危+ 汐h 2 + 去怕3 + 壶u 4 h 4 + o ( h 5 ) ( 3 - 3 ) 对以上两式，我们分别对上面的三种r - k 公式展开计算得到以下系 数对应关系。其中第一列是展式对应导数项( 有关导数的算子表示见 三种常用4 阶r - k 公式中点展开系数应该满足的关系 比较项数值经典r - k 公式k u t t a 公式g i l l 公式 f h 1 夕l + 9 2 + 夕3 + 9 49 l + 9 2 + 9 3 + 夕4吼+ 9 2 + 9 3 + 9 4 d f h 2 1 ( 夕2 + 9 3 ) + 9 4 9 2 + 3 2 - 9 3 + 9 4互1 ( 9 2 + 9 3 ) + 9 44 d 2 ， 3 l - 一1 ( 9 2 + 9 3 ) + 9 45 夕2 + 9 3 + 9 4 ( 9 2 + 9 3 ) + 9 4 1 2 d f h 3 1 j 9 3 + 9 4；9 3 + 9 4学卯+ ；9 4 2 4 d 3 ，危4 1 ( 夕2 - i - 9 3 ) + 9 41 9 2 + 嘉9 3 - t - 9 4百1 ( 9 2 + 9 3 ) + 9 43 2 九d 2 ，九4 1 ；卯+ ：肌 ；船+ ；9 4学9 3 + 9 4 9 6 d f d l , h 4 1 9 3 + 驭；驰+ 9 4学9 3 + ；9 46 4 ( 九) 2 d f h 4 l 9 4 9 4 9 41 9 2 从上表可以看出，要想使按此展开的中点公式具有4 阶精度，不 可能，因为每组方程都无法全部满足，但是我们却可以得到满足3 阶 精度的公式，因为这样只需要考虑每组方程的前面4 个而g 。，9 2 ，9 3 ，9 4 刚好有4 个自由度。同时我们也看到，在已经满足了三阶精度的情 况下，不同的g t 选择，对应h 4 的各项系数不同，这让我们觉得是不 是可以重新考察4 级龙格库塔公式，找到一组a 。，a 3 ，使此时不但节 点上有四阶精度，中点上也有或者说至少也能让中点的4 阶项系 数尽量少 下面我们就让所以系数都自由，对方程( 孓1 ) 的4 级显示龙格库 塔公式进行重新考察f 9 - 1 1 ，2 2 111 乱( + h ) t 正( t ) = u h + 去u h 2 + 去u 肿h 3 + 击u ( 4 ) 4 + o ( h 5 ) ( 3 4 ) 1 3 硕士学位论文 其中 u ( t + h ) 一u ( t ) = h ( e l 局+ e 2 k 2 - i - e 3 玛- i - e 4 k 4 )( 3 - 5 ) 引进以下算子，以方便表述参见 9 , 1 1 d = 象+ ，熹 ( 3 - 6 ) 碘淹o f = | t + lf 。，d | u = f u t + | | 。 d 2 f = 貉+ 2 f 爰羞+ ，2 羞篆= 厶+ 2 f f t u + ，2 凡u d 3 f = = ，赢+ 3 f 厶。+ 3 2 ，t u u + f 3 厶伽 下面计算函数u = u ( t ) 的对应导数( 用，表示的) 仳7 = f 畦| = t 七 u 叁d f o = | 魄+ 2 f h + f 2 | 呲+ | t l u + | | l | 浮垒d o | + a d f u ( 4 ) = 全d 3 f + 3 d f d f , , + d 2 f + ( ) 2 0 f 把上面各式代入( 3 4 ) 式就得 钍( t + h ) 一u ( t ) = f h - i - d f h 2 + d 2 f + f d f h 3 + 西1 【d 3 f + 3 d f d a + 厶d 2 ，+ ( 九) 2 0 f 1h 4 + o ( h 5 ) 上面是函数的t a y l o r 展开下面我们利用二元函数的t a y l o r 公式来 展开k ( 江1 ，2 ，3 ，4 ) 得： k 1 = f 尥= f + a 2 d f h + 互1 u 2 2 d 2 f h 2 + 石1 u 2 3 d 3 f h 3 + o ( h 4 ) k 3 = f + a 3 d f h - t - ( a 2 b s 2 d + 互1 “3 2 d 2 f ) h 2 - t - ( 虿1 2 “2 2 d 2 f + 口2 0 3 6 3 2 d ，d 丘- t - i l u 3 3 d 3 f ) h 3 + o ( h 4 ) 致= f + a 4 d f h - t - 互1u 4 2 d 2 f + ( 5 4 ：2 + b 4 3 a 3 ) f u d f h 2 + 石1 u 4 3 d 3 f + ( 互1 b 。2 + 互1 岫u 3 2 ) 厶d 2 ，+ a 4 ( b 4 2 口2 + b 4 3 a 3 ) d f d f u - t - b 4 3 b 3 2 a 2 ( f u ) 2 d f h 3 + o ( h 4 ) 将甄的展式代入( 3 - 5 ) 式得 1 4 d，副 k k +k m 配+ b k h 渺一，如以 m 一 一 k m 彬 + + + u u u 肌肌肌 2 国 q ，删m u + + + “r “r = = = = 甄鲍蚝甄 ，-_-_j、-ii-l 连续有限元求解非线性动力系统的高效算法 u ( t + h ) 一u ) = h ( e l k l + e 2 恐+ e 3 + e 4 硒) = ( e l + e 2 + e 3 + e 4 ) 厂 + ( e 2 a 2 + e 3 a 3 + e 4 a 4 ) d f h 2 + i 烈1c 2 “2 2 + e 3 a ；+ e 4 a 2 a ) d 2 f + ( e 3 a 2 b 3 2 + e 4 a 2 b 4 2 + e 4 0 3 h a ) f , , d f h 3 + l ( e 2 0 ；+ e 3 0 ；+ e 4 a s ) d a f + ( e a a ；b s 2 + e 4 a ；b 4 2 + e 4 0 錾i 鹞) 逖+ ( e 3 a 2 a 3 b 3 2 + e 4 a 4 a 2 b 4 2 + e 4 a 4 a 3 垴) d f d 凡+ e t a 2 b 4 3 b 3 2 ( 丘) 2 d ，ih 4 + j o ( h 5 ) 比较上面函数的t a y l o r 展开和r - k 法展 显示4 级4 阶r - k 法 程，1 0 个未知量。 系数应该满足的方程 开两式的系数就得下表： 组它们一共是8 个方 作比较的项函数展开式系数r - k 法展开系数 f h 1 e l + e 2 + e 3 + e 4 d f h 2 l 2 e 2 a 2 + e 3 a 3 + e 4 a 4 d 2 厂危3 1 e 2 a l + e 3 a ；+ e a a 2 4 )6 九d f h 3 1 e s a 2 b 3 2 + e 4 a 2 b 4 2 + e 4 a 3 b 4 3 6 d 3 ，危4 1 e 2 a 2 + e 3 a i + e 4 a s 4 )2 4 厶d 2 ，九4 1 l ( e 3 a 2 b 3 2 + e a a 2 b 4 2 + e 4 a b 4 3 1 2 4 d f d f 。h 4 3 e 3 a 2 a 3 b 3 2 + e 4 a 4 a 2 b 4 2 + e 4 a 4 a 3 b 4 3 2 4 ( 九) 2 d f h 4 l e 4 a 2 b 4 3 b 3 2 2 4 下面我们考虑它的求解 e 1 + e 2 十e 3 + e 4 = 1 e 2 0 2 + e 3 0 3 + e 4 n 4 = 一2 工 e 2 。2 + e 3 a 2 + e 4 a 2 百1 e 2 口3 2 + e s a i + e 4 a 3 - - _ 石1 e s 。z b 3 2 + e 4 n z 6 4 z + e t 凸s = 丢 e 3 n ：。a b 3 2 + e 4 a 4 a 2 6 4 2 + e 4 a 4 a 3 6 4 3 = 丢 e 3 n 2 2 b 3 2 + e 4 a 2 b 4 2 + e 4 a ；b 4 3 = 西1 e 4 。z 6 4 。6 3 z = 去 1 5 ( 孓7 ) ( 3 - 8 ) ( 3 - 9 ) ( 3 - 1 0 ) ( 3 - 1 1 ) ( 3 - 1 2 ) ( 3 - 1 3 ) ( 3 - 1 4 ) 硕士学位论文 由( 3 - 8 ) ( 3 - 9 ) ( 3 - 1 0 ) 式，把a ，i = 2 ，3 ，4 看作已知量，可求得岛，t = 2 ，3 ，4 如下： ，一曼竺竺! 二兰f 塑竺! ! 一 1 2 a 2 ( a 3 一0 2 ) ( 0 4 一a 2 ) 6 0 4 a 2 4 ( a 2 + a 4 ) - 4 - 3 e 3 = 2 。1 2 a 3 ( a 3 - - 。a 。2 ) ( a 4 - - a 。3 ) 6 a 2 a 3 4 ( a 2 + a 3 ) + 3 e 4 = 2 1 2 a 4 ( a 4 - a ，_ ) ( a 4 - 二a 3 ) 一 由( 3 1 1 ) 式x a 2 一( 3 - 1 3 ) 式，可得6 4 3 6 4 s = 瓦丽2 a 2 而- 1 由6 躬再结合( 3 - 1 4 ) 式可得6 3 2 一 a 3 ( a 2 一a 3 j b 3 2 5 2 2 a 2 ( 2 a 。2 - 。1 2 一) 把求得的6 4 3 ，6 3 2 代入( 3 - 1 1 ) 式，可得 蜘堕篆老端掣 把以上解得的6 4 3 ，6 3 z ，b 4 。代入( 3 - 1 2 ) 式，可得 a 4 = 1 把a 4 = 1 回代以上所解出来的表达式中，再由( 3 - 7 ) 式可求得e 。以 下为方程组的解( 把a 2 ，a a 看作自由变量) ： e 1 = 互1 + 1 - 瓦2 ( a 面2 + - a a ) e 2 = 蕊2 a 一；1 口2 - ) 1 雨 e 3 = 丽1 一- 如2 a ) 2 f 丽 e 42 a 421 5 3 2 = 6 4 2 = b 4 3 = = 连续有限元求解非线性动力系统的高效算法 a 2 ，a 。取不同的值，可以得到不同的r - k 公式。用它们计算节点值 钍( + h ) 都有4 阶精度。我们提出一个问题：这些4 阶r - k 公式中 是否有一个特殊情形，利用它的k 。，恐，尥，心的一个新的线性组合， 可以得到高精度的中点值位( + h 2 ) ? 这样就可以为二次有限元提供 好初值。 3 2 中点值的新组合公式 下面我们推导中点公式。k 。，鲍，恐，尬的表达形式不变，我们只 给它们一个新的线性组合，希望能得到中点的一个尽量精确的公式。 推导过程类似于前面节点公式函数在中点处的t a y l o r 展开和r - k 展开如下： u ( t + ) 一u ( ) = t 上7 互h + u ( ) 2 + , i t i t t ( 皇) 3 + 去u ( 4 ) ( ) 4 + o ( h 5 ) = ；，愚+ d ，危2 + 去l d 2 f + f , , d f i h 3 + 击l d 3 厂+ 3 d i d f , , + 九d 2 厂十( 九) 2 d f i h 4 + o ( h 5 ) 札( t + 鲁) 一u ( t ) = ；( 9 1 k + 9 2 恐+ 9 3 k 3 + 9 4 k 4 ) 可得在中点处应该要满足的系数对应关系： 作比较的项函数t a y l o r 展开系数r - k 法展开系数 f h 1 9 1 + 9 2 + 9 3 + 9 4 d f h 2 1 4 9 2 印+ 9 3 a 3 + g a a 4 d 2 ，忍3 1 ( 9 2 a ；+ g a a ；+ 班o i )2 4 a d f h 3 1 9 3 a 2 b 3 2 + 9 4 a 2 b 4 2 + 9 4 a z b a 3 2 4 d 3 ，危4 1 ( 9 2 a ；+ 9 3 a ；+ 9 4 a a 4 ) 1 9 2 d 2 厂 4 1 ( 9 3 a 2 b s 2 + 9 t a 2 b a 2 + 9 4 a b 4 3 )1 9 2 d f d f 。h 4 3 g a a 2 a a b 3 2 + 9 4 a 4 a 2 b 4 2 + g a a a a a b 4 3 1 9 2 ( 九) 2 d f h 4 1 9 4 a 2 b t a b 3 ： 1 9 2 对照节点方程组和上面中点处应该满足的方程组，每一组方程， 都有1 0 个末知量，8 个方程。而由节点方程组的求解过程我们知道， 每组方程中的方程都是独立的而且在节点方程组满足4 阶精度的 情况下，它只剩下两个自由度，即a 。，a 3 ( 当然也可以取其它的变量为 自由变量，只是不能取a 4 ，因为从求解过程中，我们可知他是一个常 1 7 硕士学位论文 数，a 。= 1 ) 而要保持( i = 1 ，2 ，3 ，4 ) 不变那么中点方程组的自由 度就只有6 个而方程有8 个所以不可能在保持所有k ，i = 1 ，2 ，3 ，4 不变的情况下，节点和中点都有4 阶精度。下面我们考虑适当把条 件放宽，以使中点具有尽量高的精度 方案一：保持0 = 1 ，2 ，3 ) 不变，而允许甄改变，这样我们就是 想利用计算5 个k 值得到节点和中点都有4 阶精度的公式在算 出节点值后，利用已经有的k ( i = 1 ，2 ，3 ) ，再加上专门为中点准备的 砭这样做的目的是在计算量增加不多的情况下( 多算一个函数值) ， 来增加节点公式的自由度考虑到中点处应该满足的方程组 夕1 十9 2 + 卯十9 4 9 2 口2 + 9 3 a 3 + 9 4 a 4 9 2 n ；+ 9 3 a ；+ 9 4 a i 9 2 a ；- - b9 3 a 2 + 夕4 g a a 2 b 3 2 + 9 4 a 2 b 4 2 + 9 4 a 3 6 4 3 g a a 2 a 3 b a 2 + 9 4 a 4 a 2 b 4 2 + 9 4 a 4 a 3 b 4 3 9 3 0 ；6 3 2 + 9 4 a ；b 4 2 + g a a ；b a a 9 4 a 2 5 4 3 5 3 2 由式( 3 - 1 6 ) ( 3 - 1 7 ) ( 3 - 1 8 ) ( 它们与6 3 2 无关，b 3 2 ，a 3 决定