5.2龙格库塔法
龙格-库塔方法（Runge-Kutta）
龙格-库塔⽅法（Runge-Kutta）龙格-库塔⽅法（Runge-Kutta）3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的⼀般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分⽤不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值⽅法，若⽤显式单步法(3.2.2)当，即数值求积⽤左矩形公式，它就是Euler法(3.1.2)，⽅法只有⼀阶精度，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的⼀种近似，计算时要⽤⼆个右端函数f的值，但⽅法是⼆阶精度的.若要得到更⾼阶的公式，则求积分时必须⽤更多的f值，根据数值积分公式，可将(3.2.1)右端积分表⽰为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)⼀样，⽤前⾯已算得的f值表⽰为(3.2.3)，⼀般情况可将(3.2.2)的表⽰为(3.2.4)其中这⾥均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K⽅法.它每步计算r个f值(即)，⽽k由前⾯(i-1)个已算出的表⽰，故公式是显式的.例i如当r=2时，公式可表⽰为(3.2.5) 其中.改进Euler 法(3.1.11)就是⼀个⼆级显式R-K ⽅法.参数取不同的值，可得到不同公式.3.2.2 ⼆、三级显式R-K ⽅法对r=2的显式R-K ⽅法(3.2.5)，要求选择参数，使公式的精度阶p 尽量⾼，由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6) 令，对(3.2.6)式在处按Taylor 公式展开，由于将上述结果代⼊(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯⼀.因r=2，由（3.2.5）式可知，于是有解(3.2.8)它表明使(3.2.5)具有⼆阶的⽅法很多，只要都可得到⼆阶精度R-K⽅法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)及中点公式(3.2.9)是两个常⽤的⼆级R-K⽅法，注意⼆级R-K⽅法只能达到⼆阶，⽽不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个⽅程,加上(3.2.7)的三个⽅程，共计六个⽅程求4个待定参数，验证得出是⽆解的.当然r=2,p=2的R-K⽅法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，⼀般不再给出.对r=3的情形，要计算三个k值，即其中将按⼆元函数在处按Taylor公式展开，然后代⼊局部截断误差表达式，可得可得三阶⽅法，其系数共有8个，所应满⾜的⽅程为这是8个未知数6个⽅程的⽅程组，解也是不唯⼀的，通常.⼀种常见的三级三阶R-K⽅法是下⾯的三级Kutta⽅法：(3.2.11)附：R-K 的三级Kutta ⽅法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满⾜c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K ⽅法及步长的⾃动选择利⽤⼆元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K ⽅法，其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，⽽33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定，Taylor 展开分析局部截断误差，使得精度达到四阶，即误差为5()O h 。
龙格-库塔法,求解常微分方程
隆格库塔法求解常微分方程摘要科学技术中常常需要求解常微分方程的定解问题，这里问题最简单的形式，是本章将着重考察的一阶方程的初值问题.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法求解.本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题，采用了精度较高的经典的四阶隆格库塔法,然后通过对实例运用Matlab编程进行计算求解，为了体现计算结果的精确性和方法的优越性，再采用了欧拉法和预估较正法对实例进行计算求解作为比较.通过比较三种方法的计算精度，发现四阶经典龙格-库塔方法的误差最小，预估较正法其次，欧拉方法误差则比较大.最后通过选取不同的步长，研究了不同的步长对隆格库塔法求解常微分方程初值问题的计算精度的影响.总之，本文全面分析了隆格库塔法在求解常微分方程的应用，相比与其他的数值解法，隆格库塔法计算精度较高，收敛性较好，其中四阶的隆格库塔法的效率最高，精度也最高.关键词：四阶隆格库塔法；欧拉法；预估较正法；一阶常微分方程；MatlabRunge Kutta Method For Solving Ordinary Differential EquationsABSTRACTProblem solving ordinary differential equations are often needed in science andtechnology. the problem in the simplest form is the initial value problem of first order equations in this chapter ,which will be discussed. Although there are various analytical methods for solving ordinary differential equations, the analytical method can only be used to solve some special types of equations.differential equations can be summed up the actual problems whichThis paper discusses the initial value problem of Runge Kutta Barclays by solving a differential equation, using the four order Runge Kutta method with high accuracy.for instance through classic Matlab programming calculation, the superiority in order to accurately and reflect the calculation result, then the Euler method and the prediction correction method for instance by calculation through the calculation precision. The comparison of three kinds of methods, found that the error of four order Runge Kutta method of minimum, prediction correction method secondly, Euler method error is relatively large. Finally, by selecting different step, study the affect the calculation accuracy of different step of Runge Kutta method to solve initial value problems of ordinary differential equations.In short, this paper comprehensively analyzes the application of Runge Kutta method for solving ordinary differential equations, compared with the numerical solution of other, higher accuracy Runge Kutta method, good convergence, the Runge Kutta method of order four of the highest efficiency and its precision is the highest.Key words: Four order Runge Kutta method; Euler method; prediction correction method; first order ordinary differential equations; Matlab目录1 问题的提出 (1)1.1 问题背景............................................... . (1)1.2 问题的具体内容 (1)2 问题假设 (2)3 符号系统 (2)4 问题的分析 (3)4.1 欧拉格式 (3)4.2 预估较正法 (3)4.3 四阶隆格库塔法的格式 (4)5 模型的建立与求解 (4)5.1 隆格库塔法的基本原理 (4)5.1.1 Taylor级数 (4)5.1.2 隆格库塔法的基本思想 (4)5.1.3 四阶的隆格库塔法 (5)5.2 其他求解常微分方程边值问题算法的简介 (6)5.3 模型求解 (8)5.3.1 运用MATLAB软件对模型求解结果及析 (8)6 模型的评价 (16)7 课程设计的总结与体会 (16)参考文献 (17)附录 (18)一、问题的提出1.1 问题背景：科学技术中常常需要求常微分方程的定解问题，微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：0(,)()dy f x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩ （1）其中a ，b 为常数.虽然求解此类微分方程有各种各样的解析方法，但解析方法只能用于求解一些特殊类型方程，实际问题中归结出来的微分方程主要靠数值解法求解.因为一阶常微分方程简单但又是求解其他方程的基础，所以发展了许多典型的解法.本文着重讨论一类高精度的单步法——隆格库塔法，并且运用四阶的隆格库塔格式来求解初值问题，并且通过实例运用四阶的隆格库塔格式来求解初值问题，同时与显式与隐式的Euler 格式求解出的结果进行精度比较.1.2 问题的具体内容实例一：在区间[0,1]上采用经典的四阶隆格库塔方法求解微分方程1(0)1dy y x dx y ⎧=-++⎪⎨⎪=⎩，其精确解为x y x e -=+，步长为0.5,然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.实例二：在区间[0,1]上用经典的四阶龙格库塔方法求解初值问题 (0)1x dy xe y dx y -⎧=-⎪⎨⎪=⎩， 其精确解为21(2)2xx e -+，然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，并且探究选取不同的步长对计算结果精度的影响.二、问题假设2.1 假设数值方法本身的计算是准确的.2.2 假设选取的步长趋于0时计算的结果会收敛到微分方程的准确解.2.3 假设步长的增加不会导致舍入误差的严重积累.三、符号系统3.1 符号说明符号含义h选取的步长*K平均斜率p精度的阶数∆前后两次计算结果的偏差y第n个节点的实验值n()y x第n个节点的精确值nδ实验值与精确值的绝对误差四、问题的分析问题要求运用隆格库塔算法来求解一阶微分方程的初值问题，针对前面提出的实例，本文先用经典的四阶隆格库塔法来求解上面的微分方程，为了体现隆格库塔法的优越性，同时用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，分析在选取不同的步长时，求解结果的精度变化如何.下面是欧拉法，预估校正法以及经典的四阶隆格库塔法的计算公式.4.1欧拉格式（1）显式欧拉格式1(,)n n n n y y hf x y +=+ （2） 局部截断误差22211()()()()22n n n h h y x y y y x o h ξ++''''-=≈= （3） （2）隐式欧拉格式111(,)n n n n y y hf x y +++=+ （4）局部截断误差2211()()()2n n n h y x y y x o h ++''-≈-= (5) 4.2 预估校正法预估 1(,)n n n n y y hf x y +=+ (6)校正 111[(,)(,)]2n n n n n n h y y f x y f x y +++=++ (7) 统一格式 1[(,)(,(,))]2n n n n n n n n h y y f x y f x h y hf x y +=++++ (8) 平均化格式 11(,),(,),1().2p n n n c n n p n p c y y hf x y y y hf x y y y y ++⎧⎪=+⎪=+⎨⎪⎪=+⎩ (9)4.3 四阶龙格库塔方法的格式（经典格式）112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩(10)五、模型的建立与模型求解5.1 隆格库塔法的基本原理隆格库塔法是一种高精度的单步法，这类方法与下述Taylor 级数法有着紧密的联系.5.1.1 Taylor 级数设初值问题 '00(,)()y f x y y x y ⎧=⎨=⎩有解，按泰勒展开，有2'''1()()()()....2n n n n h y x y x hy x y x +=+++； （11） 其中()y x 的各阶导数依据所给方程可以用函数f 来表达，下面引进函数序列(,)j f x y 来描叙求导过程，即(0)(1)'(0)''(1)(1)(1)'''(2),f f y f f y f f x x f f y f f x y ∂∂=≡=+≡∂∂∂∂=+≡∂∂ （12）(2)(2)()(1)j j j j f f y f f x y ---∂∂=+≡∂∂ （13） 根据上式，结果导出下面 Taylor 格式2'''()1...2!!pp n n nn n h h y y hy y y p +=++++ （14）其中一阶Taylor 格式为： '1n n n y y hy +=+提高Taylor 格式的阶数p 即可提高计算结果的精度，显然，p 阶Taylor 格式的局部截断误差为：11(1)1(1)!p p n n h y y y p ζ+++-+=+ （15） 因此它具有p 阶精度.5.1.2 隆格库塔方法的基本思想隆格库塔法实质就是间接地使用Taylor 级数法的一种方法,考察差商1()()n n y x y x h+- 根据微分中值定理，这有'1()()()n n n y x y x y x h h θ+-=+ （16） 利用所给方程 '(,)y f x y =1()()(,())n n n n y x y x hf x h y x h θθ+=+++ （17） 设 平均斜率*(,())n n K f x h y x h θθ=++，由此可见，只要对平均斜率一种算法，便相应地可以导出一种计算格式.再考察改进的Euler 格式，它可以改写成平均化的形式：1121211()2(,)(,)n n n n n n h y y K K K f x y K f x y hK ++⎧=++⎪⎪=⎨⎪=+⎪⎩（18） 这个处理过程启示我们，如果设法在1(,)n n x x +内多预测几个点的斜率值，然后将它们加权平均作为平均斜率，则有可能构造具有更高精度的计算格式，这就是隆格库塔法的基本思想.5.1.3 四阶的隆格库塔法为了方便起见，本文主要运用经典的隆格库塔算法-四阶隆格库塔格式.其格式如下：112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩（19）下面为其具体的算法流程图：5.2 其他求解常微分边值问题算法的简介5.2.1欧拉数值算法（显式）微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：(,)()dyf x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩(20) 其中a ，b 为常数.开始输入x 0,y 0,h,N x 1=x 0+hk 1=f(x 0,y 0),k 2=f(x 0+h/2,y 0+hk 1/2)k 3=f(x 0+h/2,y 0+hk 2/2),k 4=f(x 1,y 0+hk 3)y 1=y 0+h(k 1+2k 2+2k 3+k 4)/6n=1输出x 1,y 1n=N? 结束n=n+1x 0=x 1,y 0=y 1否是图5.1 龙格-库塔法流程图因为其简单但又是求解其他方程的基础，所以发展了许多典型的解法.所有算法中的f 就是代表上式中(,)f x y ，而y f 表示(,)f x y y ∂∂，x f 表示(,)f x y x∂∂. 简单欧拉法是一种单步递推算法.简单欧拉法的公式如下所示：1(,)n n n n y y hf x y +=+ (21)简单欧拉法的算法过程介绍如下:给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，计算1(,)k k k k y y hf x y +=+5.2.2欧拉数值算法（隐式）隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示：111(,)n n n n y y hf x y +++=+ (22)隐式欧拉法是一阶精度的方法，比它精度高的公式是：111[(,)(,)]2n n n n n n hy y f x y f x y +++=++ (23)隐式欧拉的算法过程介绍如下.给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，用牛顿法或其他方法求解方111(,)k k k k y y hf x y +++=+得出1k y +.5.2.3 欧拉预估-校正法改进欧拉法是一种二阶显式求解法，其计算公式如下所示：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++11(,)[(,)(,)]2n n n n n n n n t y hf x y h y y f x y f x t ++=+⎧⎪⎨=++⎪⎩ (24)四阶龙格-库塔法有多种形式，除了改进的欧拉法外还有中点法.中点法计算公式为：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++ (25)5.3 模型求解5.3.1运用MATLAB 软件对模型求解结果及分析用欧拉法、改进的欧拉格式、经典的四阶龙格库塔法来求解常微分方程的边值问题，并且比较其精度（具体的MATLAB 源程序见附录） 以下进行实例分析：实例一. 1(0)1dyy x dx y ⎧=-++⎪⎨⎪=⎩由题可知精确解为：x y x e -=+，当x=0时，y(x)=0.在这里取步长h 为0.1， 通过MATLAB 程序的计算，相应的结果如下：表5-1 步长为0.1时各方法的预测值与精确值的比较(精确到6位小数)初值 Euler 法 相对误差 预估校正法 相对误差 经典四阶库 相对误差精确值 0 -- -- -- -- -- -- 1.00000 0.1 1.00910 0.00424 1.00500 0.00016 1.00484 0.00000 1.00484 0.2 1.02646 0.00759 1.01903 0.00029 1.01873 0.00000 1.01873 0.3 1.05134 0.01011 1.04122 0.00038 1.04082 0.00000 1.04082 0.4 1.08304 0.01189 1.07080 0.00045 1.07032 0.00000 1.07032 0.5 1.12095 0.01303 1.10708 0.00049 1.10653 0.00000 1.10653 0.6 1.16451 0.01366 1.14940 0.00052 1.14881 0.00000 1.14881 0.7 1.21319 0.01388 1.19721 0.00052 1.19659 0.00000 1.19659 0.8 1.26655 0.01378 1.24998 0.00052 1.24933 0.00000 1.24933 0.9 1.32414 0.01344 1.30723 0.00050 1.30657 0.00000 1.30657 1.01.38558 0.01294 1.36854 0.00048 1.36788 0.000001.36788步长为0.1时的精确值与预测值的比较精确值欧拉法改进欧拉格式四阶龙格库塔轴Y00.10.20.30.40.50.60.70.80.91 1.1 1.2 1.3 1.4X轴图5.2 步长为0.1时各方法的预测值与精确值的比较原函数图像轴YX轴图5.3 步长为0.1时原函数图像在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-2 h=0.05时三个方法与精确值的真值表步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0 -- -- -- -- -- --1.00000 0.05 1.00250 0.00911 1.00125 0.01035 1.00123 0.01037 1.00484 0.10 1.00738 0.01711 1.00488 0.01954 1.00484 0.01958 1.01873 0.15 1.01451 0.02405 1.01076 0.02765 1.01071 0.02770 1.04082 0.20 1.02378 0.03001 1.01880 0.03473 1.01873 0.03479 1.07032 0.25 1.03509 0.03507 1.02889 0.04085 1.02880 0.04093 1.10653 0.30 1.04834 0.03930 1.04092 0.04610 1.04082 0.04619 1.14881 0.35 1.06342 0.04277 1.05480 0.05053 1.05469 0.05063 1.19659 0.40 1.08025 0.04555 1.07044 0.05422 1.07032 0.05432 1.24933 0.45 1.09874 0.04772 1.08775 0.05724 1.08763 0.05734 1.30657 0.50 1.11880 0.04933 1.10666 0.05964 1.10653 0.05975 1.36788 0.55 1.14036 0.05045 1.12709 0.06150 1.12695 0.06161 1.00484 0.60 1.16334 0.05113 1.14895 0.06286 1.14881 0.06298 1.01873 0.65 1.18768 0.05143 1.17219 0.06379 1.17205 0.06391 1.04082 0.70 1.21329 0.05139 1.19674 0.06433 1.19659 0.06445 1.07032 0.75 1.24013 0.05106 1.22252 0.06453 1.22237 0.06465 1.10653 0.80 1.26812 0.05048 1.24949 0.06443 1.24933 0.06455 1.14881 0.85 1.29722 0.04969 1.27757 0.06408 1.27742 0.06419 1.19659 0.90 1.32735 0.04871 1.30673 0.06349 1.30657 0.06361 1.249330.95 1.35849 0.04759 1.33690 0.06272 1.33674 0.06283 1.306571.00 1.39056 0.04634 1.36804 0.06178 1.36788 0.06189 1.36788欧拉格式 改进欧拉格式四阶龙格库塔精确值图5.4 步长为0.05时各方法的预测值与精确值的比较11.051.11.151.21.251.31.351.4X 轴Y 轴原函数图像图5.5 步长为0.05时原函数图像实例2 (0)1xdy xe ydx y -⎧=-⎪⎨⎪=⎩由题可知精确解为：21(2)2x x e -+当x=0时，y(x)=0.在这里取步长h为0.1，通过MATLAB程序的计算，相应的结果如下：步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.1 0.9000 0.0318 0.9096 0.0214 0.9094 0.0216 0.9295 0.2 0.8192 0.0611 0.8359 0.0420 0.8356 0.0424 0.8726 0.3 0.7544 0.0876 0.7761 0.0614 0.7757 0.0619 0.8268 0.4 0.7027 0.1109 0.7277 0.0793 0.7272 0.0799 0.7903 0.5 0.6617 0.1310 0.6886 0.0956 0.6881 0.0963 0.7615 0.6 0.6294 0.1480 0.6572 0.1103 0.6567 0.1110 0.7387 0.7 0.6040 0.1621 0.6320 0.1234 0.6315 0.1241 0.7209 0.8 0.5841 0.1737 0.6116 0.1349 0.6111 0.1355 0.70690.9 0.5686 0.1830 0.5951 0.1450 0.5946 0.1456 0.69591.0 0.5563 0.1905 0.5815 0.1538 0.5811 0.1544 0.6872原函数图像轴Y00.10.20.30.40.50.60.70.80.91X轴图5.6 步长为0.1时原函数图像各方法的预测值与精确值的比较欧拉格式改进的格式四阶龙格库塔精确值轴Y0.10.20.30.40.50.60.70.80.91X轴图5.7 步长为0.1时各方法的预测值与精确值的比较在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-4 步长为0.1时各方法的预测值与精确值的比较(精确到5位小数)步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.050.95000 0.01342 0.95245 0.01088 0.95242 0.01090 0.96292 0.100.90490 0.02650 0.90947 0.02158 0.90942 0.02163 0.92953 0.150.86428 0.03916 0.87067 0.03206 0.87061 0.03213 0.89951 0.200.82774 0.05137 0.83567 0.04228 0.83559 0.04237 0.87256 0.250.79491 0.06307 0.80414 0.05219 0.80405 0.05230 0.84842 0.300.76545 0.07423 0.77576 0.06176 0.77566 0.06188 0.82682 0.350.73904 0.08482 0.75023 0.07096 0.75013 0.07110 0.80754 0.40 0.71541 0.09482 0.72730 0.07977 0.72719 0.07991 0.79035 0.45 0.69427 0.10422 0.70672 0.08817 0.70660 0.08832 0.77505 0.50 0.67539 0.11302 0.68825 0.09615 0.68813 0.09630 0.76146 0.55 0.65855 0.12123 0.67168 0.10370 0.67156 0.10386 0.74940 0.60 0.64352 0.12886 0.65683 0.11084 0.65671 0.11100 0.73871 0.65 0.63012 0.13593 0.64351 0.11756 0.64340 0.11773 0.72925 0.70 0.61819 0.14245 0.63157 0.12388 0.63145 0.12405 0.72087 0.75 0.60754 0.14847 0.62085 0.12981 0.62073 0.12998 0.71347 0.80 0.59805 0.15400 0.61121 0.13537 0.61110 0.13553 0.70691 0.85 0.58957 0.15908 0.60254 0.14058 0.60243 0.14073 0.70109 0.90 0.58198 0.16374 0.59470 0.14545 0.59460 0.14560 0.695930.95 0.57517 0.16802 0.58761 0.15001 0.58751 0.15016 0.691321.00 0.56904 0.17194 0.58117 0.15429 0.58107 0.15443 0.687190.10.20.30.40.50.60.70.80.910.550.60.650.70.750.80.850.90.951X 轴Y 轴各方法下的预测值与精确值的比较欧拉格式改进欧拉格式四阶龙格库塔精确值图5.8 步长为0.05时各方法的预测值与精确值的比较六、模型的评价本文着重讨论了4阶的隆格库塔法来求解微分方程，并且通过两个实例验证了隆格库塔法在求解初值问题的优越性.从上面的实例比较可知，在计算精度上，四阶经典龙格-库塔方法的误差最小，改进欧拉方法其次，欧拉方法误差则比较大，所以四阶经典龙格-库塔方法得到最佳的精度.而在计算量上面，相应地，很明显的四阶经典龙格-库塔方法也是最大，改进欧拉方法其次，欧拉方法计算量最小.这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比.我们在实际运用与操作中，可以根据实际情况，选择这3种方法中的其中一种最适合的，追求精度的话，可以使用四阶经典龙格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一般情况；而欧拉方法更主要的是适用于对y的估计上，而精度则有所欠缺，以上各方法的选择，都取决于具体的情况.七、课程设计的总结与体会本文着重采用隆格库塔法运用MATLAB编程来求解微分方程，相比于欧拉法以及预估校正法，隆格库塔法在提高近似值解的精度上是非常起作用的，而且又具有计算量不大、算法组织容易.其次，每一次的课程设计总是让我学到了更多的知识，不论是C++、SPASS还是MATLAB软件，这些都让我学到了如何解决实际问题的好工具，通过这些工具，是自己能够得到突破和成长.以下是我完成此次课程设计的几点体会：（1）必须学好基础知识，在做的过程中，发现自己有很多东西都不懂，要博学必须从一点一点做起.以往训练得少只是把握的不牢靠.所以做起来感到有点吃力.所以，无论什么学科，一定要打好基础.（2）程序设计要靠多练，多见识，那样形成一种编程思维，我想对我是有很大好处的.尤其像我这种平时学得不扎实的人.（3）做事情要有恒心，遇到困难不要怕，坚决去做.如果做出来了，固然高兴，如果没有做出来也没关系，自己努力了对得起自己就好.同时，把它看做是对自己的锻炼. （4）做程序特别是做大程序是很有趣的.虽然有的问题很难，要花很多时间很多精力，但是那种解决了一个问题时的喜悦足以把付出的辛苦补偿回来.得到一种心里的慰藉.参考文献[1] 李庆杨，王能超，易大义编．数值分析（第四版）[M]：华中科技大学出版社，2006.[2] 姜启源，谢金星，叶俊编．数学模型（第三版）[M]．北京：高等教育出版社，2005[3] 刘琼荪，数学实验[M]，高等教育出版社，2004[4] 王建伟，MATLAB7.X程序设计[M]，中国水利水电出版社，2007[5] 王高雄，周之铭等编.常微分方程（第三版） [M]:高等教育出版社，2006[6]何坚勇编著. 运筹学基础（第二版）[M]. 北京：清华大学出版社，2008.附录附录一：显示欧拉法matlab程序%欧拉法clear allclcx=[];y=[];y1=[];h=0.1;x=0:h:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold ony1(1)=1;for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);迭代n次的后退的欧拉格式matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%迭代n次的后退的欧拉格式clear allclcN=input('请输入迭代次数:');x=[];y=[];y1=[];h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); T=y1(j);for k=1:NT=y1(j-1)+h*f(x(j),T);endy1(j)=T;endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y); fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)附录二：预估校正法matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%预估校正法%欧拉法clear allclcx=[];y1(1)=1;y2=[];y2(1)=1;h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:ny1(j)=y2(j-1)+h*f(x(j-1),y2(j-1));y2(j)=y2(j-1)+(h/2).*[f(x(j-1),y2(j-1))+f(x(j),y1(j))]; endY=[x;y2];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y2,'r-');figure(2)DT=y-y2;plot(x,DT)附录三：四阶龙格库塔法matlab程序%四阶龙格库塔法clear allclcn=length(x);y=[];y1=[];y1(1)=1;for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:nK1=f(x(j-1),y1(j-1));K2=f(x(j-1)+h/2,y1(j-1)+h/2*K1);K3=f(x(j-1)+h/2,y1(j-1)+h/2*K2);K4=f(x(j-1)+h,y1(j-1)+h*K3);y1(j)=y1(j-1)+(h/6)*(K1+2*K2+2*K3+K4); endY=[x;y1];fid=fopen('data1.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)。
龙格库塔算法
龙格库塔算法龙格库塔算法（Runge-Kutta method）是一种常用的数值解微分方程的方法，其基本原理是通过逐步逼近的方式，根据初始条件和微分方程的表达式，计算出方程的近似解。
该方法具有较高的精度和稳定性，在科学计算、物理模拟、工程建模等领域得到广泛应用。
龙格库塔算法的核心思想是将微分方程的解按照一定的步长进行离散化，从而将连续的求解问题转化为离散的迭代过程。
具体来说，龙格库塔算法通过计算函数在一定步长内的平均斜率，来估计下一个点的函数值。
这个平均斜率是通过多次计算函数在不同点上的导数得到的，从而提高了计算的精度。
龙格库塔算法的一般形式可以表示为：k1 = f(tn, yn)k2 = f(tn + h/2, yn + h/2 * k1)k3 = f(tn + h/2, yn + h/2 * k2)k4 = f(tn + h, yn + h * k3)yn+1 = yn + h/6 * (k1 + 2k2 + 2k3 + k4)其中，tn是当前时间点，yn是当前函数值，h是步长，f是微分方程的表达式。
通过多次迭代，可以逐渐逼近微分方程的解。
龙格库塔算法的优点在于其精确度较高，可以通过调整步长来控制计算的精度和效率。
此外，该算法具有较好的数值稳定性，可以有效处理非线性、刚性或高阶微分方程等复杂问题。
因此，在科学和工程计算中，龙格库塔算法被广泛应用于各种数值模拟和求解问题。
需要注意的是，龙格库塔算法并非万能的，对于一些特殊的问题，可能存在数值不稳定性或计算精度不够的情况。
此外，算法的步长选择也需要根据具体问题进行调整，过小的步长会增加计算量，而过大的步长可能导致精度下降。
因此，在使用龙格库塔算法时，需要根据具体问题的特点和要求来选择合适的步长和算法参数，以获得满意的计算结果。
总结起来，龙格库塔算法是一种常用的数值解微分方程的方法，具有较高的精度和稳定性。
通过离散化和迭代的方式，可以逐步逼近微分方程的解。
精品课件-龙格-库塔方法---数学课的讲解
其中 i(xi,xi1)
即
y (x i 1 ) y (x i) h y (i)
2019/5/30
8
引入记号 y(xi1)y(xi)K
K h y (i) h i,f y (i)
龙格-库塔方法---数学课的讲 解
若取前三项，可得到截断误差为O(h3)的公式
y(xi 1 )y(xi) h y(xi)h 2 2y(xi) O (h 3)
yi h(fxi,yi)
h2 2
fx(xi,yi)f(xi,yi)fy(xi,yi)O(h3)
类似地，若取前P+1项作为y(xi+1)的近似值，便得到
yi1yihyih 22 !yi h P P !yi(P) P阶泰勒方法
其中 y if,y if(x i,y i)xfx ffy
y i (fx ffy )xfx x2fx fyfy fy2fx fy (fy )2f
y(xi)h(fxi,yi)h 2 2 !fx ffy O (h3)
将 K 2 h y ( x i a 2 h ) h ( x i f a 2 h ,y i b 2 K 1 ) 1 在x=xi处进行Taylor展开：
2019/5/30
12
K 2 h f ( x i , y i ) a 2 h f x b 2 K 1 f y 1 O ( h 2 ) K1 =hf(xi, yi) h f( x i,y i) a 2 h fx b 2 h f 1 y fO ( h 3 )
确定系数 c1、c2、a2、b21 ，可得到有2阶精度的算法格式
数值计算中的龙格库塔算法
数值计算中的龙格库塔算法龙格库塔算法，又称龙格-库塔算法，是一种数值计算方法，主要用于求解微分方程。
它的好处是通过迭代得到更加精确的数值解，对于很多科学和工程问题，如天体力学、化学反应动力学、电路分析等，都有广泛的应用。
一、初识龙格库塔算法最早提出龙格库塔算法的是瑞士数学家卡尔·龙格和德国数学家马丁·库塔，他们在20世纪初期分别提出了一种求解常微分方程组的方法，后来又被合并为一种更为完善的算法，即现在我们所说的龙格库塔算法。
它的基本思想是将微分方程分解成一系列递推的步骤，通过不断迭代，逐渐逼近准确的解。
龙格库塔算法的核心是求出微分方程在某个时刻的斜率。
一般而言，我们可以使用欧拉法或者梯形法来求解，但这些方法往往会出现舍入误差，导致数值解偏离实际解。
相比之下，龙格库塔算法则将微分方程的初始值向前推进一个尽可能小的步长，通过不断缩小步长的大小进行迭代，在保证精度的同时大大提高了计算效率。
在实际应用中，我们通常会使用四阶龙格库塔算法（RK4）来求解微分方程。
具体做法是先求出微分方程在 $t$ 时刻的斜率$k_1$，然后将$t$ 向前推进一半的步长，求出此时的斜率$k_2$，再用 $k_2$ 推进一半的步长，求出此时的斜率 $k_3$，最后以$k_3$ 推进一个步长，求出微分方程在 $t+h$ 时刻的斜率 $k_4$。
最终的数值解为：$$y_{n+1} = y_n + \frac{1}{6}(k_1+2k_2+2k_3+k_4)h$$其中 $y_{n+1}$ 表示下一个时刻的函数值，$y_n$ 表示当前时刻的函数值，$h$ 表示步长。
这个公式看起来比较复杂，但实际上只是对斜率的加权平均。
通过不断迭代，我们就可以得到越来越精确的解。
二、优缺点及应用与其他数值计算方法相比，龙格库塔算法具有以下优点：1. 高精度：通过四阶跑格库塔公式，可达到高精度计算。
2. 稳定可靠：在每一步均会进行收敛性检验，确保计算结果准确无误。
龙格库塔法
一、高阶泰勒法
假设初值问题
龙格—库塔法 龙格 库塔法
dy = f (t , y ) dt y (a) = α 的解y (t)及f (t , y )足够光滑.
将y (ti +1 )在ti处作n阶泰勒展开, 得
a≤t ≤b
(1)
y′′(ti ) 2 y ( n ) (ti ) n y ( n +1) (ξ i ) n +1 y (ti +1 ) = y (ti ) + y′(ti )h + h +L+ h + h n! 2! (n + 1)! 其中, ti < ξ i < ti +1.
2
i
i
1
3
i
i
2
4
i
i
3
i +1
i
6123 Nhomakorabea4
作业 教材P198 习题3
(2)
(3)
首先将y (ti +1 )在ti处展成幂级数 h2 y (ti +1 ) = y (ti ) + hy′(ti ) + y′′(ti ) + O(h 3 ) 2 将 y′(t ) = f (t , y (t )) y′′(t ) = f t′(t , y (t )) + f y (t , y (t )) f (t , y (t )) 代入上式, 得 h2 y (ti +1 ) = y (ti ) + hf + ( f t + ff y ) + O(h 3 ) (3) 2 其中f , f t , f y′分别表示相应函数在点(ti , y (ti ))处的函数值.
龙格库塔法解轨道参数
龙格库塔法解轨道参数引言龙格库塔法（Runge-Kutta method）是一种常用的数值计算方法，用于求解常微分方程的数值解。
在天体力学中，我们经常需要通过数值方法来计算天体的轨道参数，如轨道椭圆的长短轴、离心率、倾角等。
本文将介绍龙格库塔法在解轨道参数中的应用，并详细探讨该方法的原理和实现过程。
基本原理龙格库塔法是一种迭代求解的方法，在每个时间步长内利用当前的状态来估计下一个状态。
具体而言，龙格库塔法将微分方程的求解问题转化为一个迭代的求解问题，通过逐步迭代来逼近精确解。
在解轨道参数的问题中，我们通常需要根据已知的初始条件以及天体的质量和力学模型来求解天体的轨道参数。
常用的力学模型有开普勒模型和牛顿模型。
龙格库塔法可以根据力学模型的不同进行相应的求解。
开普勒模型下的轨道参数求解步骤一：确定初始条件在使用龙格库塔法求解轨道参数之前，我们需要确定一些初始条件。
这些初始条件包括天体的质量、位置和速度。
步骤二：选择时间步长在求解过程中，我们需要选择一个合适的时间步长。
时间步长越小，计算的精度会越高，但计算的时间会增加。
步骤三：迭代求解利用龙格库塔法进行迭代求解的具体步骤如下：1.根据当前时刻的位置和速度，计算天体在该时刻的加速度。
2.根据当前时刻的位置、速度和加速度，计算下一个时刻的位置和速度。
3.更新当前时刻的位置和速度为新的位置和速度。
4.重复上述步骤，直到达到指定的终止条件。
步骤四：计算轨道参数通过迭代求解，我们可以得到天体在不同时刻的位置和速度。
根据这些位置和速度，我们可以计算出轨道参数，如离心率、倾角、长短轴等。
常用的轨道参数计算公式如下：1.离心率：e=√1+2El2μ(GM⊕)2)2.倾角：i=arccos(ℎzℎ3.长轴：a=−μT22E4.短轴：b=a√1−e2其中，E表示能量，l表示轨道角动量，μ表示标准引力参数，G表示引力常数，M⊕表示地球的质量，ℎ和ℎz分别表示轨道角动量和轨道角动量在z轴上的分量。
82第二节 龙格—库塔法
k
(1)
h h k 若令 yn1 y xn hy xn y xn y xn (2) 2! k! 则 y xn1 yn1 O hk 1
y0 k1 2 k2 hf x0 h 2, y0 k1 2
y0 k3
k4 hf x0 h, y0 k3
y0 k2 2 k3 hf x0 h 2, y0 k2 2
k
x1 x0 h y1 y0 k
数学学院 信息与计算科学系
0.1832292
0.1584376
数学学院 信息与计算科学系
接上图
0.4 0.5 0.5 0.6 0.6 1.341667 1.416026 1.412676 1.482627 1.483281 0.0745394 0.0710094 0.0708400 0.0673253
0.1416245
数学学院 信息与计算科学系
数学学院 信息与计算科学系
由表8-4可见，虽然四阶龙格-库塔方法每步要 计算四次 f 的值，但以h=0.2为步长ቤተ መጻሕፍቲ ባይዱ计算结果就
有5 位有效数字，而欧拉法与预估计-校正方法以
h=0.1为步长的计算结果才具有2 位与3 位有效数字.
如果步长 h 也取0.2，则结果的精度会更低.
即公式(2)为k 阶方法.
数学学院 信息与计算科学系
二、龙格-库塔方法(R-K方法)
R-K方法不是通过求导数的方法构造近似公式, 而是通过计算不同点上的函数值, 并对这些函数值作 线性组合, 构造近似公式, 再把近似公式与解的泰勒 展开式进行比较, 使前面的若干项相同 , 从而使近似 公式达到一定的阶数.
龙格库塔法介绍
h 0.005稳定.
2) 改进欧拉法(预测 — 校正，即二阶R K法)：
yn
1
yn
h
k1 2
k2 2
,
k1
f (xn, yn ) yn,
k2 f (xn h, yn hk1) ( yn hyn )，
即yn1 故
yn
当初值准确即e0 0时，整体误差为en O(h p ).
证明
考察单步法的收敛性归结为验证增量函数(x, y,h)是否
满足Lipschitz条件.
欧拉法 : (x, y,h) f (x, y)，L L.
改进的欧拉法：
yn1
yn
h[ 2
f
(xn, yn )
f
( xn 1,
xn
f
( x,
y ( x))dx
r
h ci
i 1
f
( xn
ih,
y ( xn
ih)).
yn1 yn h(xn, yn, h),
(3.4)
其中
r
(xn, yn, h) ciki ,
(3.5)
i 1
k1 f (xn, yn )，
欧拉法r 1, p 1.改进 欧拉法r 2, p 2.
k1)
k3）
k3 f (xn h, yn hk1 2hk2 )
称为库塔三阶方法.
阶数p和段数r（计算函数值次数）的关系
r12 p1 2
3 4 5 6 7 r≥8 3 4 4 5 6 r-2
常用的经典四阶龙格 库塔方法:
第二节 龙格-库塔方法
h( k1 + k2 ) yn +1 = yn + 2 hk1 1 1 3 3 k1 = f [ xn + ( + )h, yn + )hk2 ] +( + 2 6 4 4 6 hk2 1 3 1 3 k2 = f [ xn + ( − )h, yn + ( − )hk1 + ] 2 6 4 6 4
1 库塔三阶方法 库塔三阶方法
h h k1 = f ( xn , yn ); k2 = f ( xn + , yn + k1 ) 2 2
h yn +1 = yn + ( k1 + 4k2 + k3 ) 6
k3 = f ( xn + h, yn − hk1 + 2hk2 )
四级方法： 四级方法：N = 4 常见的2种四阶方法： 常见的 种四阶方法： 种四阶方法 经典龙格 库塔方法 龙格1经典龙格-库塔方法
1 中点法（修正的Euler法） 取 c1 = 0, c2 = 1, a2 = 1 中点法（修正的 法 2 h h
y n +1 = y n + hf ( x n +
常见的3种二级方法： 常见的 种二级方法： 种二级方法
二阶龙格库塔方法 2 二阶龙格库塔方法
h yn +1 = yn + [ f ( xn , yn ) + f ( xn + h, yn + hf ( xn , yn ))] 2
1 取 c1 = c2 = , a2 = 1 2
2
, yn +
2
f ( x n , y n ))
类似于N 的推导方法， 的推导方法 三级方法： 三级方法：N = 3 类似于 = 2的推导方法，可得到 1 c1 + c 2 + c 3 = 1; c 2 a 2 + c 3 a 3 = ; 2 O ( h4 ) 1 1 2 2 c 2 a 2 + c 3 a 3 = ; c 3 a 2 b32 = 6 3 常见的2种三阶方法： 常见的 种三阶方法： 种三阶方法
(完整版)第二节龙格-库塔方法
因为单步法是 p阶的：h0 , 0 h h0 满足| Tn1 | Chp1
| en1 || en | hL | en | Ch p1 | en |
其中 1 hL, Ch p1
en O(hp )
即取
K * 1K1 2 K2 yi1 yi h(1K1 2 K2 )
其中1 和 2 是待定常数。若取 K1 f ( xi , yi ) ，则
问题在于如何确定 xi p 处的斜率 K2 和常数 1 和 2 。
仿照改进的欧拉方法，用欧拉方法预测 y( xi p ) 的值，
yi p yi phK1
k2
f ( xn
h 2
,
yn
h 2 k1)
hh k3 f ( xn 2 , yn 2 k2 )
k4 f ( xn h, yn hk3 )
例2：用经典的龙格-库塔方法
求解下列初值问题 h 0.1
dy 。 dx
y
2x y
x (0,1)
解：经典的四阶龙格-库塔公式： y(0) 1
E(h) 1 h
绝对稳定域： 1 h 1
当 R 时， 1 h 1 2 h 0
绝对稳定区间：(2, 0)
❖经典的R-K公式：yn1
yn
h 6
(k1
2k2
2k3
k4 )
k1 f ( xn , yn ) yn
k2
f
(
yn1 yn hf ( xn , yn )
n1 yn1 yn1
n1 n h[ f ( xn , yn ) f ( xn , yn )] [1 hf y ( xn , )]n
龙格-库塔方法---数学课的讲解共34页
16、人民应该为法律而战斗，就像为 了城墙 而战斗 一样。 ——赫 拉克利 特 17、人类对于不公正的行为加以指责 ，并非 因为他 们愿意 做出这 种行为 ，而是 惟恐自 己会成 为这种 行为的 牺牲者 。—— 柏拉图 18、制定法律法令，就是为了不让强 者做什 么事都 横行霸 道。— —奥维 德 19、法律是社会的习惯和思想的结晶 。—— 托·伍·威尔逊 20、人们嘴上挂着的法律，其真实含 义是财 富。— —爱献 生
21、要知道对好事的称颂过于夸大，也会招来人们的反感轻蔑和嫉妒。——培根 22、业精于勤，荒于嬉；行成于思，毁于随。——韩愈
23、一切节省，归根到底都归结为时间的节省。——马克思 24、意志命运往往背道而驰，决心到最后会全部推倒。——莎士比亚
25、学
第二节_龙格_库塔方法
yn1 yn hf ( xn , yn )
n1 yn1 yn1
n1 n h[ f ( xn , yn ) f ( xn , yn )] [1 hf y ( xn , )]n
如果 |1 hf y | 1，则误差是不增的，故可认为是稳定的
例如：对于初值问题
y y
y(
x0
)
a
精确解为
k2
f ( xn
h 2
,
yn
h 2 k1)
hh k3 f ( xn 2 , yn 2 k2 )
k4 f ( xn h, yn hk3 )
例2：用经典的龙格-库塔方法
求解下列初值问题 h 0.1
dy 。 dx
y
2x y
x
(0,1)
解：经典的四阶龙格-库塔公式： y(0) 1
h
y
ae x x0
y y
而实际求解的初值问题为
y(
x0
)
a
a
精确解为 y (a a)e x x0在 xn 处的误差为ae xn x0
可见误差随着 xn的增加呈指数函数增长
y y
如果初值问题为
y( x0 ) a
精确解为 y ae x0 x
y y
实际求解的初值问题为
y( x0 )
2C( h)5 2
y( xn1)
y(h/ 2) n1
y( xn1)
y(h) n1
1 16
记
|
y(h/ 2) n1
y(h) n1
|
16( y( xn1)
y(h/ 2) n1
)
y( xn1)
y(h) n1
y( xn1)
y(h/ 2) n1
龙格库塔法
c1
c2
1
0,
1 2
c2
0
即常数c1, c2 , 满足条件
c1 c2 1
c2
1 2
方程组有三个未知数，但只有两个方程，因此可得到
局部截断误差为O(h3 )的计算公式.
如果取c1
c2
1 ,
2
1,递推公式为
y0
k1 f (ti , yi )
k2 f (ti h, yi hk1)
yi
代入上式, 得
yi1 yi h(c1 f c2 f ) c2h( ft ffy ) O(h2 )
在局部截断误差的前提假设yi y(ti )下,得
y(ti1)
yi1
h(c1
c2
1)
f
h2(1 2
c2 )( ft
ffy ) O(h3)
要使局部截断误差y(ti1) yi1 O(h3 ),当且仅当
§9-3 龙格—库塔法
一、高阶泰勒法
假设初值问题
dy f (t, y) a t b dt
(1)
y(a)
的解y(t)及f (t, y)足够光滑.
将y(ti1)在ti处作n阶泰勒展开 , 得
y(ti1)
y(ti )
y(ti )h
y(ti ) h2 2!
y(n) (ti ) hn n!
y(n1) (i ) hn1
f
(ti
1 2
h,
yi
1 2 hk1)
(10)
yi1 yi hk2 )
公式(8)、(9)、(10)三式是三种常见的二阶龙格—库塔公式
局部截断误差为 O(h3).
三、三、四阶龙格—库塔法
三阶龙格—库塔法
龙格-库塔(Runge-Kutta)方法
( )
dy dy y′ = = f, y′′ = f x + f y = f x + ffy = F dx dx F F dy F F y′′′ = + = +f x y dx x y F = f xx + f x f y + ffyx = f xx + f x f y + ffxy x F 2 2 f = f(fxy + f y f y + ffyy ) = ffxy + ffy + f f yy y
y ( x) = e
但
x2
∫
x
0
e dt
t2
∫
x
0
e dt 难以求积
t2
ODE数值解的基本思想和方法特点 数值解的基本思想和方法特点
1. 离散化 级数、 用Taylor级数、数值积分和差商逼近导数， 级数 数值积分和差商逼近导数， 将 ODE转化为离散的代数方程 称差分方程 。 转化为离散的代数方程(称差分方程 转化为离散的代数方程 称差分方程)。
(ha2 ) + (ha3 ) +
f 2 2! x
2 2 2
(ha f ) +
2
2
K3 f + ha3 f x + h (a3 b32 ) f + b32 K2 f y 2! 2! + h2a3 (a3 b32 ) f + b32 K2 f xy f xx + h (a3 b32 ) f + b32 K2
2
Euler法 后退 法 ym+1 = ym + hK2 + O(h2 ) K2 = f ( xm + h, ym+1 )
计算方法龙格库塔方法.PPT课件
其中i (xi , xi1)
即
y(xi1) y(xi ) hy(i )
2019/7/13
8
引入记号 y(xi1) y(xi ) K
K hy(i) hf i, y(i)
yi1 yi K
h2 2
y(xi ) O(h3)
yi hf (xi , yi )
h2 2
f x(xi , yi ) f (xi , yi ) f y (xi , yi ) O(h3 )
类似地，若取前P+1项作为y(xi+1)的近似值，便得到
பைடு நூலகம்
yi1
yi
hyi
h2 2!
y(xi1) y(xi ) hy(xi ) O(h2 )
y(xi ) hf (xi , y(xi )) O(h2 )
yi hf (xi , yi ) O(h2 )
1
2019/7/13
1
若取前三项，可得到截断误差为O(h3)的公式
y(xi1)
y(xi ) hy(xi )
R-K方法不是直接使用Taylor级数,而是利用它的思想
2019/7/13
3
Runge-Kutta 方法是一种高精度的单步法,简称R-K法 9.4.1 龙格-库塔(R-K)法的基本思想
Euler公式可改写成
Kyi
1
hf
yi (xi ,
K yi )
则yi+1的表达式与y(xi+1)的Taylor展开式的前两项 完全相同,即局部截断误差为O(h2)。
龙格库塔法
第九章 常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。
然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。
在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。
这些方法可以给出解的近似表达式，通常称为近似解析方法。
还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。
利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题的精确解记为y (x )。
数值方法的基本思想是：在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。
这些h i 可以不相等，但一般取成相等的，这时nab h -=。
在这些节点上采用离散化方法，（通常用数值积分、微分。
泰勒展开等）将上述初值问题化成关于离散变量的相应问题。
把这个相应问题的解y n 作为y (x n )的近似值。
这样求得的y n 就是上述初值问题在节点x n 上的数值解。
一般说来，不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法1.欧拉法 1．对常微分方程初始问题(9.2))((9.1)),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第五章 常微分方程的数值解法
5.2 龙格-库塔法
一、教学目标及基本要求
通过对本节课的学习，使学生掌握常微分方程、常微分方程方程组的数值解法。
二、教学内容及学时分配
本节课主要介绍常微分方程的数值解法。
具体内容如下： 讲授内容：龙格库塔方法、收敛性与稳定性
三、教学重点难点
1．教学重点：龙格库塔方法、收敛性与稳定性。
2. 教学难点：收敛性与稳定性。
四、教学中应注意的问题
多媒体课堂教学为主。
适当提问，加深学生对概念的理解。
五、正文
龙格－库塔方法
引言 龙格－库塔方法的基本思想
'11()()
()()()(,())n n n n y x y x y y x y x hf y h
ξξξ++-=⇒=+
令*(,())K f y ξξ=称为区间1[,]n n x x +上的平均斜率，只要对*K 提供一种算法，即可推导出一种计算公式。
欧拉公式只是取点n x 的斜率作为区间1[,]n n x x +的平均斜率*K ，精度自然很低。
考察改进的欧拉公式
11211
()22
n n y y h k k +=++1(,)n n k f x y =21(,)n n k f x h y hk =++
它利用n x ，1n x +两点的斜率取算术平均，1n x +处斜率通过已知信息n y 用欧拉
公式预报得到。
可以考虑设法在1[,]n n x x +多取几个点的斜率值，将它们加权平均作为区间
1[,]n n x x +的平均斜率*K 。
这就是龙格－库塔方法的基本思想。
1、二阶龙格－库塔方法
考察1[,]n n x x +内一点,01n p n x x hp p +=+<≤，希望通过,n n p x x +两个点的斜率
12,K K 加权平均得到*K ，即令112((1))n n y y h K K λλ+=+-+
取1(,)n n K f x y =，如何预报n p x +处斜率2K ？
仿照改进的欧拉公式，先用欧拉公式预测()n p y x +的值n p y +：
1n p n y y phK +=+
然后用n p y +计算2(,)n p n p K f x y ++=，从而得
112121[(1)](,)(,)
n n n n n p n y y h K K K f x y K f x y phK λλ++=+-+==+
适当选取,p λ，使上述公式具有较高得精度。
假定()n n y y x =，分别将12
,K K 泰勒展开：
'1(,)()n n n K f x y y x ==
212
'
''
2
(,)
(,)[(,)(,)(,)]()()()()
n p n n n x n n n n y n n n n K f x y pkK f x y ph f x y f x y f x y O h y x phy x O h +=+=+++=++
代入得
'2''31()()()()n n n n y y x hy x ph y x O h λ+=+++
按泰勒展开2'
''
31()()()()2
n n n n h y y x hy x y x O h +=+++ 比较得，只要1
2
p λ=
，公式截断误差为3()O h 特别，当1,1/2p λ==，就是改进的欧拉公式， 改取1,1/2p λ==，
12n n y y hk +=+，1(,)n n k f x y =，21(,)22
n n h h
k f x y k =++
称为变形的欧拉公式，亦称中点公式。
2、三阶龙格－库塔方法
为提高精度，设除n p x +外再考察一点,01n q n x x hq p q +=+<≤≤，用三个点
,,n n p n q x x x ++的斜率123,,K K K 加权平均得出*123(1)K K K K λμλμ=--++
1123((1))n n y y h K K K λμλμ+=+--++
为预测n q x +的斜率值3K ，将12,K K 加权平均得出[,]n n q x x +上的平均斜率，从而得
12((1))n q n y y qh K K αα+=+-+
然后用n q y +计算3(,)n q n q K f x y ++= 设计的计算公式具有如下形式
1123121312[(1)](,)(,)
(,((1)))
n n n n n p n n q n y y h K K K K f x y K f x y phK K f x y qh K K λμλμαα+++=+--++==+=+-+
运用泰勒展开选择参数,,,,p q λμα 下列库塔公式是其中一种
112312112
312[4]
6
(,)(,)
2(,(2))
n n n n n n n q n h
y y K K K K f x y h
K f x y K K f x y h K K +++=+++==+=+-+
3、四阶龙格－库塔方法
继续上述过程，得出四阶方法，下列是经典格式
11234121123122
413[22]
6
(,)(,)2(,)2
(,)
n n n n n n n n n n h
y y K K K K K f x y h K f x y K h K f x
y K K f x y hK ++
+
+=++++==+=+=+
流程图P105 例P105
1211121312222413332(,)2()
2(,)()22()
2
2()
2(,)()22()
22()
(,)()()
n n n n n
n n n n n n n n n n n n n n n x K f x y y y h x h h K f x y K y K h y K h x h h K f x y K y K h y K x h K f x y hK y hK y hK +++==-
+=+=+-
++=+=+-
++=+=+-
+ 4、变步长龙格－库塔方法
步长越小，截断误差越小，但步长太小，计算量大，舍入误差也大，需要数值解法选择步长。
步长选择考虑的问题： 1）如何衡量和检验计算结果的精度； 2）如何依据判定的精度来处理步长。
以经典龙格－库塔方法为例，先以步长h 求得一个近似值记为1h
n y +，经典龙格－库塔方法的截断误差为5()O h ，故有5
11
()()h n n y x y CO h ++-≈ 步长折半/2
511
()2(())2
h n n h y x y CO ++-≈ /2
/2/2/211
11111111()11()()()()1615
h h h h h h n n n n n n n n h
n n y x y y x y y y y y y x y δ++++++++++-≈⇒-≈-⇒=-- 区分两种情况处理步长：
1）δε>，继续将步长折半，直至δε<为止，取步长折半后的新值作为结果；
2）δε<，步长加倍直至δε>，取步长加倍前的老值作结果。
5、介绍： 亚当姆斯方法
单步法的特点是计算1n y +时，只用到n y 的值，此时，011,,......,,n n y y y y -的值均已算出，如果在计算1n y +时，除用n y 的值外，同时还用到11,...,n n k y y --+的值，则称此方法为多步法。
k 步亚当姆斯公式：
1
11k n n i n i i y y h f β-+-=-=+∑
当式中的10β-=时，称为亚当姆斯显式公式。
当10β-≠时，称为亚当姆斯隐式公式。
小结
本节课主要介绍了解常微分方程的龙格库塔法。
要求大家掌握低阶龙格库塔法。
作业：课后习题4-8.。