数值求解常微分方程PPT课件

方程中包含未知函数的导数或者微分的方程叫做微分方程, ,求解微分方程必须附加某种定解条件. .微分方程和定解条件一起组成定解问题, ,定解条件分为初始条件( (初值问题) )和边界条件( (边值问题) )两种. .未知函数为一元的微分方程叫做常微分方程, ,未知函数为多元函数, ,叫做偏微分方程. .微分方程中导数的最高阶叫做微分方程的阶. .本章主要讨论一阶常微分方程. .1.引言1第1页/共51页 常微分方程在物理学中具有重要意义, ,牛顿第二定律就表示为一组二阶微分方程. .前面提到的著名的三体问题就是由九个二阶常微分方程组成的方程组; ;在量子力学和电动力学中, ,由薛定谔及麦克斯韦方程分离变量得到本征值问题也是二阶常微分方程的求解问题. .2第2页/共51页现考虑上述一阶常微分方程初值问题. .在常微分方程中，我们已经掌握了一些典型方程的解析解法。但解析方法只能解一些特殊类型的方程, ,求解从实际当中得出来的微分方程只能用数值方法求近似解. .初值问题也就是由已知点(x(x0 0,y,y0 0) )求解在某些点x x上满足一定精度的y y的近似解。00(,)()d yfxyd xyxy3第3页/共51页即使对有些函数可以求出解析解. .但实际应用时, ,仍需求它在某一点的值, ,如果解析解的形式比较复杂, ,还是需要数值方法来求得. .例 y=1-2xyy=1-2xy y(0)=0 y(0)=0其解析解如下, ,要计算具体的函数值y,y,还需要积分的方法, ,不如直接用数值法求解. .xdttxy022)exp()exp(4第4页/共51页2.1 2.1 欧拉法( (折线法) ) 若将函数y(x)y(x)在点x xi i处的导数y(xy(xi i) )用( (向前) )差商来表示，即111()()()()()iiiiiiiy xy xy xxxy xy xh2.欧拉法再用yi近似地代替y(xi)，则初值问题就化为500(,)()d yfxyd xyxy第5页/共51页100( ,)( )0,1,2,iiiiyyhf x yyy xi 上式就是上式就是欧拉公式欧拉公式, ,它有很明显的几何意义它有很明显的几何意义: :我我们知道初值问题中的微分方程的解是平面上的一条们知道初值问题中的微分方程的解是平面上的一条曲线，这条曲线上任意点曲线，这条曲线上任意点(x(x，y)y)的切线的斜率的切线的斜率y(x)y(x)等于函数等于函数f(xf(x，y)y)在这点的值，而初值问题在这点的值，而初值问题的解是过点的解是过点(x(x0 0,y,y0 0) )的一条特定的曲线。的一条特定的曲线。 6第6页/共51页过点过点(x(x0 0,y,y0 0) )以以f(xf(x0 0,y,y0 0) )为斜率的方程为为斜率的方程为当当x=xx=x1 1时得下式时得下式, ,取取y1y1作为解作为解y(x1)y(x1)的近似解的近似解过过(x(x1 1,y,y1 1) )以以f(xf(x1 1,y,y1 1) )为斜率的方程为为斜率的方程为)(,(0000 xxyxfyy)(,(00001xxyxfyy)(,(1111xxyxfyy7第7页/共51页若已求得点若已求得点(x(xn n,y,yn n),),过这点过这点, ,以以f(xf(xn n,y,yn n) )为斜率的直线为斜率的直线当当x=xx=xn+1n+1时得时得: :取取y(xy(xn+1n+1)y)yn+1n+1; ;这样每个这样每个x xn n求出对应数值解求出对应数值解y yn n)(,(nnnnxxyxfyy)(,(11nnnnnnxxyxfyy8第8页/共51页还可使用泰勒展开证明欧拉公式：还可使用泰勒展开证明欧拉公式：对对y(xy(xn+1n+1) )在在x xn n处泰勒展开有处泰勒展开有舍去余项可得欧拉公式舍去余项可得欧拉公式! 2/ )( )( )()()(21nnnnnyhxhyxyhxyxy9第9页/共51页例1 1: :用欧拉法求初值问题的数值解( (取h=0.1)h=0.1)。2(0)0yxyy 解解 由已知条件可得欧拉法迭代公式由已知条件可得欧拉法迭代公式2( , )(0)0,0.1f x yxyyh2100.1 ()0,0,1,2,iiiiyyxyyi10第10页/共51页11第11页/共51页f(x,y)=x-y*y;float y=0;float h=0.1;float a=0.0;float b=1.0;for(a=0;a=b;a+h)y=y+h*f(a,y); printf(“x=%f,y=%f”,a,y); 12第12页/共51页隐式欧拉公式若将函数y(x)y(x)在点x xi i处的导数y(xy(xi i) )用( (向后) )差商来表示，即欧拉公式变为: :hxyxyxxxyxyxyiiiiiii)()()()()( 1111),(111iiiiyxhfyy13第13页/共51页两步欧拉公式 为了改善精度, ,将函数y(x)y(x)在点x xi i处的导数y(xy(xi i) )用中心差商来表示，即欧拉公式变为: :误差正比于h h3 3, ,是二阶精度hxyxyxxxyxyxyiiiiiii2)()()()()( 111111),(211iiiiyxhfyy14第14页/共51页2.2 2.2 改进的欧拉法欧拉法虽然形式简单，计算方便，但比较粗糙，精度也低。特别当y=y(x)y=y(x)的曲线曲率较大时，欧拉法的效果更差。为了达到较高精度的计算公式，对欧拉法进行改进，将在一点(x(xi i，y yi i) )的切线斜率f(xf(xi i，y yi i) )用两点的平均斜率来代替，即 111( ,) ( ,)(,)2iiiiiif x yf x yf xy15第15页/共51页 代入欧拉公式得 111 ( ,)(,)20,1,2,1iiiiiihyyf x yf xyin 这样得到的点列仍为一折线，只是用平均斜率这样得到的点列仍为一折线，只是用平均斜率. .来来代替原来一点处的斜率。称为代替原来一点处的斜率。称为改进的欧拉公式改进的欧拉公式。 16第16页/共51页 不难发现，欧拉公式是关于y yi+1i+1的显式，只要已知y yi i，经一次计算可立即得到y yi+1i+1的值; ;而改进的欧拉公式中的y yi+1i+1以隐式给出，且y yi+1i+1含在函数f(xf(xi+1i+1 ,y ,yi+1i+1) )中，因此，通常用迭代法求解。具体做法是：先用欧拉公式求出一个y y(0)(0)i+1i+1作为初始近似，然后再用改进的欧拉公式进行迭代，即17第17页/共51页(0)1(1)( )111( ,) ( ,)(,)20,1,2,iiiikkiiiiiiyyhf x yhyyf x yf xyk 直到满足 则 (1)( )11(1)11kkiikiiyyyy(为预给精度为预给精度) )18第18页/共51页2.3 2.3 预估校正法 改进的欧拉公式在实际计算时要进行多次迭代，因而计算量较大。所谓预估校正法，就是先用欧拉公式算出y yi+1i+1的预估值y y(p)(p)i+1i+1，然后再用改进的欧拉公式式进行一次迭代便得到校正值y y(c)(c)i+1i+1，即19第19页/共51页预估：预估：校正：校正：并取并取 ( )11ciiyy 虽然上式仅迭代一次，但因进行了预先估计，故精度也有较大的提高。()1( )()111( ,) ( ,)(,2piiiicpiiiiiiyyhf x yhyyf x yf xy20第20页/共51页121112(,)(,)()20,1,2,iiiiiikf xykf xh yhkhyykki在实际计算时，还常常将上式写成下列形式： 误差估计: 欧拉公式的截断误差为一阶的,误差正比于h2;改进的欧拉公式和预估校正公式的截断误差为二阶的,误差正比于h321第21页/共51页22第22页/共51页 例2 2: :求解初值问题 2(0)1,01,0.1xyyyyxh解:分别用欧拉公式和改进的欧拉公式进行计算.这里改进的欧拉公式的形式为:23课后作业：请使用课后作业：请使用C语言程序实现此求导！语言程序实现此求导！第23页/共51页121111222()()2iiiiiiiixkyyxhkyhkyhkhyykk 其解析解为 21yx24第24页/共51页25第25页/共51页26第26页/共51页27第27页/共51页 具有二阶以上精度的方法称为龙格- -库塔法，做为一种训练，我们给出一种直接推导法，再给出一种比较容易推广到高阶的推导法，希望同学们自己体会这里介绍的方法和思路！3. 龙格-库塔法28第28页/共51页29第29页/共51页30第30页/共51页31第31页/共51页32第32页/共51页33第33页/共51页34第34页/共51页由四阶龙格库塔公式计算得到的例由四阶龙格库塔公式计算得到的例2 2的结果的结果35第35页/共51页36W. C. Gear 的的Numerical Initial Value problems In Ordinary Differential Equations一文一文给出了一般龙格库塔法给出了一般龙格库塔法的推导，以及给出了更多的四阶龙格库的推导，以及给出了更多的四阶龙格库塔的形式，大家可以参考一下！塔的形式，大家可以参考一下！第36页/共51页 我们已经知道，初值问题等价于积分方程，即11()()( , ( )nnxnnxy xy xf x y x dx对积分式分别采用矩形公式和梯形公式可得到欧拉公式和改进的欧拉公式，截断误差分别为O(h2)和O(h3)。4. 阿达姆斯方法 37第37页/共51页 为此，我们自然可以想到，若用更高次的插值多项式来代替f(xf(x，y)y)，则所得公式的精度会更高。这就是线性多步法的起源思想。 本章前面介绍的方法称为单步法，因为在计算y yi+1i+1时，只用到前面yiyi的值。而对于线性多步法是要利用前面已经算出的若干个值y yi-i-k k，y yi-1i-1，y yi i来求y yi+1i+1。38第38页/共51页现用k k次多项式P Pk k(x)(x)来代替f(xf(x，y(x) y(x) 1( )iixkkxRRx dx111()()( )( )iiiixxiikkxxy xy xP x dxRx dx舍去余项 并设yi=y(xi)，而yi+1为y(xi+1)的近似值，于是可得到线性多步法的计算公式11( )iixiikxyyP k dx39第39页/共51页阿达姆斯(Adams)显式:为了易于在电子计算机上实现，常将式中的 用各点的已知函数值表示。1100( 1)iiqxmmiiixmqmimimtyydxfmyhf 特别，当特别，当q=2,3q=2,3时，有时，有: :其中其中f fi i=f(x=f(xi i,y,yi i) ) )9375955(24)51623(123211211iiiiiiiiiiiffffhyyfffhyy40第40页/共51页 阿达姆斯隐式的数值公式。 例如，当q=2q=2时1111112(58)12(9195)24iiiiiiiiiiihyyfffhyyffff当当q=3q=3时时 41第41页/共51页阿达姆斯预估校正公式 我们常把阿达姆斯显式及隐式联立使用，即构造所谓阿达姆斯预估校正公式。现以q=2q=2为例构造预估校正公式()111( )111(23165)12(58)12piiiiiciiiiihyyfffhyyfff并取 ( )11ciiyy42第42页/共51页 与同阶的龙格库塔方法相比较，阿达姆斯方法计算量小，公式简单，程序易于实现。但它的主要缺点是不能自动开始，开始的前几个值要依赖于其它方法获得。这里介绍两种计算开始值的方法。 (1) (1) 用单步法中的数值方法求出开始值。 (2) (2) 使用y(x)y(x)的泰勒展开式43第43页/共51页 设二阶线性常微分方程的边值问题为 y+p(x)y+q(x)y=f(x)y+p(x)y+q(x)y=f(x) y(a)= y(a)=，y(b)=y(b)=，axbaxb 其中p(x)p(x)，q(x)q(x)，f(x)f(x)为区间a a，b b上足够光滑的已知函数，且q(x)0q(x)0，、为已知常数。5. 二阶线性常微分方程边值问题 44第44页/共51页 在上述条件下，边值问题存在连续可微的解，且是唯一的。若采用差分方法来解边值问题，其基本步骤是： 将区间a a，b b“离散化”，即给a a，b b一个分划，此分划常考虑等距; ; 对每一个基点，将各阶导数用差商来近似表示，将微分方程转化为差分方程，进而转化为线性代数方程组; ; 解线性代数方程组，求得各基点上的近似解。45第45页/共51页现具体给出求解边值问题的方法步骤。 首先将区间a,ba,b进行等距分划，即令 x xi i= a + i h= a + i h ， i = 0 , 1 , 2 , , n ;i = 0 , 1 , 2 , , n ; 其 中 h = ( b -h = ( b -a)/n;xa)/n;x0 0=a,x=a,xn n=b=b其次在各基点x xi i上，将y,yy,y用差商来近似表示。这里要求有相同阶数的截断误差，以保证精度协调. .211()()()()2iiiyxyxyxOhh2112()2()()()()iiiiy xy xy xyxO hh46第46页/共51页代入二阶线性常微分方程的边值问题得到近似差 分 方 程 为 : : 其 中 p pi i= p ( x= p ( xi i) , q) , qi i= q ( x= q ( xi i) , ) , f fi i=f(x=f(xi i) )11112022,1,2,1iiiiiiiiinyyyyypq yfhhyyin略去上两式的截断误差略去上两式的截断误差O(hO(h2 2) )，用，用y yi i代替代替y(xy(xi i) )有有 11112()22()iiiiiiiyyyxhyyyyxh47第47页/共51页 定理: :若a ai i0 0，c ci i0 0，-b-bi iaai i+c+ci i，i=1i=1，2 2，n-1n-1，则二阶差分方程收敛 a ai iy yi-1i-1+b+bi iy yi i+c+ci iy yi+1i+1=d=di i y y0 0=，y yn n=，i=1i=1，2 2，n-1n-1 2212212,1, 2 ,1iiiiiiiihapbhqhcpdhfin 48第48页/共51页例3: 用差分法解边值问题 y-y=x; y(0)=0，y(1)=1， 0 x1，h=0.149第49页/共51页50第50页/共51页感谢您的观看！第51页/共51页