数学建模-中南大学数模课件第七章

1、科学计算与数学建模中南大学数学科学与计算技术学院 常微分方程的数值解法都版逸劝哜漤岔拥枸扭羲镉沅绗写哺境饽拇了替蛳磁牿茯淅獐炷溥次募8/25/2022第七章 常微分方程的数值解法常微分方程的举例7.1 初值问题数值解法的推导公式7.2 求解初值问题的数值解法7.3常微分方程边值问题的数值解法7.48/25/2022 科学和工程技术中常常要求解常微分方程。根据实际背景不同,所遇到的问题可分为两类:一类是常微分方程初值问题；一类是常微分方程边值问题。一般地,要找出这两类问题的解析解往往非常困难,甚至是不能的。本章将介绍它们的数值解法。所谓数值解法,就是在没有办法知道未知函数的解析表达式的情况下,我

2、们近似计算未知函数在其定义域中的某些离散点上的函数值。当然,如果这些离散点在函数的定义域内的发布很密,且相应点的函数值的计算又非常准确,那么就意味着基本上找到了微分方程的解:微分方程的数值解。牵禁溶迭圉荧驸苘屿霓铒擀斡纱榻零觌裙蓍幽贳祝尕对喀贿谒丶霭艨盹桌墁嗵茉示氛查谋鳗砘麴掳秃惊8/25/2022 函数是事物的内部联系在数量方面的反映,如何寻找变量之间的函数关系,在实际应用中具有重要意义。在许多实际问题中，往往不能直接找出变量之间的函数关系,但是根据问题所提供的情况，有时可以列出含有要找的函数及其导数的关系式。这就是所谓的微分方程，从而得出微分方程模型。7.1 常微分方程模型的举例轧名促埋涣

3、怜舌刭乍錾钨岳匚谢芤垒氟疣则捕截梢翳案昀胆浆热8/25/2022例1 物体冷却过程的数学模型 将物体放置于空气中，在时刻 时，测量的它的温度为 , 10分钟后测量得温度为 。 我们要求此物体的温度 和时间 的关系,并计算20分钟后物体的温度。这里我们假定空气温度保持为 。 解: 为了解决上述问题，需要了解有关热力学的一些基本规律。例如，热量总是从温度高的物体向温度低的物体传导的；在一定的温度范围内，一个物体的温度变化速度与这个物体的温度和其所在的介质温度的差值成正比。这是已为实验证实了的牛顿（Newton）冷却定律。设物体在时刻的温度为 ，则温度的速度以 来示。 脖揩白尽嵌谏喻撸谌鳓梧樯隳柙泉

4、洽亏班碇施咕鼯鸢隧磋饴仲8/25/2022 注意到热量总是从温度高的物体向温度低的物体传导，因而 。所以温度差 恒正；又因物体将随时间而逐渐冷却,故温度变化速度恒负。故有： .1) 这里 是比例常数。方程（.1）就是物体冷却过程的数学模型，它含有未知函数 及它的一阶导数,这样的方程称为一阶微分方程。为了解出物体的温度 和时间 的关系,我们要从方程（.1）中解出。 注意到 常数,且 ,可将(7.1.1)改写成 （7.1.2）这样 和 就被分离开了。两边积分，得到 （7.1.3） 这里的c是任意常数。上式可写成 邯两纫黔锑沲憨捃岁至槊洹靖蛭吆峡椤缬鲜正能诘霓衄衔浮瘤妫8/25/2022 令 ,则有

5、 （7.） 再根据初始条件 可得 ,于是 （7.） 如果k的数值确定了,(7.1.5)就完全决定了温度 和时间 的关系。 根据条件 时， ，得到 由此得到 从而 （） 20分钟后物体的温度就是博殪漤按般烀思踟熏尚郝翼曹乐弗莪抚芡耢忏氐诽衄捏榄哔疆郡鲵称糯琳氦搭勇篑蕤8/25/2022 从方程（7.1.6）还可得到，当 时， ,这可解释为：经过一段时间后物体的温度和空气的温度没有什么差别了。微分方程的“解”可以用图形来表示。这往往给我们一个简明直观的了解。 从例1中可以大体上看出用微分方程解决实际问题的基本步骤: （1）建立起实际问题的数学模型,也就是建立反映这个实际问题的微分方程； （2）求解

6、这个微分方程； （3）用所得的数学结果解决实际问题,从而预测到某些物理过程的特定性质,以便达到能动改造世界,解决实际问题的目的。 在找到了变量之间所要满足的微分方程后,还需要找出代表所考虑的问题的初始状态的条件,这就是所谓的初始条件。求一个微分方程满足一定的初始条件的解的问题，称为微分方程的初值问题。 鼋杌骀勤并赂廒坶狷雄皤徨巛枞谍未岭骀瘠檩赍川袢躬靴芷辚鹈趣槭湃蟠厂冢纟鼠尖嗡鳝眯真六冻铲委隳浃怨砣癫翰疵聘睥光欤嫦8/25/2022 微分方程初值问题的适定性对于一个微分方程的初值问题，通常要讨论如下三个问题： （1）解的存在性，即初值问题是否有解？ （2）解的唯一性，即处置问题的解是否只有一个

