常微分方程的数值解法PPT讲稿.ppt
常微分方程的数值解法第1页，共40页，编辑于2022年，星期六 实际中，很多问题的数学模型都是微分方程。我们可以研究它们的一些性质。但是，只有极少数特殊的方程有解析解。对于绝大部分的微分方程是没有解析解的。常微分方程作为微分方程的基本类型之一，在自然界与工程界有很广泛的应用。很多问题的数学表述都可以归结为常微分方程的定解问题。很多偏微分方程问题，也可以化为常微分方程问题来近似求解。本章讨论常微分方程的数值解法 引引 言言第2页，共40页，编辑于2022年，星期六 本章讨论一阶常微分方程的初值问题9.19.19.29.2虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，大量从实际问题当中归结出来的微分方程主要靠数值解法。第3页，共40页，编辑于2022年，星期六定义：定义：所谓数值解法数值解法，就是寻求初值问题上的近似值相邻两个节点间的距离称为步长步长。今后如不特别申明，总假定步长h为定数。第4页，共40页，编辑于2022年，星期六一、几何解释：一、几何解释：图图9 91 1 欧拉法的几何解释欧拉法的几何解释9.1 9.1 欧拉（欧拉（EulerEuler）方法）方法第5页，共40页，编辑于2022年，星期六二、计算格式：二、计算格式：1 1、公式推导：、公式推导：第6页，共40页，编辑于2022年，星期六2 2、几何意义几何意义：第7页，共40页，编辑于2022年，星期六3 3、计算格式：、计算格式：9.39.3 欧拉（欧拉（EulerEuler）法）法（也叫欧拉折线法）（也叫欧拉折线法）是最古老的一种数值解法，它体现了数值方法的基本思想，但精度很低，单独用它来作计算往往不能满足精度要求。第8页，共40页，编辑于2022年，星期六9.2 9.2 改进的欧拉方法改进的欧拉方法 同一种计算格式往往可以通过多种途径构造出来，本节与下一节就会看到这一点。一、计算格式：一、计算格式：1 1、公式推导：、公式推导：将方程（9.1）的两端同时积分，9.49.4第9页，共40页，编辑于2022年，星期六 选择不同的近似方法计算这个积分项会得到不同的计算格式。例如：用矩形公式计算积分项代入（9.4）得第10页，共40页，编辑于2022年，星期六这样建立起来的格式就是欧拉法的计算格式（9.3）。用矩形公式求积分值很粗糙，故欧拉格式精度也很低。为了改进精度，我们改用梯形法计算左端积分第11页，共40页，编辑于2022年，星期六将其代入（9.4）得9.59.5（9.5）式被称为解常微分方程的梯形法则梯形法则。格式（9.3）与（9.5）有本质上的区别，欧拉格式（9.3）是个直接的计算公式，这类格式称作显式的。第12页，共40页，编辑于2022年，星期六这个方程可以用迭代法求解（参看第五章）（参看第五章），不过计算量比较大。2 2、预报校正系统：、预报校正系统：综合使用上述两种格式，先用欧拉格式，求得一个称为预报值预报值。这样建立起来的预报校正系统称为改进的欧拉格式改进的欧拉格式。第13页，共40页，编辑于2022年，星期六3 3、改进的欧拉格式：、改进的欧拉格式：9.69.6格式（9.6）的每一步需要两次调用函数f f，它可以改写成下列形式：第14页，共40页，编辑于2022年，星期六二、算法与流程图：二、算法与流程图：1 1、算法分析：、算法分析：欧拉法每一步只需对f调用一次，而改进的欧拉法则不然，需对f调用两次，其计算量比欧拉法增加一倍，付出这种代价的目的是为了提高精度。由此可见，它比欧拉格式的截断误差提高了一倍。第15页，共40页，编辑于2022年，星期六2 2、流程图：（略）、流程图：（略）3 3、C C源程序：源程序：#include#include#define H 0.1#define N 10float f(x,y)float x,y;return(y-2*x/y);第16页，共40页，编辑于2022年，星期六main()float x0=0;float y0=1;float x1,y1;float yp,yc;float h=H;int i;for(i=1;i=N;i+)x1=x0+h;yp=y0+h*f(x0,y0);yc=y0+h*f(x1,yp);y1=(yp+yc)/2;printf(x=%f,y=%fn,x1,y1);x0=x1;y0=y1;第17页，共40页，编辑于2022年，星期六例例解解解初值问题我们分别用两种格式进行计算，这里欧拉格式的具体形式是 第18页，共40页，编辑于2022年，星期六而改进的欧拉格式是计算结果见下表：第19页，共40页，编辑于2022年，星期六同准确解比较，第二列欧拉格式的结果大致只有两位有效数字，而第三列改进的欧拉格式的结果则有三位有效数字。第20页，共40页，编辑于2022年，星期六结点欧拉法改进欧拉法准确解00.10.20.30.40.50.60.70.80.91.011.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.7177791.7847711.0959091.1840971.2662011.343361.4164021.4859561.5525141.6164751.6781661.73786711.0954451.1832161.2649111.3416411.4142141.4832401.5491931.6124521.6733201.732051第21页，共40页，编辑于2022年，星期六9.3 9.3 龙格库塔龙格库塔(Runge-kutta)(Runge-kutta)方法方法 最常用的是四阶的龙格库塔格式，但推导极为繁琐，我们以二阶为例，说明其思想方法。一、二阶龙格库塔法：一、二阶龙格库塔法：1 1、基本思想：、基本思想：设初值问题：第22页，共40页，编辑于2022年，星期六对差商9.79.7第23页，共40页，编辑于2022年，星期六平均斜率平均斜率。由此得知，只要对平均斜率提供一种算法，由（9.7）式便相应地得到一种计算格式。欧拉格式：欧拉格式：改进的欧拉格式：改进的欧拉格式：由于仅取一个点，所以精度很低。第24页，共40页，编辑于2022年，星期六 这就是龙格库塔方法的基本思想1 1、二阶龙格库塔法：、二阶龙格库塔法：第25页，共40页，编辑于2022年，星期六（同改进欧拉法改进欧拉法）这样设计出的计算格式：第26页，共40页，编辑于2022年，星期六9.89.8我们希望适当选择参数的值，第27页，共40页，编辑于2022年，星期六代入（9.8）知和二阶泰勒展开式第28页，共40页，编辑于2022年，星期六比较系数即可发现，要使（9.8）的截断误差为只要成立下列条件：9.99.9这里共有三个参数，但满足两个条件，因此有一个自由度。满足条件（9.9）的一族一族格式（9.8）统称二阶龙格库塔格二阶龙格库塔格式。式。这时龙格库塔格式称作变形的欧拉格式变形的欧拉格式，其形式是：第29页，共40页，编辑于2022年，星期六9.109.10二、三阶龙格库塔格式：二、三阶龙格库塔格式：第30页，共40页，编辑于2022年，星期六仍用（9.8）所给的形式可以使上述格式（9.10）的截断误差为这类格式统称为三阶龙格库塔格式三阶龙格库塔格式。第31页，共40页，编辑于2022年，星期六三、四阶龙格库塔法：三、四阶龙格库塔法：继续上述过程，可以进一步讨论四阶龙格库塔格式。一种最常用的经典龙格库塔格式为9.119.11第32页，共40页，编辑于2022年，星期六 若使四阶龙格库塔格式与改进欧拉格式的总体计算量相同，可以取较大的步长，但计算精度比改进欧拉法高很多。四、经典龙格库塔格式算法与流程图四、经典龙格库塔格式算法与流程图：四阶龙格库塔格式（9.11）的每一步需要四次调用函数2 2、流程图、流程图：（略）：（略）3 3、C-C-源程序源程序：1 1、算法分析：、算法分析：第33页，共40页，编辑于2022年，星期六#include#include#define H 0.2#define N 5float f(x,y)float x,y;return(y-2*x/y);main()float x0=0;float y0=1;float x1,y1,k1,k2,k3,k4;float h=H;int i;for(i=1;i=N;i+)第34页，共40页，编辑于2022年，星期六 x1=x0+h;k1=f(x0,y0);k2=f(x0+h/2,y0+(h/2)*k1);k3=f(x0+h/2,y0+(h/2)*k2);k4=f(x1,y0+h*k3);y1=y0+(h/6)*(k1+2*k2+2*k3+k4);printf(x=%f,y=%fn,x1,y1);x0=x1;y0=y1;第35页，共40页，编辑于2022年，星期六例例用四阶经典龙格库塔格式解初值问题4 4、MathematicaMathematica求解函数求解函数：NDSolveeqns,y,x,xmin,xmax函数功能：函数功能：对常微分方程或方程组eqns，求函数y关于x在xmin,xmax范围内的数值解。第36页，共40页，编辑于2022年，星期六解解由公式（9.11），有第37页，共40页，编辑于2022年，星期六例例用四阶经典龙格库塔格式解初值问题第38页，共40页，编辑于2022年，星期六解解计算结果如下表：结结 点点改改进进欧拉法欧拉法龙龙格格库库塔法塔法准确解准确解01110.21.134096 1.1832291.1832160.41.343360 1.3416671.3416410.6 1.485956 1.4832811.4832400.11.6164761.6125131.6124511.737869 1.7321401.732051第39页，共40页，编辑于2022年，星期六 一般认为，更高阶的龙格库塔格式所提高的精度已被所需增加的计算量抵消了，因此一般不采用。龙格库塔法的导出利用了泰勒展开（高阶可导），因此要求所求得解有较好的光滑性，否则可能精度反而更低。实际计算中应根据问题的具体情况选择合适的算法。第40页，共40页，编辑于2022年，星期六
