（计算机软件与理论专业论文）微分方程周期边值问题的数值求解.pdf

微分方程周期边值问题的数值求解 乔田田( 计算机软件与理论) 指导教师：李维国教授 摘要 本文的研究主要分为两个方面，包括常微分方程周期边值问题在不 同情况下的数值解法和椭圆型偏微分方程广义解的存在唯一性的论证。 首先，详细地讨论了对于常微分方程的一般周期边值问题、奇异周期边 值问题和在计算周期解时出现刚性现象这三类情况的不同求解策略和新 的方法，并通过对典型的数值实例进行编程求解，得到了相应的数值试 验结果，绘制了解曲线。结果表明，求解的方法是十分有效和实用的。 特别是对于计算周期解时出现刚性问题的求解，本文使用了变步长的 g e a r 方法，并给出了相应的理论推导，这种方法的使用大大减少了工作 量，节约了计算时间。本文所提供的方法和程序，可以作为数值求解微 分方程周期边值问题的参考。另一方面，本文对于椭圆型偏微分方程的 d i f i c h l e t 边值问题广义解的存在唯一性进行了讨论，利用了一个全局反 函数定理和i l l i n - m a x 原理的一个新的非变分形式，在共振的条件下给出 了解的存在唯一性定理及其证明的过程，同时还列举出实例说明了本文 的定理推广了已知的结果。 关键词：常微分方程，周期边值问题，数值解法，椭圆型偏微分方程广 义解，存在唯一性 n u m e r i c a ls o l u t i o n so fp e r i o d i cb o u n d a r yv a l u ep r o b l e m o fd i f f e r e n t i a le q u a t i o n s q i a ot i a n d a n ( c o m p u t e rs o f t w a r ea n dt h e o r y l d i r e c t e db yp r o f e s s o rl iw e i - g u o a b s t r a c t i nt h i st h e s i s ，t h em a i nr e s e a r c hi n c l u d e st w oa s p e c t s t h e ya r en u m e r i c a l s o l u t i o n st 0p e r i o d i cb o u n d a r yv a l u ep r o b l e m ( p b v p ) o f o r d i n a r yd i f f e r e n t i a l e q u a t i o n s ( o d e s ) i nd i f f e r e mc a s e sa n dp r o o f st oau n i q u ee x i s t e n c eo f g e n e r a l i z e ds o l u t i o no f e l l i p t i c a lp a r t i a ld i f f e r e n t i a le q u a t i o n s ( p d e s ) f i r s t l y ， i nt h eg e n e r a lc a s e s ，s i n g u l a rp b v po fo d e sa n dc a l c u l a t i n gp e r i o d i c s o l u t i o n sw i t l ls t i f fp r o b l e m s d i f f e r e n tm e t h o d sa r ed i s c u s s e d t h er e s u l t so f n u m e r i c a le x a m p l e ss h o wt h a tt h e s em e t h o d sa r ep r a c t i c a b l ea n de f f i c i e n t s p e c i a l l y , g e a r sm e t h o dw i mv a r i a b l es t e pi su s e df o rc a l c u l a t i n gp e r i o d i c s o l u t i o n sw i t hs t i f fp r o b l e m sa n dc o r r e s p o n d i n gf o r m u l a sa r ed e d u c e dw i t h u s eo ft h em e t h o dd e c r e a s i n gt h ea m o u n to fw o r ka n dc o m p u t a t i o n a lt i m e s e c o n d l y , au n i q u ee x i s t e n c eo fg e n e r a l i z e ds o l u t i o no fd i r i c h l e tb v po f e l l i p t i c a lp d e si s d i s c u s s e d u n d e rt h er e s o n a n c ec o n d i t i o n s ，an e w n o n - v a r i a t i o n a lv e r s i o no ft h er a i n m a xp r i n c i p l ea n dag l o b a li n v e r s e f u n c t i o nt h e o r e ma r eu s e d u n i q u ee x i s t e n c et h e o r e mo fs o l u t i o ni sp r o v e d t h i st h e o r e me x t e n d st h er e s u l t st h a ta r ek n o w n b yn o w k e yw o r d s ：o d e s ，p b v p , n u m e r i c a ls o l u t i o n s ，g e n e r a l i z e ds o l u t i o no f e l l i p t i c a lp d e s ，u n i q u ee x i s t e n c e i i i 独创性声明 本人声明所星交的论文是我个人在导师指导下进行的研究工作及取 得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论 文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得中国 石油大学或其它教育机构的学位或证书而使用过的材料。与我一同工作 的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了 谢意。 签名：煎望望 州年，月弓曰 关于论文使用授权的说明 本人完全了解中国石油大学有关保留、使用学位论文的规定，即： 学校有权保留送交论文的复印件及电子版，允许论文被查阅和借阅；学 校可以公布论文的全部或部分内容，可以采用影印、缩印或其他复制手 段保存论文。 ( 保密论文在解密后应遵守此规定) 学生繇竺州年，月夕日 新勰锛耻州年s 月弓日 中国石油大学( 华东) 硕士论文第1 章前言 第1 章前言 1 1 课题来源、提出背景及其研究意义 本课题来源于山东省自然科学基金资助项目。 周期现象是自然界普遍存在的现象之一，在物理、化学、生物、 电子学和工程技术等问题中经常会遇到，在数学上微分方程领域内关 于微分方程周期解的问题是一个常见而又困难的问题，至今为止这方 面的进步甚少【l l 。 由于计算机的普及和计算技术的发展，人们的研究思路和研究方 法发生了改变，常常是先计算后分析。从计算所绘出的图像中观察到 周期解可能具有的性质，然后再从理论上给予证明。此外，还有很多 周期边值问题在工程应用中并不是要求对其周期解性质的准确了解， 而更关注的是用数值算法将周期解的具体结果求得出来。所以，对于 微分方程周期解的数值计算既是可视化的需要又可以弥补理论上的先 天不足。目前已有的数值计算方法，仍对周期边值问题的求解效果不 尽如人意，而且很多算法的使用仅仅局限于具体的模型j ，而缺少一 般性的数值算法，需要进一步的创新。因此，对于微分方程周期边值 问题数值解法的研究是极为重要的。这不仅仅在实际计算中可以解决 实质性的问题，而且也可以给理论进展提供较为直观的依据。 对微分方程周期边值问题的求解，与解的存在唯一性是密不可分 的。长期以来，对于各种类型的微分方程解的存在唯一性的理论研究 是非常多的，很多专家学者从事这方面的理论研究，也取得了丰硕的 中国石油大学( 华东) 硕士论文第l 章前言 成果。众所周知，线性微分方程的研究相对来说已经较为成熟，而非 线性微分方程的研究却是非常复杂多样的，研究它的基本方法就是线 性方程的非线性扰动。因此，对于常见的重要类型的非线性微分方程 的可解性和多解的存在性问题，是讨论较多的。有关的书籍记载了使 用不同的判断工具得到的诸多的相关结论【4 1 。所以，可以通过将原有 的判断条件加以改变，使定理判断的使用范围扩大，增强理论判断的 有效性 本课题就是在这些情况下提出的，在数值算法方面，以原有方法 为基础通过对计算方法采用一些有效的策略，寻求到更多的更好的求 解微分方程周期解的计算方法，并通过数值计算结果来总结出有关周 期解的求解方法的一些重要性质。另一方面，在理论上根据以往解的 存在唯一性的理论研究，能够对某种非线性微分方程解的存在唯一性 结果进行讨论。 k 1 2 国内外研究现状分析 一直以来，国内外很多专家学者从事关于微分方程周期解方面的 研究，并取得了大量的骄人的成果，尤其对特殊微分系统的周期解的 存在性、唯一性、稳定性等解的性质的研究成果更为突出嘶】。但是， 由于周期问题本身的复杂性，给研究带来了较大的困难。因此，对周 期解性质的研究仍然存在诸多的问题尚待解决还有很多地方需要进 一步的完善，进一步的放宽条件，以求使用判断更加容易。另一方面， 对于求解微分方程周期边值问题数值算法的研究相对来说较少，大都 针对实际问题的某种情况，也有经验数据的加入，使得求解在算法和 程序实现上都缺乏通用性。沈祖和教授和李维国教授这两位学者长期 2 中国石油大学( 华东) 硕士论文第1 章前言 从事这些方面的研究并取得了显著的成果，其中的思想将对本课题的 研究有着重要的启发。 1 3 主要研究内容和目标 本课题研究的主要内容，大体可以分为两个方面。一方面，研究 常微分方程周期边值问题的数值求解，努力寻求更为有效的求得周期 解的数值计算方法；另一方面，在理论上针对一类偏微分方程解的存 在唯一性进行研究。 众所周知，数学上求微分方程周期解和周期的问题可归结为求解 下面的数学问题：寻找 ( 嘎r ) ，使得! 挚= 4 ( ，) ) 成立，其中“( f ) l 甜( f ) = u ( t + z ) 是r ”中的实值向量函数，是时间变量，a ( o ) 是一个非线性算子，r 是 最小正周期。对于一般的周期问题，求解的难点在于周期r 是未知的， 这种情况下主要利用已有的数值计算方法在求解周期解u ( t ) 的同时也 求出周期r ，但往往由于周期r 未知，而使得未知数的个数与方程的 个数不同，就需要在计算之前对某一个特定的未知数给予初值估计， 这样准确的估计一个未知数的初值就给求解带来了较大的困难；对于 较为特殊的固定周期问题，有时会由于步长和迭代次数之间的矛盾等 原因导致求解时可能出现数值不稳定的情况，针对这种数值不稳定的 情况本课题将结合刚性微分方程初值问题的求解思想来解决。这种思 想主要体现在将刚性微分方程的稳定有效的求解方法【7 j ( 例如：g e a r 方法、隐式r u n g e k u t t a 方法等等) 用来求解出现刚性现象的周期边 中国石油大学( 华东) 硕士论文第1 章前言 值问题，以达到快速准确求解的目的。近年来，刚性微分方程初值问 题的数值解法得到了蓬勃的发展嘲，这给本课题提供了较为有力的条 件。在研究过程中，首先要提到数值软件的使用。本课题将使用多种 有效的策略形成不同的数值算法，并利用大型的数值计算软件 m a t l a b 结合大量的实例程序来实现。对于本课题使用的数值计算方 法，主要研究的微分方程周期边值问题的模型为： 警刊c 删， i “( r ) = u ( t + r ) 其中l ，( f ) 是秽中的实值向量函数，是时间变量，一( ) 是一个非线性 算予，r 是最小正周期。 在理论方面，本文将讨论在共振的条件下，关于半线性椭圆型偏 微分方程广义解的存在唯一性，并在文中给出相应的定理及其证明。 证明的出发点不同于以往的一些存在性结果，不是由古典的 l e r a y s c h a u d e r 度理论确立的 9 1 ，而是利用了一个全局反函数定理和极 小极大原理的一个新的非变分形式得到的。关于解的存在唯一性的探 讨，主要研究的模型为： f a u + g ( x ，u ) = ( x ) i n q 1 = 0 d 疗弛 其中x q c r “，h e 上2 ( q ) ，g ：q r r 是一个连续函数且关于“ 是连续可微的。 本课题具体实现的内容归结如下： ( 1 ) 使用初值问题方法和边值问题方法对一般的周期边值问题 进行求解，包括周期已知和周期未知的情况，并通过数值 4 中国石油大学( 华东) 硕士论文第1 章前言 实验结果显示出这些方法的运算效果。 ( 2 ) 对奇异周期边值问题采用b v p 4 c 和s b v p 相应的有效策略 进行数值求解，同时举例说明算法的实用性。 ( 3 ) 对一般的周期边值问题，当使用通常的求解方法出现计算 不稳定现象时，结合微分方程刚性初值问题的稳定的求解 方法来解决，以使得求解的效果更好，并通过数值实例表 明方法的有效性。 ( 4 ) 对半线性椭圆型偏微分方程的d i r i c h l e t 边值问题广义解的 存在唯一性，给出相应的判定定理，并对定理给予严格的 推导证明。 本文将按照上面叙述的具体实现的内容和目标，来详细地论述得 到的相应的研究结果。 本文具体内容共分为六章。第l 章为本文的前言部分，阐述了论 文课题的主要研究内容、耳标和国内外研究现状，简单的介绍了论文 的研究背景；第2 章至第4 章是针对不同类型的微分方程周期边值问 题在数值算法上的研究和具体实现；第5 章则是在微分方程解的存在 唯一性理论上的研究结果；第6 章是全文的结论部分，是对论文的整 体内容进行的概括和总结。 5 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 第2 章一般周期边值问题的求解 本章主要讨论求解一般周期边值问题的数值解法，包括初值问题 方法和边值问题方法。在周期已知和周期未知的情况下分别进行求解， 并分析对比数值计算结果。计算数据表明，初值问题方法和边值问题 方法求解常微分方程周期边值问题是各有优势的，具体问题可以采用 具体方法，以求达到程序简洁、工作量小、准确度高的目的。 2 1 问题的提出 常微分方程是描述连续变化现象的数学模型，在数学、自然科学 和工程学中随处可见。常微分方程系统本身有很多解，一般我们感兴 趣的是所有分量在单一点x = a 都有指定值的解，这就是常微分方程的 初值问题，初值问题的求解有很多有效的途径。而对于常微分方程边 值问题指定的却不仅仅是一个石的值，这就使得边值问题不象初值问 题那样，它可以没有解，或者有有限多个解，或者有无限多个解。正 因为如此，求解边值问题的程序需要事先给出一个解的初始猜测。在 实际求解过程中，还经常地需要确定一些参数使得边值问题有解。 一类重要的边值问题就是周期边值问题，它的边界条件是不分离 的。因此，它比初值问题和其它边值问题的求解还要困难复杂。即使 是在解和未知参数有很好的猜测的情况下，也可能对求解结果很失望。 本章主要讨论的就是常微分方程周期边值问题两种有效的求解方法， 并举例说明、使用m a t l a b 程序计算数值实验结果、绘出相应的解 曲线。 6 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 考察下列厅阶常微分方程周期边值问题 j ，“( f ) + 厂( f ，y ，y ，y ”1 ) = o( 2 ，1 ) l y o ) ( = y o ) ? ) ，f i o ，l ，。n - 1 其中，l 珂e n ，是连续的非线性函数，周期丁可以已知，也可以 未知。 2 2 边值问题方法 在m a i t , a b 软件中，b v p 4 e 是求解边值问题的很有效的算法。但 是，这个算法和计算环境不适合对要求高精度和解内部变化剧烈的情 况求解。尽管b v p a e 可以接受相当一般的边值问题( 2 1 ) ，但是在非 常复杂酌模型中，仍然存在一些闯题，为了求解可能需要颈处理。在 求解周期边值问题的过程中，还包括对模型解的存在唯一性的探测， 这与求解初值问题大不相同。众所周知。打靶法可以用来求鲤边值问 题，但b v p 4 c 不是由打靶法的代码组成，而是由配置法实现的求解过 程【1 13 1 。 配置法的一般思想就是，将所求问题的近似解写成 ，。o ) = c 。仍( f ) 的形式，其中织( f ) ，f = 1 , - - - , m 是基函数， f i c 。，f = l ，m 是待定系数。然后将这个近似形式代入问题的原方程， 这样就得到了关于未知数e ，i = i ，m 的( 非) 线性方程组。这样所 求问题就转化为求解埘维( 非) 线性方程组。常用的基函数有拉格朗 日多项式、勒让德多项式等等。在m a t l a b 中，b v p 4 c 就是执行配置 7 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 法来求解问题的。假设s ( t ) 是边值问题的近似解，在区间【o ，t 】上是连 续的，并且在每个子区间【f 。，t 。】上是三次多项式，这里 0 = t o ，l ，= t 。那么，s ( t ) 在子区间的端点和中点满足 f s ( f 。) = 厂( ，s ( t d ) s ( ，。+ i ) = f ( t 。+ ，s ( ) ) 【s ( ( + f 。+ i ) 2 ) = 厂( ( f 。+ + 1 ) 2 ，s ( ( + + i ) 2 ) ) 这样就可以通过求解非线性方程组来求出每个子区间以，+ 。】上多项 式的系数，继而求出原问题的近似解。在b v p 4 c 中使用配置法求解， 能够很好的控制残量r ( t ) = s ( f 。) 一，( ，s ( ) ) 的大小，并且对不是很 好的初始猜测，程序还可以自动的稍作调节，以使得求解的近似效果 更好。 由于非自治常微分方程可以通过增加阶数而转换为自治常微分方 程f 。因此，我们只考察自治的常微分方程周期边值问题 妒( f ) + m 办，y 扣。) = 0( 2 2 ) 【y o ) ( o ) = y o ) ( r ) ，f = o ，l ，雄一l 其中，1 r n ，f 是连续的非线性函数，周期丁可以已知，也可以 未知。 2 2 1 周期已知的情况 当周期r 已知时，利用b v p 4 c 函数求解问题( 2 2 ) ，程序是十分简 洁的，只是需要构造几个简单的函数。基本的求解语句： s o l i n i t = b v p i n i t ( x i n i t ，y i n i t ) 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 s 0 1 一b v p 4 c ( o d e f u n ， b c f u n ，s o l i n i t ，o p t i o n s ) 这样得到s 0 1 y 就是问题( 2 2 ) 的满足周期r 的数值解。为了使用上 面的求解语句，我们要给出x 和y 的初始猜测值，这个初始猜测值要 求不是非常高，譬如x i n i t = l i n s p a c e ( o ，t , 2 0 ) ，y i n i t = s i n ( x ) ，c o s ( x ) 等等。 另外，还要构造两个函数o d e f u n 和b c f u n ，它们分别是高阶常微分方 程写成一阶常微分方程组的形式和方程满足的边界条件。 例2 1 使用边值问题方法，求解二阶常微分方程周期边值问题 j r + l o g ( x + 瓜) ：s i n ( t ) 4 工( 0 ) = x ( 2 x ) 石( o ) = z 7 ( 2 石) 这是一个非线性常微分方程，根据文 1 4 1 中定理4 3 知，此周期边 值问题有唯一解。求得解曲线和极限环如图2 1 所示，计算时间和误 差见表2 1 。 _ l 町期边值问题的数值求解 图2 1 ( a ) 9 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 图2 - 1 ( b ) 例2 2 使用边值问题方法，求解四阶常微分方程周期边值问题 x 4 + x 气一1 6 x x 一1 5 s i n ( 2 t ) = 0 x ( o ) = 工( 石) x ( o ) = x ( 万) 工”( o ) = x 。( 石) 工。( o ) = x 一( 石) 可以看到，这也是一个非线性常微分方程周期边值问题，能够计 算出一个近似的数值周期解。通过使用边值问题方法进行计算，求得 的解曲线和极限环分别如图2 - 2 ( a ) 和图2 - 2 所示，计算时间和误差见 表2 1 。 1 0 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 2 j 町媚边值阿题的数值求解 周勰解t = p l 图2 - 2 ( a ) 图2 - 2 ( b ) l i 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 例2 1例2 2 计算时间 1 0 0 1 0 0 0 0 0 0 0 0 0 0 00 9 4 2 0 0 0 0 0 0 0 0 0 0 0 计算误差阶 l e 1 3 l e 6 2 2 2 周期未知的情况 表2 1 当周期r 是未知参数时，求解变得复杂起来，且对初始参数的猜 测较为敏感。为了便于求解，我们对问题( 2 2 ) 作变换，令r = 去。 这样原问题( 2 2 ) 就变为( 2 3 ) 的形式。 胪( f ) + 巩m ，”- , y 扣q ) = o( 2 3 ) i y ( o ( o ) = y o ) ( 1 ) ，f = o ，l ，疗一l 其中，l 疗n ，f 是连续的非线性函数，r 是未知周期。 求解的基本语句： s o l i n i t = b v p i n i t ( x i n i t ，y i n i t , p a r a m s ) s o l = b v p 4 e ( o d e f t m ， b e f u n , s o l i n i t ，o p t i o n s ) 与周期丁已知的情况相比增加了参数项p a r a m s ，这个参数是周期r 的 初始猜测。 例2 - 3 使用边值问题方法，求解二阶常微分方程周期边值问题 f y ：= 3 ( 乃+ j ，：一；y ? 一t - 3 y ；= 一( y i 一0 7 + 0 8 y 2 ) 3 。 ij ，i ( o ) = y l ( t ) 【y 2 ( o ) = y ：( r ) 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 j 时期边值四题的歙值求解 周期解t = l f l7 1 图2 - 3 ( a ) 图2 - 3 ( b ) 例2 4 使用边值问题方法，求解三阶常微分方程周期边值问题 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 f y 。+ 。+ y - ( y ) 2 + l = 0 i y ( o ) = y ( r ) i y ( o ) = y ( 7 ) 【y 。( o ) = y ( 丁) 周期边值阿题的歙值求解 周期解t 叼p i 图2 - 4 ( a ) 极限环 图2 - 4 0 ) 1 4 ， 们 o o ， 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 例2 3 例2 4 计算时间 4 6 9 7 0 0 0 0 0 0 0 0 0 0 0o 9 2 2 0 0 0 0 0 0 0 0 0 0 0 计算误差阶 l e - 1 3l e 3 周期t l o 7 1 0 8 0 8 5 4 6 6 2 4 1 06 2 8 318 5 3 2 3 6 3 3 6 0 2 3 初值问题方法 表2 - 2 初值问题方法求解周期边值问题的主要思想，就是应用类似打靶 法的基本原理，来获得周期边值问题的数值解 1 41 5 1 。 我们要求解的是下面甩阶常微分方程周期边值问题 其中，1 力n ，是连续的非线性函数，丁是已知周期。 ( 2 4 ) 首先，我们假设口为玎维向量，r a = y ( o ( 0 ) ，f = 0 , 1 ，甩一1 。 由于，( r ，y ，y ，y 洲) 是连续的，则我们可以得到初值问题 j y ! ：( f ) + ，( r ，y ，y ，j ，枷1 ) = o ( 2 5 ) i ) ，”( o ) = 口l + i ，f = 0 , 1 ，摊一1 的解，记作：y ( t ，口) 。 我们取门维向量占( 口) ，其中g “1 ) = y o ) ( 7 ，口) ，f = 0 , 1 ，雄一1 。为 了满足原方程( 2 4 ) 的周期边值条件，y ( ，口) i = 0 , 1 ，一1 在f = 0 时和，= t 时值是应该相等的，即y o ) ( 0 ，口) = y m ( r ，口) ，i = 0 , 1 ，万一i 。 o l = 一 ) 厅 矗 ， v = ， y x弘口 以m 一， 十 = 时 o y y ，、，【 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 这样，就将周期边值问题转换为求解不动点 g ) = o f ( 2 6 ) 的问题。实际上，就是求解非线性方程组 g ( 口) ；g ( 口) - a = o ( 2 7 ) 于是，我们就可以将常微分方程周期边值问题( 2 4 ) 的求解转化为求 解非线性方程组( 2 7 ) 。 求解非线性方程组的方法很多，这里我们采用秩l 拟n e w ( o n 法。 令a ；d a c o b i ( g ) ，迭代格式为 口川= 一4 1 g k a a i2 口i “一口i a g , = g k “一g k 她= 盟掣 a k “= a k + 削t ( 2 8 ) 这样得到的解口就是问题( 2 4 ) 的初始近似值。 算法如下： s t e p0 给定初始向量a r “，一个小正数o 万 1 ； s t e p1 计算初始矩阵a o = d a c o b i ( g ) ( a ) 取定一个扰动矩阵艿r 且满足4 一= 悸? ；再确定 一组向量a ，户：，鼠，满足屈= 瓯+ 风，其中瓯是矩阵万 的第k 列向量； 1 6 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 ( b ) 使用o d e 4 5 计算对应初始向量风，屈，展，展的解并取其 在t = t 处的向量值记为9 0 ，f - o ，l ，h ； ( c ) 计算a 。= j a c o b i ( g ) = g - i ，使用差商计算导数g ， gi，：；l：蕊90j(i)-900(i)：下go,(i)-goo(i)，f，：l，疗： 6 e 母ip | u ) 一p o u ) 否 ! ”j 1 1 s t e p2 计算非线性方程组的解口。 ( a ) 取= 岛，g o = g o o - i t o ，一个小正数0 占 l ，置 k ：= 0 ： ( b ) 计算( 2 8 ) 中第一个式子，得到吼+ 。作为初值问题的初始 近似值，再使用o d e 4 5 计算( 2 5 ) 得到解y ( a k + 。) ，向量 g 川= r ( r ，吼+ i ) ，由( 2 7 ) 可得到g k + l ，再计算( 2 8 ) 其 它各式得到a a k ，a g k ，她，a ，并置k = k + 1 ，直到 慨+ ，忙； s t e p3由s t e p 2 计算得到的解口就是( 2 4 ) 的近似初始值，则 使用o d e 4 5 求解初值问题( 2 5 ) ，就可以得到原问题( 2 4 ) 的数 值解y ( ，) 。 例2 - 5 对例2 - 1 用初值问题方法进行数值求解。 采用初值问题方法进行数值求解，具体实现的数值结果见表2 - 3 ， 使用m a t l a b 软件绘制的相应的解曲线和极限环分别见图2 - 5 ( a ) 年1 1 2 5 ( b ) 。 周期边值河题的数值求解 周期解t - - 2 p i 图2 4 ( a ) 图2 5 ( b ) 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 例2 - 6 对例2 2 用初值问题方法求解。 数值计算结果见表2 3 ，相应的解曲线和极限环见图2 - 6 。 周期边值简酝的叛值求解 周期解t = p l 图2 - 6 ( a ) 图2 - 6 ( b ) 1 9 中国石油大学( 华东) 硕士论文第2 章一般周期边值问题的求解 例2 5例2 6 计算时间o 3 8 0 9 9 9 9 9 9 9 9 9 9 9 0 2 8 0 0 0 0 0 0 0 0 0 0 0 0 + 计算误差阶 l e 1 5l e 1 3 迭代次数k = 2 6k = 1 4 2 4 小结 表2 3 从上面的数值试验我们可以看到，这两种方法求解常微分方程周 期边值问题( 2 1 ) 的效果都是很好的，它们各自都有自己的特点。 ( 1 )边值问题方法可以使用参数代入来解决因周期未知而出 现的未知数比方程个数多的复杂情况，这也是初值问题 方法所不能解决的。 ( 2 )使用初值问题方法来求解边值问题，这不仅是非常好的 思想，而且也行之有效地解决了困难的周期边界条件的 满足。 ( 3 ) 一般来说，初值问题方法和边值问题方法求解的精度虽 然不是对任何问题都是非常高的，但是我们还是可以接 受的。在初值问题方法的计算过程中，使用了拟n e w t o n 的方法，采用差商代替导数从而大大减少了计算量，节 约了计算时间。 总之，使用边值问题方法和初值问题方法来求解问题( 2 1 ) ，要 结合问题本身和方法的特点灵活运用。 中国石油大学( 华东) 硕士论文第3 章奇异周期边值问题的求解 第3 章奇异周期边值问题的求解 在上一章我们讨论了一般周期边值问题的求解，给出了较为有效 的计算方法。而在实际应用中也常常会出现奇异的周期边值问题，所 以在这一章就将讨论在计算周期边值问题时出现奇异点的情况，并给 出相应求解的有效策略和具体的实现方法。 3 1 问题的提出 奇异问题是数学中常见的一类特殊的问题。在周期边值问题的求 解中，也常会出现。众所周知，周期边值问题本身的特殊性已经给求 解带来了较大的困难，这我们从上一章的论述也可以清楚地看到，再 加上奇异性的存在就使得求解更加困难。因此，我们在求解时针对问 题不同的性态，将结合不同的计算技巧，这样才能够达到有效求解的 目的。我们知道，打靶法以及多级打靶法等是具有很好的收敛性等好 的性质，但是在实现上具体的程序编辑会较为繁琐【1 6 1 。因此，本章在 求解奇异的边值问题过程中，对于程序的具体实现将使用m a t l a b 软件包s b v p 和b v p 4 c 等来进行。对于b v p 4 c 求解边值问题的计算效 果，我们通过上一章的计算实例也已经看到，而在出现奇异点的情况 下b v p 4 c 同样也可以使用，实际计算的效果仍然是很好的。软件包 s b v p 是从b v p 4 e 求解边值问题的结果得到的启发，形成的关于奇异 边值问题的数值计算方法的软件包。根据有关文献的记录，该软件包 通过反复测试，从得到的误差估计表明，其对求解奇异边值问题的效 果是比较好的f l ”s b v p 软件包也是近几年内计算数学工作者最新编 2 l 中国石油大学( 华东) 硕士论文第3 章奇异周期边值问题的求解 写组成的，在很多功能上和b v p 4 e 相比尚需进一步的完善。因此，我 们这里仅限于对周期已知的情况下的奇异周期边值问题进行求解，并 根据奇点所在的具体位置将问题分为内部和边界两种情况进行论述。 考虑如下奇异边值问题 艘n ，i 厂蟹p ) ，旭( 毋6 ) ， ( 3 1 ) i r ( y ( 口) ，y ( 6 ) ) = 0 这里微分方程的右端是含有奇异项的，即 f ( t 9 ) ) = 高m ( f ) “r ) + g ( t ( f ) ) ， ( 3 2 这里肘是一个矩阵，g 是光滑的向量值函数，a j ，s b 是问题的奇 异点。容易看到，( 3 1 ) 中的边界条件是具有一般性的，一当 月( y ( 口) ，y ( 6 ) ) = j ，( 6 ) 一y ( a ) = 0 就形成了以t = b a 为周期的边值问 题，这时问题( 3 1 ) 就是一个奇异的周期边值问题。 3 2 内部奇点的情况 在这一节中，我们将结合两点边值问题在内部出现奇点时采用的 策略，来讨论周期边值问题出现奇异情形时的处理方法。对于一般的 内部奇点的边值问题，有很多求解的策略，譬如幂级数方法等等【1 8 1 。9 1 但是，幂级数方法在使用上没有b y t e 简洁，而在求解的效果上b v p 4 e 也比较好。因此，使用b v p 4 e 是很便宜的，在使用之前我们需要在方 程函数的编写中进行一些预处理。例如，当解曲线光滑时，根据条件 ，、 l i m 兰1 2 = y 。( o ) ，我们可以将在奇异点t = 0 处的解析式直接定义；当 t - - + o f 解曲线是间断的，我t f 贝j j 可以采用左或右导数的存在来直接写出奇异 中国石油大学( 华东) 硕士论文第3 章奇异周期边值问题的求解 点的解析表达式，其它的实现语句都与第2 章的使用完全一样。这样 就使得原问题的奇异性被有效地解除了，同时也不会影响计算的准确 性。对于直接使用s b v p 我们需要编辑一个函数文件，这个函数文件 主要包括写成一阶微分方程组，及其j a c o b i n 的形式、边界条件也同 b v p 4 c 样写成残量的形式、边界条件的j a c o b i n 矩阵以及可以规定的 误差。对于这些使用方法，在下面的数值实例中将更详细的说明。 例3 1 求解一般奇异边值问题jj ，+ 手y = 一半，一1 ， l 。 【y ( 一1 ) = 一i , y ( d = l 可以看到，t = 0 是方程的内部奇点。采用直接给出奇点解析式的 办法，体现在方程函数奇点处表达式的写法上，通过下面使用的语句 可以清楚的看到。 基本的使用语句为：s o l = b v l m c ( o d e f u n ， b c f u n , s o l i n i t ，o p t i o n s ) ； 预处理的方程函数具体写为：f u n c t i o nd y d x = o d e f u n ( x , y ) i f x 一0 d y d x ( 2 ，1 ) = 1 3 ；e l s ed y d x ( 2 ，1 ) - ( 2 + y ( 2 ) x + s i n ( x ) x ) ； 使用b v p 4 c 通过程序实现，求得的数值结果见图3 1 。 图3 - 1 2 3 中国石油大学( 华东) 硕士论文第3 章奇异周期边值问题的求解 肚，+ 击朋+ ( c o s ( x - 1 ) + x - 与s i n ( x - 1 ) ) = o 【y ( o ) = y ( 2 r r ) ，y ( o ) = y ( 2 j r ) 中国石油大学( 华东) 硕士论文第3 章奇异周期边值问题的求解 o u t = s b v p s e t ( r e t o l , l e - 9 ，a b s t o l , l e 一9 ，r i ) ； o t h e r w i s e e r r o “i l i l k l l o w nf l a g ) ； e n d 求解的数值结果见图3 - 2 。 奇异j 对期边值耐磨的数值解 图3 - 2 ( a ) 图3 - 2 ( b 1 中国石油大学( 华东) 硕士论文第3 章奇异周期边值问题的求解 3 3 边界奇点的情况 对于边界奇点的情况，使用不同的工具在处理上也稍有不同。使 用b v p 4 c 解题器，由于它本身含有奇异选项s i n g u l a r t e r m ，并且这 个计算功能也是对于边界奇点的情况适用的，所以求解的程序语句也 是和b v p 4 c 其它的求解大致相同，只是添加奇异选项而己。另外，在 使用b v p 4 c 时，还可以在奇异点处添加微小的扰动，如1 0 - 5 等等。这 样可以将奇异的情况减弱，求解就象第2 章一般周期边值问题那样容 易。而直接采用s b v p 进行计算，是需要预处理的，对于边界奇异的 周期边值问题而言，需要将奇点的位寅平移到区间的内部而且不能是 特殊的节点( 例如，和奇点数值相同的t 2 点等等) ，其余的编辑都是 和内部奇点情况相同。下面就分别对一般奇异边值问题和奇异周期边 值问题分别举例说明算法的可行性和有效性。 例3 3 在b v p 4 c 中选择奇异项s i n g u l a r t e r m ，来求解奇异两点边值 r1 问题j 少+ 5 i y ( o ) = 0 ，y ( o = 4 3 4 同计算一般边值问题的使用一样，b v p 4 c 的使用是简洁有效的。 与第2 章中的基本语句有着不同的地方就是在o p t i o n s 选项中添加了 s i n g u l a r t e r m 项，求解的基本语句为： s o l i n i t = b v p i n i t ( x i n i t ，y i n i t ) o p t i o n s = b v p s e t ( ，s i n g u l a r t c r m ，s ) s o l = b v p 4 c ( o d e f u n ， b e f u n ，s o l i n i t ，o p t i o n s ) 其中的s 是奇异项的系数矩阵，要在使用前将数值赋上。另外，要注 中国石油大学( 华东) 硕士论文第3 章奇异周期边值问题的求解 意的就是o d e f u n 的编辑要除去奇异项的一阶常微分方程组的形式。由 上面的处理方法，计算得到的数值结果见图3 - 3 。 图3 3 例3 - 4 用添加微小扰动的方法，求解奇异周期边值问题【2 0 , 2 1 1 暑 j ，c x ，+ y c x ，+ ( c 。s x + 三s i n x = 。 【y ( o ) = y ( 2 u ) 这里边界奇点x = 0 ，采用b v p 4 c 添加微小扰动万= l e 一8 的处理方 法，基本的求解语句与前面一般周期边值问题的相同，只是编辑微分 方程的函数时，将奇异项分母上的自变量进行了扰动，即 x = x + 万= x + l e 一8 。 函数编写的语句为： f u n c t i o nd y d x ：= = o d e f u n ( x ，y ) d y d x ( 2 ，1 ) - , ( y ( 2 ) ( x + l e - 8 ) + c o s ( x ) + s i n ( x ) ( x + l e - 8 ) ) ； 添加这样微小的扰动，在误差允许的范围内是不影响计算的数值结果 中国石油大学( 华东) 硕士论文第3 章奇异周期边值问题的求解 的。绘制的解曲线见图3 4 ，我们可以从图形中看到，计算的实际效 果是很理想的。 图3 - 4 ( a ) 图3 - 4 ( b ) 中国石油大学( 华东) 硕士论文 第3 章奇异周期边值问题的求解 例3 - 5 采用s b v p 求解周期边值问题例3 - 4 。 由于问题是边界奇点的情况以及周期边界条件的特殊性，我们需 要将奇点的位置平移到计算区阈的内部，根据上蠢的叙述，我们要做 一些预处理，将此例中的解曲线向区间内部平移1 个单位，即将 工= x 一1 代入即可。这样，原问题的方程就变为： 卜，+ 击y + ( e o s ( x - 1 ) + 击s i n ( x - 1 ) 一- ， i 岁( o ) = y ( z x ) ，) ，( o ) = y ( 2 芹) 而边界奇点t = o 就变成了内部奇点0 s x s = l 2 x 。求解的具体情 况，可以通过内部奇点的例3 - 2 看到，故不再赘述。 3 。4 小结 对一般的奇异边值问题，内部奇点和边界奇点可以根据使用工具 本身的要求选择适用的方法，而对于奇异周期边值问题，由于在边界 出现奇点的情况加上周期边界条件的特殊性，处理上不能够直接进行， 所以要将奇点的位置平移到计算区间的内部，也就是将边界奇点转化 为内部奇点，同时由于解的周期性使得平移仍然会保持解的真实性和 准确性。这些都是可以从上面叙述的实际的数值例子当中看到，这些 处理的策略和方法得到的结果显示出，它们是有效可行的。当然，对 于周期未知的情况，较为复杂，因为在计算的过程中会需要代入未知 的参数，使用上面所述的方法计算结果还存在一定的问题，还需要进 一步的试验和探讨。 中国石油大学( 华东) 硕士论文第4 章计算周期解时出现刚性问题的求解 第4 章计算周期解时出现刚性问题的求解 从前面章节的