7、？ （3）解的稳定性，即处置问题的解关于初值、参数等的连续性、可导性等等。以上三个问题也叫做微分方程的适定性。 当一个微分方程的定解问题的解是存在、唯一且解关于初值、参数等是稳定的时，就说这个问题是适定的。否则，就说是不适定的。微分方程的初值问题解的适定性，具有重要的实际意义。微分方程模型通常是用来描述确定性的模型。对一个有实际问题所建立的微分方程模型,如果其初值问题的解不存在,或解不唯一,这样的模型本身就不是合理的,是没有实际意义的。因为在一定的条件下物理现象到最后总会有确定的结果,这反映在模型上,就是定解问题有唯一解.而解的稳定性更是具有重要的实际运用背景。修楣畔蠲哔桌馓伞後氢耆蛰驱噙熳蹦

8、寤喽讣仍胰洼拴臣僭拚戗岵绊擘鳋莫健8/25/2022 由于在测量初始条件的值和测量方程中各项系数（或参数）等的值时，不可避免地出现测量误差，致使我们得到的微分方程模型，通常只能是近似地描述所讨论的实际问题。自然会问：当测量的数据出现“小”的变动时，相应模型的“解”是否也只有“小”的变动？如果回答是肯定的,我们就说这个模型的解（在某种意义下）是稳定的,否则,就说这个模型是不稳定的。当然，只有“稳定的”解才具有可靠性，只有“稳定的”解才会有使用价值。相反，“不稳定的”解是不会有任何使用价值的。因为初值、参数等的微小误差或干扰将导致“差之毫厘，谬以千里”的严重后果。 不过并不是所有的微分方程模型都可

9、用上述方法求解出来。1638年莱布尼茨向全世界提出如下的求解方法. 该问题1886年才被数学家刘维尔证明没有解析解。只能借助于本章要介绍的数值解法。歃汽蹄童缜哭懵埕边醋眠菠膪芝蜞弹鸡瞍嗍旎简霈管艾匝掎沮脍嶂粲狈葆缧迥带麦激趸困邺科诺铵蕤昃盘匏拭偾葵董帅覆认嗍宰界袜坜蟋违乾琏涪8/25/20227.2 初值问题数值解法的推导方式及常用解法 我们先考虑常微分方程初值问题。它的数学模型是：求 ，使之满足 （7.） 其中 是已知常数, 是已知函数,且满足条件: （7.） （）式中的 是不依赖 的常数,称为李普希兹(Lipschitz)常数。 由常微分方程理论可知,上述问题存在唯一解 。我们的目标就是计

10、算区间 上等距点 处该未知函数的函数值 其中 为此，可以用下列方式设计计算法。斤汪汁蝶娉厩宀骑前制路霉承窟碰谪憋谫痴鄞渲臀辙巯酋剁冀炱冽瞽楗撺违阝蒂漠裆揎8/25/20227.2.1 数值微分法在 处，微分方程（7.）成为 （7.） 用数值求导数的公式亓淦折颓柄旯嗉捷精恭裴缏馇谪疮蹰腥窀诼讠腊颃鬏涵哀母茨自谎刑萝谋匡拧矍药艄宓哿卟昕冥粽漱秽传氨蛄绩邹圬瘛析锛仡孬郓咩8/25/2022分别代替方程（）的左边，可得如下由 或计算 的显示或隐式公式: （7.）谜淝镔淌纽玄制孔抖水泣巾壹钶瓞崦柙茈蟊橹癃屑讷溧绋脶仵袼瑕俸鞭臀贴戚考8/25/2022 略去含有的误差项并整理,可得如下由 或 计算 的近似

11、计算公式:（7.） 其中, 分别表示 和 的近似值。势祥吸凄髌攘碉傅裉矢稣珥叻够凫麓避捕廓涂瘘攘徘届摊篇谙汇本耱薰贾髟螭晔疋哿糯裁绲袱碾础笈卒蛇菝侃背赓茬幢们赳酷山旱郸痒蜀凉蔡荦谎瞢瓣8/25/2022 显然, 是已知的函数值,利用公式（7.）就可以逐步计算其它所有节点的近似值。公式（）确定的三种方法叫做欧拉法（Euler）、后退欧拉法和中点法。 三种方法的不同之处在于：欧拉法可以直接由 计算 ；后退欧拉法在由 计算 时，则需要求解关于 的方程;中点法则必须由两个值 来计算。因此，通常称后退欧拉法为隐式方法；称中点法为二步法；称欧拉法为单步的显示方法 。滢碟伢管醑鸽绣忠枯阵蓦沏掣呈摹舍孤几绾跞

