应用基础数值 14

应用数值分析第十六讲常微分方程的数值解法1问题引入

23Runge-Kutta方法4提纲内容小结欧拉方法教学目标与重难点教学目标Euler方法与Euler方法的改进。实际问题中，运用Euler方法建立与求解。

1.了解Runge-Kutta方法；2.理解Euler方法的改进；

3.掌握Euler方法。教学重点教学难点一问题引入一、问题引入

在许多应用问题中得出的数学模型，常常是一个或一组常微分方程，然而只有极少数的微分方程能够用初等的方法求出它们的解析解，多数微分方程只能用近似的方法求其近似解。近似方法有两类，一类称为近似解析法，如级数解法，逐次逼近法；另一类称之为数值解法。本章主要介绍数值解法，其基本思想是利用近似计算求出微分方程的解在一些离散点处的近似值。由于高阶的微分方程可以转化为一阶常微分方程组，而一阶常微分方程组又可写成向量形式的单个方程。所以本章主要考虑常微分方程初值问题(10.1)为了使利用数值解法得出的常微分方程的近似解具有实际意义，就必须保证常微分方程的解存在且唯一。一、问题引入定理1

设常微分方程(10.1)中的二元函数满足1)在区域上连续。2)在D上关于y满足Lipschitz条件，即存在常数L，使(10.2)其中L称为Lipschitz常数。则问题(10.1)在区间[a,b]上存在唯一连续可微的解以后我们总是假定f(x,y)在D上关于y满足Lipschitz条件。一、问题引入

数值解法的具体做法是在区间[a,b]上插入n-1个点，使记，称hi为从xi到xi+1的步长，在这些节点{xi}上通过离散化的方法，求得y(xi)的近似值yi。数值解法常分为两大类1)单步法：计算节点xi+1处的近似值yi+1时只用到了前面一个点xi处的信息。多个点2)多步法：计算节点xi+1处的近似值yi+1时用到了前面处的信息一、问题引入但至少有一个r+1次的多项式不能准确成立，则称计算公式是r阶的。记为准确值，为计算的近似值，称当公式是r阶时，有称C为渐进误差常数。为近似值yi+1的局部截断误差定义1

若对，微分方程的数值解法的计算公式均准确成立，二欧拉方法二、欧拉方法Euler方法的精度较低，因此在实际计算时，Euler方法没有多大的实用性，但由于Euler方法的推导比较简单，而且能说明一般的数值计算公式构造时的一些技巧及思想，因此我们从Euler法开始介绍常微分方程的数值方法。1.Euler方法的推导对微分方程（10.1），将[a,b]n等分，步长，节点又设方程（10.1）存在唯一充分连续可微的解y(x)，现介绍三种推导Euler公式的方法。二、欧拉方法（1）Taylor展开在点xi处将y(xi+1)进行Taylor展开有省去误差项得到近似公式局部截断误差（10.3）设yi

是y(xi)的近似值，则对方程（10.1）有差分方程二、欧拉方法(10.4)(10.4)称为求解初值问题(10.1)的Euler方法。从(10.3)(10.4)可知，Euler方法是一阶的单步显式法。推导公式(10.4)还可以用其它的方式。（2）数值微分方法其基本思想是用差商代替微商，对充分小的h有因此有近似公式二、欧拉方法从而有差分方程（10.4）。（3）数值积分法在上对积分有对右端的积分式，利用数值积分的左矩形公式有从而得方程（10.4）二、欧拉方法对差分方程（10.4），如图10.1所示，设y=y(x)是（10.1）的解，线段AB为曲线在点A处的切线，差分方程（10.4）的几何意义是利用线段AB来近似代替曲线段AD。局部截断误差为线段BD，而且BD的值与h的大小和的值有关。如图10.2所示，Euler方法的几何意义是用折线段来近似代替方程（10.1）的解二、欧拉方法3.改进的Euler法(1)中点方法

在区间上对积分有对右端的积分项使用数值积分的中矩形公式有由此有差分方程二、欧拉方法（10.5）局部截断误差差分方程（10.5）称为中点公式，中点公式（10.5）是二步二阶显式方法。(2)梯形方法在上对积分有对右端积分项使用数值积分公式中的梯形公式，有二、欧拉方法由此，有差分方程（10.6）局部截断误差为差分方程（10.6）称为梯形公式，它是单步二阶方法。由于（10.6）式的两端均含有yi+1的信息，它是关于未知量yi+1的一个隐函数方程，称之为隐式方法。二、欧拉方法对梯形公式常采用迭代法来求解yi+1的值，其迭代格式为（10.7）定理2

若f(x,y

)在区域上关于y满足Lipschitz条件（10.2），且，则迭代公式（10.7）生成的序列收敛于(3)Euler预测—校正法再用（10.7）式的梯形方法迭代一次得到的值作为yi+1的近似值,应用中常将（10.4）式Euler方法求出的解作为初值，二、欧拉方法例1设有微分方程

分别用Euler法和Euler预测－校正法作数值计算，并比较其计算结果（取h=0.1）这种方法称为Euler预测—校正法，计算公式为称（10.8）为预测式，（10.9）为校正式。二、欧拉方法解1）微分方程的解析解为2）Euler方法的计算公式为3）Euler预测－校正法的计算公式为计算结果如表10.1所示。二、欧拉方法表10.1二、欧拉方法clearall;clc%预测-校正法symsx

