常微分方程初值问题的数值解法

《常微分方程初值问题的数值解法》由会员分享，可在线阅读，更多相关《常微分方程初值问题的数值解法（20页珍藏版）》请在金锄头文库上搜索。

1、第8章 常微分方程初值问题的数值解法v基础知识v欧拉方法v龙格库塔方法18.1 基础知识 一、问题的提出很多实际问题都需要求解常微分方程。例如单摆问题。 常微分方程分为线为线 性常微分方程和非线线性常微分方程，又可以分为为一阶阶 常微分方程和高阶阶常微分方程。通过变过变 量的替换换，可以把高阶阶常微分方 程转转化为为一阶阶常微分方程再求解。对对于一阶阶常微分方程组组，可以写成向 量形式的单单个方程，求解方法与一阶阶常微分方程相似。因此本章只讨论讨论 一阶阶常微分方程的初值问题值问题 ： 目前在常微分方程理论论中，只能求出某些特殊类类型常微分方程的解析解 ，对对大部分常微分方程，用解析方法求出常

2、微分方程的精确解非常困难难 ，甚至不存在解的解析表达式。为满为满 足工程实实践的需要，常常用数值值解 法求常微分方程的近似解。 在本章中假设讨论设讨论 的一阶阶常微分方程的初值问题值问题 的解y(x)存在、唯一且 足够够光滑，方程本身是稳稳定的，即精确解y(x)连续连续 且依赖赖于初始值值及右 端函数。 28.1 基础知识二、数值解法 一阶阶常微分方程初值问题值问题 的数值值解法的主要思想，是对对区间间a,b上的节节 点 a=x0x1xnxn+1b建立y(xn)的近似值yn的某一递推格式，利用初值y0和已计算出的y1,y2, yk-1递推出yk，并且用这一方法反复递推，依次得到yk+1,yk+

3、2,yn。这一 求解方法称为步进式求解，相邻2个节点的距离称为步长，记为hi=xi+1xi 。为便于计算，常取成等距节点，称为定步长，这时把步长记为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。 38.1

4、 基础知识二、数值解法（续）另一种分类方法为： 显显式方法：递递推公式的右端都是已知量，可以直接计计算出递递推的结结 果，递递推格式为为：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)隐隐式方法的递递推公式其实实是一个方程。解方程的运算量可能较较大，为为避免 解方程，常采用预测预测 校正系统统。 预测预测 校正系统统：每一轮递轮递 推包括预测预测 和校正这这2个步骤骤。先用显显式 方法计计算出

5、yk，作为为迭代的初值值，这这一过过程称为预测为预测 ；再把隐隐式方法的 递递推公式作为为迭代公式，把预测值预测值 yk代入迭代公式右端进进行迭代，这这一 过过程称为为校正。在校正时时往往迭代1次或几次，校正值值的精度就会有大幅 提高。 一阶阶常微分方程初值问题值问题 的数值值解法一般是对连续对连续 的初值问题进值问题进 行离散 化处处理，把微分方程转转化为为代数方程来求解。常用的离散化方法有： 基于数值微分的离散化方法， 基于数值积分的离散化方法， 基于泰 勒展开的离散化方法。 48.2 欧拉方法一、显式欧拉法在一阶阶常微分方程初值问题值问题 的数值值解法中，显显式欧拉（Euler）法是最

6、简单简单 的一种。显显式欧拉法有明显显的几何含义义，缺点是精度不高。对对于一 阶阶常微分方程的初值问题值问题 ： 显式欧拉法的递推公式为：yk=yk-1hf(xk-1,yk-1)，k=1,2,3,。 显显式欧拉法每一轮递轮递 推只用到前面一轮递轮递 推的结结果，因此它是单单步法。由基于数值微分的离散化方法、基于数值积分的离散化方法、基于泰勒 展开的离散化方法都可以推导导出显显式欧拉法的递递推公式。 用显显式欧拉法求解一阶阶常微分方程初值问题值问题 的过过程，就是以已知的 (x0,y0)作为为起点，代入显显式欧拉法的递递推公式的右端，计计算出y(x)在x1处处 的近似值值y1；再以(x1,y1)