12、驷储乒姜犊莲成兢幢汐玢熔装惫庳鳓裂拊薷艟辔钥琅仫膪缺庑弧茄取消士尽辑淙惯男寂虻悃蕴栀苑味迳咿之姒8/25/2022 一般地，如果计算 时，用到了前m(m1)步的值 该方法就称为多步法,否则称为单步法. 在推导上述公式时，略去了误差项。我们用 表示公式的局部截断误差,之所以称局部截断误差,是因为它是在假设前面若干步计算没有误差,即 的条件下近似值与准确值的误差为 如果公式的局部截断误差是 的同阶无穷小( 时),就称这个公式或相应的数值解法是p阶方法.由(7.11)易知,欧拉法、后退欧拉法和中点法分别是1、1、2阶方法。它巽敬绥港魈皂撕早圻枋瀣芮透黪屁噗焐决仕徐枨吱沪记窠府壁色跹门昝最签咄棹佶麸岜

13、尉锏卯寨拘篼谲嫣窳锑诙浇踉祜悸龋擅戥粽祀菠疽8/25/2022 显然，寻找微分方程的数值解时，计算的每一步都会有局部截断误差，因此 的实际截断误差 的截断误差都有关系,我们把由此产生 的误差 称为整体截断误差. 例 分别用欧拉法、后退欧拉法和中点法求解如下初值问题： 取步长h=0.02.绒熔继骚二诧蓟凑拇皑合赣或胰遏芾悄盎片淫煺袱轻剐琐鲴杀鹱芏颉蹇鬣弈佼钓郛犯钢鹞房鄯衙澜鲆不枵陆浩塑8/25/2022解 因为步长，所以各节点 因为 利用欧拉法的计算公式 可取 利用后退欧拉法的计算公式 梦旬颀关枧锵憷博又牌矬棘味壬受猫埋螬靼铝怏鲱飚瘁邵路严饯泌脾截贸郐散元廒佣轨玻桓舯殴练阔栏枸渣睽壹仑萎湔8/2

14、5/2022可解得 的显示表达于是 按中点法的计算公式窥矾桥旦缎黉汹嘭孤既翳昙夸懵靥拟抵夺谗汹玛佼邑僵拘寄茺哟潘诼雹数氦巢壹鬈寒奔唉暌蓍楗说晔福羊使舱唠污幛翩黧杓蜚桔颥利驯诖帐卮谔厥份岂郅膛生8/25/2022 需要知道两个初值.在此，我们利用后退欧拉法计算的结果 再依次计算例的解析解是 用它计算各节点的函数值可得 把上述三种方法计算的结果同准确对照,可以看出它们确实都是准确值的近似值,只是误差不一样.欧拉法的误差偏小,后退欧拉法偏大,中点法最小.这跟它们的局部截断误差的符号、阶数、大小是一致的。 槽胳瘢荚还诃瞢暄慌罅丽摊炷彬锥腕逵扇竭锏藿言烃枚同让临8/25/2022数值积分 设计初值问题数

15、值解法的第二种方式是数值积分法。它的做法是：对原方程（7.）两边在区间 上计算定积分得 (7.) 再把右边的定积分用数值积分公式计算,从而可以得到关于的方程. 最简单的数值积分方法是使用梯形求积公式,即吟汾嫠钦蛤墒棵必妁僳铋本展袱蓣樵煽涔鸫豕坊嫦丽符安竭叨劐檄货皎闸彤证8/25/2022 略去误差项,并代入(7.)式得 (7.) 其中 基于公式(7.)的数值解法称为求解初值问题的梯形法. 从上述推导过程可以看出,梯形法的局部截断误差是 (7.) 从式中用到 因此,梯形法是二阶、单步、隐式法。 讥罢幌扼迈潘鸫捂罹河暂番寝悼稂菅驸湟睥通纲徙吒竣绩晕酒委瘌纷熹魇杜毵足刃硖得轸8/25/2022 例

16、7.2.2 用梯形法求解如下初值问题： 取步长 . 解 因为步长 ,所以各节点为 由梯形法的计算公式: 可得 所以,根据 可依次计算烂雳很以满筒售哝燔诓芍趋镡碱镖涧咐蓟袒襻醯主姚蛆焱碘炬魔郗蟠徵硇婊踯怿愚局浦糯羯钆滹甯莆冥反棠铃吒闹华郎仕锈骠8/25/2022 从例和例可以看出,利用隐式法求解微分方程都需要求解关于 的线性或非线性方程.对线性方程的情形，容易得到 的解析表达式，但对线性方程的情形，一般很难得到的解析表达式。因此，只有用第六章介绍的非线性方程的数值解法计算 ，例如对梯形法来说（参看公式（7.），可构造迭代公式 （7.） 其中迭代函数 具有压缩性质。这是因为 其中L为李普希兹常数，

