1、1,第8章 常微分方程初值问题的数值解法,基础知识欧拉方法龙格库塔方法,2,8.1 基础知识,一、问题的提出,很多实际问题都需要求解常微分方程。例如单摆问题。,常微分方程分为线性常微分方程和非线性常微分方程，又可以分为一阶常微分方程和高阶常微分方程。通过变量的替换，可以把高阶常微分方程转化为一阶常微分方程再求解。对于一阶常微分方程组，可以写成向量形式的单个方程，求解方法与一阶常微分方程相似。因此本章只讨论一阶常微分方程的初值问题：,目前在常微分方程理论中，只能求出某些特殊类型常微分方程的解析解，对大部分常微分方程，用解析方法求出常微分方程的精确解非常困难，甚至不存在解的解析表达式。为满足工程实
2、践的需要，常常用数值解法求常微分方程的近似解。,在本章中假设讨论的一阶常微分方程的初值问题的解y(x)存在、唯一且足够光滑，方程本身是稳定的，即精确解y(x)连续且依赖于初始值及右端函数。,3,8.1 基础知识,二、数值解法,一阶常微分方程初值问题的数值解法的主要思想，是对区间a,b上的节点,a=x0x1xnxn+1b,建立y(xn)的近似值yn的某一递推格式，利用初值y0和已计算出的y1,y2, yk-1递推出yk，并且用这一方法反复递推，依次得到yk+1,yk+2,yn。这一求解方法称为步进式求解，相邻2个节点的距离称为步长，记为hi=xi+1xi。为便于计算，常取成等距节点，称为定步长，
3、这时把步长记为h。,一阶常微分方程初值问题的数值解法有多种分类方法。,多步法不能自行启动，必须先用单步法计算出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 基础知识,二、数值解法（续）,另一种分类方法为：, 显式方法：递推公式的右端都是已知量，可以直接计算出递推的结果，递推格式为：yk=yk-1hT(xk
4、-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),隐式方法的递推公式其实是一个方程。解方程的运算量可能较大，为避免解方程，常采用预测校正系统。, 预测校正系统：每一轮递推包括预测和校正这2个步骤。先用显式方法计算出yk，作为迭代的初值，这一过程称为预测；再把隐式方法的递推公式作为迭代公式，把预测值yk代入迭代公式右端进行迭代，这一过程称为校正。在校正时往往迭代1次或几次，校正值的精度就会有大幅提高。,一阶常微分方程初值问题的数值
5、解法一般是对连续的初值问题进行离散化处理，把微分方程转化为代数方程来求解。常用的离散化方法有： 基于数值微分的离散化方法， 基于数值积分的离散化方法， 基于泰勒展开的离散化方法。,5,8.2 欧拉方法,一、显式欧拉法,在一阶常微分方程初值问题的数值解法中，显式欧拉（Euler）法是最简单的一种。显式欧拉法有明显的几何含义，缺点是精度不高。对于一阶常微分方程的初值问题：,显式欧拉法的递推公式为：yk=yk-1hf(xk-1,yk-1)，k=1,2,3,。,显式欧拉法每一轮递推只用到前面一轮递推的结果，因此它是单步法。,由基于数值微分的离散化方法、基于数值积分的离散化方法、基于泰勒展开的离散化方法
6、都可以推导出显式欧拉法的递推公式。,用显式欧拉法求解一阶常微分方程初值问题的过程，就是以已知的(x0,y0)作为起点，代入显式欧拉法的递推公式的右端，计算出y(x)在x1处的近似值y1；再以(x1,y1)作为起点，用显式欧拉法的递推公式计算出y(x)在x2处的近似值y2；。,6,8.2 欧拉方法,一、显式欧拉法（续）,显式欧拉法的递推公式是斜率为f(xk-1,yk-1)，经过点(xk-1,yk-1)的直线方程。,上述显式欧拉法递推过程的几何含义，就是用曲线y(x)在点(x0,y0)处的切线段代替y(x)在区间x0,x1内的曲线段；再把曲线段的终点(x1,y(x1)近似为切线段的终点(x1,y1
7、)，把曲线y(x)在点(x1,y1)处切线的斜率f(x1,y(x1)近似为f(x1,y1) ，做曲线y(x)在点(x1,y1)处的切线，并在区间x1,x2内用近似的切线段代替y(x)的曲线段；。欧拉法用一系列折线近似地代替曲线y(x)，因此它又称为欧拉折线法。,定义：假设某一步递推的起点是精确的，即yi-1=y(xi-1)，那么这一步递推的截断误差Ri=y(xi)yi称为局部截断误差。,定义：若某算法的局部截断误差为O(hp+1)，则称此算法有p阶精度。,显式欧拉法具有一阶精度。,显式欧拉法在步长过大时误差较大；在步长较小时需要多步递推，可能出现误差积累的现象。由于显式欧拉法的精度不高，因此在
8、实际应用中用的较少。,7,8.2 欧拉方法,显式欧拉法的算法,8,8.2 欧拉方法,显式欧拉法对应的程序,#include #define MAXSIZE 50double f(double x,double y);void main(void)double a,b,h,xMAXSIZE,yMAXSIZE;long i,n;printf(n请输入求解区间a,b：);scanf(%lf,%lf,double f(double x,double y)return(); /*计算并返回函数值f(x,y)*/,9,8.2 欧拉方法,二、欧拉方法的变形,方法一：隐式欧拉法（后退的欧拉法）,递推公式为：
9、yk=yk-1hf(xk,yk)，k=1,2,3,。,用隐式欧拉法求解一阶常微分方程初值问题的过程，就是以已知的(x0,y0)作为起点，计算出y(x)在x1处的近似值y1；再以(x1,y1)作为起点，计算出y(x)在x2处的近似值y2；，象这样反复地递推。,显然，隐式欧拉法是隐式方法，递推公式的右端有未知量yk。隐式欧拉法需要求解1个方程。为避免解方程，常用显式欧拉法的计算结果作为迭代的初值yk(0)，把隐式欧拉法的递推公式作为迭代公式反复迭代，得到迭代序列yk(0),yk(1),yk(2),。如果步长足够小，那么迭代序列收敛于yk。,隐式欧拉法具有一阶精度。,隐式欧拉法精度不高，计算复杂，用
10、的比较少。,隐式欧拉法由点(x0,y0)递推到点(x1,y1)的几何含义，是把曲线y(x)在点(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处的近似值y1；再以(x1,y1)作为起点，计算出y(x)在x2处的近似值y2；，象这样反复地递推。,梯形公式法具有2阶
11、精度。,显然，梯形公式法是隐式方法，需要求解方程。为避免解方程，常用显式欧拉法的计算结果作为迭代的初值yk(0)，把梯形公式法的递推公式作为迭代公式反复迭代，得到迭代序列yk(0),yk(1),yk(2),。如果步长h足够小，那么迭代序列收敛于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.
12、2 欧拉方法,二、欧拉方法的变形,方法三：中点欧拉方法（两步欧拉方法 ）,递推公式为： 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)做一直线段，在区间x0,x2内用此直线段代替y(x)的曲线段。此直线段的斜率等于曲线y(x)在点(x1,y1)处的切线斜率。类似地，由点(xk-1,yk-1)和点(xk,yk)递推出点(xk+1,yk+1
13、) ，其中k=1,2,3,。,中点欧拉方法具有2阶精度。,12,8.2 欧拉方法,二、欧拉方法的变形,上述各种方法的比较：,由拉格郎日中值定理，在区间(xk-1,xk)内必定存在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近似为区间xk-1,xk的起点xk-1；隐
14、式欧拉法把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 欧拉方法,三、改进的欧拉法,改进的欧拉法是一种预测校正方法，它的每一轮递推包括预测和校正这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),这一步称为“校正”。改进的欧拉法精度比显式欧拉法
15、高，不需要解方程，是一种更实用的方法。,14,8.2 欧拉方法,改进的欧拉法的算法,15,8.2 欧拉方法,改进的欧拉法对应的程序,#include #define MAXSIZE 50double f(double x,double y);void main(void)double a,b,h,xMAXSIZE,yMAXSIZE;long i,n;printf(n请输入求解区间a,b：);scanf(%lf,%lf,double f(double x,double y)return(); /*计算并返回函数值f(x,y)*/,16,8.3 龙格库塔方法,一、泰勒展开方法,显式欧拉法、隐式欧拉
16、法、梯形公式法、中点欧拉方法的缺陷是精度较低。本节介绍高精度的一步法。,直接用泰勒公式展开是一种高阶显式一步法。在由点(xk-1,yk-1)递推出点(xk,yk)时，把y(xk)在x=xk-1处做泰勒展开：,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,
17、，完成求解的过程。,如果y(x)在考察的区间内具有直到(n+1)阶导数，那么泰勒公式可以计算到前(n+1)项，此时泰勒展开方法具有阶精度。,从理论上讲，如果y(x)在考察的区间内足够光滑，那么泰勒展开方法可以具有任意阶精度。但是高阶泰勒展开方法计算量大，求f(x,y(x)的高阶导数困难，因此泰勒展开方法并不实用。,17,8.3 龙格库塔方法,二、龙格库塔法基本思想,龙格库塔（RungeKutta）法，简称RK方法，是一种高阶显式一步法，而且不需要计算导数。Runge首先提出了间接使用泰勒展开式的方法：用f(x,y)在一些点上函数值的线性组合代替y(x)的各阶导数，构造yn+1的表达式，比较这个
18、表达式与y(xn+1)在x=xn处泰勒展开式，确定yn+1的表达式中的系数，使yn+1的表达式与y(xn+1)泰勒展开式前面的若干项相等，从而具有对应泰勒展开式的精度阶次，这就是龙格库塔法的主要思想。,p级龙格库塔公式的一般形式为：,yn+1=yn+h ，n=1,2,3,。,其中k1=f(xn,yn)，ki= ，i=2,3,p。,上式中ci、ai、bij是与具体常微分方程初值问题以及n、h无关的常数。ki是f(x,y)在某些点处函数值，ci是在线性组合求yn+1时ki的“权”，ai和bij用来确定ki在f(x,y)上的位置。,以p级龙格库塔公式作为递推公式，以已知的(x0,y0)作为起点，依次
19、地由点(xn,yn)递推出点(xn+1,yn+1)，n=0,1,2,，完成求解的过程。,18,8.3 龙格库塔方法,二、龙格库塔法基本思想（续）,确定常数ci、ai、bij的原则是使龙格库塔公式与泰勒公式前面的尽可能多的项相等。如果p级龙格库塔公式等于泰勒公式的前(q+1)项，那么这个龙格库塔公式具有q阶精度，称此公式为p级q阶龙格库塔公式。,一级一阶龙格库塔公式就是显式欧拉公式。,二级龙格库塔公式与泰勒公式前面3项相等，具有二阶精度。这个方程组有3个方程，4个变元，因此它有1个自由参数，二级二阶龙格库塔公式有无穷多个。改进的欧拉法（预测校正法）和中点欧拉方法（两步欧拉方法）是二级龙格库塔公式
20、的特例。,标准（经典）龙格库塔公式为：,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=p1；当p=8,9时，q=p2；。四级之上的龙格库塔公式，计算量增加较大，精度增加较
21、小。一般情况下，四级龙格库塔公式已经能满足实际需求，因此很少用到四级之上的龙格库塔公式。,与显式欧拉法、隐式欧拉法、改进的欧拉法、中点欧拉法相比，龙格库塔法可以达到较高精度，但运算量较大，标准龙格库塔方法的计算量大约是改进的欧拉法计算量的2倍。龙格库塔法是基于泰勒公式的方法，它要求y(x)在考察的区间内足够光滑。如果解的光滑性较好，那么在相同步长下，标准龙格库塔方法比上述改进的欧拉法等方法的精度高得多；反之，如果解的光滑性较差，标准龙格库塔方法可能比改进的欧拉法的精度还要低。,20,本章主要介绍了一阶常微分方程初值问题的数值解法。1一阶常微分方程初值问题数值解法的分类分类方法1： 单步法； 多步法。分类方法2： 显式方法； 隐式方法； 预测校正法。分类方法3： 基于数值微分的离散化方法； 基于数值积分的离散化方法； 基于泰勒展开的离散化方法。2欧拉方法及其变形本章对下列方法给出了基本原理、几何含义和精度分析，并进行了比较。 显式欧拉法（欧拉折线法） 隐式欧拉法（后退的欧拉法） 梯形公式法 中点欧拉方法（两步欧拉方法） 改进的欧拉法（预测校正法）3龙格库塔方法本章从泰勒展开方法入手，引出了龙格库塔法的基本原理，介绍了常用的龙格库塔公式，给出了标准龙格库塔法的算法和程序。,第8章 常微分方程初值问题的数值解法 小结,
