（计算数学专业论文）一类流体混合模型的数值解法.pdf

摘要 摘要 本文主要研究了一类流体混合模型的数值解法该流体混合 模型是关于可收缩间叶细胞组织变形的模型，是由非线性的双曲 型和椭圆型方程组成的混合方程组 本文考虑一维情形，采用了三种方法对该方程组进行求解：广 义差分法，迎风格式和中心有限差分法以及广义迎风差分法本文 给出了两种不同模型的求解，并分别用数值实验表明了三种方法 的有效性 前言里主要介绍了模型的背景和求解方法，第二章到第四章 分别给出了三种方法的具体离散格式和数值实验最后一章给出 了结论和展望 关键词：流体混合模型；迎风格式；中心有限差分方法；广义差分 法：广义迎风差分法 a b s t r a c t a b s t r a c t i nt h i st h e s i s ，w es t u d yn u m e r i c a lm e t h o d sf o raf l u i dm i x t u r e m o d e l t h em o d e li sa b o u tc o n t r a c t i l em e s e n c h y m a lt i s s u ed e f o r m a t i o na n di sam i x e ds y s t e mo fn o n l i n e a rh y p e r b o l i ca n de l l i p t i cp a r t i a l d i f f e r e n t i a le q u a t i o n s w ec o n s i d e rt h eo n ed i m e n s i o nc a s ea n da d o p tt h r e em e t h o d sf o r t h es y s t e m ：ag e n e r a l i z e dd i f f e r e n c em e t h o d ；t h eu p w i n ds c h e m ea n df i n i t ed i f f e r e n c em e t h o d ；t h eg e n e r a l i z e du p w i n dd i f f e r e n c em e t h o d t w o d i f f e r e n tc a s e sa r ec o n s i d e r e di nt h i sp a p e r n u m e r i c a le x p e r i m e n t sa r e p r e s e n t e dt oi l l u s t r a t et h a to u rp r o p o s e dm e t h o d s a r ee f f i c i e n tf o rs o l v i n g t h ep r o b l e m i nt h ei n t r o d u c t i o n ，w ei n t r o d u c et h eb a c k g r o u n di n f o r m a t i o n f r o m t h es e c o n dc h a p t e rt ot h ef o u r t hc h a p t e rt h es p e c i f i cd i s c r e t ef o r m a t sa n d n u m e r i c a le x p e r i m e n t so ft h et h r e em e t h o d sw i l lb eg i v e n i nt h el a s t c h a p t e rw eg i v ec o n c l u s i o na n dp r o s p e c t s k e y w o r d s ： f l u i dm i x e dm o d e l ；u p w i n ds c h e m e ；t h ec e n t e r e df i n i t e d i f f e r e n c em e t h o d ；t h eg e n e r a l i z e dd i f f e r e n c em e t h o d ；g e n e r a l i z e du p 。 w i n dd i f f e r e n c em e t h o d 一i v 第1 章前言 第1 章前言 过去的几十多年里，人们已经很好地把连续介质理论和关于固体流体运 动变形的化学原理结合起来应用于生物学组织和细胞的研究中这些研究 主要集中在组织的生长和建模、力学、非反应混合物的传播以及组织工程学 上在这一领域内，模型建立的本身比数值计算更重要1 9 8 6 年，m d e m b o 干l e h a r l o w 基于反应的内部敏锐流这一概念，提出了关于可收缩生物聚合体网状 系统的物理模型【1 】1 9 9 6 年，d m a n o u s s a k i ，s r l u b k i n 等四人提出了体外血 管网络构成的一种机械模型【2 】近来，又有不少作者基于人类软组织变形建立 了模型【3 ，4 ，5 】，然而数值计算方面的报道却很少在【6 】中，l u b k i n 和j a c k s o n 币l j 用 了带有一步预处理校正的有限差分方法：l a x w e n d r o f f 方法，此方法被用于预 处理和校正步上来输送细胞和纤维的体积分数及新生物的细胞纤维态分数 在【7 】中，x h e 和d m e b o 在带有自适应网格的圆柱几何体上应用了有限元方法 s t a s t n a 幂u 用显式有限差分格式( 空间二阶，时间一阶) 计算了【8 】中与时间有关的 情形然而，以上文章中作者都没有给出所用方法的精度分析结果 本文研究的是一维可收缩间叶细胞组织变形的数学模型，该模型的建立 参见文献【l ，2 ，6 】另外在【9 】中给出下面三个假设： ( 1 ) n 叶组织细胞是由水成状态和细胞纤维状态两个状态组成 ( 2 ) 水成状态具有s t o k e s 流的特征 ( 3 ) 细胞纤维状态的压力可由相关时间尺度上的永久变形来消散，且也可当 做s t o k e s 流来处理在此基础上建立的一维数学模型为： 赛+ 掣_ 0 ，o z l ， ( 1 o 1 ) 杀( 等塑o z 刊_ o a z 、口 。7 ( 1 0 2 ) 嘉( 2 m 塞- p + o 妒+ t r l n ( 1 叫) _ 0 ， ( 1 o 3 ) 第1 章前言 其中口表示细胞和纤维的体积分数，u 表示细胞纤维态的速度，妒表示牵引 系数，p 是压力，砂表示收缩系数，盯表示膨胀系数，m 为细胞纤维分数的粘性系 数由此可知，0 彘 ( 1 0 9 ) 本文其他部分安排如下：第二章利用广义差分法对方程组( 1 0 5 ) ( 1 0 6 ) 进 行离散求解，试探函数空间采用一次元空间，检验函数空间采用分片常函数空 间给出了两个不同模型的计算来验证方法的有效性第三章采用迎风格式和 中心有限差分法对方程组( 1 o 5 ) ( 1 0 6 ) 进行求解给出了两种不同的模型计算， 每个模型给出两个数值实验，即初值为光滑和初值为间断的情况来验证方法 的有效性第四部分采用广义迎风差分方法对方程组( 1 0 5 ) ( 1 0 6 ) 进行求解， 同样给出了两个不同模型的数值实验来验证方法的有效性第五部分是结论 分析和展望 一4 一 第2 章广义麓分法 第2 章广义差分法 2 1 离散格式 本章我们采用广义差分法对方程组( 1 0 5 ) 和( 1 0 6 ) 进行求解首先对区 f - - j o ，纠做剖分孔，节点为： 0 = x o x l x n 一1 x n 2 l 剖分乃= 厶：厶= 【x i - - l ，兢】，i = 1 ，2 ，) ，单元厶= 【x i 一1 ，z 胡的长度为： 九t2z t z 一1 再做对偶剖分聪，节点为： 0 = x o x u 2 x 3 2 x n 一1 2 x n = l ， = _ 譬：鬈= 【耽一1 2 ，x i + l 2 ，i = 1 ，2 ，一1 ，瑶= 【x o ，x l 2 ，晡= 【x n 一1 2 ，z 】) ，其中，z 一i 2 = ( z t 一1 + z t ) ，i = 1 ，2 ，n ，x i + l 22 ( z t 十 x i + i ) ，i = 1 ，2 ，一1 ，x i + l 2 一z t 一1 2 = 女( + h i + i ) 现选取试探函数空间既为相应于剖分孔的一次元空间v h = u h c 1 ( ，) ，u h ( o ) = u h ( 厶) ，在每一个厶上是线性函数，完全由单元端点的值所 确定】r 与z t 相应的基函数为： f 1 一 九( z ) ： 1 一 h 选取检验函数空间为：= 叫h ( z ) = 为常数】- ，所以的基函数可取为： 奶( z ) = 1 ， 0 ， 一5 一 z t 一1 z z t ， z z x i + i ， 其他 0 ，z 塔或z 晴，叫h 在每一个e 上 j = 1 ，2 ，一1 警嚣 譬譬 乒 z z 笫2 章广义筹分法 因为在单元厶上， u = u i - l i l + 乱i 咖= u i - 1 【1 一九f 1x x i 一1 ) 】+ u i 1 一危f 1 ( z t z ) 】 所以，有： 乱：= u i - 1 【1 一 f 1 ( z z i 一1 ) 】7 + u t 【1 一九f 1 ( z t z ) 】7 危f 1 ( u t 一札t 一1 ) ， 即有，v x ( x j + 1 2 ) = h j + l ( v j + l 一) ，v z ( x j 一1 2 ) = h j ( v j 一一1 ) 对方程( 1 0 5 ) 两边分别乘以检验函数空间基函数奶，然后在区间 o ，驯上积 分，转化为在对偶单元譬上积分，可得到下面的离散格式： 厶警d z + 巳+ 1 2 v j + 1 2 - - 巳一- ，2 一，z 2 o ( 2 1 1 ) 对方程( 1 0 6 ) 的各项也分别乘以检验函数空间基函数奶，然后在区 f , - - o ，驯上积分，转化为在对偶单元譬上积分，具体如下： l i 吵她d x = 2 谶詈 2 2 屿一1 2 屯一v j 一1 - ( o r ) z 奶d x = ( 口砂) j v 2 一( 臼矽) j + l 2 ， 厶七l n ( 1 - 纠删z = ( 盯l n ( 1 一o ) b 一1 2 一( 盯l n ( 1 一口) ) j + 1 2 = 一1 2l n ( 1 一巳一1 2 ) 一o j + l 2l n ( 1 一巳+ l 2 ) ， 啬奶如= j + 1 2 + j v 2 2厶 一6 一 2 v 巧 卜z zd 一口v 筇= 一1 2 o r z 第2 章广义差分法 ：竺! 兰竺二! ! 卫丝! 生 2 l 一日 2 其中，对啬积分的处理，是采用先对u 取平均值，然后对p 利用数值积分公式 进行的 综合上面各式可得( 1 0 6 ) 的广义差分离散格式为： 2 ，2 鼍2 ，2 一 + b + 1 4 = 巳一1 2 奶一1 2 一o j + l 2 3 j + l 2 + 乃一1 2l n ( 1 一巳一1 2 ) 一o j - l - 1 2i n ( 1 0 j + 1 2 ) ( v j + l 2 + 啪2 ) 盏 ( 2 1 2 ) 对( 2 1 1 ) 式时间导数项采用向前差分离散，可得方程( 1 0 5 ) 的离散格式为： 矿1 = 哆一瓦a t 【n + 1 2 v j “+ l 2 0 ，v 一1 2 孕1 2 】 ( 2 1 3 ) 由格式( 2 。1 2 ) 和( 2 1 3 ) 9 1 1 ，广义差分法的离散格式中出现t o 和 在对偶 端点+ 1 2 和巧一1 2 处的值，而我们所需要求的是p 和u 在节点处的值，所以需要 对巳+ 1 2 和巳一l 2 及+ 1 2 和v j 一1 2 进行处理本文中，对易+ 1 2 和巳一1 2 采用线性 插值，即如一l 2 = ；( 如一1 + 易) ，0 j + 1 2 = ( 易+ o j + 1 ) ，刘；v j + l 2 * n v j 一1 2 采用二次 差值，l i p v j l 2 = v j l - t - ；吻一 v j + 1 ，+ i 2 = 一 一1 - i - i + i + 1 ，并设剖 分磊为等距剖分，则有下面的离散格式： 够+ 1 = 哆一瓮 ( 哆+ 咯。) ( 一孕。+ 6 哆+ 3 蠕。) 一( 呸l + 0 2 ) ( 3 v ) 生1 + 6 哆一蠕1 ) 】， ( 2 坞叫2 h 。2 上q o j o j ) v j - 1 - - ( 2 坞- - 1 1 2 - 1 - 2 坞+ t 2 + 差 2 兰) + ( 2 坞- - 1 2 一百h 2 可妒j o j ) + l = ( 2 1 4 ) 笔 ，b j - - 1 2 ( 巳一+ 巳) 一奶+ - 2 ( 岛+ 岛+ t ) 一7 一 第2 章广义差分法 + 2 0 一1 2l n ( 1 一专岛一1 一吉易) 一2 c r j + l 2l n ( 1 一言易一吉如+ 1 ) ( 2 1 5 ) 格式( 2 1 4 ) 的稳定性限制条件为：a t = m i n 鲁，：汪h 一】- ，u m 口z = m a x l v ( x ，t ) 1 因为0 0 1 ，所以易知上式中v 的系数矩阵是三对角阵，且是对角占优 的，那么椭圆型方程的解u ? 存在且唯一，而求解出的u ? 的值代入( 2 1 4 ) 式即可 求出鳄+ 1 的值，+ 1 的值代入( 2 1 5 ) 式又可以求出 ? “的值，如此循环就可以把 未知量口和v 在各节点处的值全部求出下面，我们将基于上面给出的离散格式。 给出不同的模型计算，并用数值实验来验证方法的有效性 模型一 2 2 方法验证 此模型中取m = o 5 ，t r - c l ，0 - c 2 ，妒= ，其中c 1 和c 2 为正的常数，e 为小的正 常数 由格式( 2 1 4 ) 和( 2 1 5 ) 知，此模型的离散格式为： 哆+ 1 = 0 2 一i a 6 t o n + 1 ) ( 一v ，n 一1 + 6 哆+ 3 咯1 ) 一( ，7 1 一1 + 哆) ( 3 哆l 1 + 6 哆一哆1 1 ) 】， ( 2 2 1 ) ”等e 尚啪邛申尚啪卜等尚 = 缸呸t 嘞) q - 2 c 2 i n 薯等】 ( 2 2 2 ) 因为o 50 6 2 5 x 。6 0 5 3 5 用广义差分法求解的t = 2 0 ，n = 8 0 时解的图像如下图2 2 ： 图2 2 利用广义筹分法得到的初值为间断情况的解的变化图像，其中黑点值表示初值， 红圆圈表示最末值，实线表示初末值之间的解，这里取砂= 盯= 1 0 ，= 0 ，t = 2 0 n = 8 0 这种取值下解口在中间部分衰减由上图可知，目在间断处产生了震荡由 此可见，对于初值为间断的情况，用广义差分法求解消除不了解在间断处的震 荡 模型二 此模型中取m = o 5 e 口，妒= e 口，妒- - c l ，a - - c 2 ，其中c l 和c 2 为正的常数由格 第2 章广义莠分法 式( 2 1 4 ) 和( 2 1 5 ) 知，此模型的离散格式为： 哆+ 1 = 够一j a 6 九t 0 n + 曙。) ( 一喙1 + 6 哆+ 3 咯1 ) 一( 。n 一1 + 够) ( 3 哆l 1 + 6 哆一嵋1 ) 1 ， ( 2 2 5 ) ( e 。一等高蜘七：+ 差九z 离) 哆 + e o j + 1 1 2 _ 等急 ：啪h i ( 咯j 一纠+ c 2i n 警】 ( 2 2 6 ) = 一c 1 ( 咯j 一喙 ) + 蠢】 ( 2 2 6 ) i 司模型1 ，易知上式中，u 的系数矩阵仍是三对角阵，且是对角占优的下面 的计算对如一1 2 和易+ 1 2 采用线性插值 数值实验 对于方程组( 1 0 5 ) ( 1 0 6 ) ，取c l = 1 8 ，c 2 = 1 0 ，初值为0 0 ( z ) = + ；c o s ( 2 7 i z ) ，边界为( 2 2 4 ) 式我们仍用l l ，l 2 及l o 。三种范数下的误差及如前 定义的比率来说明方程组中臼和v 的精度取空间剖分节点n = 1 2 8 0 和时间剖 分节点n = 蔫了r t 五时的解作为真解以，结果如表3 和表4 所示 表3 和表4 分别给出t = 0 0 5 时刻p 和v 的精度由下表可以看出，随着的成 倍增大，三种范数均在相应减小而所定义的范数的比率整体在增大，最后达到 数字3 ，表明广义差分法解模型2 方程组也是可行的，且使得0 和口达到了一阶精 度 图2 3 给出的是t = o 0 5 初值为0 0 ( z ) = 丢+ c o s ( 2 7 r x ) 时，用广义差分法得 到的l g ( 1 | u h 一( u ) l l o o ) 随l g ( h ) 的变化图像，其中左图中代表0 ，右图中v 代表v 由下图可知两条直线的斜率都近似为一，说明此方法使得口和v 均达到了一阶 精度 第2 章，一义差分法 表3 初值为秽o ( z ) = j 1 + 百1c o s ( 2 7 r z ) ，t = o 0 5 时亥t方程组中0 的精度 n l o 。范数 l 比率 l 2 范数l 2 比率l 1 范数l 1 比率 1 09 6 2 6 2 e 0 0 45 9 4 1 3 e 0 0 45 0 2 6 2 e 一0 0 4 2 04 9 7 8 2 e 一0 0 41 9 32 9 4 6 9 e 一0 0 42 0 22 5 6 8 6 e 一0 0 41 9 6 4 02 5 2 2le 一0 0 41 9 71 4 6 4 0 e 0 0 42 o l1 2 8 2 0 e 一0 0 42 0 0 8 01 2 3 7 7 e 0 0 42 0 47 1 3 6 1 e 0 0 52 0 56 2 519 e 0 0 52 0 5 1 6 05 81 7 6 e 0 0 52 1 33 3 4 4 2 e 0 0 52 1 32 9 2 9 8 e 一0 0 52 1 3 3 2 02 5 0l8 e 一0 0 52 3 31 4 3 6 4 e 0 0 52 3 31 2 5 8 4 e 一0 0 52 3 3 6 4 08 3 5 3 4 e 一0 0 6 3 0 04 7 9 3 3 e 0 0 63 o o4 19 9 2 e 0 0 63 0 0 表4 v j f f i y g o o ( z ) = 互1 + 吾1c o s ( 2 7 r x ) ，t = o 0 5 时亥t方程组中v 的精度 n l o o 范数l o o 比率l 2 范数l 2 比率l 1 范数 l 1 比率 1 0 0 0 0 7 30 0 0 6 00 0 0 5 9 2 00 0 0 3 81 9 20 0 0 3 02 0 00 0 0 2 92 0 3 4 00 0 0 2 0 1 9 0 o 0 0 1 5 2 0 00 0 0 1 4 2 0 7 8 09 7 5 41e 一0 0 42 0 57 3 3 6 l e 0 0 42 0 46 8 7 7 8 e 一0 0 42 0 4 1 6 04 6 0 8 1e 一0 0 4 2 1 23 4 3 6 l e 一0 0 42 1 43 2 0 9 3 e 一0 0 42 1 4 3 2 01 9 8 6 4 e 0 0 4 2 3 2 1 4 7 5 3 e 。0 0 42 3 3 1 3 7 5 4 e 0 0 4 2 3 3 6 4 06 6 4 0 1e 一0 0 5 2 9 9 4 9 2 2 3 e 一0 0 53 o o 4 5 8 4 7 e 一0 0 6 3 o o 从上面的两个表格可知，比率逐渐增大到数字3 ，表明此方法对0 和u 达到了 一阶精度对于初值为间断情况，类似模型1 ，此方法对模型2 也是无法消除解 在间断处的震荡 由上述三个数值例子结果可知，随着空间和时间步长的减小，解的精度也 得到相应的提高，由比率和l g l g 图像可知，用广义差分法来解决本文中的方 程组是有效的，且对o $ o v 均是一阶精度但是此方法对于初值为间断情况不理 想，无法消除解在间断处的震荡第三章和第四章里给出的方法则可以消除解 在间断处的震荡 第2 章广义差分法 图2 - 3 # j j f 直为o o ( z ) = + c o s ( 2 7 i z ) ，用广义差分法得到的l g ( | | 魄一u i i o o ) 随l g ( 危) 的变化 图像其中，左图中u 代表0 ，右图巾u 代表 2 1 4 第3 章迎风格式和中心有限莠分方法 第3 章迎风格式和中心有限差分方法 3 1 离散格式 如前讨论，一些传统的有限差分方法，女l ：i l a x w e n d r o f f - 方法，c r a n k n i c h o l s o n 格式和b e a m w a r m i n g 方法并不能应用到我们的方程组中的双曲型 方程上，于是我们这里对( 1 0 5 ) 采用迎风格式来离散对区f n q 0 ，l 作等距剖分， 节点为： 0 = x o x l x n 一1 x n 2l 则方程( 1 0 5 ) 离散格式为： 蹩_ 0 恻， ( 3 1 1 ) 鲣岳螋_ 0 订 o u jj 对于方程( 1 0 6 ) ，令甄= 口- 4 - i 木h ，i = 0 ，1 ，n ，h = 守( 此处n = 0 ，b = l ) 定义：x i + 1 2 = x i + 2 ，x i 一1 2 = 一鲁对( 1 0 6 ) 左边在带有半网格步长的标准 节点翰处利用中心差分离散，得： 其中，m i + 1 2 = m ( x i + l 2 ) ，咄= 矽( 既) ，o i = o ( x d 再对一阶导数项和有端项均 用中心差分离散，并用饥代替甜( z t ) 可得下面的离散格式： 2 m i + v i + 1 - 2 ( m i _ 丽+ m i + ) v i + 2 m i _ 一；v i - i 一而q o i o i 2 1 一阮1 = 一业嘎# 虹一盟坠选嚎堕必( 3 2 ) ：= 一i 1lz l 2 九2 九 。、一7 格式( 3 1 1 ) 的稳定性限制条件为：a t = m i n h ，瓦急) ，u m 口z = m o z l ( z ，t ) 1 由式( 3 1 2 ) 化简可以看出， 的系数矩阵是三对角且是对角占优的，因此解u 存 攀 在且唯一，而将解u ? 代入( 3 1 1 ) 式即可求出解p r l ，如此可以将未知量全部求 出下面针对格式( 3 1 1 ) 进行相容性和收敛性分析 设 3 2 相容性分析 假设日托和( p u ) z 霉有界先考f g v o 的情况 甲o ( z i ，t n + 1 ) 一o ( x i ，t n ) 。( o v ) ( x i ，t n ) 一! 翌! ! ! 兰兰! ! 1 2 2 a t a z o ( x l ，t n ) + a t o tx i ，t n ) + t k t 2o t t ( x i ，锯) 一o ( z i ，t n ) 其中，t f ( k ，t 礼+ 1 ) ，z c x i 一1 ，x i ) 显然，由假设，n a t _ 0 ，z _ o 时， 有t _ 0 同理，当v 0 时，有： t = 缸+ 等k n a t _ 0 ，z 一0 时，有t _ 0 故差分格式( 3 1 1 ) 是相容的 3 3 收敛性分析 令t 是时间步长， 是空间步长，g 孑= 哼一v ( x j ，t n ) 为已知量，e 芗= 凹一 o ( z j ，t n ) ，下面考虑格式( 3 1 1 ) 的收敛性问题 首先来看u o 时的情况设p ( z ，亡) 是初值问题的解，留是差分格式的解 令t ( ，t n ) 表示差分格式在点( ，t n ) 处的截断误差，则有： 丁 ) ：堕等掣+ 型坐乓掣 第3 幸迎风格式和中心有限芽分方法 此式可改写成： o ( x j ，t n + 1 ) = p ( 巧，t 亿) 一p 【( p u ) ( z j ，t n ) 一( 口 ) ( z j l ，t n ) 】+ v t ( x j ，t n ) ，( 3 3 1 ) 其中，p = 瓮= 丢差分格式写成： 哆+ 1 = 町一川( 乳) 孑一( p u ) ，】 由( 3 3 2 ) 式减去( 3 3 1 ) 式，可得 e 了+ 1 = 髟n p 【哆e 孑一弓l 1 e 务1 】一t t ( z j ，t n ) = ( 1 一p 哆) 亏+ p 孕l e ，n 一1 一t t ( x j ，t n ) ( 3 3 2 ) 若假设o p 哆 1 ，则上式等号右边第一项和第二项系数均为非负，由此可得： e 罗+ 1 i ( 1 一肛e 梦) i e 引+ l p 孕1 i i e 务1 i + t i t ( x j ，t ) 1 假设o ( z ，t ) 为原问题的解，由截断误差可知， t ( x j ，t n ) i m ( 7 + 危) ， 其中，m 为一正常数再令既= s u pi e 芗i ，所以 j 从而有 e j n + 1 i ( 1 一p 孑+ p 孕1 ) l e v i + m r ( t 十九) 磊+ m t ( t + 允) 由此不等式递推得 r + 1 r + m 7 - ( 7 + ) b e o + m n 7 ( 7 - + 九) 注意到，在初始时间层t o 上，有 鳄= p ( ，o ) = o o ( x a 1 7 第3 章迎风格式和中心有限茇分方法 所以有e ? = 0 ，因此， e o = s u pl e o l = 0 j 故有， b m n 7 - ( 丁+ 危) 而初值问题中n r t ，所以， 风m t ( t + h ) 令7 ，h _ 0 ，有风_ 0 即0 2 _ 口( ，k ) v 0 的情况证明与此类似 综上所述，迎风离散格式( 3 1 1 ) 是收敛的 模型一 3 4 方法验证 此模型g , y , m = 0 5 ，妒= c 1 ，盯= c 2 ，妒= ，其中c 1 和c 2 为正的常数，为小的正 常数 由格式( 3 1 1 ) 和( 3 1 2 ) 知，此模型的离散格式为： 嵋o ( 3 删 嵋 0 ， u ：1 一( 2 + - 2 砸0 7 ) 嵋+ 苒1 = h zi c l ( 曙l 一+ c 2 l ni 1 - 一仃o ? + 1 1 ( 3 蚴 因为0 0 1 ，所以易知上式中 的系数矩阵是三对角阵，且是对角占优 的，那么椭圆型方程的解 存在且唯一下面我们将通过数值实验来说明此方 法的有效性 数值实验一 哆嵋移一钟 ，卜 一 日 l 一 竹冲 一。 扣 孔竹件 咿p 笪畿瓦 一 一 叩卵 = = + + 野盱 ，。j、 第3 章迎风格式和中心有限茏分方法 组 对于方程组( 1 0 5 ) 和( 1 0 6 ) ，取e = 0 5 ，c 1 = 1 8 ，c 2 = 1 0 ，则有下面的方程 边界为： lo t + ( p ) z = 0 ， k 一篙一1 8 ( 1 n ( 1 叫k 。a 3 o ( o ，t ) = o ( 1 ，t ) = 0 ，v ( o ，t ) = v ( 1 ，t ) = 0 ( 3 4 4 ) 我们取初值为：o o ( z ) = + ；c o s ( 2 1 r z ) ，取空间剖分节点n = 1 2 8 0 5 1 时间剖分节 点n = = 了# 1 j 时的解作为真解以，且用第二章中定义的比率来验证o g l v 的 精度结果如表5 和表6 所示 表5 和表6 给出的分别是在t = 0 5 时n o 和u 的误差和比率由下表可以看 出，随着的成倍增大，三种范数均在相应减小而所定义的范数的比率在增大， 最后达到数字3 ，表明迎风格式和中心差分法解模型l 方程组使得0 和u 达到了 一阶精度 图3 1 给出的是t = 0 5 ，s = 0 5 时，用迎风格式和中心差分法得到 的l g ( 1 l 巩一以i i o o ) 随l g ( h ) 的变化图像，其中左图中u 代表0 ，右图中己，代表 两 条直线的斜率都近似为一，说明此方法使得0 和v 均达到了一阶精度 表5 初值为o o ( x = ；+ ic o s ( 2 7 r x ) ，t = 0 5 ，= o 5 时万栏组中曰脏 精度 n l 。o 范数l o o 比率 l 2 范数l 2 比率l 1 范数l 1 比率 1 00 0 0 2 8o o o l 90 0 0 1 8 2 0o o o l 51 8 7o 0 0 1 01 9 09 3 9 2 8 e 0 0 41 9 2 4 07 7 4 7 3 e 一0 0 41 9 45 0 5 3 8 e 0 0 41 9 84 6 0 7 8 e 0 0 42 0 4 8 0 3 8 2 6 6 e 一0 0 4 2 0 2 2 4 5 6 5 e 一0 0 42 0 62 213 9 e 0 0 42 0 8 1 6 01 8 0 l o e 一0 0 42 1 21 1 4 8 l e 0 0 42 1 4 1 0 2 9 5 e 0 0 4 2 1 5 3 2 07 7 5 4 8 e 一0 0 52 3 24 9 2 3 5 e 。0 0 52 3 34 4 0 3 2 e 0 0 52 3 3 6 4 02 5 9 0 7 e 一0 0 5 2 9 91 6 4 1 6 e 一0 0 53 o o1 4 6 5 9 e 一0 0 53 0 0 表6 初值为( z ) = j + c o s ( 2 7 r x ) ，t = 0 5 ，e = o 5 时方程组中u 的精度 第3 章迎风格式和中心有限荠分方法 n l o o 范数l o o 比率 l 2 范数l 2 比率l 1 范数l 1 比率 l o0 0 0 3 20 0 0 2 lo 0 0 1 9 2 0o o o l 71 8 80 0 0 1 02 1 09 3 0 6 7 e 一0 0 42 0 4 4 08 7 5 7 9 e 一0 0 41 9 45 2 1 1 2 e 0 0 41 9 24 5 6 4 6 e 0 0 42 0 4 8 04 3l5 6 e 0 0 42 0 32 5 319 e 0 0 42 0 6 2 2 0 3 4 e 一0 0 42 0 7 16 02 0 3 2 5 e 0 0 42 1 21 1 8 3 9 e 一0 0 42 1 41 0 2 6 9 e 0 0 42 1 5 3 2 08 7 5 1 1 e 一0 0 5 2 3 25 0 7 8 6 e 一0 0 5 2 3 34 3 9 7 7 e 0 0 52 3 4 6 4 02 9 2 3 8 e 一0 0 52 9 91 6 9 3 7 e 一0 0 53 0 01 4 6 5 3 e 0 0 53 o o 从上面的两个表格可以看出，比率逐渐增大到数字3 ，表明迎风格式和中 心差分法解方程组使得0 和u 达到了一阶精度 图3 1t = 0 5 ，e = 0 5 ，初值为p o ) = + e o s ( 2 7 r x ) ，用迎风格式和中心有限差分格式 得n l 拘l g ( 1 t u n u , i i 。o ) 随l g ( 危) 的变化图像其中，左图中u 代表0 ，右图中u 代表u 从上述两个例子可以看出，随着空间和时间步长的减小，解的精度得到 相应的提高当空间步长变为原来的 时，最大绝对误差，l 2 范数及l 1 范数 都相应地变为原来的 这说明，迎风格式和中心有限差分求解本文的方程 组( 1 0 5 1 0 6 ) 是有效的，且由比率和l g l g 图像知此方法使得0 和v 均达到一阶 精度下面看初值为间断的情况 数值实验二 第3 章迎风格式和中心有限差分方法 对于方程组( 1 0 5 ) 和( 1 0 6 ) ，取e = 0 0 5 ，c l = 1 0 ，c 2 = 1 0 初值为： c z ，= 50 6 2 5 x 。6 0 5 3 5 - c 3 4 5 ， 用迎风格式和中心有限差分法求解的t = 2 0 ，n = 8 0 时解的图像如下图3 2 ： 这种取值下解口在中间部分衰减由下图可知，臼在间断处没有震荡 图3 2 利用迎风格式和中心差分法得到的初值为间断情况的解的变化图像，其中黑点值 表示初值，红圆圈表示最末值，实线表示初末值之间的解，这里取砂= 盯= 1 0 ，= 0 0 5 ，t = 2 0 ，n = 8 0 由上述两个数值实验可知，迎风格式和中心有限差分方法解模型l 是可行 的，且使得日和v 均达到一阶精度，与理论分析一致，而对于初值为间断的情况， 此方法消除了解在间断处的震荡下面看另一模型 模型二 此模型中i 汉m = o 5 e 口，妒= ，妒= c 1 ，o - - c 2 ，其中c 1 和c 2 为正的常数由格 式( 3 1 1 ) 和( 3 1 2 ) 知，此模型的离散格式为： 叩0 ， ( 3 4 6 ) 诈 0 ，；o 嵋朋卵 卜 一 目 1 一 n 冲 一。 孔n 件 p 归 坐缸瓦 一 一 卵印 = = 盱卵 ，-，、【 第3 章迎风格式和中心有限差分方法 跏- ( e 咏h 2 十九2 啬) 嵋 1 - e o i + x 2 v 拜， = 一互h 【c 1 ( 目耳1 一p 。1 ) + c 2 l n i 1 一- 0 二h i l ( 3 4 7 ) 同样由于0 0 1 ，易知上式中，u 的系数矩阵是三对角阵，且是对 角占优的，所以椭圆型方程的解 ? 存在且唯一但同样地出现了p 在对偶端 点戤一1 2 和轨+ 1 2 处的值，故仍需对吼一1 2 和吼+ 1 2 进行处理本文中，仍采用对其 进行线性插值的方法来处理的下面我们将通过两个数值实验来说明此方法 对模型二也是可行的 数值实验一 我们取c 1 = 2 0 ，c 2 = 1 0 ，初值为如( z ) = + s i n ( 2 1 r x ) ，边界为( 3 4 4 ) 式我