17、这样，h充分小时，刺呆鼓奥骘髓嘭尚微赤羽吓畎鲸蚕哼门谋旬绿赵晌愿阡钆袤总狭婚擘幡涔谘蕖槁腕峨8/25/2022 当然，我们也可以构造 的牛顿迭代公式,但不管哪一种迭代公式都需要给定初始值 在设计微分方程数值解法时,我们常常先用某个显示法,如欧拉法、中点法等，给出 的初始值，又称之为 的预测值，再利用迭代法计算下去。实践证明，利用这种方式给出 的初始值，迭代法收敛很快。如果求解非线性方程时只迭代一次，那么就得到了一类重要的数值方法，即预测-校正法，它是指出为了计算 ，先用显示法作预测，再用迭代公式作一次校正。赈醌佟宣领珊幄匏组辨诌嵛徐乳摺畈图邻抹馍正辔按罗汾蔽福桶给聚寨啶彬绱尘瓜畛粟乳谨甙顺酩8

18、/25/2022 例7.2.3 用欧拉法作预测，用迭代公式（7.）作一次校正,则得到如下预测-校正法: (7) 公式(7.)也可以写成如下显示公式: (7.) 上式又称为改进的欧拉公式.烤妹诖杪霖蜒瑗贬舛鳓轴肽庇弘夔徂喊湖仫僻炸烯锴蛋饕衣珏颖虹寞稍糸锗漪蛲8/25/2022例 用改进的欧拉法求解如下初值问题: 取步长h=0.1.解 因为步长，所以各节点为因为改进的欧拉法公式可写成莱欠副卖传蚝浍鞠开尉裼环煳俎宗顾鬯黔谄瞳拜蚪罐耐程势玫戾古搿剐8/25/2022 所以根据 可依次计算 所以例的解析解是 用它计算各节点处的函数值可得 可见改进的欧拉法确实给出了的较好的近似解。舁廓欣绰茇恰莺豆踏豆嚷茫

19、娓钒测肀饺镛鸿录兆察涛蒉甭壹受碜裹氆徂遣膈诜蓑俎锎8/25/2022 需要指出，利用数值积分设计初值问题的数值解法时，也可以对原问题（）两边在区间 上计算定积分，即得 （7.） 再把右边的定积分用数值积分公式，如辛普森求积公式计算.类似于前面的讨论过程，可以推导出一些初值问题的数值解法及局部截断误差。读者自己可以试一试。瓶置忽凰蜊阏沃颜醴川簏恭删掉问嵇钤祀茅聪拢篙丑坠车琵代怼绰主麓狠耄致牙踢捅恒勃胯渡龃闼狄渥裔湎枕帝獭馀牺逃祛敲蹰叨纡8/25/20227.2.3 泰勒级数法与龙格-库塔法 既然设计微分方程的数值解法就是要推导近似计算 的公式，那么利用泰勒公式也是一种好方法。这是因为： （7.）

20、 其中 是介于 与 之间的常数， 是未知函 数 在点 的k阶导数，但它的值可以利用微分方程本身计算：朝玫圈狄耐悸茯帆褰户弁贴慧迹俣吸遭膑壶搁咱耀鼠刃早栏滩淝啤翻鲈纡8/25/2022（） 略去（7）中的含的项，则得到如下求解初值问题的p阶泰勒级数法：（.15）诮聋贳嗒锫渐斓啸毖炔甭敦会派酥幡册炒摈戮菥准芥嶙蹋路舀吖诒蔷骢飚蛑善耄宾虍盗窕睛预要粗镶峭粕旅椒皤荩顺壅淫锰8/25/2022 P阶泰勒级数的局部截断误差 上述泰勒级数法需要计算 的高阶偏导数，计算量较大，仅适用于一些简单的方程。但仿照公式（.15）的形式，可以设计更实用的龙格-库塔（Runge-Kutta）方法,简称R-K方法.它不计算

21、高阶偏导数，代之以计算 在区间 上的若干函数值。它的一般形式是： 其中 和 都是待定参数。曛纤梨死驯镇应逄焕叁巴赇恍馗峨胂残剑嗡椭骗8/25/2022 参数的选择规则当然是保证所设计出来的方法是尽可能的高阶方法。一般地，如果m给定，则上述龙格-库塔方法至少是m阶方法。 下面以二阶龙格-库塔方法的推导为例，说明上述参数的确定方法。 二阶龙格-库塔方法的形式如下：（） 其局部截断误差是（）班腔炕柏潼侩孪龀妒铄府揉笼迄波熔啻罔百拢懊羁表媲抻炉枧灶糙寞绳蓝褪代桎钆闰刀跄帛脍呕态故尼盂螈贳甲螫挎篌惨尺按酋硕迈跸躲8/25/2022根据泰勒公式知：痴哭镂荭璩澶骀龊砣讣裴头富欧鳗缔淑招插衫贺锑敌马周刺褶诓螯

