1、在连续系统的仿真中，主要的计算工作是求解一阶微分方程 y =f(x,y)y（ x0） = y0解析法只能用来求解一些特殊类型的方程，实际仿真问题中归结出来的微分方程主要靠数值解法。,由于实际运算只能完成有限项或有限步运算，因此要将有些需用极限或无穷过程进行的运算有限化，对无穷过程进行截断，这样产生的误差成为截断误差。 根据实际情况建立的数学模型往往难以求解。通常需要通过近似替代，将所求解的数学模型简化为易求解的数值计算问题后再进行求解。数学模型的理论解与数值计算问题的精确解之间的误差称为截断误差。这是计算方法本身带来的误差，所以也成为方法误差。,2019/3/16,3,得到高精度方法的一个直接
2、想法是利用Taylor展开假设式 y =f(x,y) (axb) 中的 f(x,y) 充分光滑,将y(xi+1)在x i点作Taylor展开，若取右端不同的有限项作为y(xi+1)的近似值，就可得到计算y(xi+1)的各种不同截断误差的数值公式。 例如：取前两项可得到,3.2 龙格库塔方法,2019/3/16,4,其中,P阶泰勒方法,若取前三项，可得到截断误差为O(h3)的公式,类似地，若取前P+1项作为y(xi+1)的近似值，便得到,2019/3/16,5,显然p=1时,y i+1=y i+hf(xi,y i)它即为我们熟悉的Euler方法。当p2时,要利用泰勒方法就需要计算f(x,y)的高
3、阶微商。这个计算量是很大的，尤其当f(x,y)较复杂时，其高阶导数会很复杂。因此，利用泰勒公式构造高阶公式是不实用的。但是泰勒级数展开法的基本思想是许多数值方法的基础。,R-K方法不是直接使用Taylor级数,而是利用它的思想,2019/3/16,6,龙格-库塔(R-K)法的基本思想,Euler公式可改写成,则yi+1的表达式与y(xi+1)的Taylor展开式的前两项完全相同,即局部截断误差为O(h2)。,Runge-Kutta 方法是一种高精度的单步法,简称R-K法,2019/3/16,7,同理，改进Euler公式可改写成,上述两组公式在形式上共同点:都是用f(x,y)在某些点上值的线性组
4、合得出y(xi+1)的近似值yi+1, 且增加计算的次数f(x,y)的次数,可提高截断误差的阶。如欧拉法:每步计算一次f(x,y)的值,为一阶方法。改进欧拉法需计算两次f(x,y)的值，为二阶方法。,局部截断误差为O(h3),2019/3/16,8,于是可考虑用函数f(x,y)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求高阶导数,又提高了计算方法精度的阶数。或者说,在xi,xi+1这一步内多计算几个点的斜率值，然后将其进行加权平均作为平均斜率，
5、则可构造出更高精度的计算格式，这就是龙格库塔（Runge-Kutta）法的基本思想。,一般龙格库塔方法的形式为,2019/3/16,9,其中ai,bij,ci为待定参数，要求上式yi+1在点(xi,yi)处作Tailor展开，通过相同项的系数确定参数。,称为P阶龙格库塔方法。,10,Runge-Kutta方法的推导思想,对于常微分方程的初值问题,的解y=y(x)，在区间xi, xi+1上使用微分中值定理，有,即,2019/3/16,11,引入记号,就可得到相应的Runge-Kutta方法,2019/3/16,12,如下图,即,则上式化为,即Euler方法,Euler方法也称为一阶Runge-K
6、utta方法,2019/3/16,二阶龙格库塔法在xi, xi+1上取两点xi和xi+a2= xi +a2h,以该两点处的斜率值K1和K2的加权平均(或称为线性组合)来求取平均斜率k*的近似值K，即,式中:K1为xi点处的切线斜率值 K1 =hf(xi, yi)=hy(xi)K2为xi +a2h点处的切线斜率值,比照改进的欧拉法,将xi+a2视为xi+1，即可得,2019/3/16,13,确定系数 c1、c2、a2、b21 ，可得到有2阶精度的算法格式,2019/3/16,14,因此,将y(xi+1)在x=xi处进行Taylor展开：,将 在x=xi处进行Taylor展开：,2019/3/16
7、,15,K1 =hf(xi, yi),2019/3/16,16,这里有 4 个未知数，3 个方程。,存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库塔格式。,令,对应项的系数相等，得到,2019/3/16,17,注意到， 就是二阶龙格 - 库塔公式，也就是改进的欧拉法。,因此，凡满足条件式有一簇形如上式的计算格式，这些格式统称为二阶龙格库塔格式。因此改进的欧拉格式是众多的二阶龙格库塔法中的一种特殊格式。,若取 ，就是另一种形式的二阶龙格 - 库塔公式。,2019/3/16,18,此计算公式称为变形的二阶龙格库塔法。式中为区间 的中点。也称中点公式。,Q:为获得更高的精度，应该如何进一步
8、推广？,2019/3/16,19,二级R-K方法是显式单步式,每前进一步需要计算两个函数值。由上面的讨论可知,适当选择四个参数c1,c2,a2, b21，可使每步计算两次函数值的二阶R-K方法达到二阶精度。能否在计算函数值次数不变的情况下,通过选择不同的参数值,使得二阶R-K方法的精度再提高呢?答案是否定的！无论四个参数怎样选择,都不能使公式的局部截断误差提高到三阶。这说明每一步计算两个函数值的二阶R-K方法最高阶为二阶。若要获得更高阶得数值方法,就必须增加计算函数值的次数。,三阶龙格库塔法,2019/3/16,20,为进一步提高精度，在区间xi, xi+1上除两点xi和xi+a2= xi +
9、a2h,以外，再增加一点xi+a3= xi +a3h ，用这三点处的斜率值K1、K2和K3的加权平均得出平均斜率K*的近似值K，这时计算格式具有形式:,2019/3/16,21,同理推导二阶公式，将y(xi+1)和yi+1在x=xi处进行Taylor展开，使局部截断误差达到O(h4)，使对应项的系数相等，得到系数方程组：,参数的选择不唯一，从而构成一类不同的三阶R-K公式，下面给出一种常用的三阶R-K公式，形似simpson公式：,2019/3/16,22,2019/3/16,23,四阶(经典)龙格库塔法,如果需要再提高精度，用类似上述的处理方法，只需在区间xi,xi+1上用四个点处的斜率加权
10、平均作为平均斜率K*的近似值，构成一系列四阶龙格库塔公式。具有四阶精度，即局部截断误差是O(h5)。推导过程与前面类似，由于过程复杂，这里从略，只介绍最常用的一种四阶经典龙格库塔公式。,2019/3/16,24,K1=hf (xi, yi)K2=hf (xi+a2h, yi+b21K1)K3=hf (xi+a3h, yi+b31K1+b32K2)K4=hf (xi+a4h, yi+b41K1+b42K2+b43K3) 其中c1、c2、c3、c4、a2、a3、a4、b21、b31、b32、b41、b42、b43均为待定系数。,这里K1、K2、K3、K4为四个不同点上的函数值，分别设其为,设,yi
11、+1=yi+c1K1+c2K2+c3K3+c4K4,2019/3/16,25,类似于前面的讨论，把K2、K3、K4分别在xi点展成h的幂级数，代入线性组合式中，将得到的公式与y(xi+1)在xi点上的泰勒展开式比较，使其两式右端直到h4的系数相等，经过较复杂的解方程过程便可得到关于ci，ai，bij的一组特解a2=a3=b21=b32=1/2b31=b41=b42=0a4=b43=1c1=c4=1/6c2=c3=1/3,26,四阶(经典)Runge-Kutta方法,2019/3/16,27,例1. 使用高阶R-K方法计算初值问题,解:,(1) 使用三阶R-K方法,2019/3/16,28,其余
12、结果如下:,(2) 如果使用四阶R-K方法,i xi k1 k2 k3 yi1.0000 0.1000 0.1000 0.1103 0.1256 1.1111 2.0000 0.2000 0.1235 0.1376 0.1595 1.24993.0000 0.3000 0.1562 0.1764 0.2092 1.42844.0000 0.4000 0.2040 0.2342 0.2866 1.66645.0000 0.5000 0.2777 0.3259 0.4163 1.9993,2019/3/16,29,其余结果如下:,i xi k1 k2 k3 k4 yi1.0000 0.1000 0
13、.1000 0.1103 0.1113 0.1235 1.1111 2.0000 0.2000 0.1235 0.1376 0.1392 0.1563 1.25003.0000 0.3000 0.1562 0.1764 0.1791 0.2042 1.42864.0000 0.4000 0.2040 0.2342 0.2389 0.2781 1.66675.0000 0.5000 0.2777 0.3259 0.3348 0.4006 2.0000,2019/3/16,【例】 已知一阶系统的微分方程为：，初始条件 ，取仿真步长h=0.1，分别用欧拉法、梯形法和龙格库塔法计算该系统仿真第一步的值
14、。,解：原方程可变为:,即,数值积分公式应用,（1）用欧拉法计算根据欧拉公式，将函数表达式及其初始值代入后，可得该系统仿真第一步的值：,数值积分公式应用,（2）用梯形法计算：根据预报校正公式，将函数表达式及其初始值代入后，可得仿真第一步的值。用预报公式求起始值：,数值积分公式应用,再用校正公式得到系统仿真第一步的值：,数值积分公式应用,二阶龙格-库塔公式,（3）用二阶龙格库塔法计算根据公式先计算出两个系数，再计算仿真第一步的值：,数值积分公式应用,则系统仿真第一步的值为：,数值积分公式应用,四阶龙格库塔（RungeKutta）法,（4）用四阶龙格库塔公式计算 根据公式先计算出4个系数，再计算仿
15、真第一步的值：,数值积分公式应用,数值积分公式应用,则系统仿真第一步的值为：,数值积分公式应用,从上述结果可以看出:对于同一个系统进行仿真计算时，其值的精度是随着数值积分公式的变化而改变的，其中欧拉法计算精度最低，其次为梯形法和二阶龙格库塔法，四阶龙格库塔法计算精度最高。,数值积分公式应用,例2：用matlab演示,龙格-库塔法的误差估计,一个高精度的仿真方法必须将步长控制作为手段。实现步长控制涉及局部误差估计和步长控制策略两方面的问题。,龙格-库塔法的误差估计,RK方法的误差估计通常是设法找一个低一阶的RK公式，将两个公式计算结果之差作为估计误差。例如Runge-Kutta-Fehlberg
16、法的计算公式是RKF1-2公式用另一个一阶公式 来估计误差,龙格-库塔法的误差估计,RKM3-4公式：误差估计式用3 阶，计算公式为4 阶。RKM3-4公式误差估计公式,龙格-库塔法的步长控制,龙格-库塔法的误差估计和步长控制策略的基本思想是：每积分一步都设法估计出本步的积分误差en，然后判断是否满足允许误差E，据此选择相应的步长控制策略。每一步的局部误差通常取以下形式 en=En/(|yn|+1) 其中|yn|是利用误差估计式计算出的本步的估计误差。当|yn|较大时， en是相对误差，当|yn|较小时， en 是绝对误差。这样作的目的是避免当y 的值很小时， en变得过大。,仿真模型的运行速
17、度与实际系统运行速度一致，称为 实时仿真。一般方法难以满足实时仿真的需要：所得 模型的执行速度较慢；机理也不符合实时仿真的需要 。 假设对一般形式的系统进行仿真：以RK-2为例进行分析，其公式为,实时仿真,假定在h/2的时间内计算机刚好计算一次右端函数f,则计算分为两步：1 在tk时刻利用当前的un、yn计算K1；2 在tn+h/2时刻计算K2,此时un+1无法得到，但实时仿真除了要满足执行速度的要求外，还要求实时接收外部输入，并实时得到输出。,此种情况下，解决的方法有两个：对un+1进行预报（增大仿真误差）或仿真延迟半个计算步距。后者的计算流程如下可见，后种方法的输出也会延迟半个计算步距，为
18、了克服这个缺陷，人们提出了如下形式的实时二阶RK法,实时RK-2公式：,其计算流程：假定在h/2的时间内计算机也刚好计算一次右端函数f,则计算也分为两步：1 在tn时刻利用当前的un、yn计算K1；2 在tn+h/2时刻计算K2,此时un+1/2可以得到，不会引入新的误差，可实时得到yn+1,2019/3/16,51,由上节分析常微分方程数值解法稳定性问题的方法，可得到各阶Runge-Kutta公式的稳定性条件：,二阶,与欧拉预估校正公式一致,三阶,四阶,龙格库塔方法的稳定性条件,2019/3/16,52,龙格库塔方法的推导基于Taylor展开方法，因而它要求所求的解具有较好的光滑性。如果解的
19、光滑性差，那么，使用四阶龙格库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在实际计算时，应当针对问题的具体特点选择合适的算法。对于光滑性不太好的解，最好采用低阶算法而将步长h 取小。,2019/3/16,53,前面已经看到,二阶、四阶R-K方法可分别达到最高阶数2阶、4阶，但是N阶R-K方法的最高阶却不一定是N阶。R-K方法的级数表示公式中计算函数值f 的次数。Butcher于1965年给出了R-K方法计算函数值f 的次数与可达到的最高精度阶数之间的关系表，如下：由表可见，四级以下R-K的方法其最高阶数与计算f 的次数一致，对m阶R-K公式，当m4，虽然计算f 的次数增加，但是方法阶数不一定增加。因此四阶R-K公式是应用最为广泛的公式。,2019/3/16,54,顺便指出，当常微分方程中的f (t,u)与u无关时，常微分方程初值问题便简化为计算定积分的问题。这时，第一节介绍的Euler方法就是求定积分的矩形公式，改进的Euler方法就是求定积分的梯形公式，而三阶Runge-Kutta公式就是计算定积分的Simpson公式，它们的精度也是一致的。,
