常微分方程的数值解法

常微分方程的数值解法第1页，共40页，2023年，2月20日，星期一※实际中，很多问题的数学模型都是微分方程。我们可以研究它们的一些性质。但是，只有极少数特殊的方程有解析解。对于绝大部分的微分方程是没有解析解的。※

常微分方程作为微分方程的基本类型之一，在自然界与工程界有很广泛的应用。很多问题的数学表述都可以归结为常微分方程的定解问题。很多偏微分方程问题，也可以化为常微分方程问题来近似求解。※本章讨论常微分方程的数值解法

引言第2页，共40页，2023年，2月20日，星期一※本章讨论一阶常微分方程的初值问题——9.1——9.2虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，大量从实际问题当中归结出来的微分方程主要靠数值解法。第3页，共40页，2023年，2月20日，星期一★定义：所谓数值解法，就是寻求初值问题上的近似值相邻两个节点间的距离称为步长。※今后如不特别申明，总假定步长h为定数。第4页，共40页，2023年，2月20日，星期一一、几何解释：图9－1欧拉法的几何解释§9.1欧拉（Euler）方法第5页，共40页，2023年，2月20日，星期一二、计算格式：1、公式推导：第6页，共40页，2023年，2月20日，星期一2、几何意义：第7页，共40页，2023年，2月20日，星期一3、计算格式：——9.3♀

欧拉（Euler）法（也叫欧拉折线法）是最古老的一种数值解法，它体现了数值方法的基本思想，但精度很低，单独用它来作计算往往不能满足精度要求。第8页，共40页，2023年，2月20日，星期一§9.2改进的欧拉方法※同一种计算格式往往可以通过多种途径构造出来，本节与下一节就会看到这一点。一、计算格式：1、公式推导：将方程（9.1）的两端同时积分，——9.4第9页，共40页，2023年，2月20日，星期一※选择不同的近似方法计算这个积分项会得到不同的计算格式。例如：用矩形公式计算积分项代入（9.4）得第10页，共40页，2023年，2月20日，星期一这样建立起来的格式就是欧拉法的计算格式（9.3）。用矩形公式求积分值很粗糙，故欧拉格式精度也很低。为了改进精度，我们改用梯形法计算左端积分第11页，共40页，2023年，2月20日，星期一将其代入（9.4）得——9.5（9.5）式被称为解常微分方程的梯形法则。※格式（9.3）与（9.5）有本质上的区别，欧拉格式（9.3）是个直接的计算公式，这类格式称作显式的。第12页，共40页，2023年，2月20日，星期一这个方程可以用迭代法求解（参看第五章），不过计算量比较大。2、预报－校正系统：综合使用上述两种格式，先用欧拉格式，求得一个称为预报值。这样建立起来的预报－校正系统称为改进的欧拉格式。第13页，共40页，2023年，2月20日，星期一3、改进的欧拉格式：——9.6格式（9.6）的每一步需要两次调用函数f，它可以改写成下列形式：第14页，共40页，2023年，2月20日，星期一二、算法与流程图：1、算法分析：欧拉法每一步只需对f调用一次，而改进的欧拉法则不然，需对f调用两次，其计算量比欧拉法增加一倍，付出这种代价的目的是为了提高精度。由此可见，它比欧拉格式的截断误差提高了一倍。第15页，共40页，2023年，2月20日，星期一2、流程图：（略）3、C－源程序：#include<stdio.h>#include<math.h>#defineH0.1#defineN10floatf(x,y)floatx,y;{return(y-2*x/y);}第16页，共40页，2023年，2月20日，星期一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=%f\n",x1,y1);x0=x1;y0=y1;}}第17页，共40页，2023年，2月20日，星期一［例］［解］解初值问题我们分别用两种格式进行计算，这里欧拉格式的具体形式是第18页，共40页，2023年，2月20日，星期一而改进的欧拉格式是计算结果见下表：第19页，共40页，2023年，2月20日，星期一同准确解比较，第二列欧拉格式的结果大致只有两位有效数字，而第三列改进的欧拉格式的结果则有三位有效数字。第20页，共40页，2023年，2月20日，星期一结点欧拉法改进欧拉法准确解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页，2023年，2月20日，星期一§9.3龙格－库塔(Runge-kutta)方法※最常用的是四阶的龙格－库塔格式，但推导极为繁琐，我们以二阶为例，说明其思想方法。一、二阶龙格－库塔法：1、基本思想：设初值问题：第22页，共40页，2023年，2月20日，星期一对差商——9.7第23页，共40页，2023年，2月20日，星期一平均斜率。由此得知，只要对平均斜率提供一种算法，由（9.7）式便相应地得到一种计算格式。♀欧拉格式：♀改进的欧拉格式：由于仅取一个点，所以精度很低。第24页，共40页，2023年，2月20日，星期一※这就是龙格－库塔方法的基本思想1、二阶龙格－库塔法：第25页，共40页，2023年，2月20日，星期一（同改进欧拉法）这样设计出的计算格式：第26页，共40页，2023年，2月20日，星期一——9.8我们希望适当选择参数的值，第27页，共40页，2023年，2月20日，星期一代入（9.8）知和二阶泰勒展开式第28页，共40页，2023年，2月20日，星期一比较系数即可发现，要使（9.8）的截断误差为只要成立下列条件：——9.9这里共有三个参数，但满足两个条件，因此有一个自由度。满足条件（9.9）的一族格式（9.8）统称二阶龙格－库塔格式。这时龙格－库塔格式称作变形的欧拉格式，其形式是：第29页，共40页，2023年，2月20日，星期一——9.10二、三阶龙格－库塔格式：第30页，共40页，2023年，2月20日，星期一仍用（9.8）所给的形式可以使上述格式（9.10）的截断误差为这类格式统称为三阶龙格－库塔格式。第31页，共40页，2023年，2月20日，星期一三、四阶龙格－库塔法：继续上述过程，可以进一步讨论四阶龙格－库塔格式。一种最常用的经典龙格－库塔格式为——9.11第32页，共40页，2023年，2月20日，星期一♀若使四阶龙格－库塔格式与改进欧拉格式的总体计算量相同，可以取较大的步长，但计算精度比改进欧拉法高很多。四、经典龙格－库塔格式算法与流程图：♀四阶龙格－库塔格式（9.11）的每一步需要四次调用函数2、流程图：（略）3、C-源程序：1、算法分析：第33页，共40页，2023年，2月20日，星期一#include<stdio.h>#include<math.h>#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++)第34页，共40页，2023年，2月20日，星期一{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=%f\n",x1,y1);x0=x1;y0=y1;}}第35页，共40页，2023年，2月20日，星期一［例］用四阶经典龙格－库塔格式解初值问题4、Mathematica求解函数：NDSolve[eqns,y,{x,xmin,xmax}]函数功能：对常微分方程或方程组eqns，求函数y关于x在[xmin,xmax]范围内的数值解。第36页，共40页，2023年，2月20日，星期一［解］由公式（9.11），有第37页，共40页，2023年，2月20日，星期一［例］用四阶经典龙格－库塔格式解初值问题第38页，共40页，2023年，2月20日，星期一［解］计算结果如下表：结点改进欧拉法龙格－库塔法准确解01110.21.1340961.1832291.1832160.41.3433601.3416671.3416410.61.4859561.4832811.4832400.11.