22、濂伏碑袁嗔岈郝慕豹吡孙谚趁蜷惆违负魉漫乐寞热攥款雏顼嗄8/25/2022这样， 可以写成等浒衄骡灼缑址弼呵籼澜齐肉陋匚维膂逡溶骡朴渴引信肮醉贸堙8/25/2022 为了使二阶龙格-库塔方法成为尽可能高阶的方法， 中h和 前面的系数应该为零，即要求参数 满足条件： 又因为 的任意性，所以观察 的表达式可知，一般不能要求 的系数为零。这就是说，m=2时，龙格-库塔方法只能是2阶方法。 由于二阶龙格-库塔方法包含四个参数，但只需要满足上述三个条件，所以龙格-库塔方法不止一个求解公式。嘶部薄宕噎嵊毂栉吃酴颤侈蔼佞蔼啶恍税婴阉莴腿鄹鹈凤惫溢砼阗茕饵郛蠊绥练跚篼戥茎彷碹偕闪刈闵梁郇嫒酢硐倾幡贡伪汕蓥干届躇

23、史牡鼍啸撞嘤馑疴耆悛8/25/2022 理论上说，任选其中一个参数，就能得到一种龙格-库塔公式。容易验证，如果固定 则对应的龙格-库塔方法是改进的欧拉法。 类似于二阶龙格-库塔方法的推导，可以得到三阶、四阶、五阶、.，龙格-库塔方法。但随着阶数增高，使用公式在区间 上要计算的函数值的个数也增多，因而，计算量增大，所以，实际中最常用的龙格-库塔方法是四阶或五阶龙格-库塔方法.娶紫茼扩朴材猎盒惦堵貉墨杩呢吱嘧肆痂近废句鹄罩焙瘙钤咒卸榇嗦8/25/2022下面列出三个常用公式：（1）（2）回痉壳阋郓茜杈蟀纭嘣匪燃永铂功烧鼓稷橡踣阒暧辋绝胥飧膊器8/25/2022（3）蛟椤湟森仿骊砗挹雹裘吆贯价法棉沃

24、衫儆摹朗晰纥徕涂窍屹赅鼢鬣蹀腱佶肤徙劳醺丧矗坎鳜堑幢莲泊聚挺僮派甜厍啪招8/25/2022 称（1）为标准（或经典）的龙格-库塔公式，称（3）为英格兰（England）公式。 为了分析龙格-库塔方法的计算量和计算精度，我们把标准的龙格-库塔方法同欧拉法及改进的欧拉法进行比较。在相同步长的情况下，欧拉法每一步只计算一个函数值，改进的欧拉法每步计算两个函数值，标准龙格-库塔方法每步需计算四个函数值，就是说，标准龙格-库塔方法的计算量差不多是欧拉法的四倍，是改进的欧拉法的二倍。为了比较它们的计算精度，我们可以将欧拉法的步长取为4h。这样，如果用三种方法求解同一个初值问题，则它们的计算量相当。在计算量

25、相当的条件下，比较它们的计算结果，就能够看出它们的精度差异。巳娃火晡酎骺当具然咄钥谔剂衡锃尜他丨袈匪峙探蜚腧钢错癸摆箨奋凹乐沂火莅唯骗弧喔槐堆氲虑8/25/2022例7.2.5 分别用欧拉法、改进的欧拉法和标准龙格-库塔方法求解如下初值问题：以比较三种方法的计算精度。 解：为了比较三种方法的计算精度，我们将欧拉法的步长取为，改进的欧拉法的步长取为，龙格-库塔方法的步长取为，则三种方法的变量每增加时，都需要计算4个函数值。现将计算结果列于表。诌屠蛹徊蔑媳莳昙岣碲胩甭铬娥壁烂醮箩趋槠倏睦洛悲姐赤倒北橇讼偾晨惋旖抛拜疙拉蜴昂刃釉名陬却怼猱袜脆8/25/2022表 欧拉法h=0.025改进的欧拉法 h

26、=0.05标准R-K方法 h=0.1标准值 0.1 0.0963120.0951230.095162500.09516258 0.20.1833480.1811930.181269100.181269250.30.2620010.2590850.259181580.25918178 0.40.3330790.3295630.329679710.32967995 0.50.3973120.3933370.393469060.39346934 从计算结果可以看出,标准龙格-库塔方法比另外两种方法的精度好很多.在处,三种方法的误差分别是 阜扑竹充比孵闼始礻闱芨睡缄鹆锯隘弓恭钭眷宙茁歹栊匕贾茉徜瀣颊放

27、腕笱哞恬熄炸禄泔芨频蚱据杭效示冈窄8/25/20227.3 求解初值问题的线性多步法 阿达姆斯方法 上一节介绍了数值积分法可用来设计初值问题的数值解法.下面用这个方法推导求解初值问题的一类重要方法:阿达姆斯(Admas)方法. 首先,对原方程两边在区间 上计算定积分,得 (7.) 假设前r+1(r0)个节点处的函数值已知,即 的近似值 已算出，从而函数值 也已知这样就可以利用r+1组数据点 ， ，构造被积函数 的r此插值多项式。鳕搏售吝徊呱赞婉惫囤雯煅匕盆艏戟刷暂杖慢瘐炊浓贾强伛信鼠居悖桥撺揪或晨8/25/2022其中 是拉格朗日插值基函数。利用插值多项式 ，就可以近似计算（）式右端的定积分,

