《数值求解常微分方程》由会员分享，可在线阅读，更多相关《数值求解常微分方程（51页珍藏版）》请在金锄头文库上搜索。
1、 数数值求解常微分方程求解常微分方程1.1.引言引言2.2.欧拉法欧拉法3.3.龙格龙格- -库塔法库塔法4.4.线性多步法线性多步法5.5.高阶微分方程高阶微分方程 常微分方程数值解法常微分方程数值解法1 数数值求解常微分方程求解常微分方程 方方程程中中包包含含未未知知函函数数的的导导数数或或者者微微分分的的方方程程叫叫做做微微分分方方程程, ,求求解解微微分分方方程程必必须须附附加加某某种种定定解解条条件件. .微微分分方方程程和和定定解解条条件件一一起起组组成成定定解解问问题题, ,定定解解条条件件分分为为初初始始条条件件( (初初值值问问题题) )和和边边界界条条件件( (边边值值问问
2、题题) )两两种种. .未未知知函函数数为为一一元元的的微微分分方方程程叫叫做做常常微微分分方方程程, ,未未知知函函数数为为多多元元函函数数, ,叫叫做做偏偏微微分分方方程程. .微微分分方方程程中中导导数数的的最最高高阶阶叫叫做做微微分分方方程程的的阶阶. .本本章主要讨论章主要讨论一阶常微分方程一阶常微分方程. .1.1.引言引言2 数数值求解常微分方程求解常微分方程 常常微微分分方方程程在在物物理理学学中中具具有有重重要要意意义义, ,牛牛顿顿第第二二定定律律就就表表示示为为一一组组二二阶阶微微分分方方程程. .前前面面提提到到的的著著名名的的三三体体问问题题就就是是由由九九个个二二阶
3、阶常常微微分分方方程程组组成成的的方方程程组组; ;在在量量子子力力学学和和电电动动力力学学中中, ,由由薛薛定定谔谔及及麦麦克克斯斯韦韦方方程程分分离离变变量量得得到到本本征征值值问问题题也也是是二二阶常微分方程的求解问题阶常微分方程的求解问题. .3 数数值求解常微分方程求解常微分方程现现考考虑虑上上述述一一阶阶常常微微分分方方程程初初值值问问题题. .在在常常微微分分方方程程中中，我我们们已已经经掌掌握握了了一一些些典典型型方方程程的的解解析析解解法法。但但解解析析方方法法只只能能解解一一些些特特殊殊类类型型的的方方程程, ,求求解解从从实实际际当当中中得得出出来来的的微微分分方方程程只
4、只能能用用数数值值方方法法求求近近似似解解. .初初值值问问题题也也就就是是由由已已知知点点(x(x0 0,y,y0 0) )求求解解在在某某些些点点x x上满足一定精度的上满足一定精度的y y的近似解。的近似解。4 数数值求解常微分方程求解常微分方程即即使使对对有有些些函函数数可可以以求求出出解解析析解解. .但但实实际际应应用用时时, ,仍仍需需求求它它在在某某一一点点的的值值, ,如如果果解解析析解解的的形形式式比比较较复复杂杂, ,还是需要数值方法来求得还是需要数值方法来求得. .例例 y=1-2xyy=1-2xy y(0)=0 y(0)=0其其解解析析解解如如下下, ,要要计计算算具
5、具体体的的函函数数值值y,y,还还需需要要积积分分的方法的方法, ,不如直接用数值法求解不如直接用数值法求解. .5 数数值求解常微分方程求解常微分方程2.1 2.1 欧拉法欧拉法( (折线法折线法) ) 若若将将函函数数y(x)y(x)在在点点x xi i处处的的导导数数y(xy(xi i) )用用( (向向前前) )差商来表示，即差商来表示，即2.2.欧拉法欧拉法再用再用y yi i近似地代替近似地代替y(xy(xi i) )，则初值问题就化为，则初值问题就化为6 数数值求解常微分方程求解常微分方程 上上式式就就是是欧欧拉拉公公式式, ,它它有有很很明明显显的的几几何何意意义义: :我我们
6、们知知道道初初值值问问题题中中的的微微分分方方程程的的解解是是平平面面上上的的一一条条曲曲线线，这这条条曲曲线线上上任任意意点点(x(x，y)y)的的切切线线的的斜斜率率y(x)y(x)等等于于函函数数f(xf(x，y)y)在在这这点点的的值值，而而初初值值问问题题的解是过点的解是过点(x(x0 0,y,y0 0) )的一条特定的曲线。的一条特定的曲线。 7 数数值求解常微分方程求解常微分方程过点过点(x(x0 0,y,y0 0) )以以f(xf(x0 0,y,y0 0) )为斜率的方程为为斜率的方程为当当x=xx=x1 1时得下式时得下式, ,取取y1y1作为解作为解y(x1)y(x1)的近
7、似解的近似解过过(x(x1 1,y,y1 1) )以以f(xf(x1 1,y,y1 1) )为斜率的方程为为斜率的方程为8 数数值求解常微分方程求解常微分方程若已求得点若已求得点(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 n9 数数值求解常微分方程求解常微分方程还可使用泰勒展开证明欧拉公式：还可使用泰勒展开证明欧拉公式：对对y(xy(xn+1n+1) )在在x
8、 xn n处泰勒展开有处泰勒展开有舍去余项可得欧拉公式舍去余项可得欧拉公式10 数数值求解常微分方程求解常微分方程例例1 1: :用欧拉法求初值问题的数值解用欧拉法求初值问题的数值解( (取取h=0.1)h=0.1)。 解解 由已知条件可得欧拉法迭代公式由已知条件可得欧拉法迭代公式11 数数值求解常微分方程求解常微分方程12 数数值求解常微分方程求解常微分方程f(x,y)=x-y*yf(x,y)=x-y*y; ;float float y=0y=0; ;float float h=0.1h=0.1; ;float float a=0.0a=0.0; ;float bfloat b=1.0=1.
9、0; ;for(a=0;a=b;a+h)for(a=0;a=b;a+h) y=y+h*f(a,y);y=y+h*f(a,y); printf(“x=%f,y=%f”,a,yprintf(“x=%f,y=%f”,a,y);); 13 数数值求解常微分方程求解常微分方程隐式欧拉公式隐式欧拉公式若若将将函函数数y(x)y(x)在在点点x xi i处处的的导导数数y(xy(xi i) )用用( (向向后后) )差商来表示，即差商来表示，即欧拉公式变为欧拉公式变为: :14 数数值求解常微分方程求解常微分方程两步欧拉公式两步欧拉公式 为为了了改改善善精精度度, ,将将函函数数y(x)y(x)在在点点x
10、xi i处处的的导导数数y(xy(xi i) )用用中中心心差差商商来来表表示示，即即欧欧拉拉公公式式变变为为: :误误差差正正比比于于h h3 3, ,是是二二阶阶精精度度15 数数值求解常微分方程求解常微分方程2.2 2.2 改进的欧拉法改进的欧拉法欧拉法虽然形式简单，计算方便，但比较粗糙，欧拉法虽然形式简单，计算方便，但比较粗糙，精度也低。特别当精度也低。特别当y=y(x)y=y(x)的曲线曲率较大时，欧的曲线曲率较大时，欧拉法的效果更差。为了达到较高精度的计算公式，拉法的效果更差。为了达到较高精度的计算公式，对欧拉法进行改进，将在一点对欧拉法进行改进，将在一点(x(xi i，y yi
11、i) )的切线斜的切线斜率率f(xf(xi i，y yi i) )用两点的平均斜率来代替，即用两点的平均斜率来代替，即 16 数数值求解常微分方程求解常微分方程 代入欧拉公式得代入欧拉公式得 这这样样得得到到的的点点列列仍仍为为一一折折线线，只只是是用用平平均均斜斜率率. .来来代替原来一点处的斜率。称为代替原来一点处的斜率。称为改进的欧拉公式改进的欧拉公式。 17 数数值求解常微分方程求解常微分方程 不不难难发发现现，欧欧拉拉公公式式是是关关于于y yi+1i+1的的显显式式，只只要要已已知知y yi i，经经一一次次计计算算可可立立即即得得到到y yi+1i+1的的值值; ;而而改改进进的
12、的欧欧拉拉公公式式中中的的y yi+1i+1以以隐隐式式给给出出，且且y yi+1i+1含含在在函函数数f(xf(xi+1i+1 ,y,yi+1i+1) )中中，因因此此，通通常常用用迭迭代代法法求求解解。具具体体做做法法是是：先先用用欧欧拉拉公公式式求求出出一一个个y y(0)(0)i+1i+1作作为为初初始始近近似似，然然后后再再用用改改进进的的欧欧拉拉公公式式进行迭代，即进行迭代，即18 数数值求解常微分方程求解常微分方程 直到满足直到满足 则则 (为预给精度为预给精度) )19 数数值求解常微分方程求解常微分方程2.3 2.3 预估校正法预估校正法 改改进进的的欧欧拉拉公公式式在在实实
13、际际计计算算时时要要进进行行多多次次迭迭代代，因因而而计计算算量量较较大大。所所谓谓预预估估校校正正法法，就就是是先先用用欧欧拉拉公公式式算算出出y yi+1i+1的的预预估估值值y y(p)(p)i+1i+1，然然后后再再用用改改进进的的欧欧拉拉公公式式式式进进行行一一次次迭迭代代便便得到校正值得到校正值y y(c)(c)i+1i+1，即，即20 数数值求解常微分方程求解常微分方程预估：预估：校正：校正：并取并取 虽然上式仅迭代一次，但因进行了预虽然上式仅迭代一次，但因进行了预先估计，故精度也有较大的提高。先估计，故精度也有较大的提高。21 数数值求解常微分方程求解常微分方程在实际计算时，还
14、常常将上式写成下列形式：在实际计算时，还常常将上式写成下列形式： 误差估计误差估计: : 欧拉公式的截断误差为一阶的欧拉公式的截断误差为一阶的, ,误差正比于误差正比于h h2 2; ;改进的欧拉公式和预估校正公式的截断误差为二阶的改进的欧拉公式和预估校正公式的截断误差为二阶的, ,误误差正比于差正比于h h3 322 数数值求解常微分方程求解常微分方程23 数数值求解常微分方程求解常微分方程 例例2 2: :求解初值问题求解初值问题 解解: :分别用欧拉公式和改进的欧拉公式进行计算分别用欧拉公式和改进的欧拉公式进行计算. .这里改进的欧拉公式的形式为这里改进的欧拉公式的形式为: :24课后作
15、业：请使用课后作业：请使用C语言程序实现此求导！语言程序实现此求导！ 数数值求解常微分方程求解常微分方程 其解析解为 25 数数值求解常微分方程求解常微分方程26 数数值求解常微分方程求解常微分方程27 数数值求解常微分方程求解常微分方程28 数数值求解常微分方程求解常微分方程 具具有有二二阶阶以以上上精精度度的的方方法法称称为为龙龙格格- -库库塔塔法，法，做做为为一一种种训训练练，我我们们给给出出一一种种直直接接推推导导法法，再再给给出出一一种种比比较较容容易易推推广广到到高高阶阶的的推推导导法法，希希望望同学们自己体会这里介绍的方法和思路！同学们自己体会这里介绍的方法和思路！3. 3.
16、龙格龙格- -库塔法库塔法29 数数值求解常微分方程求解常微分方程30 数数值求解常微分方程求解常微分方程31 数数值求解常微分方程求解常微分方程32 数数值求解常微分方程求解常微分方程33 数数值求解常微分方程求解常微分方程34 数数值求解常微分方程求解常微分方程35 数数值求解常微分方程求解常微分方程由四阶龙格库塔公式计算得到的例由四阶龙格库塔公式计算得到的例2 2的结果的结果36 数数值求解常微分方程求解常微分方程37W. C. Gear 的的Numerical Initial Value problems In Ordinary Differential Equations一文一文给出
17、了一般龙格库塔法给出了一般龙格库塔法的推导，以及给出了更多的四阶龙格库的推导，以及给出了更多的四阶龙格库塔的形式，大家可以参考一下！塔的形式，大家可以参考一下！ 数数值求解常微分方程求解常微分方程 我们已经知道，初值问题等价于积分方程，即我们已经知道，初值问题等价于积分方程，即对积分式分别采用矩形公式和梯形公式可得到对积分式分别采用矩形公式和梯形公式可得到欧拉公式和改进的欧拉公式，截断误差分别为欧拉公式和改进的欧拉公式，截断误差分别为O O(h(h2 2) )和和O(hO(h3 3) )。4. 4. 阿达姆斯方法阿达姆斯方法 38 数数值求解常微分方程求解常微分方程 为为此此，我我们们自自然然
18、可可以以想想到到，若若用用更更高高次次的的插插值值多多项项式式来来代代替替f(xf(x，y)y)，则则所所得得公公式式的的精精度度会会更高。这就是线性多步法的起源思想。更高。这就是线性多步法的起源思想。 本本章章前前面面介介绍绍的的方方法法称称为为单单步步法法，因因为为在在计计算算y yi+1i+1时时，只只用用到到前前面面yiyi的的值值。而而对对于于线线性性多多步步法法是是要要利利用用前前面面已已经经算算出出的的若若干干个个值值y yi-ki-k，y yi-1i-1，y yi i来求来求y yi+1i+1。39 数数值求解常微分方程求解常微分方程现现用用k k次次多多项项式式P Pk k(
19、x)(x)来来代替代替f(xf(x，y(x) y(x) 舍去余项舍去余项 并设并设y yi i=y(x=y(xi i) )，而，而y yi+1i+1为为y(xy(xi+1i+1) )的近似的近似值，于是可得到线性多步法的计算公式值，于是可得到线性多步法的计算公式40 数数值求解常微分方程求解常微分方程阿达姆斯阿达姆斯(Adams)(Adams)显式显式: :为了易于在电子计算机上实为了易于在电子计算机上实现，常将式中的现，常将式中的 用各点的已知函数值表示。用各点的已知函数值表示。特别，当特别，当q=2,3q=2,3时，有时，有: :其中其中f fi i=f(x=f(xi i,y,yi i)
20、) 41 数数值求解常微分方程求解常微分方程 阿达姆斯隐式的数值公式阿达姆斯隐式的数值公式。 例如，当例如，当q=2q=2时时当当q=3q=3时时 42 数数值求解常微分方程求解常微分方程阿达姆斯预估校正公式阿达姆斯预估校正公式 我我们们常常把把阿阿达达姆姆斯斯显显式式及及隐隐式式联联立立使使用用，即即构构造造所所谓谓阿阿达达姆姆斯斯预预估估校校正正公公式式。现现以以q=2q=2为为例例构造预估校正公式构造预估校正公式并取 43 数数值求解常微分方程求解常微分方程 与与同同阶阶的的龙龙格格库库塔塔方方法法相相比比较较，阿阿达达姆姆斯斯方方法法计计算算量量小小，公公式式简简单单，程程序序易易于于
21、实实现现。但但它它的的主主要要缺缺点点是是不不能能自自动动开开始始，开开始始的的前前几几个个值值要要依依赖赖于于其其它它方方法法获获得得。这这里里介介绍绍两两种种计计算算开始值的方法。开始值的方法。 (1) (1) 用单步法中的数值方法求出开始值。用单步法中的数值方法求出开始值。 (2) (2) 使用使用y(x)y(x)的泰勒展开式的泰勒展开式44 数数值求解常微分方程求解常微分方程 设二阶线性常微分方程的边值问题为设二阶线性常微分方程的边值问题为 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(
22、x)p(x)，q(x)q(x)，f(x)f(x)为为区区间间a a，b b上上足足够够光光滑滑的的已已知知函函数数，且且q(x)0q(x)0，、为为已已知知常数。常数。5. 5. 二阶线性常微分方程边值问题二阶线性常微分方程边值问题 45 数数值求解常微分方程求解常微分方程 在在上上述述条条件件下下，边边值值问问题题存存在在连连续续可可微微的的解解，且且是是唯唯一一的的。若若采采用用差差分分方方法法来来解解边边值值问问题题，其基本步骤是：其基本步骤是：将将区区间间a a，b b“离离散散化化”，即即给给a a，b b一一个分划，此分划常考虑等距个分划，此分划常考虑等距; ;对对每每一一个个基基
23、点点，将将各各阶阶导导数数用用差差商商来来近近似似表表示示，将将微微分分方方程程转转化化为为差差分分方方程程，进进而而转转化化为为线线性性代数方程组代数方程组; ; 解线性代数方程组，求得各基点上的近似解。解线性代数方程组，求得各基点上的近似解。46 数数值求解常微分方程求解常微分方程现具体给出求解边值问题的方法步骤。现具体给出求解边值问题的方法步骤。 首首先先将将区区间间a,ba,b进进行行等等距距分分划划，即即令令 x xi i=a+ih=a+ih，i=0,1,2,n;i=0,1,2,n;其中其中h=(b-a)/n;xh=(b-a)/n;x0 0=a,x=a,xn n=b=b其其次次在在各
24、各基基点点x xi i上上，将将y,yy,y用用差差商商来来近近似似表表示示。这里要求有相同阶数的截断误差，以保证精度协调这里要求有相同阶数的截断误差，以保证精度协调. .47 数数值求解常微分方程求解常微分方程代代入入二二阶阶线线性性常常微微分分方方程程的的边边值值问问题题得得到到近近似似差差分方程为分方程为: :其中其中 p pi i=p(x=p(xi i),q),qi i=q(x=q(xi i), f), fi i=f(x=f(xi i) )略去上两式的截断误差略去上两式的截断误差O(hO(h2 2) )，用，用y yi i代替代替y(xy(xi i) )有有 48 数数值求解常微分方程
25、求解常微分方程 定定理理: :若若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 49 数数值求解常微分方程求解常微分方程例例3 3: : 用差分法解边值问题用差分法解边值问题 y-y=x; y-y=x; y(0)=0 y(0)=0，y(1)=1y(1)=1， 0x10x1，h=0.1h=0.150 数数值求解常微分方程求解常微分方程51
