第10章_常微分方程数值解.

计算方法课件1结束第8章常微分方程的数值解科学研究和工程技术中的问题往往归结为求某个常微分方程的定解问题.常微分方程的理论指出，很多方程的定解虽然存在，但可能十分复杂难于计算，也可能不能用简单的初等函数表示,因此常求其能满足精度要求的近似解.常微分方程的数值解法常用来求近似解，由于它提供的算法能通过计算机便捷地实现，因此近年来得到迅速的发展和广泛的应用常微分方程数值解法的特点是：对求解区间进行剖分，然后把微分方程离散成在节点上的近似公式或近似方程，最后结合定解条件求出近似解.因此数值解法得到的近似解是一个离散的函数表.§8.1引言计算方法课件2结束8.1.1基本知识复习一阶常微方程的初值问题：)2()()1(),(0yaybxayxfy其中f(x,y)是已知函数，(2)是定解条件，常微分方程的理论指出：当f(x,y)定义在区域G=(a≤x≤b,｜y｜＜∞)，若存在正的常数L，使：｜f(x,y1)–f(x,y2)｜≤L｜y1-y2|对任意的x∈[a,b]和任意的y1,y2成立，则称f(x,y)对y满足李普希兹(Lipschitz)条件,L称李普希兹常数.若f(x,y)在区域G连续，关于y满足李普希兹条件，则上述一阶常微分方程的初值问题的解存在且唯一.我们以下的讨论，都在满足上述条件下进行.计算方法课件3结束8.1.2一阶常微分方程组和高阶常微分方程一阶常微分方程组常表述为：)()()(),,,(),,,(111111bxaayayyyxfyyyxfymmmmmm写成向量形式为Tmyyaybxayxfy),,(,)(),(100计算方法课件4结束还有高阶常微分方程定解问题，如二阶定解问题：)()(),,(ayayyyxfy可化作一阶常微方程组的定解问题：本章着重讨论一阶常微方程初值问题的各种数值解法.这些解法都可以写成向量形式而用于一阶常微分方程组的初值问题.这样，也就解决了高阶方程的定解问题.)(,)(),,(),(auayuyxfuyxuy计算方法课件5结束§8.2欧拉(Euler)方法对问题)2.8()()1.8(),(0yaybxayxfy最简单而直观的方法是欧拉方法.欧拉方法在精度要求不高时，仍不失为一实用方法.同时通过欧拉方法的讨论，容易弄清常微方程初值问题数值解法的一些基本概念和构造方法的思路.8.2.1欧拉方法的导出把区间［a,b］分为n个小区间，取步长h=(b-a)/n，节点xi=x0+ih,i=0,1,2,…，n，其中x0=a，又设y(x)为上述问题的解.计算方法课件6结束把(8.1)式两端积分，积分区间取［xi,xi+1］)3.8()()(),(111iixxxxxyxydxyxfdxyiiii记g(x)=f(x,y(x))对(8.3)用左矩形公式(7.3)2)(2)()()()(abgagabdxxgba得到：iiiiiRxyxhfxyxy))(,()()(1.1,,1,0),(2)(12nixxhyRiii舍去Ri，并令yi=y(xi),得到：)4.8(1,,1,0),(1niyxhfyyiiii这就是欧拉公式.计算方法课件7结束欧拉公式还可由泰勒展开得到.假定y(x)二阶连续可导，把y(xi+1)在xi点展开：21!2)()()()(hyhxyxyxyiii舍去h2项，并令yi=y(xi)，注意到y’(xi)=f(xi,yi),也得到)4.8(1,,1,0),(1niyxhfyyiiii欧拉公式有明显的几何意义，如图8.1，过点(x0,y0)的曲线是解y(x).欧拉方法是在(x0,y0)作y(x)的切线，它与直线x=x1交于(x1,y1)，过(x1,y1)作过此点的积分曲线的切线，又与x=x2交于(x2,y2)，…如此下去，得到一条折线，欧拉方法就是用这条折线近似地代替曲线y(x)，故欧拉方法有时也称欧拉折线法.8.2.2欧拉隐式公式和欧拉中点公式对(8.3)式的右端用右矩形公式代替积分，可得到计算方法课件8结束),()(2))(,()()(12111iiiiiixxyhxyxhfxyxy)5.8(,1,,1,0),(111niyxhfyyiiii推出：不难看出，(8.5)式中公式两端都含有yi+1,一般情况不能由yi的值计算yi+1,而需要解方程，故称为欧拉隐式公式.),()(3))(,(2)()(1311iiiiiixxyhxyxhfxyxy)6.8(,1,,2,1),(211niyxhfyyiiii推出：称为欧拉中点公式.容易看出,中点公式计算yi+1时,不仅需要yi的值，还需要yi-1的值.把(8.3)式积分区间改为[xi-1,xi+1]，并对右端用中矩形公式代替积分，可得到：计算方法课件9结束凡是从已知(或已算出)的y0,y1,…,yi能直接从公式算出yi+1的公式称显式，否则称隐式公式(8.4),(8.6)都是显式，(8.5)是隐式在计算yi+1时，只需要yi的值，则称公式为单步法；若除yi之外还需要以前的yi-1等多个值，则称多步法公式.(8.4)、(8.5)是单步法，(8.6)是多步法，确切说是二步法.8.2.3局部截断误差和方法的阶以上三个公式都有误差.假定yi,yi-1的值没有误差，用公式计算yi+1时产生的误差称为局部截断误差，也称为单步截断误差,假定yi+1以前的y值都没有误差，局部截断误差表示为：11)(iiyxyR对欧拉公式：)(2)(221hOhyR计算方法课件10结束对欧拉隐式公式：)(2)(222hOhyR对欧拉中点公式的局部截断误差：)(3)(333hOhyR凡是局部截断误差为O(hp+1)的方法称为p阶方法，这样欧拉方法和欧拉隐式方法称为一阶方法，而欧拉中点方法称为二阶方法.一般讲高阶方法比低阶方法要精确.计算方法课件11结束8.2.4梯形公式及其预估校正法对积分式(8.3)的右端用梯形积分公式近似代替，得把y(xi+1)记为yi+1，并略去余项-h3y(ξ)/12得到梯形公式：),()(12))(,())(,(2)()(13111iiiiiiiixxyhxyxfxyxfhxyxy)8.8(),(),(2111iiiiiiyxfyxfhyy它是单步隐式方法，也是一个二阶方法.又是隐式方法，一般求yi+1需要解方程，常用迭代法解此方程，初值y(0)i+1由显式的欧拉公式给出,于是完整的公式写为:),()0(1iiiiyxhfyy)9.8(1,,1,0,,2,1,0),(),(2)(11)1(1nikyxfyxfhyykiiiiiki计算方法课件12结束它称为梯形公式的预估校正法，用低精度的显式欧拉公式给出下一步yi+1的初值y(0)i+1，再用较高精度的隐式方法梯形公式进行迭代,一般迭代一两次就能达到精度要求.迭代法的收敛条件：||2|),(),(|2||)1(1)(1)1(11)(11)(1)1(1kikikiikiikikiyyLhyxfyxfhyy其中L为Lipschitz常数。当120Lh时迭代过程收敛。即迭代过程中误差没有扩大。计算方法课件13结束例1用欧拉方法，隐式欧拉方法和欧拉中点公式计算1)0(10yxyy的近似解，取步长h=0.1，并与精确值比较解:欧拉方法的算式为：1.9,,1,0)1(01yihyyii欧拉隐式方法在本题可解出方程，不必迭代，公式为：1.9,,1,0101yihyyii欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式1.9,,2,12011yihyyyiii计算方法课件14结束本题精确解为y=e-x,计算结果（略）例2用欧拉公式和梯形公式的预估校正法计算：的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较.方程的解析解为:1)0(102yxyxyyxy21解:本例中欧拉公式为：1.9,,1,02.01.121.001yiyxyyxyyyiiiiiiii计算方法课件15结束梯形公式只校正一次的格式为12205.02.01.10)0(11)0(11)0(1yyxyyxyyyyxyyiiiiiiiiiiii结果（略）计算方法课件16结束§8.3龙格-库塔(Runge-Kutta)法欧拉方法是显式的一步法,使用方便,但精度较低.本节将构造出高精度的显式一步法：龙格-库塔法,简称R-K法.8.3.1二阶R-K法欧拉法的公式为：yi+1=yi+hf(xi,yi)i=0,1,2,…,n-1决定其精度的是函数f(xi,yi).如能改进这个函数,就可能提高公式的精度.为此把公式改写成：yi+1=yi+h(xi,yi,h)i=0,1,2,…,n-1(8.10)选择函数(xi,yi,h),一种方法是用若干个点的函数值的线性组合代替(xi,yi,h),如：计算方法课件17结束1111,,3,2),(),(),,(jlljlijijiipjjjiipjKbhyhaxfKyxfKKchyx其中cj,aj,bjl是待定参数,aj和bjl满足pjbajljlj,,3,211以上方法称为p级R-K法,选择cj,aj和bjl,可能使以上方法为p阶方法.显然欧拉法就是一阶R-K法.计算方法课件18结束二级R-K法的形式是:),(),(1222122111hKayhaxfKyxfKKcKchyyiiiiii此时)11.8())()(,()()))(,()(,())(,(),,(22212221iiiiiiiiiiiixyhaxyhaxfcxycxyxhfaxyhaxfcxyxfchyx由二元函数的泰勒展开：)()())(,())()(,(22222hOfxyhafhaxyxfxyhaxyhaxfyixiiiii其中所有的偏导数都是它们在点(xi,y(xi))的值,下同计算方法课件19结束又由于：)()())(,()(xyffxyxyxfxyyx所以)()()())()(,(2222hOxyhaxyxyhaxyhaxfiiiii代入(8.11))()()(),,(2221hOxyhaxycchyxiiii代入(8.10))()()()(3222211hOxyhcaxyhccxyyiiii而Taylor展开式)()(2)()()(321hOxyhxyhxyxyiiii计算方法课件20结束二式相减,得局部截断误差)()(21)(1)(322221111hOxyhcaxyhccyxyRiiiii:021,012221得令cacc2112221cacc只要c1,c2,a2满足以上方程,就得到一个二阶的R-K法.这是一个不定方程,有无穷多解.比如：计算方法课件21结束(1)取c1=c2=1/2,a2=1得)12.8(),(),(2121211hKyhxfKyxfKKKhyyiiiiii这实际上是(8.9)公式,即梯形公式的预估-校正公式只迭代一次的形式,通常称为改进的欧拉法.(2)取c1=0,c2=1,a2=1/2得:)13.8(2,2)