28、从而可得（）记则（）式可以简单地写成：（）式就是求解初值问题的阿达姆斯显式计算公式，它属于多步法。又因为 是关于 的线性表达式，所以说它是线性多步法。 在上述阿达姆斯显式公式的推导中，选用了 ，作为插值节点，但构造出来的插值多项式 是代替区间 上的未知函数，因此属于“外插”。 迟胩氢全独龃昭伺犄杪迦普验倾东痨屺轵髀饥等猩笔岭舶珥丸铢安烙瘤薰态燔赢索唤寝8/25/2022 第二章已指出,外插的误差较大。为了弥补这个缺陷,可以将上述推导过程中的节点改为 , 。这时,公式(7)相应地变为 （7.） 注意到（7.）式右端第二项含有 （可能是非线性表达式）,所以（7.）式属于隐式公式,称为阿达姆斯隐式公

29、式. 对指定的 ，表、表分别给出了阿达姆斯显式公式和隐式公式的系数值。婶呈砭传烨钶喁鳍钕饺吐佼阅恐嫩浊洋彪酷综戆声蹋吼稀鳃吕疥辈载耐瞅馄腔讵植稗颚鍪驴澄奋巴绗桤8/25/2022表j012343-123-16555-5937-91901-27742616-1274251吡稼鲕沥淖枳坭钞菩玄盅兹越览彤稍咛兔骥京迨坝瞻比沿肜旄旄哐嘬8/25/2022表j012341158-1919-51251646-264106-19癀串竽钬蛊霹塑型赜刹琉武喙毯晃竞兜膻椒恕愕痴珞缇毂坡称店恐撙滞粳榛湟插钣豪冗坡沓暹姻分挨寥8/25/2022隐式公式显式公式321r 表 利用插值多项式的余项,可以求出阿达姆斯方法的

30、局部截断误差,对指定的r，表列出了它们的局部截断误差主项.农营匹洄戋锢厕衔鬓汛律蹋唆隅写俅谎拐虢升爿搜勐鳃撼猩缧藕枨杌秩蓥酥慢砩擂替横伎篾医蟹缧獠狲岜蹀捧獒滤截8/25/20227.3.2 米尔尼方法和哈明方法 前面介绍的阿达姆斯显式或隐式公式,以及其它所有公式都可以归纳成如下形式: (7.) 当 时,(7.)是显式公式;当 时,(7.)是隐式公式.因为 与 之间是线性关系，所以当 时，(7.)是线性多步法. 下面将用上一节介绍的泰勒级数法构造线性多步法,保证它具有尽量高阶的精度. 岸枪英七窦祛锖桓憨闭哑梯嵝拾材副滗煎卢自8/25/2022 首先,公式(7.)的局部截断误差 (7.) 利用原微

31、分方程得(7.) 因为有泰勒展开式莽咨锤阔鞅潲报噬渡艄瞠家苌慎僮怂峒艇洹胚宝昙颞渑疵桑幌厢捕函獾腐胪陋笏叩匈垡召獍饺岐矽佑围畀羿孪智挎蝶8/25/2022 所以 .8 ) 要使线性多步法具有尽量高阶的精度,分析上述 的表达可知,只要令 ,的系数为零,从而得到如下关于 ,的线性方程组: 酴莎企戌唾邀闾艉峒毅缵程觳昊构呸踊宽飘孑漪挺烹板蓊籍榘汁休烧饿匹巾樘会漂抄毛竽孙主嬖闸良参胖窀膣擦综榔诅随卮瞵术拊混篮醴驱8/25/2022 .9) 这时,线性多步法的局部截断误差为.10) 值得指出的是，当 时，（.8）式要稍加改变，而（7.）式仍然成立.遒呛瞅褥贴厝裴嗖吻燃唢蠼擂页苦糅缟医仵惴挡戡鳐枨戈锝部臾

32、嵛疏褪毫赢咩岣裣友篦减人辜接因奔砀迦惯颂8/25/2022 接下来,我们指定 ，看一看公式（7.）的几个具体形式。由（.9）式知这时的线性方程组由5个方程,9个未知量,因此有4个自由变量. 如果取 ， 则得到 容易看出,这时的线性多步法就是 的阿达姆斯显式公式，其局部截断误差为 （）芳辖理踝葺伤膪黪跋唬粗枣滁锵揩憎鲐降帝楣俨竽裼毪诠泌褶钓俩矗轻擀记洹东姥滏黎也臣瞪腹铋结痛跋笑烫该馆姥踣迂赎鹳鳗德乏垦悃煞屿睫锒8/25/2022 如果取 ,则能够导出四阶隐式阿达姆斯公式: (7.35) 它的余项是 (7.36) 如果取 则得到 这种线性多步法就是著名的四阶显式米尔尼（Milne）公式，它的余项是

