7微分方程数值解课件-11

120第6章常微分方程初值问题数值解法本章探讨常微分方程特解的常用数值方法的构造和原理，主要介绍求常微分方程初值问题的常用方法和有关知识。重点论述Euler方法、Runge-Kutta方法和线性多步法的原理、构造、局部截断误差和稳定性等内容。1216.1实际案例工程技术里某些振动问题可以表示为单摆的运动，其运动规律的微分方程为：sin0(0)30(0)0gl怎样求出其特解()t？该微分方程不能用通常的解析方法来求解！怎样解不能用解析方法求解的微分方程特解问题是本章要解决的问题。1226.2问题的描述和基本概念1、常微分方程初值问题一般形式0(,)()yfxyyay()axb（6.1）式中(,)fxy已知，0()yay称为初值条件。初值问题的数值方法和数值解求函数()yyx在若干离散点kx上的近似值(0,1,,)kykn的方法称为初值问题的数值方法，而称(0,1,,)kykn为初值问题的数值解。1232.建立数值解法的思想与方法微分方程初值问题的数值解法是用离散化方法将初值问题化为差分方程后再求解的方式。设区间[a,b]上的一组节点为01naxxxb距离1kkkhxx称为步长。求数值解一般是从0y开使逐次顺序求出12,,yy。初值问题的解法有单步法和多步法两种：单步法：计算1ky时只用到ky一个值；多步法：计算1ky时要用1,,,kkklyyy多个值。数值解法还有显格式和隐格式之分。124微分方程离散化方法主要有数值微分法，数值积分法和Taylor展开法。初值问题化为差分方程的方法1.用离散方法去掉方程中的导数y得到近似离散化方程；2.在近似离散化方程中用ky代替()kyx；3.在近似离散化方程中将近似号“”用等号“=”代替。1251)数值微分法由初值问题（6.1）有'()(,())kkkyxfxyx，用数值微分的2点前差公式代替'()kyx，得近似离散化方程11()()'()(,())kkkkkkkyxyxyxfxyxxx记1kkhxx，做()kkyyx，“”，得差分方程1(,)kkkkyyfxyh写容易计算的形式1(,)kkkkyyhfxy(0,1,2,,1)kn（6.2）（Euler公式）由初值条件0()yya及式（6.2）可求出（6.1）的数值解12,,,nyyy。公式（6.2）是显式单步法。1262）数值积分法在1[,]kkxx上对'(,)yfxy两边取定积分，得111()()'(,())kkkkxxkkxxyxyxydxfxyxdx对右端积分采用梯形公式（数值积分公式）得近似离散化方程：111()()[(,())(,())]2kkkkkkhyxyxfxyxfxyx于是得到求初值问题（6.1）的梯形方法111[(,)(,)]2kkkkkkhyyfxyfxy(0,1,,1)kn该公式是隐式单步法。1273）Taylor展开法因为初值问题中函数(,)fxy是已知函数，由(,)yfxy，可以计算''y，'''y，…，于是有函数()yyx在kx处的Taylor展式212()()'()''()2!()(,())[(,())]2!kkkkkkkkxxhyxyxhyxyxhdyxhfxyxfxyxdx取上式右端前若干项，得近似离散化方程。例如取前两项有1()()(,())kkkkyxyxhfxyx于是又得到Euler公式：1(,)kkkkyyhfxy1286.3数值解法的误差、阶与绝对稳定性单步法数学描述为显式：1(,,)kkkkyyhxyh（6.4）隐式：11(,,,)kkkkkyyhxyyh（6.5）其中称为增量函数，与(,)fxy有关。129显式单步法的一些概念定义6.1设()yx是问题（6.1）的解，ky是经过式（6.4）求出的的计算解，则称111()kkkeyxy为单步法（6.4）在节点1kx的整体截断误差，而称11()()(,(),)kkkkkTyxyxhxyxh（6.6）为在1kx点的局部截断误差。()kyx表示解()yx在kx的值，是准确值，没有误差；ky表示由数值解公式得出()kyx的近似值，是数值解，有截断误差；ky表示用计算机计算ky给出的计算解，有舍入误差。130局部截断误差1kT的理解假设在计算()kyx时没有误差（()kkyyx）下，由式（6.4）计算出的1ky（1()(,(),)kkkkyyxhxyxh）与1()kyx的误差111kkkTyxy。整体截断误差还要加上()kyx与ky的误差。考察初值问题解法的优劣，引入阶的概念。定义6.2如果初值问题（6.1）对某种数值解法的局部截断误差为11()PkTOh则称该方法具有p阶精度或该方法是p阶方法。方法的阶越高，方法越好。131主局部截断误差或局部截断误差的主项如果某方法是p阶方法，11()PkTOh按h可展为1121()(,())()PPPkkkTOhgxyxhOh则称1(,())Pkkgxyxh为主局部截断误差。在同阶方法中，主局部截断误差越小，方法越好。例6.1常微分方程初值问题00(,)()yfxyyxy的单步法为1111[(,)2(,)],3nnnnnnkkhyyfxyfxyhxx试求其局部截断误差主项并回答它是几阶精度的？132解该单步公式的局部截断误差是1111()()[(,())2(,())]3nnnnnnnhTyxyxfxyxfxyx11()()()2()3nnnnhyxyxyxyx234(4)()()()()()2!3!4!nnnnnhhhyxhyxyxyxyx22()()()()()332!nnnnnhhyxyxhyxhyxyx231212(1)()()()()3323nnhyxhyxOh2321()()()6nhyxOhOh故局部截断误差主项是21()6nhyx，方法是一阶的。133求阶p的另一方法因为1()()(,(),)kkkkkTyxhyxhxyxh去掉下标，有()()()(,(),)Txyxhyxhxyxh若将()yxh在x点展开有1()()PTxOh则知该方法的阶是p。134例如，对Euler方法1(,)kkkkyyhfxy，有(,(),)(,())xyxhfxyx那么()()()(,())Txyxhyxhfxyx将()yxh在x点展开，有2()()'()''()2!hyxhyxhyxyx2()(,())''()2!hyxhfxyxyx故有22()''()()2hTxyxOh因此，Euler方法是一阶方法。135定义6.3设用某种数值方法求初值问题（6.1）在任意节点kx的数值解ky时，满足1kk，则称该数值方法是绝对稳定的。这里k是计算机计算ky时得出的计算解ky的舍入误差，kkkyy。通常用试验方程'yy（为复数）来讨论求解初值问题的数值方法绝对稳定性；对具体初值问题，可取(,)kkxyfy。稳定性常与步长h有关。136定义6.4某方法在试验方程中绝对稳定的复平面h范围称为该方法的绝对稳定域，它与复平面实轴的交称为该方法的绝对稳定区间。绝对稳定域包含复平面左半平面|Re0hh的方法称为是A-稳定的。绝对稳定域越大，对应的方法绝对稳定性越好。1376.4Euler方法的有关问题用Euler公式1(,)kkkkyyhfxy，0()yay求解初值问题（6.1）数值解的方法称为Euler方法。1)Euler方法的几何意义Euler方法常称为折线法。2)Euler方法的误差设1ky为1ky的计算解，(,)fxy满足(,())(,)()kkkkkkfxyxfxyLyxy其中0(0,1,2,,)Lkn。则有Euler方法的局部截断误差2()kTOh；138Euler方法的总体截断误差1111111()()kkkkkkkeyxyyxyyy由1,kkkkyyxhfxyx，1(,)kkkkyyhfxy,有11(1)kkkeThLe2111(1)(1)(1)kkkkThLThLThLT因为对任意m都有2()mTOh，可得21100(1)()(1)kkmmkkmmmehLTOhhL12(1)1()()11khLOhOhhL由k的任意性，可知Euler方法的总体截断误差为()keOh说明当h足够小时，由Euler方法计算所得数值解可以很好地逼近准确解，从而Euler方法是收敛的。1393)Euler方法稳定性将Euler公式用于试验方程'yy，得到1(1)kkkkyyhyhy（6.8）设计算ky时有舍入误差,0,1,2,kk，则有11(1)()kkkkyhy得1(1)kkh要想1kk，只须11h，因此Euler方法在11h时是绝对稳定的，其绝对稳定域为复平面h上以-1为中心的单位圆盘。由图形可知Euler方法的绝对稳定区间为2Re0h。140若指定是负实数，则有当步长h满足20h时，可保证Euler方法的计算绝对稳定。hOO-1-2Oh图6.3图6.4141例6.2'100,00.1,01yyxy，取0.025h用Euler方法求其数值解并与其准确解作比较。解本题准确解为100xye，Euler法计算公式为101001100,1kkkkyyhyhyy直接计算结果与对应准确值的误差有如下表xk0.0250.050.0750.1yk-1.52.25-3.3755.0625y(xk)-yk1.58-2.243.37-5.06可见计算结果波动大，不稳定，计算结果不可用。原因：100，h的稳定性范围02/0.02h，而本题的h=0.025〉0.02，故计算结果不稳定。1424)改进的Euler方法易证梯形方法（6.3）的绝对稳定域为复平面h上1122hh的部分，它是复平面h的整个左半平面（图6.4）。因此是A-稳定的。其绝对稳定域比Euler方法大，是二阶方法，故梯形方法比Euler方法好，但美中不足的是梯形方法是隐式方法。为避免迭代，实用中采取如下两次计算格式由ky来计算1ky（改进的Euler方法）：预测1(,)kkkkyyhfxy（6.9）校正111[(,)(,)]2kkkkkkhyyfxyfxy（6.10）预测—校正公式。易证它也是二阶方法。1436.5Runge-Kutta方法理论上，由Taylor展开法可以构造出解初值问题（6.1）的高阶数值方法，但这涉及到要计算(,)fxy的高阶导数，因此很不方便。本节的Runge-Kutta方法从函数(,)fxy本身着手，避开计算(,)fxy的高阶导数，来构造高阶数值方法。基本思想将微分方程初值问题转化为积分方程问题，再对积分方程中的定积分使用待定的m点插值型求积公式构造高阶的函数展开模式以获得高阶数值方法。1441.构造原理将初值问题的方程'(,)yfxy在1[,]kkxx上积分11()()(,())kkxkkxyxyxfxyxdx对定积分用m个点的插值型求积公式11()()(,())mkkiiiiyxyxhcfy11[,],ikkkkxxhxx取ikixah，并对(,())iify中的()ik