1、1一、一、常微分方程的定解问题与应用常微分方程的定解问题与应用应用：自然科学领域，如物理；应用：自然科学领域，如物理；工程技术问题，如石油勘探。工程技术问题，如石油勘探。常微分方程的定解问题主要有常微分方程的定解问题主要有初值问题初值问题和和边值问题边值问题两两大类，我们仅考虑大类，我们仅考虑初值问题初值问题。例：例：马尔萨斯人口模型：假设某特定区域在马尔萨斯人口模型：假设某特定区域在t0时刻的人口时刻的人口p(t0)=p0为已知的，该区域人口的自然增长率为为已知的，该区域人口的自然增长率为。人口的增长与人口的总数成正比，所以人口的增长与人口的总数成正比，所以 t 时刻的人口时刻的人口总数总数
2、 p(t)满足如下的微分方程：满足如下的微分方程：00,ptp tp tp一般地，称这样的方程为一般地，称这样的方程为模型方程模型方程。2二、二、常微分方程的解法：常微分方程的解法：解析法解析法：给出精确解析解。只适合少数简单情况。：给出精确解析解。只适合少数简单情况。近似解法近似解法：给出解的近似表达式。如级数法：给出解的近似表达式。如级数法,逐步逼近法逐步逼近法。数值方法数值方法：给出方程在离散点上的近似解。它适合计算机：给出方程在离散点上的近似解。它适合计算机求解求解,应用广泛应用广泛,具有理论应用价值。具有理论应用价值。三、常微分方程初值问题的数值方法三、常微分方程初值问题的数值方法单
3、步法单步法Euler方法方法Taylor方法和方法和Runge-Kutta方法方法多步法多步法Adams方法和一般线性多部法方法和一般线性多部法线性多部法的收敛性与稳定性线性多部法的收敛性与稳定性3一、常微分方程初值问题的一般提法一、常微分方程初值问题的一般提法问题问题:求函数求函数(),y xaxb满足满足(,),()dyf x yaxbdxy a其中其中:f(x,y)为已知函数为已知函数,是已知值是已知值.(可能是观察值或实验值可能是观察值或实验值)(2)(1)基本条件基本条件:设设(,),Dx y axb y f(x,y)在在D上连续；上连续；f(x,y)在在D上关于变量上关于变量y满足
4、满足Lipschitz连续条件：连续条件：121212.(,)(,),f x yf x yL yyaxby y其中其中L为为Lipschitz常数。常数。(3)4若若f(x,y)在在D上满足基本条件上满足基本条件,一阶常微分方程初值问题一阶常微分方程初值问题(1),(2)对任意给定的对任意给定的存在唯一解存在唯一解且在且在a,b上连续可微上连续可微.定理定理1关于解关于解y(x)的适定性：的适定性：定义定义1 方程方程(1),(2)的解的解y(x)称为称为适定的适定的,若存在常数若存在常数0，K 0,对任意满足条件对任意满足条件(),(),xx 及的常微分方程初值问题：常微分方程初值问题：(,
5、)(),dzf x zx axbdx存在唯一解存在唯一解z(x),且有且有()()y xz xK()z x摄动摄动(扰动扰动)误差误差(4)(1),(2)上各加一个上各加一个摄动摄动(扰动扰动)项项.(1),(2)上各加一个上各加一个摄动摄动(扰动扰动)项项.5定理定理2 适定问题的解适定问题的解y(x)连续依赖于连续依赖于(1)式右端的式右端的f(x,y)和初值和初值.或者说解或者说解y(x)关于关于(1)式右端的式右端的f(x,y)和初值稳定和初值稳定.假设假设f(x,y)在在D上满足基本条件上满足基本条件,从而方程从而方程(1),(2)的解的解y(x)存在且适定存在且适定.注注:若若f(
6、x,y)在在D上满足基本条件上满足基本条件,则微分方程则微分方程(1),(2)的解的解y(x)是是适定适定的的.6假设初值问题假设初值问题(1)的解的解y=y(x)唯一存在且足够光滑唯一存在且足够光滑,对求解区域对求解区域a,b做剖分做剖分 1 构造数值解法的基本思想构造数值解法的基本思想在区间在区间xk,xk+1上对方程上对方程(1)做积分做积分,则有则有二、二、初值问题数值解的基本概念初值问题数值解的基本概念00121,nkknkkkxaxxxxb hxx 常用等步长常用等步长:()hb a n,则有则有(0,).kxakh kn(1)，(2)的准确解记为的准确解记为y(x),(,).kk
7、kf x yf记称为称为步长步长。(1,2,)kn()ky x的近似解记为的近似解记为,ky11()()(,()kkxkkxy xy xf x y xdx7因此，建立节点处近似值因此，建立节点处近似值yn满足的差分公式满足的差分公式称为称为Euler公式公式.称为称为梯形公式梯形公式.10(,),0,1,2,1kkkkyyhfxyykn1110(,)(,)2,0,1,2,1kkkkkkhyyf xyf xyykn若对若对*式右边的积分应用梯形求积公式，则可导出差分公式式右边的积分应用梯形求积公式，则可导出差分公式对右边的积分应用左矩形公式，则有对右边的积分应用左矩形公式，则有1()()(,()
8、kkkky xy xhfxy x8称为称为Euler中点公式中点公式或称或称双步双步Euler公式公式.若在区间若在区间xk-1,xk+1上对方程上对方程(1)做积分做积分,则有则有1111()()(,()kkxkkxy xy xf x y x dx对右边的积分应用中矩形求积公式，则得差分公式对右边的积分应用中矩形求积公式，则得差分公式1102(,),0,1,2,1kkkkyyhfxyykn因此，求初值问题数值解的基本方法是因此，求初值问题数值解的基本方法是步进法步进法.即逐个节点计即逐个节点计算算,由由 yi(i k)计算计算yk+1步进法步进法单步法单步法:多步法多步法:由由yk-j(j=
9、0,1,l-1)计算计算yk+1共用到共用到l个个值值.即即yk+1，yk-1,yk-l+1称为称为l步法步法。仅由仅由yk计算计算yk+19单步法与多步法的区别单步法与多步法的区别:(1)计算方面计算方面:l步方法只用于步方法只用于,kykl的计算的计算,即即y0，y1,yl-1的计算要用其它方法。的计算要用其它方法。(2)理论分析理论分析:单步法比单步法比 l 1的多步法容易分析的多步法容易分析(稳定性稳定性).(3)选步长方面选步长方面:单步法容易改变步长单步法容易改变步长.(4)精度精度:多步法精度高一些多步法精度高一些.单步法与多步法均有单步法与多步法均有显式显式和和隐式方法隐式方法
10、之分之分.在在Euler公式和公式和Euler中点公式中中点公式中,需计算需计算yk+1的已显式表示的已显式表示,称为称为显式公式显式公式,而梯而梯形公式中形公式中,需要计算的需要计算的yk+1隐含在等式两侧隐含在等式两侧,称其为称其为隐式公式隐式公式.计计算公式依次可抽象写成算公式依次可抽象写成:显式单步法显式单步法:1(,)kkkkyyhx y h 隐式单步法隐式单步法:11(,)kkkkkyyhxyyh(5)(6)10该式右端项含有该式右端项含有1,ky因此若求因此若求1,ky需要解方程。需要解方程。注注:显式多步法显式多步法:111(,)kkkkkk lyyhx y yyh 隐式多步法
11、隐式多步法:111(,),1kkkkkk lyyhx yyyh kl 线性多步法线性多步法:11101,1llkk ik iiiiiyyhfkl 注注:(9)关于关于yk-i,fk-i都是线性的。都是线性的。其中其中011,101,ll 是独立于是独立于k和和f 的常数。的常数。10,若则(9)是显式是显式(右端不含有右端不含有yk+1)；10,若则则(9)是隐式的。是隐式的。(7)(8)(9)三、微分方程初值问题的数值解法讨论的问题三、微分方程初值问题的数值解法讨论的问题:1.方法构造方法构造2.误差分析误差分析3.稳定性稳定性11考虑考虑一、显式一、显式Euler方法方法(折线法折线法)1
12、.显式欧拉公式显式欧拉公式其中其中fk=f(xk,xk+1).(,),()dyf x yaxbdxy a01,0,1,1kkkkyknyyh f(1)(2)(10)设节点为设节点为:a=x0 x1 0,kkbahhxakhn00maxmax()(1)pkkkknkney xychp 常数常数c 独立于独立于h(与与h无关无关),称称(17)是是 p 阶方法，即阶方法，即O(hp+1)判断某种方法的阶数往往通过局部截断误差的阶数来确判断某种方法的阶数往往通过局部截断误差的阶数来确定，而局部截断误差的阶容易由公式来确定。定，而局部截断误差的阶容易由公式来确定。说明：说明：(19)设设y(x)是方程
13、是方程(1)(2)的准确解，的准确解，yk,k=0,1,n是单步法是单步法(17)的数值解，称的数值解，称 ek=y(xk)-yk,是单步法是单步法(18)在在xk点的点的整体截断误差整体截断误差。成立成立定义定义4302.局部截断误差与整体截断误差之间的关系局部截断误差与整体截断误差之间的关系若单步法若单步法(17)的局部截断误差是的局部截断误差是p+1阶的，即阶的，即111,(),0,1,1,pkkchhh kn(,),x u v h axb u0,0vhh c1独立于独立于h(不依赖于不依赖于h或与或与h无关无关)，且函数且函数(x,u,v,h)在区域在区域定理定理4上关于上关于u,v
14、满足条件满足条件12MuMv 则单步法则单步法(17)是是 p 阶阶方法。方法。(20)31推论推论:当当f 在在D上满足基本条件时上满足基本条件时,单步法单步法(17)的阶由的阶由局部截断误局部截断误差的阶差的阶来确定。来确定。结论：结论：一般情况下，若局部截断误差是一般情况下，若局部截断误差是 p+1阶的，则单步法是阶的，则单步法是p 阶方法。阶方法。121(,()()ppkkkx y xhO h说明：说明：可用可用Taylor展开法估计展开法估计1k的阶，即的阶，即1k 的主项的主项高阶项高阶项则对应的单步法是则对应的单步法是p 阶方法。阶方法。111()()(,(),(),),kkkk
15、kky xy xhxy xy xh 由于由于(21)323.定理的应用（讨论各种方法的阶数）定理的应用（讨论各种方法的阶数）显式显式Euler方法是一阶方法方法是一阶方法事实上事实上,当当(1),(2)的解的解y(x)二阶连续可导时二阶连续可导时,f(x,y)关于关于y满满足足Lipschitz条件，条件，Lipschitz常数为常数为L，其局部截断误差为，其局部截断误差为:11()()(,()kkkkky xy xhf x y x2111(),2!kkkkkyhxx将将y(xk+1)在在xk点展开可得点展开可得),(1kkkkyxhfyy 1()(,()kkkkyy xhf xy x局部截断
16、误差有界，存在局部截断误差有界，存在M，使得，使得221()22khhyM记记于是于是3311111111112()()()(,()(,)()(,()(,)(1)()(1)2kkkkkkkkkkkkkkkkkkkkkkkkky xyy xyyyy xhfxy xyhfxyy xyhfxy xfxyhhLy xyMhL221122110(1)(1)22(1)1 (1)(1)12(1)12kkknkkkhh MMhLhLhh MhLhMMhLhLhLL341()()1(1)(1)(1)1 2banL bahL baknhbahLhLehMeL35n 隐式隐式Euler方法是一阶方法方法是一阶方法1
17、111()()(,()kkkkky xy xhf xy x事实上，局部截断误差为事实上，局部截断误差为:),(111 kkkkyxhfyy将将 f(x,y)在在y用微分中值定理，有用微分中值定理，有1111111(,)(,()(,)()kkkkykkkf xyf xy xf xyy x231111()()()(,)()()kkkkykkkyy xhy xh y xhf xyy xO h则则将方程的解作将方程的解作Taylor展开展开231()()()()()2kkkkhy xy xhyxyxO h362311111()()(,)()()2kkkykkkhy xyy xhfxyy xO h 因此
18、因此23111231312()()()2()()2()()211(,)1(,).kkkkynynnnhfxhfhy xyyxO hhyxO hhyxO hx37n 梯形方法是二阶方法梯形方法是二阶方法将将 f(xk+1,yk+1)在在xkTaylor展开，有展开，有2314111()()()()24(,)()()2kkkkkykkkhhyy xhy xyxyxhfxyy xO h则则将方程的解作将方程的解作Taylor展开展开),(),(2111 kkkkkkkyxfyxfhyy2341()()()()()()26kkkkkhhy xy xhy xyxyxO h23111(,()()()()(
19、)()2kkkkkkhf xy xy xy xhy xyxO h1111111(,)(,()(,)()kkkkykkkf xyf xy xf xyy x将将 f(x,y)在在y用微分中值定理，有用微分中值定理，有383411111()()(,)()()122kknynnnhhy xyyxfxyy xO h 因此因此11134341()()()12(11(,)(12)/2)ykkknnky xyhyxO hhyxO hhfx39n 预估校正预估校正Euler方法是二阶方法方法是二阶方法事实上，局部截断误差为：事实上，局部截断误差为：111111111111111()(),(),(),()2(,(
20、)(,()(,()2(,)()()(,()2()1(,)2122kkkkkkkkkkkkkkkkkykkkkkkkkykkkkyhy xy xf xy xf xy xh f xy xhf xy xf xy xhf xy xhfxy xy xhf xy xfhfxf 梯梯梯欧31311(,)()2kkkxyhc h),(,(),(211kkkkkkkkyxhfyxfyxfhyy .)(12131hykk 梯梯.)(!2121hykk 欧欧4020max.kknech二阶方法二阶方法即得整体截断误差即得整体截断误差:结论：结论：一般情况下，若局部截断误差是一般情况下，若局部截断误差是 p+1阶的，
21、则单步法阶的，则单步法是是p 阶方法。阶方法。方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低,计算量大计算量大精度提高精度提高计算量大计算量大41五、五、单步法解的唯一性、收敛性和稳定性单步法解的唯一性、收敛性和稳定性u 唯一性唯一性11()()(,)(,)22kkkkhLg yg yf xyf xyhyy若若f(x,y)在在D上满足基本条件，上满足基本条件，f(x,y)关于关于y的的Lipschitz常数常数为为L时，只要时，只要hkL1，则由定理，则由定理1可知可知(13)确定了唯一的确定了唯一的yk+1。同理，只要同理，只
22、要hkL2，(14)确定了唯一的确定了唯一的yk+1.u 收敛性收敛性(以以(14)式为例说明式为例说明)当当2kh L 时时,以以y 为变量的函数为变量的函数:1()(,)(,)2kkkkkhg yyf x yf xy同样，从任意初始值出发同样，从任意初始值出发,迭代迭代(13)、(14)，都收敛到，都收敛到 1ky事实上，事实上，42在在 上关于上关于y 满足满足Lipschitz条件条件,且且L常数为常数为 ,()y12kLh由由压缩不动点定理压缩不动点定理得方程得方程1(,)(,)2kkkkkhyyf x yf xy有唯一不动点有唯一不动点 ,1ky而且从而且从 出发出发,迭代迭代(0
23、)1ky(1)()111(,)(,),2iikkkkkkkhyyf x yf xy0,1,1in(15)都收敛到都收敛到1.ky见非线性方程数见非线性方程数值解的定理值解的定理343u 稳定性稳定性例：例：考察初值问题考察初值问题 在区间在区间0,0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。()30()(0)1y xy xy 30 xye44定义定义若某算法在计算过程中任一步产生的误差在以后若某算法在计算过程中任一步产生的误差在以后的计算中都的计算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的。
24、讨论数值方法的稳定性讨论数值方法的稳定性,通常仅限于典型的试验方程通常仅限于典型的试验方程 其中其中 是复数。在复平面上是复数。在复平面上,当方法稳定时要求变量当方法稳定时要求变量 h的取的取值范围称为方法的值范围称为方法的绝对稳定域绝对稳定域,它与实轴的交集称为它与实轴的交集称为绝对稳绝对稳定区间定区间.yy45将将Euler方法应用于方程方法应用于方程y =y,得到得到 设在计算设在计算yn时产生误差时产生误差 n,计算值计算值 yn=yn+n,则则 n将对以后各将对以后各节点值计算产生影响节点值计算产生影响.记记 ym=ym+m,m n,由上式可知误差由上式可知误差 m满足方程满足方程
25、可见可见,若要若要|m|n|,必须且只须必须且只须|1+h|1,因此因此Euler法的绝对稳定域为法的绝对稳定域为|1+h|1,绝对稳定绝对稳定区间是区间是-2Re()h0.1(1)nnyh y1(1).(1),m nmmnhhmn 0-1-2ReImg对于前述例题，因为对于前述例题，因为 h=-3不属于不属于Euler方法的绝对稳定区方法的绝对稳定区间，故间，故Euler方法计算出现严重错误。方法计算出现严重错误。461211211nnhyyh类似前面分析类似前面分析,可知绝对稳定区域为可知绝对稳定区域为1212111hh由于由于Re()0,所以此不等式对任意步长所以此不等式对任意步长h恒成立恒成立,这是隐这是隐式公式的优点式公式的优点.解出解出yn+1得得 隐式单步方法可类似讨论隐式单步方法可类似讨论.如对梯形公式使用试验方程如对梯形公式使用试验方程,有有 11()2nnnnhyyyy47结论 单步显式方法的稳定性与步长密切相关,在一种步长下是稳定的差分公式,取大一点步长就可能是不稳定的.收敛性是反映差分公式本身的截断误差对数值解的影响;稳定性是反映计算过程中舍入误差对数值解的影响.只有即收敛又稳定的差分公式才有实用价值.一些常用方法的绝对稳定区间一些常用方法的绝对稳定区间48作业作业1思考题思考题1中的（中的（a）（f）2习题习题3
