数值分析-第五章-常微分方程数值解法

图５畅２令珔h＝hλ，则yn＋１＝１＋珔h＋１２珔h２＋１６珔h３＋１２４珔h４yn．由此可知，绝对稳定性区域在珔h＝hλ复平面上满足｜１＋珔h＋１２珔h２＋１６珔h３＋１２４珔h４｜≤１的区域，也就是由曲线１＋珔h＋１２珔h２＋１６珔h３＋１２４珔h４＝eiθ所围成的区域．如图５畅２所示．例22　用Euler法求解y′＝－５y＋x，y（x０）＝y０，　x０≤x≤X．从绝对稳定性考虑，对步长h有何限制？解　对于模型方程y′＝λy（λ＜０为实数）这里λ＝抄f抄y＝－５．由｜１＋hλ｜＝｜１－５h｜＜１得到对h的限制为：０＜h＜０畅４．四、习题１畅取步长h＝０畅２，用Euler法解初值问题y′＝－y－xy２，y（０）＝１．　（０≤x≤０畅６），２畅用梯形公式解初值问题y′＝８－３y，y（１）＝２，　（１≤x≤２），２１２取步长h＝０畅２，小数点后至少保留５位．３畅用改进的Euler公式计算初值问题y′＝１xy－１xy２，y（１）＝０畅５，　１＜x＜１畅５，取步长h＝０畅１，并与精确解y（x）＝x１＋x比较．４畅写出用梯形格式的迭代算法求解初值问题y′＋y＝０，y（０）＝１的计算公式，取步长h＝０畅１，并求y（０畅２）的近似值，要求迭代误差不超过１０－５．５畅写出用四阶经典Runge唱Kutta法求解初值问题y′＝８－３y，y（０）＝２的计算公式，取步长h＝０畅２，并计算y（０畅４）的近似值，小数点后至少保留４位．６畅证明公式yn＋１＝yn＋h９（２K１＋３K２＋４K３）．K１＝f（xn，yn），K２＝fxn＋h２，yn＋h２K１，K３＝fxn＋３４h，yn＋３４hK２，至少是三阶方法．７畅试构造形如yn＋１＝α（yn＋yn－１）＋h（β０fn＋β１fn－１）的差分格式，使其为二阶方法，并求出其局部截断误差的首项．３１２８畅导出具有下列形式的三阶方法：yn＋１＝a０yn＋a１yn－１＋a２yn－２＋h（b０y′n＋b１y′n－１＋b２y′n－２）．９畅证明yn＋１＝１３（４yn－yn－１）＋２３hy′n＋１是二阶公式．１０畅就初值问题y′＝ax＋b，y（０）＝０导出改进Euler方法的近似解的表达式，并与准确解y＝１２ax２＋bx相比较．１１畅用Adams显式公式作预估公式，Adams隐式公式作校正公式，求解初值问题y′＝－y－２sinx，　０＜x＜２，y（０）＝１，取h＝０畅２，已知y１＝０畅７８１１３８，y２＝０畅５３１３０７，y３＝０畅２６０４３５．１２畅已知三阶Adams显式公式及截断误差yn＋１＝yn＋h１２（２３fn－１６fn－１＋５fn－２），T１＝y（xn＋１）－yn＋１＝３８h４y（４）（ξ１），　xn－２＜ξ１＜xn＋１和三阶Adams隐式公式及截断误差yn＋１＝yn＋h１２（５fn＋１＋８fn－fn－１），T２＝y（xn＋１）－yn＋１＝－１２４h４y（４）（ξ２），　xn－１＜ξ２＜xn＋１，试求其预估—改进—校正系统．１３畅设有常微分方程初值问题y′＝f（x，y），y（x０）＝y０的单步法yn＋１＝yn＋h３［f（xn，yn）＋２f（xn＋１，yn＋１）］，４１２证明该方法是无条件稳定的．１４畅对任意选取的步长h，试证明改进Euler法yn＋１＝yn＋h２［f（xn，yn）＋f（xn＋１，yn＋hf（xn，yn））］，中点公式yn＋１＝yn＋hfxn＋１２h，yn＋１２hK１（其中K１＝f（xn，yn））及Heun公式yn＋１＝yn＋h４f（xn，yn）＋３fxn＋２３h，yn＋２３hf（xn，yn），分别对初值问题y′＝－２y－４x，y（０）＝２给出相同的逼近，为什么？１５畅证明解初值问题y′＝f（x，y），y（x０）＝y０的二步法yn＋１＝１２（yn＋yn－１）＋h４（４fn＋１－fn＋３fn－１）（fi＝f（xi，yi））是二阶的，并求其局部截断误差主项．１６畅证明线性多步法yn＋１＋α（yn－yn－１）－yn－２＝１２（３＋α）h（fn＋fn－１）存在α的一个值，使方法是四阶的．１７畅证明线性二步法yn＋１＋（b－１）yn－byn－１＝１４h［（b＋３）fn＋１＋（３b＋１）fn－１］，当b≠－１时方法为二阶的，当b＝－１时方法为三阶的．１８畅用经典四阶R唱K方法计算初值问题y′＝－２０y　（０≤x≤１），y（０）＝１．试分析当步长分别取h＝０．１及０畅２时，计算的稳定性．５１２１９畅讨论求解初值问题y′＝λy，　λ＜０，y（０）＝a的二阶中点公式yn＋１＝yn＋hfxn＋h２，yn＋h２f（xn，yn）的稳定性．２０畅对于初值问题y′＝－１０００（y－g（x））＋g′（x），y（０）＝g（０），其中g（x）为已知函数，其解y（x）＝g（x）．（１）若用Euler法求解，从稳定性考虑步长应在什么范围内选取？（２）若用隐式Euler法解题，从稳定性考虑，步长有没有限制？（３）若g（x）为不超过一次的多项式，用Euler法求解此问题时，从精确阶考虑，步长的选取有无限制？２１畅填空题（１）解初值问题y′＝f（x，y），y（x０）＝y０的梯形格式yn＋１＝yn＋h２［f（xn，yn）＋f（xn＋１，yn＋１）］是阶方法．（２）解初值问题y′＝１０y－xy，１≤x≤２，若用梯形法求解，要使迭代法y（m＋１）n＋１＝yn＋h２［f（xn，yn）＋f（xn＋１，y（m）n＋１）］，m＝０，１，…收敛，步长h＜？（３）用隐式Euler法求解y′＝f（x，y），y（x０）＝y０，要使数值计算是稳定的，步长０＜h＜．（４）显式Euler格式的绝对稳定区间为．五、习题解答１畅解　f（x，y）＝－y－xy２，Euler格式为６１２yn＋１＝yn＋hf（xn，yn）＝yn＋０畅２（－yn－xny２n）＝０畅８yn－０畅２xny２n．由y０＝１计算得y（０畅２）≈y１＝０畅８，　　y（０畅４）≈y２＝０畅５８８８，y（０畅６）≈y３＝０畅３８９５．２畅解　f（x，y）＝８－３y，梯形公式为yn＋１＝yn＋h２［f（xn，yn）＋f（xn＋１，yn＋１）］＝yn＋０畅２２［８－３yn＋８－３yn＋１］．整理得显格式为yn＋１＝７１３yn＋１６１３．由y（１）＝y０＝２计算得y（１畅２）≈y１＝２畅３０７６９，y（１畅４）≈y２＝２畅４７３３７，y（１畅６）≈y３＝２畅５６２５８，y（１畅８）≈y４＝２畅６１０６２，y（２畅０）≈y５＝２畅６３６４９．３畅解　f（x，y）＝１xy－１xy２，改进的Euler公式为y０＝０．５，珔yn＋１＝yn＋hf（xn，yn），Δyn＝h２［f（xn，yn）＋f（xn＋１，珔yn＋１）］，yn＋１＝yn＋Δyn，　i＝０，１，２，３，４．计算结果见表５畅５．７１２表5畅5nxnynf（xn，yn）珔yn＋１f（xn＋１，珔yn＋１）Δyny（xn）０１０畅５０畅２５０畅５２５０畅２２６７０４０畅０２３８３５０畅５１１畅１０畅５２３８３５０畅２２６７５６０畅５４６５１１０畅２０６５３１０畅０２１６６４０畅５２３８０９２１畅２０畅５４５４９９０畅２０６６０８０畅５６６１６００畅１８８９４１０畅０１９７７７０畅５４５４５５３１畅３０畅５６５２７６０畅１８９０３００畅５８４１７９０畅１７３５１００畅０１８１２７０畅５６５２１６４１畅４０畅５８３４０３０畅１７３６０３０畅６００７８３０畅１５９８９８０畅０１６６７５０畅５８３３３５１畅５０畅６０００７８０畅６０００００　　４畅解　梯形格式的迭代算法为y（０）n＋１＝yn＋hf（xn，yn），y（k＋１）n＋１＝yn＋h２［f（xn，yn）＋f（xn＋１，y（k）n＋１）］，　　k＝０，１，２，…，n＝０，１，２，…．于是取f（x，y）＝－y，有y（０）n＋１＝０畅９xn，y（k＋１）n＋１＝０畅９５yn－０畅０５y（k）n＋１．由y（０）＝y０＝１，经计算有y（０）１＝０畅９，　y（１）１＝０畅９０５，　y（２）１＝０畅９０４７５，y（３）１＝０畅９０４７６２５，　y（４）１＝０畅９０４７６１８７５．因｜y（４）１－y（３）１｜＝６畅２５×１０－７＜１０－５，于是取y（０畅１）≈y１＝y（４）１＝０畅９０４７６１８７５，则y（０）２＝０畅８１４２８６，　y（１）２＝０畅８１８８０９，y（３）２＝０畅８１８５８３，　y（４）２＝０畅８１８５９５，y（５）２＝０畅８１８５９４．因｜y（５）２－y（４）２｜＝１０－６＜１０－５，故得y（０畅２）≈y２＝y（５）２＝０畅８１８５９４．５畅解　f（x，y）＝８－３y，四阶经典Runge唱Kutta格式为８１２yn＋１＝yn＋１６（K１＋２K２＋２K３＋K４），K１＝hf（xn，yn）＝１畅６－０畅６yn，K２＝hf（xn＋h２，yn＋K１２）＝１畅１２－０畅４２yn，K３＝hf（xn＋h２，yn＋K２２）＝１畅２６４－０畅４７４yn，K４＝hf（xn＋h，yn＋K３）＝１畅２２０８－０畅４５７８yn．故yn＋１＝１畅２６４８＋０畅５２５７yn．由于y（０）＝y０＝２，故y（０畅２）≈y１＝２畅３１６２，y（０畅４）≈y２＝２畅４８２４．６畅证　只需证明公式的局部截断误差为O（h４）即可．　　为求局部截断误差，假设yn＝y（xn），因K１＝f（xn，yn）＝f（xn，y（xn））＝y′（xn），K２＝y′（xn）＋h２［fx（xn，y（xn））＋fy（xn，y（xn））y′（xn）］＋１２！h２４［fxx（xn，y（xn））＋２fxy（xn，y（xn））y′（xn）＋fyy（xn，y（xn））（y′（xn））２］＋O（h３）＝y′（xn）＋h２y″（xn）＋h２８［fxx（xn，y（xn））＋２fxy（xn，y（xn））y′（xn）＋fyy（xn，y（xn））（y′（xn））２］＋O（h３）＝y′（xn）＋h２y″（xn）＋h２８y碶（xn）＋O（h３），K３＝y′（xn）＋３４hy″（xn）＋３４h２２y″（xn）fy（xn，y（xn））＋１２！９１６h２［fxx（xn，y（xn））＋２fxy（xn，y（xn））y′（xn）＋fyy（xn，y（xn））（y′（xn））２］＋O（h３）．９１２由此得yn＋１＝yn＋h９（２K１＋３K２＋４K３）＝y（xn）＋hy′（xn）＋h９３２＋３y″（xn）＋h３９［３８（fxx（xn，y（xn））＋２fxy（xn，y（xn））y′（xn）＋fyy（xn，y（xn））（y′（xn））２）＋３２y″（xn）fy（xn，y（xn））＋９８（fxx（xn，y（xn））＋２fxy（xn，y（xn））y′（xn）＋fyy（xn，y（xn））（y′（xn））２］＋O（h４）＝y（xn）＋y′（xn）h＋h２２！y″（xn）＋h３３！y碶（xn）＋O（h４）．（５畅３２）而y（xn＋１）＝y（xn）＋y′（xn）h＋h２２！y″（xn）＋h３３！y碶（xn）＋O（h４）．（５畅３３）比较（５畅３２）、（５畅３３）两式可知，公式的局部截断误差至少是四阶的．因此该公式至少是三阶方法．７畅解　依题意只需公式的局部截断误差为O（h３），从而确定出待定系数α、β０、β１．为求得局部截断误差，假设yn－１＝y（xn－１），yn＝y（xn）．将yn＋１在xn处Taylor展开，即分别将差分格式右端的每一项在xn处展开，yn＋１＝α（y（xn）＋y（xn－１））＋h（β０y′（xn）＋β１y′（xn－１）），y（xn－１）＝y（xn）－y′（xn）h＋h２２y″（xn）－h３３！y碶（xn）＋O（h４），y′（xn－１）＝y′（xn）－y″（xn）h＋h２２y碶（xn）－h３３！y（４）（xn）＋O（h４）．０２２代入后得　yn＋１＝α（y（xn）＋y（xn）－y′（xn）h＋h２２y″（xn）－h３３！y碶（xn）＋O（h４））＋β０hy′（xn）＋β１hy′（xn）－β１y″（xn）h２＋β１２y碶（xn）h３－β１h４３！y（４）（xn）＋O（h５）＝２αy（xn）＋（β０＋β１－α）