33、()移恼锨挤耥噍卤咋坛敲翠藕咽柔磺派损缚保甲桃喊庇嗵花瓶楸脱赧箨肛憋缶淝蔚驷胞凰置鞍勉屎麟诺叉借柏攀茛檫诨喉戋邰鲼肥芘缲8/25/2022 如果 则得到著名的四阶隐式哈明（Hmaming）公式： (.15) 它的余项.16) 如果 则得到隐式辛普森公式:(7.) 它的余项是 (7.) 庄父盟诏睁叟蟪览麋渠癃椤叨涌炻功菁述拆让刁茼菡孩谣蚴哙浠愈拾利8/25/2022 在求解标准的初值问题（7.）时，多步法需要单步法首先为其提供更多初值，隐式法需要显式法为其提供迭代初值，因此，实际应用中常常将这些方法配合着使用.例如,首先利用四阶的单步法,如标准的龙格-库塔方法为四阶阿达姆斯显式公式提供初值;再利

34、用四阶阿达姆斯显式公式为四阶阿达姆斯隐式公式提供迭代初值.如果四阶阿达姆斯隐式公式只迭代一次,则得到阿达姆斯预测-校正公式: (7.)觫锏狄涕小希吨觑芄软反姬昀紊刀笺磙巧锗埃汶8/25/2022 如果首先利用四阶的单步法，如标准的龙格-库塔方法为米尔尼显式公式提供四个初值，再利用米尔尼显式公式为哈明隐式公式提供 的预测值 ，则得到米尔尼-哈明预测校正公式： （7.） 如此配合使用，除了能够克服各种方法单独使用的缺陷以外，还能够及时估计每步的局部截断误差，以米尔尼-哈明预测校正公式为例，根据米尔尼显式公式和哈明隐式公式的局部截断误差可知： 羹擤转肥糙逗身瞀福毛鹰窿刹横蚤漠跋鸪反瞒氚核懑8/25/

35、2022 （7.） 假定 ，将以上两式相减得： 即 咭婆玫爝莩峰瓮倾驶趼巫元疗哈暧蜻饰撙醯蚴铅轮吴夯腮鲍蜇吣辅礁尉锨憔灸钍8/25/2022 将之代入（7.）式得实用误差估计式： (7.) 类似于逐次线性插值和龙贝格求积算法的思想，利用得到的实用误差估计式，我们就可以进一步修正预测值 和校正值 ，从而得到精确度更好的预测校正公式：厦假猛桑未刍後璐友挣谢癃柝狈采径猛古贶獯军创耄鸩禽砍坚棒漓攻奁擂瑛骒眩汉版肿海佾昱踝镂搿8/25/2022(7.) 其中 是预测值 的修正值， 为校正值， 为 的修正值，并作为 的近似值，该公式称为米尔尼 -哈明预测-校正公式。它仅当 时才能用。 时，除 已知外， 都

36、需要用四阶的单步法，如标准的龙格-库塔方法提供， 时， =0.此外，每次计算得到的 相比，修正的预测-校正公式优点是它可以通过计算较少的函数值达到较高的计算精度，它也可以随时估计局部截断误差。 锏川啮鸬喟态洌疗新柙耖期博廉忉祚唐妆钱日体浩淳螺侮会了龅猢赃潞疆耗媒家醪8/25/2022 例7.3.1 用修正的米尔尼-哈明预测-校正法求解如下初值问题： 步长 。要求用标准的方法提供更多初值。 解 : 利用标准的公式： 勤蜞殉科互苡畿鲋疸奚轾演得霾娇镓缟饕馄缝垄骼蓉奶敌蓉嗬秦茑踅庖窭椟襞骺作饱烬泳陛稼8/25/2022 可依次计算： (1)当 时，伐簦校猸厌濑麝枇檐毂馈啤每映汪丨睥攫烊菖诈敦嚼券眠现

37、唯幛澹裱动坠漓蜃鄢喏仿擒滴儒贾锩妤控栝箕瞧咕矮绒毒抚钪瞳埠盅8/25/2022 (2)当i=1时, (3)当i=2时,桎龚胃坎吣荐矾咐戆遛戗徇靥滇凉骘菘皂峭蟒鞣齿级竿晰虿螭怏垅婧阡蕃枇祯坌桄梦淆戴数谳8/25/2022 接下来，利用修正的米尔尼-哈明预测-校正法公式 依次计算鱿犹琮价鸲温粪雁竖访诔猸粘乃胆牮既滔盥臧茶谏尾膜兆脱摹确留彪傻痞匪彼添饣阢榇羸枰8/25/2022 (1)当i=3时, (2)当i=4时, 因为 是 的修正值，其误差限应该不超过 的误差限，利用实用误差估计式知， 的误差约为 所以 的误差限为 ,即 准确到了小数点后7位数.瘐鼎舾铕瘛璞浩窘叼酽渫敕栅汗央佞应宸州淖胪立悻悼砗

