（计算数学专业论文）同时满足kdv方程两个守恒律的数值算法.pdf

摘要 k o r t e w e g - d ev r i e s ( k d v ) 方程是人们在研究一些物理问题时得到的非线性波 动方程，其解满足无穷多个守恒律。本文为该方程设计了一种守恒型差分格式， 其采用的是有限体积法，它是g o d u n o v 型的但与传统的守恒型差分格式不同的 是，它的数值解同时满足k d v 方程前两个相关的守恒律这样可以更好地保持解 的物理上的守恒性质，从而格式本身有更好的稳定性数值例子表明这一算法是 有效的 关键词：k d v 方程，守恒律，网格平均，函数重构 a b s t r a c t k o r t e w e g - d ev r i e s ( k d v ) e q u a t i o ni s an o n l i n e a rw a v ee q u a t i o na r i s i n gf r o m p h y s i c ss t u d y ，w h o s es o l u t i o n ss a t i s f i ei n f i n i t e l ym a n yc o n s e r v a t i o nl a w s t h i sp a p e r d e s i g n sad i f f e r e n c es c h e m ef o rt h ee q u a t i o n w h i c hi s o ft h ef i n i t e - v o l u m et y p e ，i t a c o d u n o vt y p e d i f f e r e n tf r o mt h et r a d i t i o n a lf i n i t e - v o l u m es c h e m e sh o w e v e r j t h e s c h e m es a t i s f i e st w or e l a t e dc o n s e r v a t i o nl a w so ft h ek d ve q u a t i o nd e s i g n e di n s u c haw 日孔t h es c h e m ep r e s e r v e sb e t t e rt h ep h y s i c a lc o n s e r v a t i o np r o p e r t i e so ft h e k d ve a u a t i o a n dt h e ni th a sb e t t e rs t a b i l i t yt h en u m e r i c a le x a m p l e ss h o wt h e e f f i c i e n c yo ft h es c h e m e k e y w o r d s ：k d ve q u a t i o n ，c o n s e r v a t i o nl a w s ，c e l l a v e r a g e ，r e c o n s t r u c t i o n i i 原创性声明 本人声明：所呈交的论文是本人在导师指导下进行的研究工作除了文中特别 加以标注和致谢的地方外，论文中不包含其他人已发表或撰写过的研究成果参与 同工作的其他同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示 了谢意 弥苹粉啉m 7 7 本论文使用授权说明 本人完全了解上海大学有关保留、使用学位论文的规定，即：学校有权保留论 文及送交论文复印件，允许论文被查阅和借阅；学校可以公布论文的全部或部分内 容 ( 保密的论文在解密后应遵守此规定) 酶榔导师始撇铱嘲妒5 。厶f 7 2 0 0 5 上海大学硕圭学位论文 第一章引言 孤立子现象是在1 8 3 4 年由英国数学家s c o t tr u s s e l 首先观察到的，并于1 8 4 4 年在英国科学促进协会第十四届会议报告，发表了”论波动”一文，详细阐述了 他所观察到的孤立波现象，但他始终没有给出令人信服的论断，没有从流体力方 程中推得孤波的存在性直到1 8 9 5 年，荷兰著名数学家k o r t e w e g 和他的学生d e v r i e s 在对孤波进行全面分析后指出这种波可近似为小振幅得长波，并建立了非线 性浅水波运动方程，这就是著名的k o r t e w e g - d ev r i e s ( k d v ) 方程，并从中求出了 孤立子解，从理论上证明了孤立子的存在，孤立子也称孤立波 但由于工具的限制孤立子的研究一直处于停顿状态，许多孤立子的基本问题 得不到满意的答案直到6 0 年代随着计算机技术的发展，孤立子的研究才出现了 进展，尤其是1 9 6 5 年，美国物理学家k r u s k a l 和z a b u s k y 通过数值计算详细研究 了k d v 方程两波相互作用的全过程，他们对作用前后所得的数据进行分析后发现 孤波的形状和速度保持不变而具有弹性散射的性质，从而开辟了孤立子研究的新 时代 具有如下形式的偏微分方程； 札t + 6 u u z + u z 。z = 0 称为k o r t e w e g - d ev r i e s ( k d v ) 方程( 见【2 8 ， 2 9 ) ，是人们在研究物理问题时得到 的非线性波动方程，其中“是关于z 和t 的函数，是为方程的解，毗是“对时 间t 的导数，“。是u 对空间x 的一阶导数，“是“对。的三阶导数 k d v 方程具有无穷多个守恒律，孤立子碰撞以后形状与波速保持不变，就是 无穷多个守恒律的体现其前两个守恒律为： “+ ( 3 u 2 + “。k = 0 ， ( 1 2 ) ( 2 ) + ( 4 u 3 一“：+ 2 u u 。k = 0 ， ( 1 3 ) 2 0 0 5 上海大学硕士学位论文2 第一个守恒律表示动量守恒，第二个守恒律表示能量守恒 记能量u 2 为u ( u ) ，n s 寸i g 量和能量流函数分别为： ，( “。) = 3 u 2 + “。 和 f ( u ，札。，“。) = 4 u 3 一“：+ 2 u u 。 在以下的讨论中我们 ( 14 ) ( 15 ) k d v 方程是一类非线性波动方程，其解满足无穷多个守恒律近年来这类 非线性波动方程的数值解法也迅猛发展起来，目前所应用的数值方法大致可以归 纳为三大类第一类方法是差分方法， k r u s k a l 和z a b u s k y 就是应用二阶精度的 l e a p - f r o g 格式计算k d v 方程的( 见【2 7 ) ，以后，又分别提出了二次守恒格式， h o p s c o t c h 格式和耗散格式( 见f 3 3 j ，f 3 j ，f 2 5 j ) ，最近叉发展起来很多新静差分算法 ( 见吼 2 1 】， 1 a 1 ) 第二类方法是有限元方法，w a h l b i n 提出了计算方程的耗散有 限元格式( 见【2 6 1 ) 第三类方法是谱方法和拟谱方法( 见 3 4 ，【1 】) 高精度的差分格式不一定给出好的数值解，所以人们经常从物理定律出发 构造合理的差分格式，使其尽可能的保持原问题的物理性质对于k o r t e w e g - d e v r i e s ( k d v ) 方程，其解满足无穷多个守恒律，孤子运动的稳定性与守恒律有关， 所以，人们希望设计一种差分格式，使得它的解满足相应的离散守恒律然而， 一般很难模拟全部守恒律，故只能选择其中若干个 在本文中设计了一种新的差分格式，这种差分格式同时模拟了k d v 方程的 前两个相关的守恒律，这样设计的格式能更好地保持解的物理守恒性质，从而格 式本身有更好的稳定性其思想来源于求解双曲守恒型方程时所用的守恒型差分 方法在求解守恒型方程时，一个好的格式应该是守恒的，非线性稳定的其中 g o d u n o v 型格式就是一类典型的守恒型格式，我们的格式也属于g o d u n o v 型格 式近年来发展起来的众多高分辨g o d u n o v 型格式都是满足这些条件的( 见 卧【1 7 ，【1 8 ，【2 0 ) 原始的g o d u n o v 格式只有一阶精度，为了提高精度，近年来 采用了高次多项式重构而发展了所谓的g o d u n o v 型格式然而采用二次及以上的 重构以后，格式就会遇到非线性稳定性的问题，为此提出了由此提出了t v d 和 e n o ( 见 6 ， 7 i ，【2 j ) 的概念及想园相应的实施技术，近二十多年来发展起来的所 谓g o d u n o v 型格式就是基于这一思想的 2 0 0 5 上海大学硕士学位论文3 对于单个守恒型方程来说，传统的守恒型格式都只模拟一个守恒律近年来， 李红霞和茅德康设计了一种新的守恒型差分格式，它是g o d u n o v 型的但和传统 格式不同的是，该格式同时模拟了两个相关的守恒律，关于解的和关于熵的。其 做法是格式定义和计算了两个相关的数值量，数值解和数值熵，并用数值熵来帮 助重构数值解，见 1 1 】， 1 2 3 0 】和 3 l 】文特别是当将这一格式用于线性传输 方程计算光滑解时，由于保持了两个守恒律，其有意想不倒的稳定性质对于像 w a v e - p a c k e t 这样高难度的问题，极其长时问的计算效果都非常好，没有任何数值 耗散，见 3 1 恰是这一现象启发我们是否可以将这一设计思想用于设计k d v 方 程的差分格式，使其同时满足多个守恒律的以提高稳定性使得数值解的长时间 的计算效果都非常好 本文将【1 1 ， 1 2 】， 3 0 和 3 1 】文中的设计方法应用于k d v 方程作为一个初 试阶段的工作，我们仅设计了一个同时满足两个相关守恒律的差分格式数值实 验表明，我们所设计的格式有非常好的稳定性，可以在长的时间内保持数值解的 物理我们的这种设计方法可以推广到同时满足多个守恒律的情况，这将是我们 以后的工作 本文的具体结构如下，第一章为引言；第二章作为后面的准备，回顾了有限 体积法和g o d u ! l o v 型格式的一些基本概念和理论；第三章对算法进行丁详细地描 述，第一节给出了算法的基本思想，第二节给出了算法的详细构造；第四章给出 了一些数值算例来检验我们的算法；第五章是结论 2 0 0 5 上海大学硕士学位论文4 第二章基本的概念和理论 我们所设计的算法，其思想来源于双益型守恒型方程的差分格式，作为对以 后几章的准备，在本章我们给出了求解双曲型守恒型方程( 组) 有限体积法中的一 些基本概念和理论，并介绍了g o d u n o v 型方法 2 1 守恒型方程的有限体积法 具有如下形式的一维空间的偏微分方程 ( 2 1 1 称为双曲型守恒型方程，简称守恒型方程，( 2 ，1 ) 的形式称为守恒形式，其中“是 关于x 和t 的函数，是方程( 2 1 ) 的解，f ( u ) 称为流函数。显然k d v 方程也是守 恒型方程，本章主要考虑对这种守恒型方程的数值模拟该守恒型方程是由物理 定律在任意两点x z 和z 2 之间的如下积分形式 一r f o 云“( z ，t ) d x = ，( “( t ) ) 一，( u ( t ) ) ( 2 2 ) u 1 得到的( 2 2 ) 表示在区间k ，z 2 中总的流体的量( 如动量，能量等) 的变化仅仅 与两端点处的流有关，这就是守恒的基础，其中，( n ( z t ，) ) 和，( n ( 逝，t ) ) 分别表 示在z l 和z 2 点的流人流出的量 首先对( z ，t ) 平面进行网格剖分，即在( x ，t ) 平面内作两组平行线，由这两组 乎行线构成的长方形网格覆盖了整个( z ，t ) 平面，其中网格线的交点称为网格节 点，见图2 1 在空间方向上，第j 个网格记为( q a 2 ，x j + l 2 ) ，分别以x j _ l 2 和 巧+ 1 2 为左、右网格边界，网格中心为q ，网格长度记为b = x j + u 2 一q 一1 2 在 本文中我们采用的是正规网格划分，即网格长度就是空间步长，h = h j = a s ， 则x j l 2 = q 士；= 0 土；) ，时间步长r = k + l t 2 0 0 5 上海大学硕士学位论文5 圈2l 由网格边界的流来更新网格平均畸的有限体积法，显示的是( z ，) 平面 对于一般的差分方法来说，所求得的数值解畸是对原方程在节点( 巧，t 。) 处 的精确解u ( q ，k ) 的逼近而在有限体积法中，数值解叼是在时刻t 。对精确解 在阿格( 有限体积) ( q 1 ，2 ，q + 1 2 ) 上的积分的平均( 网格平均) 的逼近，即 哼= m j 。心，缸 ( 2 3 ) 采用网格平均的意义是很重要的，不但可以很容易地运用守恒型方程( 组) 的一些 重要的性质，并可保证得到的数值方法是守恒的 由守恒律的积分形式( 2 2 ) 得 爰e 2 州- ，( “( z j - u 2 , t ) ) _ m ( 。，m ( 2 4 ) 将( 2 4 ) 从t 。到t 。+ 1 积分可得 ；z ：2 “( 置如上i ) 如= ；z 二二2 “( 薯u 如 一；”1 ，锄+ l 。泓弦一f ”1 ，( “( 勺1 1 。) ，蜘叼 ( 2 5 j 一般来说我们不能精确地对( 2 5 ) 式的右端的时间积分进行赋值，因为u ( 。2 1t ) 在每一网格边界匕随着时间变化侣我们可研究具有以下形式的数值方法 矿1 = u ? 一a ( 魂 一茁 ) ， ( 2 6 ) 其中尊称为数值流函数，是在网格边界+ z 2 处对精确解的流函数，( u ) 的网 格边界平均的逼近 魂 = m 删，“) = ；1 ，( u ( ，啪d t ， ( 2 7 ) 2 0 0 5 上海大学硕士学位论文 6 而a = r 是网格步长比。 ( 2 6 ) 称为守恒型格式，因为它模拟了精确解的性 质如果我们在所有网格上求 u r l 的加和，可得 nn “= 哆一r ( 扁+ l 。一髓啪) ， ( 2 8 ) t = ji = j 由此可知，除了在边界处的流外在整个区域上我们得到了精确的守恒性 2 2 g o d u n o v 格式及g o d u n o v 型格式 在求解方程( 2 1 ) 时，个好的格式应该是守恒的、非线性稳定的近年来发 展起来的众多的g o d u n o v 型格式( 见【2 2 】【2 3 ， 2 4 ) 就是这样的格式。g o d u n o v 格式可称为r e a ( r e c o n s t r u c t i o n - e v o l u t i o n - a v e r a g i n g ) 算法下面我们来简述这一 格式 数值解如( 2 3 ) 式所示，是对网格平均的逼近假定已求得t = 。时的数值解 n ? ，下面求解t = t 。+ ，时的数值解u n ，“g o d u n o v 型格式是按照重构( r ) ，发 展( e ) 和求网格平均( a ) 三个步骤来进行 重构：在如层上由网格平均u7 ，1 , 对解进行重构，得到分片多项式函数月如；u “) ， 可根据精度要求采用分片常数函数，分片线性函数或者更高阶的重构最简单的 重构是如下形式的分片常数函数重构 r ( z ；u ”) = 哆z ( 一j ，吩+ ) ，j = 1 2 - ( 29 ) 发展：考虑方程( 2 1 ) 的初值问题 ! 也t + ，( 也) z = o n 2 n + 1( 2 1 0 ) f 砬( z ，t 。) = r ；) ， 显然也( z ，t ) 是由一系列以( q + ；1 ，t n ) 为中心的局部r i e m a n n 问题组成利用c f l 条件的限制，我们可以保证在t = t l 时，相邻的两个局部r i e m a n n 问题所产生 的波还没有相互作用，这样就可翰由这些局部r i e m a n n 问题的解求得此初值问题 在，t n + 1 上的精确解 2 0 0 5 上海大学硕士学位论文7 求网格平均：t = k + - 时的数值解由铲( 。，佗+ - ) 的网格平均给出： 嵋“= ；r 5 照x , t n + l 廊 ( 2 1 1 ) 这样就求得了 哼“ ，然后重复上述过程求得( 哼帕 、 哼“，、 在实际计算中，因( 2 1 1 ) 式可由守恒律的积分形式算出而使上述算法得到简 化，即 矿17 哼一a ( 嚣 一孕 ) ， ( 2 1 2 ) 其中，翌l 如( 27 ) 式所示 由于求解r i e m a n n 问题是一件很费时的工作，因此在实际计算中人们往往用 近似方法求解r i e m a n n 问题，最常用的是r o e 的近似方法( 觅 15 f 1 6 】) ，故数值 流函数只能是对( 2 1 2 ) 中的流函数的逼近，记对网格边界上流通量的一个逼近 曩 = ；e ”懈( ) ) o ( h q ( 2 1 3 ) c o d u n o v 型格式中的一个重要的步骤是在各时间层上构造一分片多项式函数 r ( z ；矿) ，它在每一网格上就等于该网格上的网格平均，即由网格平均重新构造该 时间层上的数值解函数，以此作为初值来求下一时间层上的数值解，这一步骤称 为重构，很显然这样的重构满足 ；f 翟州) 虻哼， ( 2 a ) 称为重构保持守恒性。 显然( 2 9 1 式的分片常数重构是不够理想的，因为它使得格式只有一阶精度 为了提高精度，可做分片多项式重构，最简单的为分片线性重构 r ( z ；矿) = 哼+ 譬扛一巧) 2 ( 吩一；，+ j ) ，j = i ，2 ， ( 2 1 5 ) 这儿斜率彤n 如满足 s j n = u z ( x j ，t n ) + 0 ( ) ， ( 2 1 6 ) 贝重构函数r ( 矿) 就是对“t 。) 二阶精度的逼近斜率哼可由数值解的差商 来算得，这样就构造出了l a x - w e n d r o f f 和l a x - w e n d r o f f 类型的格式然而，众所 2 0 0 5 上海大学硕士学位论文 8 周知，这类格式有非线性不稳定问题，数值解在间断附近会有非物理振荡，由此， 诸如t v d 和e n o 等技巧就发展起来以克服这类非线性不稳定性，如见f 1 4 和 1 9 所有这些格式都只模拟了一个守恒律 近年来李红霞和茅德康在 1 1 ， 1 2 l 【3 0 】和 3 1 】文中发展了同时模拟了两个 相关的守恒律有限体积差分格式。下面我们来简单地介绍以下这个格式众所周 知，对任意一个关于“的非线性函数u ( ) ，对( 2 1 ) 方程的光滑解也成立 u ( “) t + f m ) 。= 0 ，( 2 1 7 ) 这儿 f ( u ) = ，7 ( u ) u 7 ( u ) ( 21 8 ) 如果u ( u ) 是“的凸函数，那麽它就是一个熵 基于这样一个事实，李和茅的格式中除了计算( 2 3 ) 意义的数值解以外，还计 算一个数值熵叼，它是对解的网格平均的逼近 叼= ；r 瓤蝴圳也 ( 2 1 9 ) 此时重构是分片线性的，如( 2 ，5 ) ，其中斜率s n ，由下面方程解出， 1 f 夸z j + 女6 u ( r ( 薯= ；鬈? u ( 哼+ 8 3 ( 。一巧) ) d z 2 叼； ( 22 。) 即重构函数的熵的网格平均应等于该网格上的数值熵可以证明如此解得的毋满 足( 2 1 6 ) 在重构了数值解之后，下一时间层上的数值解还是按发展和求网格平 均步骤来进行，由此两个根据( 2 1 7 ) 可导出叼“的守恒的计算格式如此设计 的格式，由于同时满足两个相关的守恒律，因此计算光滑解时非常稳定 2 0 0 5 上海大学硕= 学位论文 9 第三章算法的描述 本章的第一节将对算法的基本思想进行简单的描述。第二节将对算法进行详 细描述，读者将会看到本算法和传统的守恒型方法在形式上有很大的差异，特别 是本算法涉及了第二个守恒律，同时模拟了k d v 方程的前两个守恒律，这样可以 更好的保持解的性质 3 1 算法的基本思想 本文对于k o r t e w e g - d e 、f r i e s ( k d v ) 方程所设计守恒型差分格式，沿用丁李红 霞和茅德康在【1 l 】， 1 2 1 ，【3 0 】和 3 1 文中的设计思想，同时模拟了k d v 方程的前 两个相关的守恒律。然而，李和茅的格式是为双曲守恒型方程所设计的，而k d v 方程和双曲守恒型方程有很大的差别，因此在贯彻这一思想是会遇到不同的问题 和用到不同的解决手段。为此我们有必要在本章的开始来介绍一些基本的思想 1 如引言中所述，k d v 方程的解满足两个守恒律( 12 ) 和( 1 3 ) 。但和双曲 守恒型方程不同的是其对应的两个流函数( 14 ) 和( 1 5 ) 不仅和有关，而且还和 它的一阶和二阶导数n 。和u 。有关。下面的讨论中会看到，这给格式的设计带来 一定的复杂性。 2 ) 我们所设计的格式还是属于g o d u n o v 型的，因此其数值解仍是如同( 2 3 ) 那样的阿格平均的逼近。如同李和茅的格式那样，我们的格式还计算了数值动量 u ( u ) ，其是对动量u ( u ) = “2 的网格平均的逼近，如( 21 9 ) 。格式仍然按前述的重 构，发展和求网格平均来进行。 3 ) 在重构步骤中，我们须在每个网格中重构数值解。然而( 21 5 ) 那样的线性 重构是不够的，这是因为流函数中有解的一阶和二阶导数为了对它们达到二阶 逼近。重构函数在每个网格中必须是如下形式的三阶多项式 r 扛；u “) = n ? o + 哆1 扛一x j ) 十畸1 2 扛一) 2 + 喀3 0 一q ) 3 ， ( 3 1 ) r ( z ：u “) = n 芦+ 哼1 扛一q ) 十畸1 2 如一z ，) 2 + 喀3 0 一q ) 3 ， ( 3 1 ) 2 0 0 5 上海大学硕士学位论文9 第三章算法的描述 本章的第一节将对算法的基本思想进行简单的描述。第二节将对算法进行详 细描述，读者将会看到本算法和传统的守恒型方法在形式上有很大的差异，特别 是本算法涉及了第二个守恒律，同时模拟了k d v 方程的前两个守恒律，这样可以 更好的保持解的性质 3 1 算法的基本思想 本文对于k o r t e w e g - d ev r i e s ( k d v ) 方程所设计守恒型差分格式，沿用了李红 霞和茅德康在【1 1 ， 1 2 l ，【3 0 和 3 l 】文中的设计思想，同时模拟了k d v 方程的前 两个相关的守恒律然而，李和茅的格式是为双曲守恒型方程所设计的，而k d v 方程和双曲守恒型方程有很大的差别，因此在贯彻这一思想是会遇到不同的问题 和用到不同的解决手段为此我们有必要在本章的开始来介绍一些基本的思想 1 ) 如引言中所述，k d v 方程的解满足两个守恒律( 1 2 ) 和( 1 3 ) 。但和双曲 守恒型方程不同的是其对应的两个流函数( 1 4 ) 和( 1 5 ) 不仅和u 有关，而且还和 它的一阶和二阶导数u 。和u 。有关下面的讨论中会看到，这给格式的设计带来 一定的复杂性 2 ) 我们所设计的格式还是属于g o d u n o v 型的，因此其数值解仍是如同( 2 3 ) 那样的阿格平均的逼近如同李和茅的格式那样，我们的格式还计算了数值动量 u ( “) ，其是对动量u ( u ) = “2 的网格平均的逼近，如( 21 9 ) 格式仍然按前述的重 构，发展和求网格平均来进行 3 ) 在重构步骤中，我们须在每个网格中重构数值解然而( 2 1 5 ) 那样的线性 重构是不够的，这是因为流函数中有解的一阶和二阶导数为了对它们达到二阶 逼近，重构函数在每个网格中必须是如下形式的三阶多项式 r ( x ；u ”) = n ；o + 喀1 扛一x a 十畸1 2 扛一q ) 2 + 吩n ”x x j ) 3 ， ( 3 1 ) 2 0 0 5 上海大学硕士学位论文1 0 其中系数n ”哼”，n 。和吗n 1 3 是分别对“( ，。) ，u 。( 町：t n ) u 。( 岛，t ，) 和? 2 x x v ( t 。) 的逼近如同一般的有限体积法一样，这些系数应由解的阿格平均计算得。这可 由原函数重构法来得到，具体可见下面的讨论。 4 ) 如同李和茅的格式那样，我们仍要求( 2 1 9 ) 式成立，但此时重构函数是如 同( 3 1 ) 式那样的为了满足( 21 9 ) 我们让a 驴不定，然后由( 2 1 9 ) 式解出它来 5 ) 如同g o d u n o v 型的格式一样，在发展和求网格平均步骤中解( 11 ) ( 或( 1 2 ) ) 的初值问题，然后计算类似( 2 1 1 ) 那样的积分以计算数值流函数和数值动量流 然而和双曲守恒型方程不同的是，这儿初值问题的解并非是由一系列的r i e m a m 问题拼接而成，因为k d v 方程本质上并不是一个双曲型方程然而，k d v 方程 的解是光滑的，由此我们能用一定的数值方法计算数值流和数值动量流 3 2 格式的建立 首先对( x , t ) 平面进行网格剖分，其对应于( x t ) 平面内的两组平行线：q = j h ，j = 0 ，1 ，士2 ，和t 。= n r ，n = 0 ，1 ，2 ，其中h 和r 分别为空间步长和时 间步长 我们的格式是g o d u n o v 型的，其数值解u ? 是对精确解的网格平均的近似： 哼= ；r 7 瓤圳虹 ( 3 z ) 其中( 。，t 。) 表示方程的精确解，我们的算法将涉及到第二个守恒量( u ( 。，t ) ) 的 数值逼近数值函数w ，称为数值能量，是对精确解的能量u ( “( z ，t ) ) = u 2 ( 。，t 。) 的网格平均的近似： 叼= ；毋m ) 如- _ 1 f - 。识州圳z ( 3 3 ) 以下详细描述如何由t 。层上的数值动量喈和数值能量叼来计算i n + ，层上 的数值动量呼“和数值能量叼”如同一般的g o d u n o v 型格式一样( 见f 1 ) ，我 们的格式也按重构，发展和求网格平均三个步骤来进行 2 0 0 5 上海大学硕圭学位论文 重构：在t 。层上对的数值解进行重构对于k d v 方程而言，动量和能量流 函数( 14 ) 和( 1 5 ) 中分别涉及到了“( z ，t ) 的一阶导数和二阶导数，为了保证差分 格式具有二阶精度，重构函数目( z ；矿) 在每一网格( 一 ，马+ ) 中是z 的三次函 数，它是对该冈格中精确解t c ( e ) 的一个三次多项式逼近，从而使其对精确解具 有四阶精度的逼近，求两次导数以后也具有二阶精度，从而保证了差分格式的精 度。我们的重构将分两步进行 1 ) 我们这里已知的是函数“( 五t 。) 网格平均的值“? ，所以要用原函数插值 的方法( 见 5 】) 来进行重构我们记“( z ，t 。) 的原函数为 r ”( z ) = “( g ，t 。) 白，( 3 4 ) j z 0 其中。o 是一任意常数，它并不影响最后的重构结果由( 3 9 ) ，( 31 1 ) 式可以得 到原函数在插值节点的逼近函数值为； 。州u ( z ，t 。) d x 叫；g 巾- q 出+ 扳州枷+ 1f 1 - = 地? ， ( 3 5 ) 其中i o 是和。o 对应的网格指标首先用l a g r a n g e 插值方法( 见 3 5 ) 构造出原函 数“( z ) 的四次插值多项式( 曩w “) ，然后取 冗( z ；矿) = 五d 日( z ；w ”) ， ( 3 ，6 ) 则r ( z ；矿) 是u ( x ，t 。) 的三次多项式逼近 在网格( 勺一 ，。j + ) 中构造日( z ；w “) 时插值节点可有多种选取方式- 由此产 生的一个问题是应如何合理选取插值节点为此我们作如下的分析对于k d v 方 程而言，非线性项6 u u 。是关于z 的低阶导数项，所以我们暂时忽略非线性项，考 虑线性方程 u t + u = 0( 3 7 ) 2 0 0 5 上海大学硕士学位论文1 2 的差分方法稳定性，从中来得到启示该如何选取适当的插值节点 对于线性化方程( 3 , 7 ) ，用差商代替微商来离散方程，即地 嵋“) ，u 雨il n + 。一3 t 盘- + 3 哼一呼t ) ，将其代入到方程中 式： 丐n “= u ；一r ( 曙2 3 味- + 3 u ? 一u t ) 。；( 矿，一 得到差分格 ( 3 8 ) 其中r = 而t 。 下面我们对差分格式( 3 8 ) 进行稳定性分析，在这里我们用传播因子法( 见 3 2 ) 来进行稳定性分析 取咛= 旷e - ，代入格式( 3 8 ) ， v ”+ 1 e m ，= v “( e m 唧一r ( e t 。+ 2 3 e 妇q + 1 + 3 e m 一e ，一1 ) )( 3 9 ) 将( 3 9 ) 式两边同时除以e ，得到 v “+ 1 = v “( 1 一r ( e “一3 e ”“+ 3 一e - i a h ) )( 3 1 0 ) 所以放大因子为： g ( “，h ) i i1 一r ( e 4 。”一3 e 8 “+ 3 一e “) ：j1 + 4 s i n 2 ( 芸) ( c o s ( n ) 一1 ) + 4 i as i n 2 ( - 譬) s i n ( a h ) = 1 ( 1 + 4 as i n e ( - 警) ( c o s ( n ) 一1 ) ) 2 + 1 6 舻s i n 4 ( 警) s i n 2 ( a ) = 1 _ 1 6 矧n 2 ( 警) + 6 4 撬i n 4 ( 警) s i n 2 ( q h ) ( 3 1 1 ) 不难得到当lg ( a ， ) l 1 ，即o r i 时，差分格式条件稳定我们曾尝试了 用其它节点来离散三阶微商，但得到的格式均无法保证稳定性这表明设计我们 的格式时应采用插值节点x j 一；，+ ，一 ，+ ，和。，+ 如此选取插值节点， 所构造的差分格式能线性部分是稳定的 在如此选定的插值节点上可计算得在网格( q 一1 2 ，巧+ 1 2 ) 上月( z ；矿) 为 r ( xq u ”) = q n ，o + 曰1 1 扛一码) + 喀1 2 扛一q ) 2 + q n x q ) 3 ， ( 3 1 2 ) 2 0 0 5 上海大学硕圭学位论文 1 3 其中 a j “。2 畸一寺( u 知，一2 嵋+ 呼) ， 哼，1 = ；( “；一丐n ，) + 甄1 ( 知i - 2 u ；+ t z ，) 赤( u 知。一3 吼，+ 3 哼一u h n ) ， n 2 2 赤( 一2 u ；+ “ ，) ， 2 赤( “肿n 一乩知- + 舛一呼) 由( 3 5 ) 和( 36 ) 式可以看出函数r ( z ；u n ) 满足： ( 3 1 3 ) f 3 1 4 1 ( 3 1 5 ) r 3 1 6 1 ；z ? ：5 只( 。；u ”) d 。5 瓦1 日( 弓+ ；w ，”) 一日( 一 ；h ，”) = “? ， ( 3 1 7 ) 即r ( z ；铲) 的网格平均和原来的网格平均相同，重构函数r ( z ；u “) 满足第一个守 恒律。 2 ) 我们将r ( z ；u n ) 中最高次项系数吩n 3 用一个自由变量巧来代替重构函 数9 ( z ；矿) 在网格( 巧一1 2 ，z ，+ 1 2 ) 上取为： g ( z ；矿) = n ，n ，o + 回1 ( z z ，) + 哼2 ( z 一) 2 + s ? 一q ) 3 ， ( 3 1 8 ) 其中自由变量s ? 由第二个矸但律采碉足： 1f x j + ( z 出2 叼， ( 3 也即重构函数的能量的网格平均应和该网格上的数值能量相等将( 31 8 ) 代入 ( 3 1 9 ) 式，可得： ；f i 5 ( 掣嘲出 ( 3 2 0 ) z 哆 + n 町 州乇 瓶， = | | 2 0 0 5 上海大学硕三二学位论文 1 4 将等式( 3 2 0 ) 左边积分后得到一个关于譬的一元二次方程： n ( s ；) 2 枷s ；+ c = 0 ， ( 3 2 1 ) 肌= 鼢学肌( + 学+ 掣+ 半a ；o n2 2 婶 当b 2 4 a c 0 时，一元二次方程( 3 2 1 ) 可解出两个根： 一尘罢坐 ( 3 2 2 ) 显然譬应取使j 譬一哼3 较小的那个值。当b 2 4 a c 0 时，方程( 3 2 1 ) 无实根， 也即重构函数无法使( 31 9 ) 满足此时我们取 巧n = 一云 ( 32 3 ) 如此选取重构函数虽无法满足( 3 2 1 ) ，但可使得。( s ? ) 2 + b s ? + c 取得最小值，这 样就可能使( 3 2 1 ) 式左右两端的差最小 由于( x 一巧) 3 项的网格平均为零，因此重构函数口( z ；u “) 仍满足第一守恒律 ；f i 5 出；州扣n f 3 2 4 1 如此重构的函数夕( z i “) 是一分片三次多项式，其在各半节点q + ，2 处间断，其 左右两侧的值不相同 发展：解方程( 1 1 ) 的初值问题： 仇十他池) 。= o “骗+ 1 ( 32 5 ) i 削( z ，t 。) = 9 ( 。；“) 一o o 茁 o o 如果所解的偏微分方程是一个一阶的双曲守恒型方程，显然”( z ，t ) 是由一系列以 ( 弓+ ，t n ) 为中心的局部广义r i e m a n n 问题拼接而成对于k d v 方程的初值问题 ( 32 5 ) 情况则要复杂的多然而，将g r e e n 公式应用于方程( 12 ) 和( 1 3 ) ，如同 守恒型方程时那样，我们仍毹得到如下的积分守恒关系： x 5 + 1 2 叱，“肛层沁，u 如 一 m ( 胁巩嘣啪2 啪斑一1 m ( ，2 1 哦嘣2 嘞斑) 2 0 0 5 上海大学硕士学位论文， 一一 ! ! 和 f + l 2 ，f z j + 1 2 局叫2 u ( ”( ) ) 如2 上二，：矿( ”( z ，) ) 如 一 _ ：p 。t - f x , t n + l ( x j + 1 2 ，t ) t ( 啊1 仁t ) ，( x j 叫2 ，t ) ) 出( 3 2 7 ) 一上。f ( ，f j ( z j 一1 2 ，) ，v x ( x j - l 2 ，) ，。- - 1 2 ，t ) ) d t ) 求网格平均：n + l 层上的数值解计算为”( 。，t 。+ 。) 的网格平均 在( 32 8 ) 式两边除以h 并注意到( 3 2 4 ) 和( 3 2 5 ) ，我们不难得到 招2 8 1 丐n “= 哼一；( 嚣 一孕) ，9 ) 其中数值动量流函数盘 为 儡 2 ；。r o ( 毛+ ，) ，t k 。( 勺+ ，t ) ) d t + 1 层上的数值能量计算为( ”( z ，如+ ，) ) 的网格平均 叼= ；户5 ( 刈川) ) 出 如果( 3 2 1 ) 有实根，则( 3 1 9 ) 成立，此时按同样方式可得 矿12 叼一无t l ，j n + 一孽 ) 其中 1r t 曙 2 ；z 。f p ( 弓+ ，t ) ，( + ，z ) 。( 畅+ ；，) ) 出 ( 33 3 ) 如果( 32 1 ) 无实根，则( 3 1 9 ) 不成立但为了保证叼的守恒性，我们仍按( 3 3 2 ) 方式计算数值能量 由于初值问题( 3 - 2 5 ) 的精确解u ( 。，t ) 很难求得，因此数值动量流函数( 3 3 0 1 和数值能量流函数( 3 3 3 ) 只能由数值积分公式来计算我们用左矩形公式计算积 分注意到”( 。，t n ) = 9 ( 蜀“) ，同时注意到9 ( z ，u n ) 在各半节点上是间断的，其 出 +k o口 畦； f p 1 一h i | 哼 删 叫 跚 仁 2 0 0 5 上海大学硕士学位论文 1 6 左右两侧的值不一样，数值动量流的计算应同时照顾到两边的值，因此可具体计 算如下： 毋1 22 ； ，( 9 k + i 1 一o ；u n ) ，g x z ( x j + 一o ；n ”) ) + f ( g ( x i + + o ；u ”) ，胁( + i 1 + 0 ；u “) ) ) = ； 3 扩( 吩+ 一o ；矿) + 啦z ( x j + 一o ；u n ) + 3 9 2 ( z ，+ + o ；扩) + 9 。( 吩+ i + o ；扩) 同理，数值能量流可计算为： 魂l 2 = ； f ( g ( q + i o ：u n ) ，啦( 勺+ 一o ；u 8 ) ，( 吼。( + 一o ；矿) ) + f ( g ( 勺+ + o ；”) ，9 。( z j + + 0 ； i t ”) ，9 r 。( 吩+ + o ；缸“) ) ) = ； 4 ( 9 3 ( 茹j + 一o ；钍“) 一9 ：( j + ；一o ；u “) + 2 ( g ( z ，+ j o ；u ”) ( 9 z z ( z j + 一o ：u “) ) + 4 ( 9 3 ( 弓+ ；+ 0 ；u “) 一建( 码+ i 1 + 0 ；u “) + 2 囟( q + + 0 ；u ”) ( ( 巧+ + o ：”) ) ) ( 3 3 5 ) 对于k d v 方程，本文所设计的这种守恒型差分格式在形式上和传统的守恒型 格式存在很大的不同，该格式同时模拟了k d v 方程前鼹个相关的守恒律，格式中 不仅用到了数值动量，还用到了数值能量，从而控制非物理振荡，以保证非线性 稳定性，后面的数值算例表明该格式可以在更长的时间内保持解的性质 对于一般的差分格式，只模拟了一个守恒律现在只考虑第一个守恒律方程， 差分格式的构造如上所述，不同的是如果在第一步重构中，重构函数为( 3 ，1 2 ) 式 中的r ( z ；“n ) ，即v ( x ，t 。) = r ( z ；矿) ，同样兄( z ；“) 在各半节点上是间断的数 值动量流函数( 3 3 0 ) 用不同的积分公式，就可以构造出不同的差分格式如果用 左矩形公式计算积分，我们就得到一个显格式： 丐n “2u ? 一i ( 承 一孕 ) ，( 3 3 6 ) 其中数值动量流函数为； 盘l ，2 = ； ，( q + j o ；? n ) ，冠m ( 弓+ o ；u “) ) + ，( r ( 巧+ ；+ o ；札“) ，r 牡( 勺+ + o ；u n ) ) 2 0 0 5 上海大学硕士学位论文 1 7 = ； 3 n 2 ( x j + ；- 0 ；u “) + 皿如+ _ o ；“”) + 3 冗2 ( z j + + 。；札n ) + r 。( 茁，+ + 。：“n ) ) 【33 7 ) 这种差分格式只屉模拟了一个守恒律方程的显格式，通过数值例子可班验证 这种格式是不稳定的，计算无法进行下去。如果对于数值动量流函数( 3 3 0 ) 式， 用中矩形公式来计算积分，就可以构造出一个隐格式，差分格式如( 33 6 ) ，数值动 量流函数为： f j = + 1 2 = 寺( 嚣l 2 4 - 群j k ) ( 3 3 8 ) 其中 嚣- ：2 渺鹏+ 一o r t + ；一w ” 1 + f ( r ( x j + + o ；矿) ，见z ( x j + + o ；u n ) ) ) 耥z2 扣r ( 弓+ - - o ；u n + 1 ) ，r z ( 吩+ 一o ；u ”1 ) ) 眦) + ，( r ( z ，+ + o ；u n + 1 ) ，r ( + 1 + o ；u n + 1 ) ) ) 数值结果表明这样构造出的隐格式是稳定的，但是这种隐格式的计算结果和我们 的格式相比较，从程序实现上，我们所设计的同时保证两个相关守恒律的显格式 更容易实现一些，而且从第四章的算例1 中就可以看出，这种显格式可以在更长 的时间内保持解的形状，具有更好的稳定性 2 0 0 5 上海大学硕士学位论文 第四章：数值算例 算例1 ：考虑k d v 方程的单孤波解 取参数分别为 “扛，t ) = a s e c h 2 ( w o v t 一知)( 4 1 ) a = 2 k 2 ，u = 4 一， k = 1 ，。o = 0 ( 4 2 ) 表41 给出了采用我们的格式计算t = 1 时刻的l 1 误差和误差的阶以及l 。 误差和误差的阶其中误差的阶是按如下方式计算的，在区间( 一1 0 ，1 0 ) 上，当取n 个 点进行网格剖分的时候，网格步长记为h 。，l 1 误差记为e 。= 警，i 一奶l 工。 误差记为嚷= m a xz j