（工程力学专业论文）高维FPK方程的数值解法.pdf

浙江大学硕卜学位论文 南维f p k 方程的数值解法 摘要 本文研究高维f p k 方程的数值解法，由于只有在严格的限制条件下j j 饿得 到高维f p k 方程的精确平稳解，通常情况下只能寻求其数值解，虽然已提出了 不少数值求解高维f p k 方程的方法，但随着f p k 方程维数的增加，数据存储量 及计算速度的问题更为突出，用有限差分法对f p k 方程离散化后得到的代数方 程组的带宽是最小的，且是非常稀疏的矩阵，非常适用于迭代法进行求解，使 需要存储的数掘量最小，为了进一步加速迭代的收敛性，本文应用多重网格法， 即首先在较粗糙的网格上迭代计算得到较粗糙的近似解，然后进行插值得到细 网格的初值，再在细网格上迭代得到较精确的最终近似解，用多重网格法可使 计算时间大大缩短，并保持高的计算精度，且方程维数越高，用多重网格法的 计算效率也越高，本文用二维、三维、四维4 个系统的算例详细讨论了超松弛 迭代因子、差分格式阶次等因素对解的收敛性的影响，说明用多重网格法结合 差分法求解高维f p k 方程是非常有效的算法，该方法还可进一步推广到求解高 维f p k 方程的瞬态解及用并行算法提高计算速度。 关键词：差分法，超松弛迭代法，多重网格法，高维f p k 方程 a b s t r a c t t h ee x a c t s t a t i o n a r y s o l u t i o n s f o r h i g h d i m e n s i o n a l f o k k e r p l a n c k k o l m o g o r o v ( f p k ) e q u a t i o na r eu s u a l l yo b t a i n e du n d e r s t r i c tc o n d i t i o n s t h ef p ke q u a t i o ni s u s u a l l ys o l v e dn u m e r i c a l l y s e v e r a lm e t h o d ss u c ha sp a t hin t e g r a lm e t h o d ，f in it edif f e r e n c em e t h o d e t a l ，h a v eb e e ne m p l o y e dt on u m e r i c a l l ys o l v et h ef p ke q u a t i o n i ti s s t i l lat a s kw o r kt on u m e r i c a l l ys o l v et h eh i g h d i m e n s i o n a lf p ke q u a t i o n d u et ot h ev e r yh i g hd e g r e e o f f r e e d o m t h ef i n i t ed i f f e r e n c em e t h o da n d s u c c e s s i v eo v e r r e l a x a t i o ni n t e r a t i o na r e e m p l o y e dt os o l v et h e h i g h d i m e n s i o n a lf p ke q u a ti o ni nt h ep r e s e n tp a p e r t h ee f f i c i e n c ya n d a c c u r a c yo ft h ep r o p o s e dp r o c e d u r e sf o rd i f f e r e n to r d e ro fd i f f e r e n c e f o r ma n dt h ef a c t o ro fs u c c e s s i v eo v e r r e l a x a t i o na r ed i s s c u s s e d t h e 浙江大学颂 学位论文 岛维f p k 方程的数值解法 m u l t i g r i dm e t h o da r ea l s oe m p l o y e dt or e d u c et h ec p ut i m eo ft h e i l l e r a t i o np r o c e d u r e s f o u re x a m p l e sa r eg i f e nt oi l l u s t r a t e t h e a p p l i c a t i o no fp r o p o s e dp r o c e d u r e s i ti sp o i n t e dt h a tt h ec o m b i n a t i o n o ft h ef i n i t ed i f f e r e n c em e t h o da n dm u l t i g r i di l l e r a t i o nm e t h o di sa v e r yg o o dp r o c e d u r ef o rn u m e r i c a ls o l u ti o no fh i g hd i m e n s i o n a lf p k e q u a t i o na n d t h em e t h o dc a nb ef u t h e ra p p ll e dt op a r a l l e lp r o c e s s o r c o m p u t e r k e yw o r d s ： f i n i t ed i f f e r e n c em e t h o d ，s u c c e s s i v eo v e r r e l a x a t i o n ， m u l t i p l eg r i d ，h i g hd i m e n s i o n a lf p ke q u a t i o n s 独创性说明 稳态f p k 方程只有在严格条件下爿能得到其精确平稳解，一般情况下只能 得到其数值解。随着乃程维数的增加，未知量的数量急剧增加，求解的难度也大 大增加，因此寻求高效快速的f p k 方程数值解法一直是随机振动理论的难点之 一。 论文用不矧阶的差分法对原f p k 方程离散化，这是冈为差分法可使方程的 带宽最小，数据存储量最小。本文用超松弛迭代法求解离散化的代数方程组，讨 论了差分阶次，网格密度及超松弛迭代因子刘解精度及效率的影响。为进一步提 高求解效率，还用多重网格法加速其求解过程，算例表明，多重网格结合超松弛 迭代法特别适用于高维稳态f p k 方程的求解。 哥乞 扩耐0 1 卿焉 弘? 。爪w j 浙江大学硕十学位论文 商维f p k 方程的数值解法 第一章绪论 在自然及社会领域中存在着大量的随机扰动，如地震、风、浪波、喷气噪声 等等，乃至股市的波动及社会的动荡等等都是非常常见的随机扰动。关于随机动 力学的研究最早可上塑至上世纪初爱因斯坦关于白朗运动的研究，由于上世纪五 六十年代航空航天等方面的需求，诞生了随机振动这一学科“”。 很快就建立了较为完善的线性随机振动理论，而对非线性随机振动理论的研 究经过几十年的研究，虽然已发展了许多方法，但总的柬说尚很不成熟。对于一 扩散的马尔可夫过程，其响应的转移概率密度可用福克一普朗克一柯尔莫哥洛夫 ( f p k ) 方程描述， 设r l 为时不变随机动力系统，由以下随机微分方程描述 d x , 2 q ( 1 ) 西+ o - , n ! 蛾( ( 1 ) i = 1 ，”：= 1 ， 其中口f 为漂移系数，屯= o - , 。巳。为扩散系数，战( ，) 为独立并具有单位强度维纳过 程。 与肺随机微分方程( 1 ) 相应扩散过程的转移概率密度满足如下的f p k 方程 粤：一掣+ 昙掣 c z ) 瓦一百+ j 面 皑 其中p = p ( x ，t x o ，。) 为条件转移概率密度，并满足如下初始条件 p = p ( 五f l 而，f o ) | ，- - 8 ( x 一) ( 3 ) 转移概率密度还需满足如下归一化条件 p ( x ，l ，t 。) d x = 1 ( 4 ) 求该f p k 方程的解就得到该系统响应的统计分布，迄今为止对该f p k 的求解 主要集中在它的平稳解，即方程( 2 ) 中的概率密度与初值及时间无关p = p ( x 、， 其相应的稳态f p k 方程为 一旦掣+ 三掣：o ( 5 ) a t2 缸叙 浙江大学硕l ：学位论文商维f p k 方程的数值解法 及归一化条件 e j p ( 耵- 吒地丸= 1 ( 6 ) 已提出了诸如详细平衡等方法求解高维非线性随机系统的响应，近年来朱位 秋院士等应用哈密顿系统的可积性与共振性，通过将哈密顿系统分成不可积、完 全可积非共振、完全可积共振、部分可积非共振、部分可积共振五类，对每一类 分别得到了其精确平稳解结构的泛函形式及求解方法，这是迄今为止得到的最广 泛的几类多自由度强 线性系统的精确平稳解，且陔精确平稳解是能量非等分的 “1 。另外林家浩教授的随机振动的虚拟激励法在实际工程中有许多应用“”1 ，张 森文教授曾用小波分析法求解了f p k 方程的数值解“。 但精确平稳解只有在非常严格的条件下才能得到，一般情况下很难得到其精 确平稳解，虽然对某几类多自由度强非线性系统，我们可以用等效非线性系统法 及随机平均法等求得其近似解以扩大其应用范围”1 ，但也需要在一定的条件下才 能得到其近似平稳解，在通常情形下需数值求解该f p k 方程，而对高阶f p k 方程 的数值解法是随机振动理论走向实际应用必须要解决的重大课题。现己发展了一 些数值求解方法，如蒙特卡罗模拟法、随机步行法、路径积分法、有限元法、有 限差分法等，下面分别作简要介绍。 浙江大学帧十学位论文 商维f p k 方程的数值解法 1 1 蒙特卡罗模拟法“”】 蒙特卡罗模拟法又称数字模拟或随机模拟，该方法的基本思想是概率统计， 例如可用抛硬币的方法得到随机样本。在随机振动中，蒙特卡罗模拟法利用 计算机产生的随机样本来近似求解系统响应的统计特性，近似程度随随机样本数 的增加而提高，蒙特卡罗模拟法有如下三个基本步骤 1 、由于计算机真正能产生的随机数只有0 、1 两个数，即高电平和低电平， 要根据具体的系统激励过程把计算机原始的随机数转化为所需要的激励样本，由 算法转化得到的随机激励样本要与系统激励的概率统计特性尽量保持一致，算法 要尽量简练，这样j 能提高蒙特卡罗模拟法的计算速度。 2 、将随机激励样本代入随机运动方程求解其响应样本。 3 、从大量响应样本求解其统计特性，如概率密度、谱密度等。 蒙特卡罗模拟法是在随机振动中应用最广的数值计算方法，它是从原方程 ( 1 ) 出发进行统计分析的，适用于高维的线性与非线性随机系统，所以在随机 振动中是应用最广的数值计算方法，而且不受计算机存储误差的干扰，其精度也 比较高。蒙特卡罗模拟法的另一个优点是所产生的各个随机样本及求解得到的响 应样本相互之间没有联系，非常适合于用并行算法进行计算。随着计算机性能的 不断提高，用蒙特卡罗模拟法结合并行算法对高维非线性随机系统的响应计算范 围将不断扩大。当然该方法也有缺点，其主要缺点是提高精度困难，蒙特卡罗模 拟法所得到的结果精度受激励样本的约束而带有随机性，且结果统计的误差随样 本数( m ) 的增大而减小，其误差so cl , m ，因此要增多一位有效数字样数 本要增大1 0 0 倍。为提高样本的统计精度，特别是小概率事件的统计精度，又提 出不少的改进统计方法，如重要性样本方法等。 浙江大学硕十学位论文 商维f p k 方程的数值解法 1 2 随机步行法恤_ 2 帕 t o l a n d 。”，r o b e r t s 1 等人曾经用随机步行求解过f p k 方程。如图l 所示， 随机步行是指一个人在离散时刻r ，= j a t ( j = 0 ，1 ，2 ) 上沿y 轴随机地向左或 右走一步，人所在的位置只能是离散值儿= t 缈( 七= 0 ，1 ，+ - + - 2 ) ，这是必然的 任何时刻人所在的所有位置不可能离开儿。假设时间t i = f 时人在儿= k a y 上，当t j “= ( ，+ 1 ) a t 时刻，人向儿“= ( t + d a y 走的概率设为乩，折人向 乩一= ( 七一1 ) 妙走的概率为卜玩。随机步行是一种概率扩散过程，其一二阶增 量矩阵分别是 ( 儿，a t ) = ( 如一) 缈 ( 雎，j ，出) ：i ( 缈) ：+ 吼( 缈) ：( 缈) z 钾 为了与一维f p k 方程的概率进化过程相匹配，要使 ( 靠一q 自) 缈= q ( y k ，f ，) a t ( i ) ：x 舰一皿。 其中既+ = 1 ，口l ，6 1 为漂移与扩散系数，可由方程( 1 ) 给出。在给定缈后。 由漂移系数与扩散系数q ( 几，0 ) ，6 1 ( 儿，) 得到，和a t 。 概觏 概率q 灯 夕 图1 随机步行法的示意图 4 浙江大学顾十学位论文赢维f p k 方程的数值解法 用随机步行法求解f p k 方程的学者不是很多，其原因是该算法不是无条件稳 定的，当所设的时间步长大于它的条件稳定的最大步长时将使该随机步行法变得 不稳定。因此，用此方法求解时要用较小得时间步长，且该法很难推广于高维 f p k 方程得求解。 浙江大学硕十学位论文 岛维f p k 方程的数值解法 i 3 路径积分法“嘲 路径积分法( p a t h i n t e g r a lm e t h o d ) 是一种可用于求解较高维f p k 方程 的数值方法，它的基本思想是将f p k 方程之解表示为如下形式的泛函积分 p ( y ，r ) = d u ( y ) e x p 一p 【刑) ，y ( f ) 协7 ) p ( y 0 ，) ( 9 ) 托h 式中d u ( y ) 是一个积分测度，孝称为o n s a g e r m a c h l u p 泛函。在空日j 与时间上离 散化，以路径和代替路径积分，则 p ( y f ) = 罂垂肿慨) e x p - r 荟n - i 妣l y ，黼) “7 + 响 7 1 0 ) 一l = 卜柚l i m 兀 p ，g ( y “，y t f ) ) p ( y ) 式中g ( y ，y 1 f ) = u ke x p 【一( y 。，y 。r ) 】通常称为短时传播子。 不同的离散化规则导致不同的g 的形式，唯一的要求是它满足f p k 方程 o c r 2 ) 。一种特别简单的短时传播子是 式中 屯= 。k ( y ) + d b w ( y ) ，玩2 阳- ( y ) + d a r ( y7 ) v ，u = 1 ，2 ，聍，c + d = l ，0 c ，d i a v 与k 分别为f p k 方程的漂移系数与扩散系数，i 【瓦】i 表示 瓦】的行列式。 式( 1 1 ) 的物理意义为：由于f 很小，在短时日j 内从到y 的转移概率密度近似 为高斯分布。 短时传播子有许多不同的类型，w o j t k i c w i c z 等讨论了扩散矩阵奇异时的短 时传播子并求解t - 维f p k 方程的解汹1 ，h u a n g 等用路径积分法研究了二维平均 f p k 方程的解啪1 ，用路径积分法求解高维f p k 方程遇到的最大困难是随着求解维 数的提高，转移矩阵的维数急剧增加，数据的存储量非常大，因此该方程不适用 于求解3 维以上的f p k 方程的数值解。 6 ，绷 芝， 矿 蝴 盼 刊 甥 咿桫蛊， 乜 协 协嗽鳓 浙江大学硕十学位论文 商维f p k 方程的数值解法 1 4 有限元法嘲、 限元法是在现代科学与工程中得到了最为广泛的一种数值解法。在实际工程 中，对于许多力学问题或场问题，人们已经得到了它们应当遵循常微分方程或偏 微分方程和相应的边界条件。一般的固体力学问题必须满足几何条件、物理条件 和平衡方程。能用解析法求出精确解的是少数，而c l o u g h 在一篇关于平面特性 问题的论文中提出了有限元法的概念( f i n i t ee l e m e n tm e t h o d ) 。4 0 多年来有 限元法的理论也非常完善了，应用得到了飞速的发展，几乎遍及所有的工程技术 领域。有限元能够迅速发展成为现代工业与工程技术密不可分的一个组成部分， 除了依赖于现代计算机迅速发展的大环境之外，有限元本身具有的许多优点也吸 引了大量的理论研究人员和应用工程技术人员。有限元法的主要优点有： l 、概念浅显，容易掌握，可以在不同理论层面上建立起对有限元的理解， 既可以通过非常直观的物理解释来理解，也可以建立严格的数学理论分析来理 解。 2 、有很强的适用性，应用范围及其广泛。它不仅能成功地处理线特性力学 问题、非均匀材料、各项异性材料、非线性等问题，而且随着其理论和方法地逐 步完善于改进，能成功地用来求解热传导、流体力学、电磁场等领域的各类线性、 非线性问题。 3 、有限元采用矩阵形式表达，便于编制计算机软件，使其求解的问题方法 规范化、软件通用化，为有限元的推广和应用奠定良好的基础。目| j i 的通用有限 元软件主要有a b a q u s 、a n s y s 、m s c 系列。 用通用有限元法解非线性随机动力学问题的基本思想是将高维f p k 方程 写成它相应的弱形式或直接用加权残数法得到以节点未知量为基本变量的常微 分方程组。 应用加权残数法求解方程( 2 ) 的基本方程为 詹c 一喜静卅圭喜喜去( h y p ) 忡胁 ， 其中p 是n ( x ，k ，。) 的近似值，w a x ) 为权函数，可用下式表示 卢。( x ，f ) = 虬( x ) 玩( f ) 7 ( 1 3 ) 浙江大学硕t 学位论文 岛维f p k 方程的数值解法 其中n 。c x ) 是形函数，虎( f ) 是第m 个节点的未知量，共有个未知量，一般取 n o ( x ) = w a x ) ，将式( 1 3 ) 代入式( 1 2 ) 并完成积分，并取其极小值可得关于 节点变量得常微分方程组为 m 掣+ k f i ( ，) ：0 ( 1 4 ) 其中m = 【而f ，】，k = 【】 西( r ) = 惦( ，) ，氏( f ) 毛= l 兀。( x ) d x 啊。k n , ( x ) n j ( x ) d x 兀序刚帅瘩掣扣，掣 ， + 喜帅，掣掣 求解常微分方程组( 1 4 ) 即可得到原系统( 1 ) 的数值解，当仅需求解方程 ( 1 ) 的稳态解时，p ( t ) 与时间无关，此时方程( 1 4 ) 退化为代数方程组。求解 该常微分方程组即可得方程( 1 ) 的近似平稳解。文 3 3 还给出了有限元法结合 并行算法求该系统响应得方法。 8 浙江大学硕十学位论文 赢维f p k 方程的数值解法 1 5 有限差分法n ”“” 有限差分法是求解偏微分方程数值解的基本方法，它一般分为隐式格式与显 式格式二大类，对于一维系统，其中显式差分格式导出的方程一般可表示为： 差分法分为显示差分与隐式差分两种 + l2 “m + 矽( , u m ) ( 1 6 ) u o = u ( t o 、= a o 其中u 。表示时刻的值，+ ，表示0 + 时刻的值，g l o 为初值，厂为差分表达式， h 为时问步长，由式( 1 6 ) 知，其等式右端不含“。，只要通过递推就可得到“。 对于一维系统其隐式差分格式导出的差分方程一般可表示为： l = + 尝( ，( 。，) + f ( t m 1 ) )( 1 7 ) u o = u ( t o ) = a o 式( 1 7 ) 中各项的意义与式( 1 6 ) 相同，由式( 1 6 ) 知隐式差分格式等式两边都 含有只能通过解方程组的办法解出。 现已发展了多种显示与隐式差分格式，具体可参考文献 4 7 ，一般情况下隐 式差分格式的数值稳定较好或无条件稳定的，因此可取较大的时间步长，而显式 差分格式的数值稳定性较差，因此求解时时间步长较短。 f p k 方程( 2 ) 是典型的抛物型方程，而稳态f p k 方程( 5 ) 是典型的椭圆型 方程，本文只求稳念f p k 方程( 5 ) 的数值解，并考虑不同阶中心差分格式对解 精度的影响。 一维函数厂( t ) 的二阶、四阶、六阶中心差分格式分别为 二阶中心差分格式： 望：! ! 二岛 钒 2 z x , 乎ff i 。七f7 一2 f i 醚瞄 四阶中心差分格式： 9 ( 1 8 a ) 浙江人学俩t 学位论文 离维f p k 方程的数值解法 旦：! f ：! 二：! ! 二! ：! 二：! ! 瓠|1 2 缸。 a 2 一1 6 ( 丘l + 丘i ) 一( 厶2 + ，一2 ) - 3 0 f , 一= = - h - - _ - - _ - - - - - - _ - - - - 二- - - - - - - h = 二一 研 1 2 ( 缸) 2 六阶中心差分格式： 可一4 5 ( 厶一正。) 一9 ( 兀：一，一：) + ，+ ，一正， a 6 0 & x 旦z ：! ! ! l 五：! ：! ! 二! ! f 二! ：! ! ! ! ：! ：! ! 二! ! ! a 旷 1 8 0 6 x j 其中丘：，z + z ，兀，厶：的位置见示意图2 。 l t l 。l l ? 。 f l 。 图2 置方向的差分网格示意图 ( 1 8 b ) ( 1 8 c ) 利用方程( 1 8 ) 的差分格式将平稳f p k 方程( 5 ) 离散化，可得线性代数方 程组为 k p = 0 ( 1 9 ) 其中系数矩阵【k 】由不同差分格式得到，p = ( n ，p o ) 7 为空间离散化后各节点 的稳态概率密度组成的矢量。 由于系数矩阵【k 】一般是非正定的矩阵，且要保证p 满足归一化条件( 6 ) ，直 接求解代数方程组( 1 9 ) 需加入附加约束，且直接求解方程组( 1 9 ) 需要较大的 存储空b j ，本文我们用超松弛迭代法求解方程( 1 9 ) ，即用如下迭代公式“” o 一一，：# 一，+ 口f 一芝巧_ t w ，一圭巧一t 一，1 乞 ( z o ) i 5 lj 5 | 其中上标n 表示第n 步迭代的结果，口为超松弛迭代因子，为系数矩阵【k 】的 元素，为保证( 2 0 ) 式的收敛性，口取值范围一般为( 0 ，2 ) ，其最佳值取决于矩 1 0 浙江大学硕士学位论文1 l i i 维f p k 方程的数值解法 阵【k 】的性质，可按下式取值。7 1 2 赢 其中e ( k ) 为矩阵k 的条件数。 实际迭代时，初值取解空间均匀分布，为保证p 的非负性及迭代的收敛性， 将迭代过程中将可能产生的负只赋为零，且每经1 0 步迭代进行方程( 1 8 ) 的归 一化处理。有限差分法的主要缺点是对边界条件较敏感，如果边界不规则且处理 不当可使数值计算结果很差。 浙江大学硕七学位论文 离维f p k 方程的数值解法 1 6 本章小结 本章对高维f p k 方程的数值解法作了简要说明，由以上讨论知，蒙特卡罗模 拟法适用范围最广，但由于它需要大量样本进行统计分析，因此非常耗时，计算 效率较低；随机步行法由于适用性较窄已很少采用；路径积分法构思非常巧妙， 有明确物理意义，对求解低维随机系统非常有效，是较为成熟的一种方法，但是 对高维随机系统，由于其存储量的急剧扩大，已很难适用，一般该方法适用于3 维以下的随机系统的数值解法；有限元法可利用已有的通用结构分析软件的求解 器等模块，是非常有前途的数值方法，但现在研究尚较少，用该方法求解高维非 线性随机系统的解，最大的困难仍是随着求解维数的增加，需要巨大的数据存储 量。有限差分法是求解高维f p k 方程常用的方法之一，它的最大优点是数据存储 量最小且精度较高。由于差分法得到的以节点自由度为基本变量的代数方程组， 其带宽内仅有非常少的非零元，用传统的直接求解法这些带宽内零元素都需要存 储单元，当求解的系统维数大于4 时，数据存储量非常大，例如对4 维随机系统 划分为4 0 4 0 4 0 x 4 0 的网格，则该系统共有4 0 1 自由度，对于二阶差分格式 的半带宽约为4 0 ，但在4 0 半带宽内实际的非零元素仅为8 个，因此用迭代法 求解该方程组对存储量的要求大大降低，使解题规模得到最大限度的扩大，这正 是本论文选用差分法结合迭代法求高维f p k 方程稳念解得主要原因。 用差分法求解高维f p k 方程的另一个优点是f p k 方程的解往往是某一区域内 的局部分布，因此可较自由选取适当的计算边界，可取边界上的概率分佃为零， 这样克服了差分法对边界条件较为敏感的缺陷，因此用差分法结合超松弛迭代法 求解高维f p k 方程的解是很有前途的。为进一步提高迭代法求解的收敛速度，本 文将进一步用多重网格法加速其迭代过程，并讨论各种参数对解收敛性及精度得 影响。 浙江大学硕十学位论文 高维f p k 方程的数值解法 第二章多重网格法“5 ，叫 多重网格法( m u l t i p l eg r i dm e t h o d ，简称m g 方法) 的基本思想由前苏联 f e d o r e n k o 于1 9 6 2 年提出，并就舰则区域上的p o s s i o n 方程研究了方法的收敛 性。1 9 7 0 年后多重网格法发展十分迅速，它和有限元法及有限差分法结合，已 成为求解线性和非线性方程最有效的方法之一。1 9 7 7 年以来，多网格方法己广 泛地应用于解椭圆型边值问题，除此以外，还可用于解抛物型问题和其他依赖于 时间的问题，本征值问题、分歧问题和非线性问题，多网格方法可以有效地用于 求解积分方程。在向量机或在并行计算机上使用多网格方法更为有效。所以多网 格法是一种通用的迭代方法。 用直接方法或迭代方法求解选定网格上的代数方程组，通常我们不可能预 先知道解域上函数的性质或函数的行为，所以不可能预先确定一种合格的网格剖 分。剖分过密会导致方程组过大，而且用迭代方法求解时，亦不可能预先给出较 好的初值，因而造成计算时间过多。多重网格法就是克服以上弊病的“近乎最优” 的一种迭代方法。它达到预定精确度所需要的计算工作量仅与未知数的个数n 成正比，即0 ( n ) 。对于多重网格法最简单的一种双网格法，在细网格上作松弛 能使误差的高频分量迅速衰减，而低频分量却减少得相当缓慢，即误差的光滑效 应好，但收敛却很慢。在裉网格上由于未知量少，计算工作量比在细网格上小得 多，当然精确度会降低一些，但在细网格上高频分量显示不出来，低频分量却很 容易近似。细网格松弛和粗网格修f 通过限制和插值结合起来组成双网格法。双 网格法把在细网格上松弛使误差光滑的优点和在粗网格上低频分量容易收敛的 优点结合起来了。 浙江大学硕t 。学位论文 商维f p k 方程的数值解法 2 1 多重网格法的基本原理 我们要求解的连续问题是 l u = 其中工是一个微分算子或积分算子或泛函求极值算予。例如二维p o i s s o n 方 程第一边值问题 l u = 厂( x ，j ，) ( 工，力q u = g ( x ，y ) ( x ，y ) 硷 ( 2 1 ) 其中扛一c 昙+ 。 不论用有限差分法还是用有限元方法求解( 2 1 ) ，首先对平面区域进行剖分 若用有限差分法，则作分别与x 轴和y 轴平行的两族平行线，将q 剖分成矩形或 正方形子域。我们设a x = a y = h ，将这种剖分得到的网格点的集合称为q 。若 用有限元方法，则将q 作三角形或矩形划分，所得网格点的集合亦称为q 。 在瓯内的每一网格上可建立方程( 2 1 ) 的差分方程，用u q 表示“( ，玎) 的 近似值。我们用五点差分格式得到与微分方程( 2 1 ) 相应的差分方程 趣，= 一监丝半手堕鱼= ( 脚。 ( 2 2 ) q j = 蜀j ( ，y ，) 锄6 其中m 。为边界上网格点的集合。 用矩阵形式表示，得 l h i i = 瓦 ( 2 3 ) 其中冠。和毛为具有n 个分量的列向量，l h 为n x n 阶得大型稀疏矩阵。对于方程 ( 2 3 ) ，传统的迭代解法是在确定得一种网格瓯上采取某种迭代解法，例如g a u s s - - s e i d e l 迭代法，逐次超松弛法( s o r ) ，隐式交替方向方法( a d i ) 以及它们的 改进迭代技术。而多重网格法恰恰不是在一种确定得网格上求解，而采用不同等 级得网格剖分。假定有一组网格q ，q ，q 称为网格序列 4 浙江大学硕t 学位论文 岛维f p k 方程的数值解法 q = 0 ，l ，m ) ，简写成哦( 女= o ，l ，肘) 。随着七的增加，差分网格越来越细， 这些网格都用同一种方式剖分，都逼近同一个区域q 。相应的有 网格步长序列 魂 ( = 0 ，1 ，m ) 差分算子序列 厶 ( k = o ，l ，m ) 我们感兴趣的是要求解q 。i - - f 约差分方程 k面m=k(24) 下面我们称在g 上解厶吐。= 瓦的问题是q 。问题。在q 。上求解( 2 4 ) 称为q 。的 问题。多重网格方法的基本思想是用较粗网格q 。( t m ) 上的问题作为q 。问题 的近似。也就是说，首先近似地求解q 。问题，因为q 问题的代数方程组比r 问 题的代数方程组要小得多，解起来比q ，问题容易得多。再将q 。上的函数插值转 移到q 。上作为迭代解f 问题的初始近似值。同样，在求解q 。问题时，用比瓯 更粗一级的网格上的近似解作为q 。问题的迭代初值。 浙江大学硕p 学位论文离维f p k 方程的数值解法 2 2 双网格法慨“ 双网格法是多重网格法的基础，我们首先叙述附加一种网格的双网格法。 双网格法用两种网格，一种是桐网格，用q 。表示；另一种是细网格，用q 表示。网格步长分别取为日和h ，一股h = 2 h 。在q 和q 。上的差分算子分别 为厶和“，并设f 1 和一。存在。 设我们用粗细网格q 。、q 。上待求解的方程组分别为 k u h = d ( 2 5 a ) 厶毗=厶(25b) 其中l n 、厶为用粗细网格得到的矩阵，d 。、面。分别为粗细网格节点自由度组 成的向量，a 。、矗分别为粗细网格上的广义载荷向量。 厅司、厅习 、f 7习 bd 2 2 矗b2 4 2 j bd 、 、， 、 、i ，”、2 0 厅司7、厅习厂、旷司 bd 、 b纠 l ，坠 刳 1 21 31 4 、厂、， 、r 7 3 9 1 0 司、厅司， 、可 b纠坠刳b d 【 23 4 j 图3 双网格示意图 粗细两种网格上的值需要保持一致，下面定义由粗网到细网及由细网到粗网 上网格函数的转移算子，定义砰为限制算子，表示由细网到粗网上值的转移。 下面举例说明彰的含义，当h = 2 h 时，粗网节点和细网节点的分确如图3 所示。 其中o 表示细网节点，口表示粗网节点。可以用各种方式进行转移。以下仅举 j 6 浙江人学硕l 学位论文 岛维f p k 方程的数值解法 1 印- 【1 】：，称为直接映射算子，即 ( ，“ ) ( x ，y ) = “ ( x ，y ) ，( 工，y ) q 其中u h ( x ，y ) 为细网上的函数值。事实上，与细网点重合的点的粗网值就取 2 露称为完全加权转移算子。 定义t 为插值算子，表示由粗网到细网上值的转移。下面举例束说明， 例如，用双线性插值：露= 去 i 4 司例如用双线性插值：露= 去l ；：；i 用双网格方法迭代求解方程( 2 5 a 、2 5 b ) 时，由已知d 计算d ，“的一个迭 1 给定初值d ，在细网q 上对( 2 5 a 、2 5 b ) 作b 次迭代( 后面称松弛) ，一般 取h = 1 或2 。这样得到近似值面妒： 卮7 ：= r e l a x ( 面，厶，矗) 其中记号- 表示定义的意思。上式右端表示对算子厶、方程( 2 5 b ) 右端项以初 值i 松弛h 次。 ( 1 ) 计算细网亏损量：a r ) - 毛一厶d ? ” ( 2 ) 限制亏损量( 从细网到粗网转移亏损量) ：酩= 彤酚 ( 3 ) 在粗网q 。上求l m d 妒= a 妒的精确解面： ( 4 ) 插值( 从粗网到细网转移修j 下量) 面，= h - u ( 。n ( 5 ) 修i f 瓯上的节点函数值面紫。= 讧，+ 蔚， 1 7 浙江夫学顾p 学位论文 南维f p k 方程的数值解法 3 以匠r 。为初始近似值，在细网q 上对( 2 5 a 、2 5 b ) 松弛吃次，般吃= 1 或 2 ，松弛结果是d “：= r e 砌一( 面? ”，厶，最) 。 以上过程可以用图4 表示，这样用符号表示更加简单明了。图中，o 口、i 夕 分别表示松弛，求精确解，限制亏损量和插值修正量并对细网上的解作修正。 q q 图4线性( h ，) 双网格法迭代过程示意图 上述过程的算子表达式为 讧= 卵讧+ s h e i h ( 2 6 ) 其中群= 黠( 1 h 一，：吲彤厶) 跗为双网格迭代算子。最【毛】= ( 厶一辫) 1 最。 厶为单位算子。酚和辨为在细网q 上的迭代算子。 理论上可以证明：当v 2 = 0 ，v ，寸时，双网格方法对一般问题是收敛的， 且与h 取何值无关。 还有一种双网格法的思想就是先精确求解粗网格上节点值，然后经过多项式 插值把粗网格上节点值转化到细网格上，再在细网格上用超松弛迭代法多步迭 代，那计算精度就可与直接在细网格上用超松弛迭代法求解有相同的计算精度。 其过程为： l 、精确求解粗网格上节点值 2 、高阶样条插值把粗网格上节点值转化到细网格上 3 、在细网格上多步迭代 从双网格法的计算过程看出，它由二个要素组成：其一是在细网上松弛，其 二是粗网修正。在细网上松弛的作用是起到光滑误差的作用。细网松弛和粗网修 正通过限制和插值结合起来组成双网格法。双网格法把在细网上松弛使误差光滑 的优点和在粗网上低频分量容易收敛的优点结合起来了。 本文采用的双网格法是后者，即先在粗网格精确求解，然后用三阶样条插值 得到细网格的粗值，在细网格上进行迭代以得到较精确的数值解。 1 8 浙江大学顾 学位论文扁维f p k 方程的数值解法 第三章4 个算例分析 本章利用不同的差分格式并用超松弛迭代法及多重网格法分别求解高维f p k 方程，用4 个具体算例分析所提方法的有效性及精度，讨论了不同差分阶次、超 松弛迭代因子等对解的精度及效率的影响。 3 1 算例1 考虑受高斯白噪声激励的杜芬振子，其运动方程为 膏+ p 量- c o o x + s x l = w ( ，) ( 2 7 ) 其中p 为线性阻尼系数，为线性刚度系数，万为非线性刚度系数、w ( f ) 为 强度2 d 的高斯白噪声。 设y = j ，与方程( 2 7 ) 式相应的肺随机随机微分方程为 d x = 西 方：( 一p y + x 一西+ q 匾d d b ( f ) ( 2 8 ) 其中b ( r ) 为单位维纳过程。 与方程( 2 7 ) 相应的平稳f p k 方程为 掣一型一业坐二型二型生。 ( 2 9 ) 却2苏西 其中p ；p ( x ，y 1 为位移与速度联合稳念概率密度，平稳f p k 方程( 2 9 ) 式的精确 平稳解为： p c 五y ，= c e x p 一告( ；一警x 2 + 鲁x 4 c s 。， 其中c 为归一化常数。 当p = o 1 ，d = o 0 1 ，( o o = l ，占= 4 时，由( 2 9 ) 式得到的位移与速度联合稳态概 率密度精确平稳解及由方程( 2 9 ) 用4 阶差分法及超松弛迭代得到的数值解分别 见图5 ( a ) ( b ) ，本算例数值解用8 0 8 0 网格密度、超松弛迭代因子a = o 0 5 ， 经1 2 0 0 0 步迭代得到。 1 9 浙江大学顾十学位论文 岛维f p k 方程的数值解法 - 1 5- 1 5 ( a ) 精确解 x ，y 1 5- 1 5 ( b ) 数值解 图5 位移与速度的联合稳态概率密度 ， 浙江大学顾十学位论文曲j 维f p k 方程的数值解法 图6 ( a ) ( b ) 分别给出了用不同阶的差分格式及网格密度得到的边际概率密 度p ( x ) 与精确解的差的绝对值，由图6 知，差分格式的阶次越高，数值解精度 越高。 ( a ) 边际概率密度p ( x ) 及绝对误差 ( b ) 边际概率密度p ( y ) e 绷e 图6 边际概率密度p ( x ) 、p ( y ) 的精确平稳解及各阶差分格式近似解的误差 2 浙江大学硕i 擘位论文 岛维f p k 方程的数值解法 图7 ( a ) ( b ) 分别给出了由( 3 0 ) 式z 一i - - 概率密度p ( x ) ，p ( y ) 及由 三种不同的差分网格密度( 4 0 4 0 ，8 0 x 8 0 ，1 6 0 1 6 0 ) 得到的近似解的边际概率密 度与精确平稳解的差的绝对值，由图7 可知，差分网格越细，数值解精度越高， 但需迭代的步数也越多。 x ( a ) 边际概率密度p ( z ) 及误差 ( b ) 边际概率密度p ( y ) 及误差 图7 边际概率密度p ( x ) 、p ( y ) 的精确平稳解及用4 阶差分格式 浙江大学坝l 。学位论文岛维f p k 方程的数值解法 由于方程( 2 7 ) 对应的扩散矩阵是奇异的，因此( 2 8 ) 式由差分法得到的线 性代数方程组( 1 9 ) 的矩阵k 并不是对角占优的，条件数p ( k 1 较大，因此收敛 的超松弛因子口很小，本例中口的取值范围一般为1 2 = o 0 2 ，0 0 8 ，当盯大于0 0 8 时超松弛迭代不收敛，当口很小的时收敛很慢，在可能的超松弛因子口的取值范 围内，口取值越大，收敛速度越快，因此岱应尽可能取得大一些。 图8 给出了用4 阶差分法得到的近似解与精确平稳解误差绝对值的最大值 与迭代步数与网格密度的关系，由图8 知，随着迭代步数的增加，最大误差的绝 对值越小，但对给定的网格密度，其精度是有限的，网格越粗，近似解精度越差， 但收敛速度越快：网格越细，近似解精度越高，但收敛速度越慢。 由图8 可知，为了提高超松弛迭代法的收敛速度并保持高的近似解精度， 一种可行的方法是先用粗网格进行迭代，达到精度极限后，用插值方法得到细网 格迭代的初值，再用细网格进行迭代得到更高精度的近似解。多重网格法已在其 他偏微分方程求解中得到广泛应用，在本文中应用多重网格法特别适合用于求解 高维平稳f p k 方程的数值解。 i t e r a t i v en u m b e r 图8 用4 阶差分格式得到的近似解与精确平稳解误差的 最大值与网格密度与迭代步数的关系 浙江大学硕 。学位论文商维f p k 方程的数值解法 由于f p k 方程可以用其它的数值方法求解，如蒙特卡罗数字模拟法、路径 积分法等等。为了比较用差分法与其它方法的差别，图9 给出了差分法与直接由 方程( 2 7 ) 用蒙特卡罗模拟及用路径积分法得到的数值解与精确解误差的比较， 其中蒙特卡罗模拟时用1 0 0 0 个样本，每个样本用有效的1 0 0 0 0 点进行统计，路 径积分法及4 阶差分法结合超松弛迭代法都使用8 0 8 0 网格，经1 2 0 0 0 步迭代 得到，初值都为均匀分布，出图9 知，4 阶差分法结合超松弛迭代法得到的近似 解精度最高，且在相同计算机上用的时问最少。 x 图9 差分法与蒙特卡罗模拟及路径积分法结果的比较 由于本例是二维问题，用差分法结合迭代法所需c p u 很少，一般只要几分钟 即可，因此用多重网格迭代法加速收敛显得不迫切，因此本算例没有进一步用多 重网格法进行计算。 浙江大学硕t 学位论文 商维f p k 方程的数值解法 3 2 算例2 考虑如下式所示的谐和与高斯白噪声作用下的单自由度线性系统， 王+ i + 彩2 x = e c o s f 2 t + w ( r ) ( 3 1 ) 其中为线性阻尼，国为线性振子的频率，w ( ，) 为强度为2 d 的高斯白噪声。应 用标准随机平均法可导得以幅值a 及相位差r 为基本变量的伊藤随机微分方程“1 谢= 心爿+ 扣r + 南卜- - - - 历a s , ( ，) = 叭咿业1 7 0 酬一 ( 3 2 ) 订= ( 刍o s f + q - o 击+ 鱼c o a 娼( ，) 吨( 盯肌鲁础) 其中e ( r ) 为独立的单位维纳过程，r 为相位差 口：t g 一三 其中0 为相位 与方程( 3 2 ) 相应的稳态f p k 方程为 旦垒+ 旦塑c 3 ( m t p ) 猁：o c p 2a a 2 m 2 日2a ，2 a 口 a y 方程( 3 4 ) 的稳态概率密度p ( 口，) 为“” 出力p p + ( 3 3 ) ( 3 4 ) ( 3 5 ) 其中c 为归一化常数，s = 2 d ( q 一 o ) i p o ：2 。 图1 0 ( a ) 为= o 1 ，讲= q = 1 0 ，e = o 1 ，d = 0 0 0 4 时由( 3 5 ) 式得到的稳态 概率密度，图1 0 ( b ) 为用四阶差分格式( 3 b ) 、8 0 8 0 的网格密度、超松弛因 子口= 1 3 时，经3 0 0 0 步迭代得到的结果。 p ( a oo 精确解 oo 数值解 图1 0位移与速度的联合稳态概率密度 图1 1 ( a ) 给出了由( 3 5 ) 式得