计算方法第8章 常微分方程初值问题的数值解法

1、1 第8章 常微分方程初值问 题的数值解法 v基础知识基础知识 v欧拉方法欧拉方法 v龙格龙格库塔方法库塔方法 2 8.1 基础知识基础知识 一、问题的提出一、问题的提出 很多实际问题都需要求解常微分方程。例如单摆问题。很多实际问题都需要求解常微分方程。例如单摆问题。 常微分方程分为线性常微分方程和非线性常微分方程，又可以分为一阶常微分方程分为线性常微分方程和非线性常微分方程，又可以分为一阶 常微分方程和高阶常微分方程。通过变量的替换，可以把高阶常微分方常微分方程和高阶常微分方程。通过变量的替换，可以把高阶常微分方 程转化为一阶常微分方程再求解。对于一阶常微分方程组，可以写成向程转化为一阶常微

2、分方程再求解。对于一阶常微分方程组，可以写成向 量形式的单个方程，求解方法与一阶常微分方程相似。因此本章只讨论量形式的单个方程，求解方法与一阶常微分方程相似。因此本章只讨论 一阶常微分方程的初值问题：一阶常微分方程的初值问题： 00) ( , ),( yxy bxayxf dx dy 目前在常微分方程理论中，只能求出某些特殊类型常微分方程的解析解，目前在常微分方程理论中，只能求出某些特殊类型常微分方程的解析解， 对大部分常微分方程，用解析方法求出常微分方程的精确解非常困难，对大部分常微分方程，用解析方法求出常微分方程的精确解非常困难， 甚至不存在解的解析表达式。为满足工程实践的需要，常常用数值

3、解法甚至不存在解的解析表达式。为满足工程实践的需要，常常用数值解法 求常微分方程的近似解。求常微分方程的近似解。 在本章中假设讨论的一阶常微分方程的初值问题的解在本章中假设讨论的一阶常微分方程的初值问题的解y(x)存在、唯一且存在、唯一且 足够光滑，方程本身是稳定的，即精确解足够光滑，方程本身是稳定的，即精确解y(x)连续且依赖于初始值及右连续且依赖于初始值及右 端函数。端函数。 3 8.1 基础知识基础知识 二、数值解法二、数值解法 一阶常微分方程初值问题的数值解法的主要思想，是对区间一阶常微分方程初值问题的数值解法的主要思想，是对区间a,b上的节点上的节点 a=x0 x1xnxn+1b 建

4、立建立y(xn)的近似值的近似值yn的某一递推格式，利用初值的某一递推格式，利用初值y0和已计算出的和已计算出的y1,y2, yk-1递推出递推出yk，并且用这一方法反复递推，依次得到，并且用这一方法反复递推，依次得到yk+1,yk+2,yn。这一。这一 求解方法称为步进式求解，相邻求解方法称为步进式求解，相邻2个节点的距离称为步长，记为个节点的距离称为步长，记为hi=xi+1xi。 为便于计算，常取成等距节点，称为定步长，这时把步长记为为便于计算，常取成等距节点，称为定步长，这时把步长记为h。 一阶常微分方程初值问题的数值解法有多种分类方法。一阶常微分方程初值问题的数值解法有多种分类方法。

5、多步法不能自行启动，必须先用单步法计算出多步法不能自行启动，必须先用单步法计算出y1,y2,yk-1，才能启动一，才能启动一 个个r步的多步法。步的多步法。 一种分类方法为：一种分类方法为： 单步法：每一轮递推只用到前面一轮的递推结果，递推格式为：单步法：每一轮递推只用到前面一轮的递推结果，递推格式为： yk=yk-1hT(xk-1,yk-1) 多步法：每一轮递推要用到前面多轮递推的结果，递推格式为：多步法：每一轮递推要用到前面多轮递推的结果，递推格式为： yk=yk-1hT(xk-r,yk-r,xk-r+1,yk-r+1, xk-1,yk-1)，其中，其中r1。 4 8.1 基础知识基础知识

6、 二、数值解法（续）二、数值解法（续） 另一种分类方法为：另一种分类方法为： 显式方法：递推公式的右端都是已知量，可以直接计算出递推的结显式方法：递推公式的右端都是已知量，可以直接计算出递推的结 果，递推格式为：果，递推格式为：yk=yk-1hT(xk-r,yk-r,xk-r+1,yk-r+1, xk-1,yk-1) 隐式方法：递推公式左端的未知量也出现在公式的右端，递推格式隐式方法：递推公式左端的未知量也出现在公式的右端，递推格式 为：为： yk=yk-1hT(xk-r,yk-r,xk-r+1,yk-r+1, xk,yk) 隐式方法的递推公式其实是一个方程。解方程的运算量可能较大，为避免隐式

7、方法的递推公式其实是一个方程。解方程的运算量可能较大，为避免 解方程，常采用预测解方程，常采用预测校正系统。校正系统。 预测预测校正系统：每一轮递推包括预测和校正这校正系统：每一轮递推包括预测和校正这2个步骤。先用显式方个步骤。先用显式方 法计算出法计算出yk，作为迭代的初值，这一过程称为预测；再把隐式方法的递推，作为迭代的初值，这一过程称为预测；再把隐式方法的递推 公式作为迭代公式，把预测值公式作为迭代公式，把预测值yk代入迭代公式右端进行迭代，这一过程称代入迭代公式右端进行迭代，这一过程称 为校正。在校正时往往迭代为校正。在校正时往往迭代1次或几次，校正值的精度就会有大幅提高。次或几次，校

8、正值的精度就会有大幅提高。 一阶常微分方程初值问题的数值解法一般是对连续的初值问题进行离散化一阶常微分方程初值问题的数值解法一般是对连续的初值问题进行离散化 处理，把微分方程转化为代数方程来求解。常用的离散化方法有：处理，把微分方程转化为代数方程来求解。常用的离散化方法有： 基基 于数值微分的离散化方法，于数值微分的离散化方法， 基于数值积分的离散化方法，基于数值积分的离散化方法， 基于泰勒基于泰勒 展开的离散化方法。展开的离散化方法。 5 8.2 欧拉方法欧拉方法 一、显式欧拉法一、显式欧拉法 在一阶常微分方程初值问题的数值解法中，显式欧拉（在一阶常微分方程初值问题的数值解法中，显式欧拉（E

9、uler）法是最简）法是最简 单的一种。显式欧拉法有明显的几何含义，缺点是精度不高。对于一阶单的一种。显式欧拉法有明显的几何含义，缺点是精度不高。对于一阶 常微分方程的初值问题：常微分方程的初值问题： 00) ( , ),( yxy bxayxf dx dy 显式欧拉法的递推公式为：显式欧拉法的递推公式为：yk=yk-1hf(xk-1,yk-1)，k=1,2,3,。 显式欧拉法每一轮递推只用到前面一轮递推的结果，因此它是单步法。显式欧拉法每一轮递推只用到前面一轮递推的结果，因此它是单步法。 由由基于数值微分的离散化方法、基于数值积分的离散化方法、基于泰勒基于数值微分的离散化方法、基于数值积分的

10、离散化方法、基于泰勒 展开的离散化方法展开的离散化方法都可以推导出显式欧拉法的递推公式。都可以推导出显式欧拉法的递推公式。 用显式欧拉法求解一阶常微分方程初值问题的过程，就是以已知的用显式欧拉法求解一阶常微分方程初值问题的过程，就是以已知的(x0,y0) 作为起点，代入显式欧拉法的递推公式的右端，计算出作为起点，代入显式欧拉法的递推公式的右端，计算出y(x)在在x1处的近似处的近似 值值y1；再以；再以(x1,y1)作为起点，用显式欧拉法的递推公式计算出作为起点，用显式欧拉法的递推公式计算出y(x)在在x2处处 的近似值的近似值y2；。 6 8.2 欧拉方法欧拉方法 一、显式欧拉法（续）一、显

11、式欧拉法（续） 显式欧拉法的递推公式是斜率为显式欧拉法的递推公式是斜率为f(xk-1,yk-1)，经过点，经过点(xk-1,yk-1)的直线方程。的直线方程。 上述上述显式欧拉法递推过程的几何含义，就是用曲线显式欧拉法递推过程的几何含义，就是用曲线y(x)在点在点(x0,y0)处的切处的切 线段代替线段代替y(x)在区间在区间x0,x1内的曲线段；再把曲线段的终点内的曲线段；再把曲线段的终点(x1,y(x1)近似近似 为切线段的终点为切线段的终点(x1,y1)，把曲线，把曲线y(x)在点在点(x1,y1)处切线的斜率处切线的斜率f(x1,y(x1)近近 似为似为f(x1,y1) ，做曲线，做曲

12、线y(x)在点在点(x1,y1)处的切线，并在区间处的切线，并在区间x1,x2内用近似内用近似 的切线段代替的切线段代替y(x)的曲线段；的曲线段；。欧拉法用一系列折线近似地代替曲线。欧拉法用一系列折线近似地代替曲线 y(x)，因此它又称为欧拉折线法。，因此它又称为欧拉折线法。 定义：假设某一步递推的起点是精确的，即定义：假设某一步递推的起点是精确的，即yi-1=y(xi-1)，那么这一步递推，那么这一步递推 的截断误差的截断误差Ri=y(xi)yi称为局部截断误差。称为局部截断误差。 定义：若某算法的局部截断误差为定义：若某算法的局部截断误差为O(hp+1)，则称此算法有，则称此算法有p阶精

13、度。阶精度。 显式欧拉法具有一阶精度显式欧拉法具有一阶精度。 显式欧拉法在步长过大时误差较大；在步长较小时需要多步递推，可能显式欧拉法在步长过大时误差较大；在步长较小时需要多步递推，可能 出现误差积累的现象。由于显式欧拉法的精度不高，因此在实际应用中出现误差积累的现象。由于显式欧拉法的精度不高，因此在实际应用中 用的较少。用的较少。 7 8.2 欧拉方法欧拉方法 显式欧拉法的算法显式欧拉法的算法 输入求解区间的边界a,b，步长h，起点x0=a处的纵坐标y0。 递推的次数n=(b-a)/h。 for(i=0;in;i+) /循环1次完成1轮递推 xi+1=xi+h; yi+1=yi+h*f(xi

14、,yi); for(i=0;i=n;i+) /输出计算结果 输出xi,yi。 8 8.2 欧拉方法欧拉方法显式欧拉法对应的程序显式欧拉法对应的程序 #include #define MAXSIZE 50 double f(double x,double y); void main(void) double a,b,h,xMAXSIZE,yMAXSIZE; long i,n; printf(n请输入求解区间请输入求解区间a,b：);scanf(%lf,%lf, printf(n请输入步长请输入步长h：);scanf(%lf, n=(long)(b-a)/h); x0=a; printf(n请输入

15、起点请输入起点x0=%lf处的纵坐标处的纵坐标y0：,x0); scanf(%lf, for(i=0;in;i+) xi+1=xi+h; yi+1=yi+h*f(xi,yi); printf(n计算结果为：计算结果为：); for(i=0;i=n;i+) printf(nx%ld=%lf,y%ld=%lf,i,xi,i,yi); double f(double x,double y) return(); /*计算并返回函数值计算并返回函数值f(x,y)*/ 9 8.2 欧拉方法欧拉方法 二、欧拉方法的变形二、欧拉方法的变形 方法一：隐式欧拉法（方法一：隐式欧拉法（后退的欧拉法后退的欧拉法） 递

16、推公式为：递推公式为： yk=yk-1hf(xk,yk)，k=1,2,3,。 用隐式欧拉法求解一阶常微分方程初值问题的过程，就是以已知的用隐式欧拉法求解一阶常微分方程初值问题的过程，就是以已知的 (x0,y0)作为起点，计算出作为起点，计算出y(x)在在x1处的近似值处的近似值y1；再以；再以(x1,y1)作为起点，作为起点， 计算出计算出y(x)在在x2处的近似值处的近似值y2；，象这样反复地递推。，象这样反复地递推。 显然，隐式欧拉法是隐式方法，递推公式的右端有未知量显然，隐式欧拉法是隐式方法，递推公式的右端有未知量yk。隐式欧拉。隐式欧拉 法需要求解法需要求解1个方程。为避免解方程，常用

17、显式欧拉法的计算结果作为迭个方程。为避免解方程，常用显式欧拉法的计算结果作为迭 代的初值代的初值yk(0)，把隐式欧拉法的递推公式作为迭代公式反复迭代，得到，把隐式欧拉法的递推公式作为迭代公式反复迭代，得到 迭代序列迭代序列yk(0),yk(1),yk(2),。如果步长足够小，那么迭代序列收敛于。如果步长足够小，那么迭代序列收敛于yk。 隐式欧拉法具有一阶精度。隐式欧拉法具有一阶精度。 隐式欧拉法精度不高，计算复杂，用的比较少。隐式欧拉法精度不高，计算复杂，用的比较少。 隐式欧拉法由点隐式欧拉法由点(x0,y0)递推到递推到点点(x1,y1)的几何含义，是把曲线的几何含义，是把曲线y(x)在点

18、在点 (x1,y1)处的切线平行移动，移动到经过点处的切线平行移动，移动到经过点(x0,y0) ，在区间，在区间x0,x1内用此内用此 直线段代替直线段代替y(x)的曲线段。的曲线段。 10 8.2 欧拉方法欧拉方法 二、欧拉方法的变形二、欧拉方法的变形 方法二：梯形公式法方法二：梯形公式法 递推公式为：递推公式为： yk=yk-1(h/2)*(f(xk-1,yk-1)+f(xk,yk)，k=1,2,3,。 用梯形公式法求解一阶常微分方程初值问题的过程，就是以已知的用梯形公式法求解一阶常微分方程初值问题的过程，就是以已知的 (x0,y0)作为起点，计算出作为起点，计算出y(x)在在x1处的近似

19、值处的近似值y1；再以；再以(x1,y1)作为起点，作为起点， 计算出计算出y(x)在在x2处的近似值处的近似值y2；，象这样反复地递推。，象这样反复地递推。 梯形公式法具有梯形公式法具有2阶精度阶精度。 显然，梯形公式法是隐式方法，需要求解方程。为避免解方程，常用显然，梯形公式法是隐式方法，需要求解方程。为避免解方程，常用 显式欧拉法的计算结果作为迭代的初值显式欧拉法的计算结果作为迭代的初值yk(0)，把梯形公式法的递推公，把梯形公式法的递推公 式作为迭代公式反复迭代，得到迭代序列式作为迭代公式反复迭代，得到迭代序列yk(0),yk(1),yk(2),。如果步。如果步 长长h足够小，那么迭代

20、序列收敛于足够小，那么迭代序列收敛于yk。 梯形公式法由点梯形公式法由点(x0,y0)递推到递推到点点(x1,y1)的几何含义，是经过点的几何含义，是经过点(x0,y0)做一做一 直线段，在区间直线段，在区间x0,x1内用此直线段代替内用此直线段代替y(x)的曲线段。此直线段的斜的曲线段。此直线段的斜 率等于曲线率等于曲线y(x)在点在点(x0,y0)处的切线斜率和曲线处的切线斜率和曲线y(x)在点在点(x1,y1)处的切线处的切线 斜率的平均值。类似地，由点斜率的平均值。类似地，由点(xk-1,yk-1)递推出递推出点点(xk,yk)，k=1,2,3,。 11 8.2 欧拉方法欧拉方法 二、

21、欧拉方法的变形二、欧拉方法的变形 方法三：中点欧拉方法（两步欧拉方法方法三：中点欧拉方法（两步欧拉方法 ） 递推公式为：递推公式为： yk+1=yk-12hf(xk,yk)，k=1,2,3,。 中点欧拉方法是双步法，需要中点欧拉方法是双步法，需要2个初值个初值y0和和y1才能启动递推过程。一般才能启动递推过程。一般 先用单步法由点先用单步法由点(x0,y0)计算出计算出(x1,y1)，再用中点欧拉方法反复地递推。，再用中点欧拉方法反复地递推。 中点欧拉方法由点中点欧拉方法由点(x0,y0)和和(x1,y1)递推出递推出点点(x2,y2)的几何含义，是经过的几何含义，是经过 点点(x0,y0)做

22、一直线段，在区间做一直线段，在区间x0,x2内用此直线段代替内用此直线段代替y(x)的曲线段。的曲线段。 此直线段的斜率等于曲线此直线段的斜率等于曲线y(x)在点在点(x1,y1)处的切线斜率。类似地，由点处的切线斜率。类似地，由点 (xk-1,yk-1)和和点点(xk,yk)递推出递推出点点(xk+1,yk+1) ，其中，其中k=1,2,3,。 中点欧拉方法具有中点欧拉方法具有2阶精度阶精度。 12 8.2 欧拉方法欧拉方法 二、欧拉方法的变形二、欧拉方法的变形 上述各种方法的比较：上述各种方法的比较： 由拉格郎日中值定理，在区间由拉格郎日中值定理，在区间(xk-1,xk)内必定存在内必定存

23、在k，使得：，使得： y(xk) y(xk-1)=y(k)(xkxk-1) y(xk)=y(xk-1)+hy(k)=y(xk-1)+hf(k,y(k)， 如果知道如果知道k，代入，代入这个递推公式这个递推公式，那么递推过程得到的序列，那么递推过程得到的序列y0,y1, y2,没有误差。求没有误差。求k往往很困难，因此常用一个易求的值近似地代替往往很困难，因此常用一个易求的值近似地代替 k 。显式欧拉法、隐式欧拉法、梯形公式法、中点欧拉方法的区别是。显式欧拉法、隐式欧拉法、梯形公式法、中点欧拉方法的区别是 对对k 、y(k)的近似方法不同。的近似方法不同。 显式欧拉法把显式欧拉法把k近似为区间近

24、似为区间xk-1,xk的起点的起点xk-1； 隐式欧拉法把隐式欧拉法把k近似为区间近似为区间xk-1,xk的终点的终点xk； 梯形公式法把处梯形公式法把处k的导数的导数y(k)近似为区间近似为区间xk-1,xk的起点和终点导数的起点和终点导数 的平均值；的平均值； 中点欧拉方法考察的区间为中点欧拉方法考察的区间为xk-1,xk+1，k(xk-1,xk+1)，k被近似为区被近似为区 间间xk-1,xk+1的中点的中点xk。 13 8.2 欧拉方法欧拉方法 三、改进的欧拉法三、改进的欧拉法 改进的欧拉法是一种预测改进的欧拉法是一种预测校正方法，它的每一轮递推包括预测和校正校正方法，它的每一轮递推包

25、括预测和校正 这这2个步骤：先用显式欧拉公式计算出个步骤：先用显式欧拉公式计算出y*k， y*k=yk-1hf(xk-1,yk-1) 这一步称为这一步称为“预测预测”；再用梯形公式迭代一次，；再用梯形公式迭代一次， yk=yk-1(h/2)*(f(xk-1,yk-1)+f(xk,y*k) 这一步称为这一步称为“校正校正”。改进的欧拉法精度比显式欧拉法高，不需要解方程，。改进的欧拉法精度比显式欧拉法高，不需要解方程， 是一种更实用的方法。是一种更实用的方法。 14 8.2 欧拉方法欧拉方法 改进的欧拉法的算法改进的欧拉法的算法 输入求解区间的边界a,b，步长h，起点x0=a处的纵坐标y0。 递推

26、的次数n=(b-a)/h。 for(i=0;in;i+) /循环1次完成1轮递推 xi+1=xi+h; yi+1=yi+h*f(xi,yi); /预测 yi+1=yi+h*(f(xi,yi)+f(xi+1,yi+1)/2; /校正 for(i=0;i=n;i+) /输出计算结果 输出xi,yi。 15 8.2 欧拉方法欧拉方法改进的欧拉法对应的程序改进的欧拉法对应的程序 #include #define MAXSIZE 50 double f(double x,double y); void main(void) double a,b,h,xMAXSIZE,yMAXSIZE; long i,n

27、; printf(n请输入求解区间请输入求解区间a,b：);scanf(%lf,%lf, printf(n请输入步长请输入步长h：);scanf(%lf, n=(long)(b-a)/h);x0=a; printf(n请输入起点请输入起点x0=%lf处的纵坐标处的纵坐标y0：,x0); scanf(%lf, for(i=0;in;i+) xi+1=xi+h; yi+1=yi+h*f(xi,yi); yi+1=yi+h*(f(xi,yi)+f(xi+1,yi+1)/2; printf(n计算结果为：计算结果为：); for(i=0;i=n;i+)printf(nx%ld=%lf,y%ld=%lf

28、,i,xi,i,yi); double f(double x,double y) return(); /*计算并返回函数值计算并返回函数值f(x,y)*/ 16 8.3 龙格龙格库塔方法库塔方法 一、泰勒展开方法一、泰勒展开方法 显式欧拉法、隐式欧拉法、梯形公式法、中点欧拉方法的缺陷是精度较显式欧拉法、隐式欧拉法、梯形公式法、中点欧拉方法的缺陷是精度较 低。本节介绍高精度的一步法。低。本节介绍高精度的一步法。 直接用泰勒公式展开是一种高阶显式一步法。在由点直接用泰勒公式展开是一种高阶显式一步法。在由点(xk-1,yk-1)递推出递推出点点 (xk,yk)时，把时，把y(xk)在在x=xk-1处

29、做处做泰勒展开：泰勒展开： y(xk)=y(xk-1)+hy(xk-1)+(h2/2!)y(xk-1)+(hn/n!)y(n)(xk-1)+Rn(x) 其中余项其中余项Rn(x)=(hn+1/(n+1)!)y(n+1)()，(xk-1,xk) 泰勒公式中的各阶导数依次为：泰勒公式中的各阶导数依次为：y(x)=f，y(x)=fx+fy f， 以已知的以已知的(x0,y0)作为起点，依次地用泰勒公式由点作为起点，依次地用泰勒公式由点(xk-1,yk-1)递推出递推出点点 (xk,yk)，k=1,2,3,，完成求解的过程。，完成求解的过程。 如果如果y(x)在考察的区间内具有直到在考察的区间内具有直

30、到(n+1)阶导数，那么泰勒公式可以阶导数，那么泰勒公式可以 计算到前计算到前(n+1)项，项，此时泰勒展开方法具有阶精度。此时泰勒展开方法具有阶精度。 从理论上讲，如果从理论上讲，如果y(x)在考察的区间内足够光滑，那么泰勒展开方在考察的区间内足够光滑，那么泰勒展开方 法可以具有任意阶精度。但是高阶泰勒展开方法计算量大，求法可以具有任意阶精度。但是高阶泰勒展开方法计算量大，求 f(x,y(x)的高阶导数困难，因此泰勒展开方法并不实用。的高阶导数困难，因此泰勒展开方法并不实用。 17 8.3 龙格龙格库塔方法库塔方法 二、龙格二、龙格库塔法基本思想库塔法基本思想 龙格龙格库塔（库塔（Runge

31、Kutta）法，简称）法，简称RK方法，是一种高阶显式一步方法，是一种高阶显式一步 法，而且不需要计算导数。法，而且不需要计算导数。Runge首先提出了间接使用泰勒展开式的方法：首先提出了间接使用泰勒展开式的方法： 用用f(x,y)在一些点上函数值的线性组合代替在一些点上函数值的线性组合代替y(x)的各阶导数，构造的各阶导数，构造yn+1的表的表 达式，比较这个表达式与达式，比较这个表达式与y(xn+1)在在x=xn处处泰勒展开式，确定泰勒展开式，确定yn+1的的表达式表达式中中 的系数，使的系数，使yn+1的的表达表达式与式与y(xn+1)泰勒展开式前面的若干项相等，从而具有泰勒展开式前面的

32、若干项相等，从而具有 对应泰勒展开式的精度阶次，这就是龙格对应泰勒展开式的精度阶次，这就是龙格库塔法的主要思想。库塔法的主要思想。 p级龙格级龙格库塔公式的一般形式为：库塔公式的一般形式为： yn+1=yn+h ，n=1,2,3,。 p i iik c 1 其中其中k1=f(xn,yn)，ki= ，i=2,3,p。 ),( 1 1 i j jijnin kbhyhaxf 上式中上式中ci、ai、bij是与具体常微分方程初值问题以及是与具体常微分方程初值问题以及n、h无关的常数。无关的常数。ki 是是f(x,y)在某些点处函数值，在某些点处函数值，ci是在线性组合求是在线性组合求yn+1时时ki

33、的的“权权”，ai和和bij用用 来确定来确定ki在在f(x,y)上的位置。上的位置。 以以p级龙格级龙格库塔公式作为递推公式，以已知的库塔公式作为递推公式，以已知的(x0,y0)作为起点，依次地作为起点，依次地 由点由点(xn,yn)递推出递推出点点(xn+1,yn+1)，n=0,1,2,，完成求解的过程。，完成求解的过程。 18 8.3 龙格龙格库塔方法库塔方法 二、龙格二、龙格库塔法基本思想（续）库塔法基本思想（续） 确定常数确定常数ci、ai、bij的原则是使龙格的原则是使龙格库塔公式与泰勒公式前面的尽可能多库塔公式与泰勒公式前面的尽可能多 的项相等。如果的项相等。如果p级龙格级龙格库

34、塔公式等于泰勒公式的前库塔公式等于泰勒公式的前(q+1)项，那么这个项，那么这个 龙格龙格库塔公式具有库塔公式具有q阶精度，称此公式为阶精度，称此公式为p级级q阶龙格阶龙格库塔公式。库塔公式。 一级一阶龙格一级一阶龙格库塔公式就是显式欧拉公式。库塔公式就是显式欧拉公式。 二级龙格二级龙格库塔公式与泰勒公式前面库塔公式与泰勒公式前面3项相等，具有二阶精度。这个方项相等，具有二阶精度。这个方 程组有程组有3个方程，个方程，4个变元，因此它有个变元，因此它有1个自由参数，二级二阶龙格个自由参数，二级二阶龙格库库 塔公式有无穷多个。改进的欧拉法（预测塔公式有无穷多个。改进的欧拉法（预测校正法）和中点欧

35、拉方法校正法）和中点欧拉方法 （两步欧拉方法）是二级龙格（两步欧拉方法）是二级龙格库塔公式的特例。库塔公式的特例。 标准（经典）龙格标准（经典）龙格库塔公式为：库塔公式为： yn+1=yn+(h/6)(k1+2k2+2k3+k4)， 其中其中k1=f(xn,yn)， k2=f(xn+h/2,yn+k1h/2)， k3=f(xn+h/2,yn+k2h/2)， k4=f(xn+h,yn+k3h)。 此公式是一种四级四阶龙格此公式是一种四级四阶龙格库塔公式库塔公式，是最常用的龙格最常用的龙格库塔公式。库塔公式。 19 8.3 龙格龙格库塔方法库塔方法 二、龙格二、龙格库塔法基本思想（续库塔法基本思想（续2） 龙格龙格库塔公式的级与阶并不总是相等。设库塔公式的级与阶并不总是相等。设p级龙格级龙格库塔公式所能达库塔公式所能达 到的最高精度阶数为到的最高精度阶数为q阶，那么当阶，那么当p=1,2,3,4时，时，q=p；当；当p=5,6,7时，时，q=p 1