第七章 微分方程

第七章微分方程第一页，共三十七页，编辑于2023年，星期四§1引言常微分方程复习：一般的，凡表示未知函数、未知函数的导数与自变量之间关系的方程叫做微分方程。微分方程的解就是一个未知函数，将该未知函数代入微分方程中，能使方程成为恒等。

一阶常微分方程的初值问题Initial-ValueProblem该微分方程的特解y(x)为一条曲线，称为微分方程的积分曲线。第二页，共三十七页，编辑于2023年，星期四

考虑一阶常微分方程的初值问题/*Initial-ValueProblem*/:只要f(x,y)在[a,b]R1上连续，且关于y满足Lipschitz

条件，即存在与x,y无关的常数L使对任意定义在[a,b]上的(x,y1)和(x,y2)都成立，则上述IVP存在唯一解y=y(x)。要计算出解函数y(x)在一系列节点a=x0<x1<…<xn=b

处的近似值节点间距为步长，通常采用等距节点，即取hi=h

(常数)。第三页，共三十七页，编辑于2023年，星期四数值解法的基本特点（思路）：求解过程顺着节点排列次序一步步向前推进，即按递推公式由已知的y0,y1,…yi，求出yi+1。yox0y(x)p0x1p1p2x2xxnPn+1Pnxn+1第四页，共三十七页，编辑于2023年，星期四§2欧拉方法

/*Euler’sMethod*/将两端在上积分已知依次求出y1,y2,…yn第五页，共三十七页，编辑于2023年，星期四y(x)xyox0p0x1p1Pi+1xi+1xnPn+1Pnxn+1欧拉公式:yi+1

=yi+hf(xi,yi)欧拉公式的几何意义近似代替y(x)，用该直线与x=xi+1的交点Pi＋1(xi＋1,yi＋1)的纵坐标yi+1作为y(xi＋1)的近似值。xiPi(xi,yi)在区间用过点Pi(xi,yi)，以f(xi,yi)为斜率的直线y

=yi+f(xi,yi

)(x-xi)第六页，共三十七页，编辑于2023年，星期四

欧拉公式：x0x1向前差商近似导数记为定义在假设yi=y(xi)，即第i步计算是精确的前提下，考虑的截断误差Ri=y(xi+1)

yi+1称为局部截断误差/*localtruncationerror*/。定义若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。

欧拉法的局部截断误差：欧拉法具有1阶精度。Ri的主项/*leadingterm*/第七页，共三十七页，编辑于2023年，星期四