7、作为为起点，用显显式欧拉法的递递推公式计计算出y(x) 在x2处处的近似值值y2；。 58.2 欧拉方法一、显式欧拉法（续）显式欧拉法的递推公式是斜率为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) ，做曲线线y(x)在点(x1,y1)处处的切线线，并在区间间x

8、1,x2内用近似 的切线线段代替y(x)的曲线线段；。欧拉法用一系列折线线近似地代替曲 线线y(x)，因此它又称为为欧拉折线线法。 定义：假设某一步递推的起点是精确的，即yi-1=y(xi-1)，那么这一步递推 的截断误差Ri=y(xi)yi称为局部截断误差。 定义：若某算法的局部截断误差为O(hp+1)，则称此算法有p阶精度。 显式欧拉法具有一阶精度。显式欧拉法在步长过大时误差较大；在步长较小时需要多步递推，可能 出现误差积累的现象。由于显式欧拉法的精度不高，因此在实际应用中 用的较少。 68.2 欧拉方法显式欧拉法的算法输入求解区间的边界a,b，步长h，起点x0=a处的纵坐标y0。 递推的

9、次数n=(b-a)/h。 for(i=0;i #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请输入起点x0=%lf处的纵坐标y0：“,x0); scanf(“%lf“, for(i=0;i #define MAXSIZE 5

10、0 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请输入起点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(

11、xi+1,yi+1)/2; 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)*/ 158.3 龙格库塔方法一、泰勒展开方法显显式欧拉法、隐隐式欧拉法、梯形公式法、中点欧拉方法的缺陷是精度较较 低。本节节介绍绍高精度的一步法。 直接用泰勒公式展开是一种高阶显阶显 式一步法。在由点(xk-1,yk-1)递推出点 (xk,yk)时，把y(xk)在x=xk-1处做泰勒展开： y(xk)=y(xk-1

12、)+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)在考察的区间间内具有直到(n+1)阶导阶导 数，那么泰勒公式可以 计计算到前(n+1)项项，此时泰勒展开方法具有阶精度。 从理论论上讲讲，如果y(x)在考察的区间间内足够够光滑，那么泰勒展开方

13、法可以具有任意阶阶精度。但是高阶阶泰勒展开方法计计算量大，求 f(x,y(x)的高阶导阶导 数困难难，因此泰勒展开方法并不实实用。 168.3 龙格库塔方法 二、龙格库塔法基本思想龙龙格库库塔（RungeKutta）法，简简称RK方法，是一种高阶显阶显 式一步 法，而且不需要计计算导导数。Runge首先提出了间间接使用泰勒展开式的方法 ：用f(x,y)在一些点上函数值值的线线性组组合代替y(x)的各阶导阶导 数，构造yn+1的 表达式，比较这个表达式与y(xn+1)在x=xn处泰勒展开式，确定yn+1的表达式 中的系数，使yn+1的表达式与y(xn+1)泰勒展开式前面的若干项项相等，从而具 有

14、对应对应 泰勒展开式的精度阶阶次，这这就是龙龙格库库塔法的主要思想。 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)作为为起点，依次地 由点(xn,yn)递推出点(xn+1,yn+1)，n=0,1,2,，完成求解的过过程。

15、178.3 龙格库塔方法 二、龙格库塔法基本思想（续）确定常数ci、ai、bij的原则是使龙格库塔公式与泰勒公式前面的尽可能多 的项相等。如果p级龙格库塔公式等于泰勒公式的前(q+1)项，那么这个 龙格库塔公式具有q阶精度，称此公式为p级q阶龙格库塔公式。 一级级一阶龙阶龙 格库库塔公式就是显显式欧拉公式。 二级龙格库塔公式与泰勒公式前面3项相等，具有二阶精度。这个方 程组有3个方程，4个变元，因此它有1个自由参数，二级二阶龙格库 塔公式有无穷多个。改进的欧拉法（预测校正法）和中点欧拉方法 （两步欧拉方法）是二级龙格库塔公式的特例。 标准（经典）龙格库塔公式为： 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)。此公式是一种四级四阶龙格库塔公式，是最常用的龙格库塔公式。 188.3 龙格库塔方法 二、龙格库塔法基本思想（续2）龙龙格库库塔公式的级级与阶阶并不总总是相等。设设p级龙级龙 格库库塔公式所能达 到的最高精度阶阶数为为q阶