（计算数学专业论文）hamiltonjacobi方程的一个数值解法及其应用.pdf

中山大学硕士论文 h a m i l t o n j a c o b i 方程的一个数值解法 及其应用 专业 学位申请人 计算数学 杨胜龙 指导教师 ： 朱庆勇教授 摘要 本文在文i l 】【2 j 之上，分析了迎风紧致群速度格式( u p w i n dc o m p a c tf i n i t e d i f f e r e n c es c h e m ew i t hg r o u pv e l o c i t yc o n t r o l 简称u c g v c 格式) 。首先分析了数 值解行为特征，由此提出群速度概念，从物理的角度分析了它在控制数值格式稳 定性中的应用。利用h a m i l t o n - j a c o b i 方程与双曲守恒律方程的关系，将u c g v c 格式用于求解h a m i l t o n - j a c o b i 方程，并推广到多维及方程组情形。最后利用 u c g v c 格式，对图像的平移、旋转、变形等运动进行了数值模拟和计算。计算 结果表明，u c g v c 格式是令人满意的。 关键词： h a m i l t o n j a c o b i 方程，双曲守恒律，迎风紧致方法，u c g v c 格式 中山大学硕士论文 a nn u m e r i c a la p p r o a c ht o h a m i l t o n j a c o b i e q u a t i o n sa n d i t sa p p l i c a t i o n m a j o r ： c o m p u t a t i o n a lm a t h e m a t i c s n a m e ： y a n gs h e n g l o n g s u p e r v i s o r ：p r o f e s s o rz h u q i n g y o n g a b s t r a c t b a s e do np a p e r s 【”w e i n v e s t i g a t e d a n u p w i n dc o m p a c tf i n i t e d i f f e r e n c e s c h e m ew i t l lg r o u pv e l o c i t yc o n t r o l ( u c g v c ) f i r s t l y , w ed i s c u s s e dt h ec h a r a c t e r i s t i c b e h a v i o r so fn u m e r i c a ls o l u t i o n sa n dt h ea p p l i c a t i o nt ot h en u m e r i c a ls t a b i l i t y n e x t w ee m p l o y e du c g v ct os o l v et h eh a m i l t o n j a c o b ie q u a t i o n s a c c o r d i n gt ot h e r e l a t i o n s h i pb e t w e e nh a m i l t o n j a c o b ie q u a t i o n sa n dh y p e r b o l i cc o n s e r v a t i o nl a w s ， w es o l v e do n ed i m e n s i o n a lh a m i l t o n j a c o b i p r o b l e m s a n de x t e n d e dt o h i 曲 d i m e n s i o n a lc a s e s t h em o v e m e n to fi m a g ew a ss i m u l a t e db yu c g v cs c h e m e n u m e r i c a lr e s u l t si n d i c a t et h i ss c h e m ei ss a t i s f a c t o r y k e yw o r d s ： h a m i l t o n - j a c o b ie q u a t i o n s ，h y p e r b o l i cc o n s e r v a t i o nl a w s ，u p w i n dc o m p a c t m e t h o d ，u c g v cs c h e m e i i i 中山大学硕士论文 1 1 引言 第一章：前言 h a m i l t o n - j a c o b i 方程( 简称h - j 方程) ，如一维情形： 粤哆y 矿) = x e d o ( 1 1 1 ) 【矿( 工，o ) = 丸( x ) 它来源与晟优控制理论、微分方程、流体力学等学科1 4 i i ”，目前，h a m i l t o n j a c o b i 方程在光 学、计算流体、控制系统、以及网格生成等方面有非常重要的应用，并产生了像用于计算流 体界面的水平集方法【6 】等新理论。 n 维空间下一般形式的h a m i l t o n j a e o b i 方程： j 谚+ 日( _ ，恐吒，屯，噍破) 。o( 1 - 1 - 2 ) 【庐( 为，叠吒，o ) = 纯( 葺，毪矗) 方程( 1 1 2 ) 解的性质很复杂，在满足熵条件下是唯一的。一方面，h j 方程的弱解存在但 不是唯一：另一方面，即使初值和h a m i l t o n 函数充分光滑时，它的解也可以是不光滑的， 可能会在某一时刻出现间断情况( 类似于双曲守恒律问题中的激波) 导致解曲面( 线) 出 现尖点或扭结等现象。h a m i l t o n - j a c o b i 方程解的这些性质和双曲型守恒律方程解的性质相 似，实际上，他们有着紧密的联系。在一维情况下，h - j 方程与双曲守恒律方程是等价的。 对于一维情况下，由一维h a m i l t o n j a c o b i 方程： j + 磐假o “即刈( 1 1 3 ) 【庐( x ，o ) = 丸( x ) 令 i “( x ，t ) = 丸如f ) i u o ( x ) = 以( x ，0 ) 对h a m i l t o n j a c o b i 方程两边关于x 求导，推出与之等价的双曲型守恒律方程初边值问题： q + 参( “，) = o x r 伊o ( 1 1 4 ) i u ( x ，o ) = u o ( x ) 二者关系是：“= 丸，妒( x ) = c + r “o ) d s ，r o m e o ( 工) w 驷，且 以o ，0 ) = u o ( x ) 口矿僻) n r ( r ) ，b v ( r ) 表示有界变差函数空间，f ( u ) 为l i p s c h i t z 连 续的h a m i l t o n 函数。 h a m i l t o n - j a c o b i 方程的一个数值解法及其应用 对于多维问题。虽然没有上述一样的转化，但可以从形式上认为h a m i l t o n j a c o b i 方程 是双曲型守恒律的积分形式，有多维h - j 方程与弱守恒律方程组等价5 】【7 j 口1 ，求 h a m i l t o n j a c o b i 方程可以等价求解弱守恒律方程组。所以参考双曲守恒律方程组的解法可 以建立相应的多维h a m i l t o n j a c o b i 方程的数值逼近格式。 通常确定h - j 方程的某一弱解，需要额外增加一些条件。比较常见的具有物理意义的解 是粘性极限解( 简称粘性解) ，是通过将陛消去法得到的。文1 7 1 从粘性解的角度分析了双曲 型守恒律与h a m i l t o n j a c o b i 方程的关系。由于在一维情况下，h j 方程与守恒律方程有着等 价关系，o s h e r 和s h u1 1 0 l 借用守恒律方程的数值解方法推导出h j 方程的e n o 格式，而 c o r r i a s ，f a l c o n e 和n a t a l i n ii j ”反过来利用h j 方程推导出守恒律方程的差分格式。 在计算资源有限的条件下，采用高精度格式是有益的i ”1 【”1 。为了求得双曲守恒律方程 的高精度数值解以及有较好的分辨率，人们提出了许多很好的数值方法。采用现有的一些格 式在计算激波时数值解中将有数值振荡产生。为克服数值振荡。已发展了许多有效的方法， 其中包括t v d ，n n d ，e n o , w e n o 等类型的格式。这些格式已经成功地用于求解实际问题。 然而现有的一些高分辨率格式多是从数学观点出发而构造的，这些方法没有很好的研究数值 解产生振荡的根本原因，也没有针对其原因对格式进行改造。群速度控制方法是近年来出现 的求解双曲守恒律方程的一种新方法2 ”】，它是从物理角度分析非物理振荡产生 的原因，并提出改进数值解的方法。该方法构造简便，物理意义清楚，利于实现，分辨率较 高，没有非物理数值振荡产生。 h j 方程的数值方法主要分为：传统的有限差分法，其文献相对较多一些 7 1 1 ”m 1 ，另 外还有有限体积法【”，李祥贵等人用有限元方法也得到很好的结果i i i9 】i ，用迎风紧致格 式求解h j 方程的文献较少i i i ”。而迎风紧致群速度控制法有很多优点，迎风紧致格式有 着精度高、网格基架少的优点，利用数值解的群速度特性和迎风紧致格式本身的耗散性基本 抑制了数值振荡的产生，使计算易于稳定，可用于求解非定常多尺度的复杂流动问题，研究 其细微结构和机理已经成功地应用于求解双曲守恒律方程。 本文的思想是将迎风紧致群速度格式用于求解h j 方程，并借助双曲守恒律方程已有的 结果分析格式的性质。根据文l q 2 l 将u c g v c 格式应用于求解h a m i l t o n - j a c o b i 方程。对于一 维情况等价地利用了双曲守恒律已有的结果对格式的性质进行了分析，将这种新的u c g v c 格式推广到了多维和方程组得情况。采用本文所述方法，对于一些界面追踪问题进行了数值 模拟，取得了较好的效果。 中山大学硕士论文 1 , 2 古典解、弱解、熵条件的概念 1 , 2 1 古典解 例如考察如下一阶拟线性双曲型方程的初值问题 导+ 吴( 乌：0 2a t。x 、。 “i - o = p ( x ) ( f o ，- - - o o 0 。定义平移算子瓦为 ( 0 v ) ( 砷= v + 功，则，e ，可定义如下： 等量。一i ，；i 一互j ， ：；( + ) ，： ( 2 2 1 ) 9 - h a m i l t o n j a e o b i 方程的一个数值解法及其应用 其中1 为恒等算子，下标x 表明算于作用在x 方向上。基于f o u r i e r 变换，可得算子， ，霹的象征。 例如嚣= 瓦一，有： u ( x ) = ( 瓦一1 ) u ( x ) = p “舯a ( k ) a k p “站卉( ) 班 = p 朋( e “札”- 1 ) d ( k ) d k 所以有： a ( x ，k ) = ( e “一1 ) = c o s ( a x ，七) 一1 + i s i n ( a x ，七) 即 = c o s ( a x ，女) 一i + f - s i n ( a x ，七) 在一维情况中，有：毒= c o s a 一l + i s i n a 以= 女x ) 同理有：乏= i - c o s a + i s i n a 母= ；( 母+ 毒) = i s i n 口 鸷= 越置= 2 c o s a 一2 其中，一万口万，i 2 = - 1 。显然，向前算子e 的f o u r i e r 象征的实部和虚部分别为： 5 r = c o s d 一1 ，和函= s i n a 。同样的，对向后算子有西= 1 - c o s a 和西= s i n t z ，中 心算子有舌，= o 和舌。= s i n a ，窭有孑，= 2 c o s a 一2 和占，= 0 。 2 3 半离散方程 对应于式( 2 1 1 ) 的半离散逼近式( 只对空间导数进行离散) 为： 挈+ c 盟：o ( 2 3 1 ) a t x 誓为( 孰的差分逼近式。 现在讨论差分方程的准确解。在连续函数空间设有函数“( x ，t ) = ( t ) e x p ( i k x ) ，其x 方 向的阶导数为：u x ( x ，r ) = i k f i ( t ) e x p ( i k x ) 。可将其改为a x u x = d ；f i ( t ) e x p ( i o c ) ，其中 j = 口i ( 口= t x ) 。 对于离散空间可取初值为： 中山大学硕士论文 u ( x j ，o ) = e x p ( i k x j ) ( 2 3 2 ) 将u ( x j ，f ) = 五( t ) e x p ( i c c j ) ，正叶= 6 ( t ) e x p ( i ) 代入式( 2 3 1 ) ，得准确差分解： “c 力= e x p 一c 七吾r e x p ( 衍c - 一c 吾， c z a ， 舌，+ f 孑，= 舌，毒为舌实部，耷为舌虚部。 对照式( 2 i 7 ) 和式( 2 3 3 ) 可知。对取定波数k ( 或口) ，数值解逼近精确解要求 生斗o ，尘斗1 2 4 数值解的耗散性 文章指出，可依据爱对格式进行分类。 对c o 情况： ( i ) 毋 o 时，格式是耗散型的，耗散型格式是稳定的。 文章【4 9 l 同时指出，还可以利用修正方程中的余项分析格式的耗散特性，并由此对不稳定的 差分格式修改，使之成为稳定的差分格式。 例如当c o 时，e 是耗散格式，而是反耗散格式，则是非耗散格式。反之。c 譬 是耗散格式，是反耗散格式。 2 5 群速度效应及其在微分方程数值解中的应用 群速度是流体力学的概念，它是指一系列单个波的叠加而形成的波群的传播速度。 2 5 1 数值解的群速度 针对离散方程( 2 3 1 ) 的解式( 2 3 3 ) 定义群速度为 d ( o t ) = 丢( 孔嗍 ( 2 5 1 1 ) 显然，对应于微分方程( 2 1 1 ) 的准确解式( 2 1 7 ) 有参：口，d ( 口) = 1 。文章”1 ，依据 各曲线奢( 口) 随口的变化率将格式分为如下类型 h a m i l t o n - j a e o b i 方程的一个数值解法及其应用 ( 1 ) 快型格式( f a s t ，简称f s t ) d ( 口) l ，( 0 口丌) ( 2 ) 慢型格式( s l o w ，简称s l w ) d ( 口) l ，( 0 l ，( o 口蔓 ，r ) d ( c t ) 1 ，( 1 6 是，算子l = ( 霹一耐群+ o 5 耐) 是混合型算子，对于低 频波段为快算子。 证明：考虑到= j a x ，l = x j + a x ，将吩= ( f ) e x p ( j 呜) 代入则有： l = i s i nc z i 0 - s i n a ( 2 c o s a t 一2 ) + o 5 0 - ( 2 c o s a 一2 ) 2 z ：( t ) e x p ( i b c j ) 所以 万= i s i nc z - i 0 - s i n a ( 2 c o s a - 2 ) + o 5 0 - ( 2 c o s z - 2 ) 2 】 = 2 0 - ( e o s 口t - 1 ) 2 + ( s i n c r + 2 0 - s i n a 一2 0 - c o s a s i n a ) l 口i 丌 舌f ：s i n a i + 2 盯s i nc z 一2 0 c o s a s i n a t 舌，：2 0 - ( c o s a n 2 ；( 孑j ( 口”= c o s 口+ 2 仃c o s 口+ 2 盯s i n 2 口- - 2 盯c o s 2 口 = ( 1 + 2 0 - ) c o s a t + 2 0 一4 0 c o s 2 口 1 j 4 0 - c o s 2 口一( 1 + 2 0 - ) c o sc r 一2 0 - + 1 0 因为4 0 c o s 2 a - ( 1 + 2 仃) c o s a - 2 0 - + 1 = 0 有根1 ，1 4 仃一1 2 ，存在 o ，使得i _ c o s o r c o s = 石i 一；成立；否则就不存在 o ，使得不等式成立。也就是说，当i 口i a r c c 。s ( 去一圭) 时，芝( 磊( 口) ) 1 。当口 落在该区间外时，芝( 舌- ) ) 1 。所以盯 i 1 时，该算子是混合型算子。由于d 口o 西 ) = 2 0 - ( c o s a - 1 ) 2 0 ，v i 口喀7 ，这里口 0 ，该算子是耗散型的，它对各种波长 的误差都起衰减作用。l 口i 越大，衰减越快。对于低频波段为快算子。证毕 口 定理2 5 2 l - ( + 碱0 q 2 + o 5 碱2 巳2 ) 叶，j 12 盯o 是慢算子。 1 3 h a m i l t o n - j a c o b i 方程的一个数值解法及其应用 证明：类似t h2 5 1 ，略。 2 5 2 群速度控制法简述 在讨论之前先作一定义，设波的传播速度c o 。则称激波的左侧为波后，激波右侧为波 前。当c - - 0 ，= l = 0 并且要求满足c f l 条件 姚r a t ? i k na 瓦, 1 & o 在计算中，采用的是二阶或三阶形式。 二阶t v dr u n g e - k u t t a 格式为： u 1 = u “+ a t l ( u “1 u n + l = ! u n + ! u ( 1 ) + 三，f u ( 1 1 222 、 三阶t v dr u n g e - k u t t a 格式为： “1 = u 。+ f l ( u “1 u = 4 u a + 1 4 u ( 1 ) + ；r 上( u ( 1 ) “一+ 1 ) ：；u n + 1 2 - - u ( 2 ) + 2 1a t 上( u ( 4 ) jjj 同时它童1 5 0 1 中怕据羽i 了四阶、再阶糟席的时闻离散格式。 ( 3 1 3 ) ( 3 1 4 ) h a m i l t o n - j a c o b i 方程的一个数值解法及其应用 3 2 紧致算子象征 我们在这里提出这篇文章中用迎风紧致算子来估计空间微商它的符号形式为： 。；( 1 2 鳄+ ：1 吒2 ) _ 1 ( 一2 彰) ( 3 2 1 ) 由上面讨论，可知它的f o u r i e r 象征为： 妒= 丁- 4 e e o s o ! r + 4 e 一+ i s i n a ( 一万口万)( 3 2 2 ) 二c o s 搿+ = 一i 2 e s i n 口 相应地：( 3 2 3 ) 8 p 2s i n t 2 - 8 s 2s i n t t t c o s o r + ! s i n 口c o s 口+ 一2s i n 口 西= 1 i c o si 2 l c o s i 五4 i cs i n 产。2 4 一 。 + 一口+ 一+ 口 y，y 并且一； s 0 i 庐( x ，0 ) = 蔬( x ) 等价的双曲守恒律方程是： 卜+ 六( “) 。o x r ，t o ( 3 5 1 1 ) 【“( 五o ) = u 0 ( x ) 1 1 1 【2 l 考虑其半离散形式，解该方程的一个紧致差分格式可以表示为： 一d u j ：一坚丝! ( 3 ，m ) d t血 k j 荨 i - 2 斛+ 了l 巳2 o 中山大学硕士论文 誓一去【一髟( 也】- 砌) ( 3 5 1 - 3 ) 出缸 。、。4、。 这里a ，b 均为算子，其中 一= ，一2 鹾+ ：，口= 霹一2 4 ，( 一j 1 s 争 由上面的定义，可得 彳叶= ( 丢+ 占) 叶一+ ；叶t ( 一占) 叶+ 。 b - = 一( 圭+ 2 占) 叶一，+ 4 s 叶+ ( 圭一2 占) 叶。 令 弓= a u j 则格式( 3 5 1 3 ) 可以表示为： 堕=一1(aodta x 一2 ) m ) 、7 j 为数值振荡，同时保持格式的精度，使用如下开关函数”1 矿i ， 1 1 坐! ! 二! 生竺：! ! 七 纺2 卜竺兰筮笺笠! 型： 七 【o 5 ( i 叶“一叶i + i 叶一u j 1 1 ) + o 5 ( i 叶“+ 2 u j + u j 一1 1 ) 这里k 称为“阀值”，计算中可取k = o 9 5 在光滑区域纺= o ，在激波附近纺2 1 。开关函数能够保证在光滑区域具有d ( 一) 的精 度。 构造如下群速度控制的迎风紧致差分格式( 盯= 妄) ： 鲁一去m ，叫m 卜胴埘+ 去群呐吲 ( 3 ， 。) 一乏i 【伤口巧】+ 篆【， ) ，一竹( ，( “l 一，( 刃埘 其中s ，( “) 。，一j 1 g 0 ，在激波右侧有：s s j f f j 0 计算出巧后，还需要从巧恢复到叶由式( 3 5 1 6 ) 的计算公式为： u j = a - l f f _ f = ( ，- 2 鹾+ 1 2 ) 一1 巧 2 5 坚! 竺坐! ! ：! 竺! ! ! 互矍堕二尘塑堡竖鎏丝基生旦 在计算时，把3 5 1 4 式及3 5 1 5 式代入3 5 1 3 式，由3 5 1 2 式，可形成一个求解 f ( u j ) 的线性代数方程组，系数矩阵是三对角的，当取定行从左到右三个元素分别是 i 1 + ，j 2 ，i 1 一占。当给定边界条件，可用追赶法求解方程组，为了保证方程组有唯一解， 要求系数矩阵是严格对角占优的，即要求i ：一引+ l i l + s i ；成立，可求解出一； s o 时，前三项变为 一乏i 霹霹+ 乏i 霹) l 厂( 厅o 】，由定理2 5 _ 1 知道 该算子为混合型算子，在此，可以将其视为快算子。 哆 一 1 2 ， 到 。 孺 一 l 回 弧 阻d 一讲栅51 ，理己j育可接苴明证 h a m i l t o n - j a c o b l 方程的一个数值解法及其应用 即a 强对角占优矩阵。那么，彳_ 1 的l ，范数是有界的。即 i a - 1 1 0 南 证明。略。 由式( 3 5 1 4 ) 可以得出系数矩阵为： a ； o 占o 1 一一占一 6 o三+ 占 6 oo 由引理3 5 1 5 ，可得出l | a - i l l l 是有界的- o0 o0 21 一一一占 36 o三+ 6 定理3 5 1 6 半离散格式( 3 5 1 8 ) 是t v b 的。即 t v ( u 。) = i 略一衅l 蔓c j 证明由定理3 5 1 4 和引理3 5 1 5 ，有 t v ( u ”) = 略，一q i j = i ( a 。矿) 川- ( a 。矿) ， j - - - i i a 。忆略。一霹j ， s c 其中c 2