y;a=0;b=1;h=0.1;y0=1;formatlongfun=y-2*x/y;Hfun=@ImprovedEuler;[s,S]=feval(Hfun,fun,a,b,h,y0)x=0:0.1:1;f=sqrt(1+2.*x);plot(x,f,'r');holdonplot(x,s,'b--',x,S,'b+');运行结果：Ivalue=1至4列

1.0000000000000001.0959090909090911.1840965692429971.2662013608757765至8列

1.3433601514839991.4164019285369091.4859556024156691.5525140913261469至11列

1.6164747827520581.6781663636751861.737867401035414莱昂哈德·欧拉（LeonhardEuler，1707.4.15～1783.9.18）数学家、自然科学家，生于瑞士的巴塞尔。13岁时入读巴塞尔大学，15岁大学毕业，16岁获得硕士学位。1731年接替丹尼尔·伯努利成为物理教授。18世纪数学界最杰出的人物之一，他不但为数学界作出贡献，更把整个数学推至物理的领域。是数学史上最多产的数学家，平均每年写出八百多页的论文。主要成就：创立函数的符号；创立分析力学；解决了柯尼斯堡七桥问题；给出各种欧拉公式；主要作品：《无穷分析引论》《微分学原理》《积分学原理》二、欧拉方法三Runge-Kutta方法三、Runge-Kutta方法Runge-Kutta方法的基本思想是利用函数f(x,y

)在某些点上函数值的线性组合来计算y(xi+1)处的近似值yi+1。构造近似公式时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前几项重合，从而使近似公式达到所需的阶数，这样既避免了计算函数f(x,y

)的偏导数，又提高了精度。1.Runge-Kutta方法的构造Runge-Kutta方法的一般形式(10.10)三、Runge-Kutta方法其中是待定参数。由（10.10）得一般每步需要计算p次f(x,y)的值，故称之为p级Runge-Kutta方法。假设y=y(xi)是准确的，（10.11）若则称方法是m阶的，现以p=2为例来推导Runge-Kutta方法，当p=2时有因此对给定的，则（10.10）式为求解微分方程（10.1）的显式单步法，三、Runge-Kutta方法（10.12）将（10.12）中的K2在(xi

,yi

)处进行Taylor展开。对y(xi+1)在xi处进行Taylor展开有三、Runge-Kutta方法使（10.13）式中右端的h和h2的系数恒为0。为使局部截断误差Ti

的阶数尽量地高，应选择适当的参数即取（10.14）

（10.15）以为自由参数有此时局部截断误差为三、Runge-Kutta方法对(10.15)式，选取不同的值可得到相应的二级二阶Runge-Kutta

公式，取，得中点公式取，得Heun（赫恩）公式取，得Euler预测－校正公式三、Runge-Kutta方法2.四阶经典Runge-Kutta公式（10.17）（10.17）也称为标准四阶Runge-Kutta方法，它常用来作线性多步法的启动值计算，（10.17）式也可写成三、Runge-Kutta方法例2

用经典四阶Runge-Kutta方法求解三、Runge-Kutta方法解：取步长h=0.2，计算公式为三、Runge-Kutta方法计算结果如下表所示。与例1的Euler方法计算结果相比较，尽管Runge-Kutta方法的步长放大了一倍，而其数值解的精度还是比Euler法高。三、Runge-Kutta方法clc;clearall;formatlongsymsx

y;fun=y-2*x/y;%fun=x^(-2)*(x*y-y^2)a=0;b=1;h=0.2;y0=1;n=[(b-a)/h]+1;x1=a:h:b;Hfun=@FourRungeKutta;Ivalue=feval(Hfun,fun,a,b,h,y0);运行结果：s=1至4列

1.0000000000000001.1832292874453071.3416669298526071.4832814583502625至6列

1.6125140416775271.732141882691193三、Runge-Kutta方法Runge-Kutta法简练易于编制程序，而且是单步法。其计算也具有数值稳定的优点。但每一次的计算量都比较大，而且公式的局部截断误差也难以求得。所以从计算量来看，一般来说低精度问题常采用低阶Runge-Kutta公式，高精度的问题采用高阶的Runge-Kutta公式，但四阶以上的Runge-Kutta公式很少被使用。3.步长的选取步长的选取在数值解法中非常重要，步长过大，每步计算产生的局部截断误差也较大。若步长取得较小，虽然每步计算的截断误差较小，但在一定的求解范围内需要完成的计算步骤就较多。这不仅增加了计算量，而且还会造成舍入误差的累积。应用中的一种有效措施是在计算的过程中自动调整步长，即变步长技巧。现利用Richardson外推法来构造变步长的技巧。三、Runge-Kutta方法设计算公式是p阶的，从yi出发先取步长为h，经过一步计算得出的数值解记为，局部截断误差记为然后将步长折半，取步长为，经过两步计算，从yi

出发计算得出的数值解记为，其局部截断误差为上的值有关，但我们可以近似认为其中的渐近误差常数与在三、Runge-Kutta方法故有以