第十三讲常微分方程初值问题数值解法.ppt

它比欧拉法的?(xn,yn,h)=f(xn,yn),增加了计算一个右函数f的值，可望p=2.若要使得到的公式阶数p更大,?就必须包含更多的f值.实际上从方程(1.1)等价的积分形式(2.4)，即若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积节点，为此可将(3.3)的右端用求积公式表示为第31页，共59页，星期日，2025年，2月5日一般说来，点数r越多，精度越高，上式右端相当于增量函数?(x,y,h)，为得到便于计算的显式方法,可类似于改进欧拉法(3.1),(3.2)，将公式表示为其中这里ci,?i,?ij均为常数.(3.4)和(3.5)称为r级显式龙格-库塔(Runge-Kutta)法,简称R-K方法.第32页，共59页，星期日，2025年，2月5日当r=1,?(xn,yn,h)=f(xn,yn)时，就是欧拉法，此时方法的阶为p=1.当r=2时，改进欧拉法(3.1)式就是其中的一种,下面将证明其阶p=2.要使公式(3.4),(3.5)具有更高的阶p，就要增加点数r.下面我们只就r=2推导R-K方法.并给出r=3,4时的常用公式，其推导方法与r=2时类似，只是推导计算较复杂.第33页，共59页，星期日，2025年，2月5日9.3.2二阶显式R-K方法对r=2的R-K方法，由(3.4),(3.5)式可得如下计算公式这里c1,c2,?2,?21均为待定常数，我们希望适当选取这些系数，使公式阶数p尽量高.根据局部截断误差定义，推导出(3.6)的局部截断误差为第34页，共59页，星期日，2025年，2月5日其中这里yn=y(xn),yn+1=y(xn+1).为得到Tn+1的阶p，要将上式各项在(xn,yn)处做泰勒展开，由于f(x,y)是二元函数，故要用二元泰勒展开，各项展开式为第35页，共59页，星期日，2025年，2月5日将以上结果代入(3.7)，则有第36页，共59页，星期日，2025年，2月5日要使公式(3.6)具有p=2阶，必须使即(3.9)的解是不唯一的.可令c2=a≠0，则得这样得到的公式称为二阶R-K方法.第37页，共59页，星期日，2025年，2月5日则由此可以看出在改进的欧拉公式中相当于取(xn,yn),(xn+1,yn+1)两点处斜率的平均值，近似代替平均斜率，其精度比欧拉公式提高了.如取a=1/2，则c1=c2=1/2,?2=?21=1.这就是改进的欧拉公式(3.1)式.第38页，共59页，星期日，2025年，2月5日称为中点公式(变形的欧拉公式)，相当于数值积分的中矩形公式.也可以表示为如取a=1，则c1=0,c2=1,?2=?21=1/2.得计算公式第39页，共59页，星期日，2025年，2月5日对r=2的R-K公式(3.6)能否使局部误差提高到O(h4)?为此需把K2多展开一项，从(3.8)的看到展开式中的项是不能通过选择参数消掉的，实际上要使h3的项为零，需增加3个方程，要确定4个参数c1,c2,?2及?21，这是不可能的.故r=2的显式R-K方法的阶只能是p=2，而不能得到三阶公式.第40页，共59页，星期日，2025年，2月5日9.3.3三阶与四阶显式R-K方法要得到三阶显式R-K方法，必须r=3.此时计算(3.4),(3.5)的公式表示为其中c1,c2,c3及?2,?21,?3,?31,?32均为待定常数，公式(3.11)的局部截断误差为第41页，共59页，星期日，2025年，2月5日只要K1,K2将按二元泰勒展开，使Tn+1＝O(h4)，可得待定参数满足方程第42页，共59页，星期日，2025年，2月5日上页下页第十三讲常微分方程初值问题数值解法第1页，共59页，星期日，2025年，2月5日9.1引言科学技术中很多问题都可用微分方程的定解问题来描述，主要有初值问题与边值问题两大类，本章只考虑初值问题.常微分方程初值问题中最简单的例子是人口模型，设某特定区域在t0时刻人口为y(t0)=y0已知的，该区域的人口自然增长率为?，人口增长与人口总数成正比，所以t时刻的人口总数y(t)满足以下微分方程第2页，共59页，星期日，2025年，2月5日很多物理系统与时间有关，从卫星运行轨道到单摆运动，从化学反应到物种竞争都是随时间的延续而不断