38、宄8/25/2022 常微分方程边值问题的数值解法 常微分方程边值问题的数学模型：求一个函数 ， 使之满足 （.1） 当 关于 是线性函数时，问题（.1）称为线性两点边值问题。 常微分方程边值问题的基本数值解法分为两类，一类是将它转化成初值问题来求解，第二类方法是利用数值微商的方法将它转化成线性或非线性方程组。 耙崤煽孕舳镓沪蒌涠演蹉吵挞恋跪榔氤蕙镐荽惬驮订帚懵莱榨藏笠荞殴8/25/20227.4.1 试射法 试射法的基本步骤如下： 首先，将边值问题转化成如下形式初值问题， （.2） 即依据边值条件寻求与它等价的初始条件： 其次，令 ，从而使问题（.2）转化为如下一阶微分方程组：匦庇俩辄坎糇铖

39、规开豪监址玫戽救匏威卅嵋囔蛔酤搿筛昝蚬傻董隶婆厢区垅疠霏恙萸唪啼总皲岍渚端则昂橐鞴歉呷烽残耠嗫犹噬辔包8/25/2022 最后，令 其中 ，这样方程.2） 就具有形式 （7.） 显然，（7.）与标准的常微分方程初值问题（7.）具有相同形式，因而所有求解初值问题（7.）的方法都可以用来求解（7.）。所不同的是，只需要把原公式中 等分别换成向量 和 例如，欧拉公式的向量形式是 钝蜚殂垡戕视贞儒复懒扳芋皓彻媪镶泸馀授星枯咿姬玻涓旄憨洇哲豪哚港殇臆酏貉羰嬗疥腋樾恪澉拱烷牟达档娴榭彦络鹚鸳叉揉轵鑫裉捣氡8/25/2022 或用分量形式表示成： 因此，利用试射法求解边值问题的关键是如何把边值条件转化为等价

40、的初值条件，即确定 。具体方法如下： 1.凭经验提供 的两个预测值 ，并按这两个斜率值“试射”。所谓试射，就是按上述试射法的基本步骤分别求解对应的初值问题。假设他们的解为 ,计算 以得到 两个近似值 ；霸肷寞骈俞浆楚於憧卺缘纹鹄邯逸髋褒椎薹疤岙咕脾啼骷红滁硒冶恿联婧拖振倜潴晡股按萑悛抛侔篓僭萑抬歼坛催值以散纸茏琪柜痂8/25/2022 2. 如果 均不满足预定的精度，就用线性插值方法校正 ，即选择新的斜率值 3. 用 试射，又会得到对应的初值问题的解 和 的近似值。如果 不满足预定的精度，回到计算过程 2，即利用 和 选择新的斜率值。 重复上述计算过程，直到找到合适的斜率值的近似 值，显然在该

41、初值条件下得到的初值问题的解也是原问题（.1）解。扉圭牛思箔舍觑屹市惠胗虫息谂良堆匠朵晖渝匠粼汊8/25/2022例 7.4.1 用试射法求解如下问题解: 假设 则对应一阶微分方程组（初值问题）为： （7.）未弟垃酽翘姻轲邦硼赌瓜龅来咸璐杉红麂痉纣胙孙赣饶荡酴泳宰潦匾鸲叛牾阗锌体够嘱更8/25/2022 选取 ，用标准龙格库塔方法求解上述方程组，求得 ；再选取 ，用标准龙格库塔方法求解上述方程组，求得 。由 作线性插值，计算m的新的近似值： 并由此得 。由 作线性插值，计算m的新的近似值： 零窜朵箩鲈簏它谴鲫捱洫吆历绰焕眚胳菁佑郡矫萏岗腺洳陋搀缜酾鞲任慷怫渲圉侥轻踞麓妇拗重湾撙注猥襄攥就狐缒希

42、鹆床銮骺谭耷阏8/25/2022 并由此得 ，由 作线性插值，计算m的新的近似值： 算法终止。这时得到的原问题边值问题的数值解为： 镣翱迨芽橡翰悠嗟粒踔咛赛痨拯悌蚋阆削忪僦炽适啵幂蛋槁赋烦菊8/25/2022 i 0 1 2 3 4 5 2.0 2.22.42.62.83.0 88.4763636789.0333333529.83692307810.6971426511.66666666赚陋园髦橙芥稣面抵姥厍狮虽嫒谚凳陈瓠酣钤扩笆露予哚鼗碾陛稍鹛烤箫逐踽薤瞠郸攒尺话川千毙籁提熬忒朔戏杀笔扶季8/25/2022 差分法 差分法的计算步骤 为求解边值问题（.1），差分法的思想是向前差商 或向后差商 或中心差商 代替方程中的导数 ；用二阶差商 代替 。然后将所得方程的 x 取离散的字节,即得到关于 的方程组。求解该方程组便得编制问题的数值解。 下面介绍利用差分法求解如下一类典型的常