第九章常微分方程的数值解法

第九章常微分方程的数值解法,9.1欧拉（Euler）方法9.2改进的欧拉方法9.3龙格库塔（Runge-kutta）方法9.4步长的自动选择,实际中，很多问题的数学模型都是微分方程。我们可以研究它们的一些性质。但是，只有极少数特殊的方程有解析解。对于绝大部分的微分方程是没有解析解的。,常微分方程作为微分方程的基本类型之一，在自然界与工程界有很广泛的应用。很多问题的数学表述都可以归结为常微分方程的定解问题。很多偏微分方程问题，也可以化为常微分方程问题来近似求解。,本章讨论常微分方程的数值解法,引言,本章讨论一阶常微分方程的初值问题,9.1,9.2,虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，大量从实际问题当中归结出来的微分方程主要靠数值解法。,定义：所谓数值解法，就是寻求初值问题,上的近似值,相邻两个节点间的距离,称为步长。,今后如不特别申明，总假定步长h为定数。,一、几何解释：,图91欧拉法的几何解释,9.1欧拉（Euler）方法,二、计算格式：,1、公式推导：,2、几何意义：,3、计算格式：,9.3,欧拉（Euler）法（也叫欧拉折线法）是最古老的一种数值解法，它体现了数值方法的基本思想，但精度很低，单独用它来作计算往往不能满足精度要求。,9.2改进的欧拉方法,同一种计算格式往往可以通过多种途径构造出来，本节与下一节就会看到这一点。,一、计算格式：,1、公式推导：,将方程（9.1）的两端同时积分，,9.4,选择不同的近似方法计算这个积分项会得到不同的计算格式。,例如：用矩形公式计算积分项,代入（9.4）得,这样建立起来的格式就是欧拉法的计算格式（9.3）。,用矩形公式求积分值很粗糙，故欧拉格式精度也很低。,为了改进精度，我们改用梯形法计算左端积分,将其代入（9.4）得,9.5,（9.5）式被称为解常微分方程的梯形法则。,格式（9.3）与（9.5）有本质上的区别，欧拉格式（9.3）是个直接的计算公式，这类格式称作显式的。,这个方程可以用迭代法求解（参看第五章），不过计算量比较大。,2、预报校正系统：,综合使用上述两种格式，先用欧拉格式，求得一个,称为预报值。,这样建立起来的预报校正系统称为改进的欧拉格式。,3、改进的欧拉格式：,9.6,格式（9.6）的每一步需要两次调用函数f，它可以改写成下列形式：,二、算法与流程图：,1、算法分析：,欧拉法每一步只需对f调用一次，而改进的欧拉法则不然，需对f调用两次，其计算量比欧拉法增加一倍，付出这种代价的目的是为了提高精度。,由此可见，它比欧拉格式的截断误差提高了一倍。,2、流程图：（略）,3、C源程序：,#include#include#defineH0.1#defineN10floatf(x,y)floatx,y;return(y-2*x/y);,main()floatx0=0;floaty0=1;floatx1,y1;floatyp,yc;floath=H;inti;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;,例,解,解初值问题,我们分别用两种格式进行计算，这里欧拉格式的具体形式是,而改进的欧拉格式是,计算结果见下表：,同准确解比较，第二列欧拉格式的结果大致只有两位有效数字，而第三列改进的欧拉格式的结果则有三位有效数字。,9.3龙格库塔(Runge-kutta)方法,最常用的是四阶的龙格库塔格式，但推导极为繁琐，我们以二阶为例，说明其思想方法。,一、二阶龙格库塔法：,1、基本思想：,设初值问题：,对差商,9.7,平均斜率。,由此得知，只要对平均斜率提供一种算法，由（9.7）式便相应地得到一种计算格式。,欧拉格式：,改进的欧拉格式：,由于仅取一个点，所以精度很低。,这就是龙格库塔方法的基本思想,1、二阶龙格库塔法：,（同改进欧拉法）,这样设计出的计算格式：,9.8,我们希望适当选择参数的值，,代入（9.8）知,和二阶泰勒展开式,比较系数即可发现，要使（9.8）的截断误差为,只要成立下列条件：,9.9,这里共有三个参数，但满足两个条件，因此有一个自由度。满足条件（9.9）的一族格式（9.8）统称二阶龙格库塔格式。,这时龙格库塔格式称作变形的欧拉格式，其形式是：,9.10,二、三阶龙格库塔格式：,仍用（9.8）所给的形式,可以使上述格式（9.10）的截断误差为,这类格式统称为三阶龙格库塔格式。,三、四阶龙格库塔法：,继续上述过程，可以进一步讨论四阶龙格库塔格式。一种最常用的经典龙格库塔格式为,9.11,若使四阶龙格库塔格式与改进欧拉格式的总体计算量相同，可以取较大的步长，但计算精度比改进欧拉法高很多。,四、经典龙格库塔格式算法与流程图：,四阶龙格库塔格式（9.11）的每一步需要四次调用函数,2、流程图：（略）,3、C-源程序：,1、算法分析：,#include#include#defineH0.2#defineN5floatf(x,y)floatx,y;return(y-2*x/y);main()floatx0=0;floaty0=1;floatx1,y1,k1,k2,k3,k4;floath=H;inti;for(i=1;i=N;i+),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;,例,用四阶经典龙格库塔格式解初值问题,4、Mathematica求解函数：,NDSolveeqns,y,x,xmin,xmax,函数功能：对常微分方程或方程组eqns，求函数y关于x在xmin,xmax范围内的数值解。,解,由公式（9.11），有,例,用四阶经典龙格库塔格式解初