欧拉公式的改进：隐式欧拉法向后差商近似导数x0x1))(,()(1101xyxfhyxy+)1,...,0(),(111-=+=+++niyxfhyyiiii由于未知数yi+1

同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/

欧拉公式，而前者称为显式/*explicit*/欧拉公式。一般先用显式计算一个初值，再迭代求解。隐式欧拉法的局部截断误差：即隐式欧拉公式具有1阶精度。第八页，共三十七页，编辑于2023年，星期四梯形公式—显、隐式两种算法的平均注：局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。中点欧拉公式/*midpointformula*/中心差商近似导数x0x2x1假设，则可以导出即中点公式具有2阶精度。需要2个初值y0和y1来启动递推过程，这样的算法称为双步法/*double-stepmethod*/，而前面的三种算法都是单步法/*single-stepmethod*/。第九页，共三十七页，编辑于2023年，星期四方法显式欧拉隐式欧拉梯形公式中点公式简单精度低稳定精度低,计算量大精度提高计算量大精度提高,显式多一个初值,可能影响精度第十页，共三十七页，编辑于2023年，星期四例题：用Euler公式和改进的Euler公式分别求下列初值问题的数值解(取步长h=0.1计算到y3)：y´=-2xy2

y(0)=1解：由欧拉公式

yn+1=yn+hf(xn,yn)=yn-2hxnyn2计算如下y1=y0-2hx0y02=1-2·0·0·12=1y2=y1-2hx1y12=1-2·0.1·0.1·12=0.98y3=y2-2hx2y22=1-2·0.1·0.2·0.982=0.9416第十一页，共三十七页，编辑于2023年，星期四改进欧拉法的预-校公式计算如下y1=0.99;y2=0.9614;y3=0.9173精确解y(0.1)=0.99,y(0.2)=0.9614;y(0.3)=0.9173可见改进欧拉公式比欧拉公式精度高第十二页，共三十七页，编辑于2023年，星期四改进欧拉法/*modifiedEuler’smethod*/Step1:先用显式欧拉公式作预测，算出),(1iiiiyxfhyy+=+Step2:再将代入隐式梯形公式的右边作校正，得到1+iy)],(),([2111+++++=iiiiiiyxfyxfhyy注：此法亦称为预测-校正法/*predictor-correctormethod*/。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。第十三页，共三十七页，编辑于2023年，星期四为便于编制程序上机计算，改进欧拉公式可以改写为如下形式：第十四页，共三十七页，编辑于2023年，星期四§3

Runge-Kutta方法[一]基本思想Runge-Kutta方法是一种高精度的单步法,简称R-K法．得到高精度方法的一个直接想法是利用Taylor展开．假设式y'=f(x,y)(a≤x≤b)

中的f(x,y)

充分光滑,将y(xn+1)在xn点作Taylor展开:y(xn+1)=y(xn)+hy'(xn)+(h2/2!)y''(xn)+......+(hp/p!)y(p)(xn)+..... 其中 y'(x)=f(x,y(x)) y''(x)=[f(x,y(x))]'x=fx+f·fy ............................... y(p)(x)=[f(x,y(x))](p)x第十五页，共三十七页，编辑于2023年，星期四对照标准形式yn+1=yn+hф(xn,yn;h).若取ф(x,y;h)=y'(x)+(h/2!)y''(x)+......+(hp-1/p!)y(p)(x)并以yn代替y(xn),则得到一个p阶近似公式

yn+1=yn+hф(xn,yn;h)(n=0,1,2,......)(**)显然p=1时,式(**)就是 yn+1=yn+hf(xn,yn)它即为我们熟悉的Euler方法.当p≥2时,要利用公式(**)就需要计算f(x,y)的高阶微商.这个计算量是很大的.因此，利用式(**)构造高阶公式是不实用的.第十六页，共三十七页，编辑于2023年，星期四R-K方法不是直接使用Taylor级数,而是利用它的思想,即计算f(x,y)在不同结点的函数值,然后作这些函数值的线性组合,构造近似公式,式中有一些可供选择的参数.将近似公式与Taylor展开式相比较,使前面的若干项密合,从而使近似公式达到一定的精度.下面以二级二阶R-K方法为例说明这一方法的基本思想.第十七页，共三十七页，编辑于2023年，星期四[二]二级二阶R-K方法在[xn,xn+1]上,取f(x,y)在两个点的函数值作线性组合,记得到二级R-K方法:yn+1=yn+h(c1K1+c2K2)K1=f(xn,yn) （*） K2=f(xn+a2h,yn+b21hK1)其中c1,c2,a2,b21为待定参数.对照式(**)有ф(x,y,h)=c1f(x,y)+c2f(x+a2h,y+b21hf(x,y))若要求式(**)达到二阶精度，则只要局部截断误差Tn+1=O(h3).将y(xn+1)及ф(xn,y(xn),h)在xn作Taylor展开, f(xn,y(xn))=y'(xn) f(xn+a2h,y(xn)+b21hK1)=f(xn,y(xn)+a2hfx+b21hf·fy+O(h2))其中 fx=fx(xn,y(xn)),fy=fy(xn,y(xn)),f=f(xn,y(xn))由此得第十八页，共三十七页，编辑于2023年，星期四

ф(xn,y(xn);h)=(c1+c2)y'(xn)+c2(a2hfx+b21hfyf)+O(h2) 因为y(xn+1)在xn处的Taylor展开为 y(xn+1)=y(xn)+hy'(xn)+(h2/2!)y''(xn)+O(h3)由显式单步法在xn+1的局部截断误差定义Tn+1=y(xn+1)-y(xn)-hф(xn,y(xn),h)=h(1-c1-c2)y'(xn)+h2[(1/2-a2c2)fx+(1/2-c2b21)fyf]+O(h3)显然，若要求Tn+1=O(h3),则应有c1+c2=1c2a2=1/2,c2b21=1/2上方程组含有3个方程,4个未知数,其解是不唯一的.若取c2=α为自由参数,则得其的一族解c1=1-α,c2=α,a2=b21=1/（2α） （***） 满足条件（***)的(*)式称为二级二阶R-K方法.特别当α=1/2时,公式即是前面介绍的改进的Euler方法.第十九页，共三十七页，编辑于2023年，星期四

当α

=1时,c1=0,c2=1，得yn+1=yn+hK2n=0,1,….N-1K1=f(xn,yn)K2=f(xn+h/2,yn+hK1/2)这就是变形的欧拉方法或中点方法。当α

=3/4时，c1=1/4c2=3/4，得到Heun方法。yn+1=yn+h(K1+3K2)/4,n=0,1,….N-1K1=f(xn,yn)K2=f(xn+2h/3,yn+2hK1/3)

第二十页，共三十七页，编辑于2023年，星期四二级R-K方法是显示单步式,每前进一步需要计算两个函数值.由上面的讨论可知,适当选择四个参数c1,c2,a2,b21,可使每步计算两次函数值的二阶R-K方法达到二阶精度.能否在计算函数值次数不变的情况下,通过选择四个参数,使得二阶R-K方法的精度再提高呢?我们说,答案是否定的.无论四个参数怎样选择,都不能使公式(**)提高到三阶.这说明每一步计算两个函数值的二阶R-K方法最高阶为二阶.若要获得更高阶得数值方法,就必须增加计算函数值的次数.第二十一页，共三十七页，编辑于2023年，星期四

m级显式Runge-Kutta方法

仿照二级R-K方法,在[xn,xn+1]上,取f在m个点的函数值做线性组合,即得到m级R-K方法

m

yn+1=yn+h∑crKr

r=1 K1=f(xn,y

n) r-1 Kr=f(xn+har,yn+h∑brsKs)，r=2,3,....,m

s=1

第二十二页，共三十七页，编辑于2023年，星期四使用不同的方法确定参数cr,ar,brs可使上式成为不同阶的R-K方法.在m级R-K方法中,最著名的是经典R-K方法:

yn+1=yn+(K1+2K2+2K3+K4)h/6 K1=f(xn,yn) K2=f(xn+h/2,yn+hK1/2） K3=f(xn+h/2,yn+hK2/2) K4=f(xn+h,yn+hK3)

由于它每前进一步需要计算四个点的函数值,因此称为四级公式.按定义可直接验证它的局部截断误差为O(h5),故它是四阶方法.第二十三页，共三十七页，编辑于2023年，星期四

前面已经看到,二阶、四阶R-K方法可分别达到最高阶数二阶、四阶，但是N阶R-K方法的最高阶却不一定是N阶。R-K方法的级数表示公式中计算函数值f的次数。Butcher给出了R-K方法计算函数值f的次数与阶数之间的关系表，如下：

计算f的次数 1 23 4567 方法的最高阶数12 3 4 4 5 6由表可见，四级以下R-K的方法其最高阶数与计算f的次数一致，对m阶R-K公式，当m>4，虽然计算f的次数增加，但是方法阶数不一定增加。因此四级四阶R-K公式是应用最为广泛的公式。 第二十四页，共三十七页，编辑于2023年，星期四

绝对稳定性与绝对稳定域求解初值问题的数值方法，当给定不同步长计算时结果的舍入误差影响差别很大，如果舍入误差不增长算法就是数值稳定的，若舍入误差增长很快算法就不稳定。定义：用一个数值方法求解微分方程y´=y是复数(1.5)对给定的步长h,在计算yn时引起的误差n，若这个误差在计算后面的yn+k中所引起的误差n+k满足:|n+k|≤|n|(k=1,2,…)就说这个数值方法对步长h和复数是绝对稳定的,使得数值方法是绝对稳定的H=h在复平面上的允许范围称为数值方法的绝对稳定域.实验方程第二十五页，共三十七页，编辑于2023年，星期四特拉法加尔（Trafalgar）海战和纳尔森（Nelson）秘诀19世纪中叶，法国拿破伦统帅大军要与英国争夺海上霸主地位，而实施这一战略的最主要的关键是消灭英国的舰队。英国海军统帅、海军中将纳尔森亲自制定了周密的战术方案。第二十六页，共三十七页，编辑于2023年，星期四1805年10月21日，这场海上大战爆发了。英国是纳尔森亲自统帅的地中海舰队，由27艘战舰组成；另外一方是由费伦纽夫（Villenuve）率领的法国——西班牙联合舰队，共有33艘战舰。Trafalgar大海战的概况是：费伦纽夫（Villenuve）率领的法国——西班牙联合舰队采用常规的一字横列，以利炮火充分展开，而纳尔森的战术使费伦纽夫大出意外。第二十七页，共三十七页，编辑于2023年，星期四英国的舰队分成两个纵列：前卫上风纵列由12艘战舰组成，由纳尔森亲自指挥，拦腰将法国——西班牙联合舰队切为两段；后卫下风纵列由英国海军中将科林伍德(Collingwood)指挥，由15艘战舰组成。在一场海战后，法国——西班牙联合舰队以惨败告终：联合舰队司令费伦纽夫连同12艘战舰被俘，8艘沉没，仅13艘逃走，人员伤亡7000人。而英国战舰没有沉没，人员伤亡1663人，但是，作为统帅的纳尔森阵亡。第二十八页，共三十七页，编辑于2023年，星期四秘密备忘录中的纳尔森（Nelson）秘诀：预期参加战斗的英国舰队：40艘。法国—西班牙联合舰队：46艘。预计联合舰队战斗队形一字横列。英国舰队的战斗队形与任务：分成两个主纵列及一个小纵列。第二十九页，共三十七页，编辑于2023年，星期四主纵列1：16艘，由纳尔森亲自指挥，拦腰将法国——西班牙联合舰队切为两段，并攻击联合舰队的中间部分。主纵列2：16艘，由英国海军中将科林伍德指挥，从联合舰队后半部再切断，分割并攻击后部12艘。小纵列：8艘，在中心部分附近攻击其先头部分的3-4艘