数值分析9-3 龙格-库塔方法
龙格-库塔方法（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 。
龙格库塔算法
龙格库塔算法龙格库塔算法（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是微分方程的表达式。
通过多次迭代，可以逐渐逼近微分方程的解。
龙格库塔算法的优点在于其精确度较高，可以通过调整步长来控制计算的精度和效率。
此外，该算法具有较好的数值稳定性，可以有效处理非线性、刚性或高阶微分方程等复杂问题。
因此，在科学和工程计算中，龙格库塔算法被广泛应用于各种数值模拟和求解问题。
需要注意的是，龙格库塔算法并非万能的，对于一些特殊的问题，可能存在数值不稳定性或计算精度不够的情况。
此外，算法的步长选择也需要根据具体问题进行调整，过小的步长会增加计算量，而过大的步长可能导致精度下降。
因此，在使用龙格库塔算法时，需要根据具体问题的特点和要求来选择合适的步长和算法参数，以获得满意的计算结果。
总结起来，龙格库塔算法是一种常用的数值解微分方程的方法，具有较高的精度和稳定性。
通过离散化和迭代的方式，可以逐步逼近微分方程的解。
数值计算中的龙格库塔算法
数值计算中的龙格库塔算法龙格库塔算法，又称龙格-库塔算法，是一种数值计算方法，主要用于求解微分方程。
它的好处是通过迭代得到更加精确的数值解，对于很多科学和工程问题，如天体力学、化学反应动力学、电路分析等，都有广泛的应用。
一、初识龙格库塔算法最早提出龙格库塔算法的是瑞士数学家卡尔·龙格和德国数学家马丁·库塔，他们在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. 稳定可靠：在每一步均会进行收敛性检验，确保计算结果准确无误。
龙格库塔法解轨道参数
龙格库塔法解轨道参数引言龙格库塔法（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轴上的分量。
龙格-库塔法
龙格－库塔法（Runge-Kutta）数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
经典四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。
令初值问题表述如下。
则，对于该问题的RK4由如下方程给出：其中这样，下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。
该斜率是以下斜率的加权平均：k1是时间段开始时的斜率；k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值；k3也是中点的斜率，但是这次采用斜率k2决定y值；k4是时间段终点的斜率，其y值用k3决定。
当四个斜率取平均时，中点的斜率有更大的权值：RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。
注意上述公式对于标量或者向量函数(y可以是向量)都适用。
显式龙格库塔法显示龙格－库塔法是上述RK4法的一个推广。
它由下式给出其中(注意：上述方程在不同著述中由不同但却等价的定义)。
要给定一个特定的方法，必须提供整数s (阶段数)，以及系数aij (对于1 ≤ j < i ≤ s), bi (对于i = 1, 2, ..., s)和ci (对于i = 2, 3, ..., s)。
这些数据通常排列在一个助记工具中，称为龙格库塔表：c2 a21c3 a31 a32cs as1 as2 as,s ? 1b1 b2 bs ? 1 bs龙格库塔法是自洽的，如果如果要求方法有精度p则还有相应的条件，也就是要求舍入误差为O(hp+1)时的条件。
这些可以从舍入误差本身的定义中导出。
例如，一个2阶精度的2段方法要求b1 + b2 = 1, b2c2 = 1/2, 以及b2a21 = 1/2。
第三部分龙格-库塔方法
内江师范学院数学与信息科学学院 吴开腾 制作
于是有
其中
y ( xn +1 ) − y ( xn ) = y '(ξ ), ξ ∈ ( xn , xn +1 ) h y ( xn +1 ) = y ( xn ) + hf (ξ , y (ξ ))
k * = f (ξ , y (ξ )) 称作区间 [ xn , xn +1 ] 上的平均斜率。 上的平均斜率 平均斜率。 问题：计算近似值y ( xn +1 ) 的关键是如何选择算法确定平均斜率 k *
（15）
f ( xn +1 , yn + h ( − k1 + 2 k 2 ))
内江师范学院数学与信息科学学院 吴开腾 制作
注释1 可以用Taylor展示证明格式（14） 注释1：可以用Taylor展示证明格式（14）具有三阶精 展示证明格式
度，并且还可以用类似的方法得到四阶及其以上的更高 阶精度的Runge-Kutta格式 阶精度的Runge-Kutta格式。 Runge 格式。
内江师范学院数学与信息科学学院 吴开腾 制作
h yn + （ k1 + 2 k 2 + 2 k3 +k 4） 6 f ( xn , y n ) h f ( x 1 , yn + k1 ) n+ 2 2 h f ( x 1 , yn + k 2 ) n+ 2 2 f ( xn +1 , yn + hk3 ) （16）
四阶龙格－ 四阶龙格－库塔格式计算结果
xn
0.1 0.2 0.3 0.4 0.5
yn
欧拉格式计算结果 xn yn y ( xn )
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
该算法是构建在数学支持的基础之上的。
龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。
如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。
一阶常微分方程可以写作：y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'）Yn+1=Yn+h*f(Xn,Yn)另外根据微分中值定理，存在0<t<1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th))这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。
利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式：K1＝f(Xn,Yn);K2=f(Xn+h/2,Yn+(h/2)*K1);K3=f(Xn+h/2,Yn+(h/2)*K2);K4=f(Xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6);所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。
仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。
想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。
编写的定步长的龙格库塔计算函数：function [x,y]=runge_kutta1(ufunc,y0,h,a,b,Vg)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点,发酵罐体积（参数形式参考了ode45函数）n=floor((b-a)/h);%求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值，可以是向量，但是要注意维数for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii),Vg);k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2,Vg);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2,Vg);k4=ufunc(x(ii)+h,y(:,ii)+h*k3,Vg);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end调用的子函数以及其调用语句：function dy=test_fun(x,y)dy = zeros(3,1);%初始化列向量dy(1) = y(2) * y(3);dy(2) = -y(1) + y(3);dy(3) = -0.51 * y(1) * y(2);对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下：[T,F] = ode45(@test_fun,[0 15],[1 1 3]);subplot(121)plot(T,F)%Matlab自带的ode45函数效果title('ode45函数效果')[T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数，别忘记改变初始值的维数subplot(122)plot(T1,F1)%自编的龙格库塔函数效果title('自编的龙格库塔函数')运行结果如下：。
数值分析9-3 龙格-库塔方法
最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ：
y i 1 y i h ( K1 2K 2 2K 3 K 4 ) 6 K1 K2 K3 K4 f ( xi , yi )
h f ( xi h , y K1 ) i 2 2 h f ( xi h , y K2 ) i 2 2
一、再看Taylor方法
h y ( x h) y ( x) hy' ( x) y" ( x) 2 h3 h 4 ( 4) y ( x) y ( x) 6 24 一般可取公式为如下形式
2
y n 1
h h ( p) y n hy' n y"n yn 2! p!
yi 1 K1 K2 yi h [1 K 1 2 K 2 ] f ( x i , yi ) f ( xi ph, yi phK 1 )
d f ( x, y) dx 首先希望能确定系数 1、2、p，使得到的算法格式有 2阶 dy 精度，即在 yi y( xi ) 的前提假设下，使得 f x ( x, y) f y ( x, y) dx Ri y( xi 1 ) yi 1 O( h3 ) f x ( x, y) f y ( x, y) f ( x, y) y( x )
考察改进的欧拉法，可以将其改写为： 斜率 一定取K1 K2 的平均值吗？
y i 1 K1 K2
1 1 yi h K 1 K 2 2 2 f ( xi , yi ) f ( x i h, y i hK 1 )
步长一定是一个h 吗？
将改进欧拉法推广为：
龙格库塔积分算法
龙格库塔法龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。
龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。
其中4 阶龙格库塔法是最常用的一种方法。
因为它相当精确、稳定、容易编程。
在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。
如果需要较高的精度, 可采取减小步长的方法即可。
4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。
1、初值问题对于一阶常微分方程的初值问题根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。
2、离散化取步长h=(b-a)/n，将区间[a , b]分成n个子区间：a=<=b在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲线上至少有一点，它的斜率与整段曲线的平均斜率相同，得=y’() (0<<1)其中，=可以将上式改写成y()=y()+h*K (2.1)其中K为平均斜率，K=f()公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。
欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。
3、Euler法欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙格库塔法的基础。
首先，令、为y() 及y()的近似值，并且令平均斜率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1）4、改进的欧拉法此种方法是取、两点的斜率的平均值作为平均斜率K，即K= ,其中、均为y()以及y()的近似值，就得到改进后的欧拉公式（4.1）其中、分别为、两点的斜率值，即= ,=在上面的（4.1）式中，k2是未知的，采用一种叫预报法的方法来求解。
龙格库塔法
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
该算法是构建在数学支持的基础之上的。
对于一阶精度的拉格朗日中值定理有：
对于微分方程：y'=f(x,y)
y(i+1)=y(i)+h*K1
K1=f(xi,yi)
当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进拉格朗日中值定理：
y(i+1)=y(i)+[h*( K1+ K2)/2]
K1=f(xi,yi)
K2=f(x(i)+h,y(i)+h*K1)
依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。
经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：
y(i+1)=y(i)+h*( K1+ 2*K2 +2*K3+ K4)/6
K1=f(x(i),y(i))
K2=f(x(i)+h/2,y(i)+h*K1/2)
K3=f(x(i)+h/2,y(i)+h*K2/2)
K4=f(x(i)+h,y(i)+h*K3)
通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式。
第3章03龙格-库塔方法
这样有计算公式：
yn1 yn h[(1 )k1 k2 )
k1 f (xn , yn ), k2 f (xn ph, yn phk1)
y(h) n1
1 16
y( xn1)
(h)
y2 n1
1 15
(h)
y2 n1
y(h) n1
误差事后估计
可以通过检查步长折半前后两次计算结果的偏差
y y ( h ) 2
(h)
n1
n1
来判断所选取的步长是否合适.
1) 对于给定的精度ε,如果δ> ε,则反复将步长折半 进行计算,直到δ< ε为止,这时取步长折半后的 新值作为结果;
k4
k3 f ( xn 2h, yn 21hk1 22hk2 ) f ( xn 3h, yn 31hk1 32hk2 33hk3)
四阶龙格－库塔公式每步要四次计算函数值，具有四阶精 度，局部截断误差是O(h5) .
四阶龙格－库塔格式
下列经典格式是其中常用的一种:
yn1 k2
yk
1.000000 1.005000 1.019025 1.041218 1.070800 1.107076 1.149404 1.197210 1.249975 1.307228 1.368514
y( xk ) yk
0.0 1.6×10－4 2.9×10－4 4.0×10－4 4.8×10－4 5.5×10－4 5.9×10－4 6.2×10－4 6.5×10－4 6.6×10－4 6.6×10－4
调用 2 次函数。有 2 阶精度。
龙格库塔法
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).
三、三、四阶龙格—库塔法
三阶龙格—库塔法
龙格-库塔方法
8.2 龙格-库塔方法8.2.1 二阶龙格-库塔方法常微分方程初值问题：做在点的泰勒展开：这里。
取，就有(8.11) 截断可得到近似值的计算公式，即欧拉公式：若取，式（8.11）可写成：或（8.12）截断可得到近似值的计算公式：或上式为二阶方法，一般优于一阶的欧拉公式（8.2），但是在计算时，需要计算在点的值，因此，此法不可取。
龙格-库塔设想用在点和值的线性组合逼近式（8.12）的主体，即用（8.13）逼近得到数值公式：（8.14）或更一般地写成对式（8.13）在点泰勒展开得到：将上式与式（8.12）比较，知当满足时有最好的逼近效果，此时式（8.13）－式（8.14）。
这是4个未知数的3个方程，显然方程组有无数组解。
若取，则有二阶龙格-库塔公式，也称为改进欧拉公式：（8.15）若取，则得另一种形式的二阶龙格-库塔公式，也称中点公式：（8.16）从公式建立过程中可看到，二阶龙格-库塔公式的局部截断误差仍为，是二阶精度的计算公式。
类似地，可建立高阶的龙格-库塔公式，同时可知四阶龙格-库塔公式的局部截断误差为，是四阶精度的计算公式。
欧拉法是低精度的方法，适合于方程的解或其导数有间断的情况以及精度要求不高的情况，当解需要高精度时，必须用高阶的龙格-库塔等方法。
四阶龙格-库塔方法应用面较广，具有自动起步和便于改变步长的优点，但计算量比一般方法略大。
为了保证方法的收敛性，有时需要步长取得较小，因此，不适于解病态方程。
8.2.2 四阶龙格-库塔公式下面列出常用的三阶、四阶龙格-库塔计算公式。
三阶龙格-库塔公式（1）（8.17）（2）（8.18）(3)（8.19）四阶龙格-库塔公式（1）（8.20）（2）（8.21）例8.3用四阶龙格-库塔公式（8.20）解初值问题：解：取步长，计算公式为：计算结果列表8.3中。
表8.3 计算结果8.2.3 步长的自适应欧拉方法和龙格-库塔方法在计算时仅用到前一步的值，我们称这样的方法为单步法。
龙格-库塔(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 )
数值分析9-3(龙格-库塔方法)
总结词
除了Python和MATLAB，还有许多其他编 程语言可以用于实现龙格-库塔方法。
详细描述
例如C、Java和R等编程语言也提供了相应 的数值计算库或框架，可以实现龙格-库塔 方法。使用这些语言实现龙格-库塔方法需 要一定的编程基础和对相应语言的数值计算 库的了解。
龙格-库塔方法可以用于求解偏微分方程的数值解，通过将偏微分方程转化为常微分方程组，利用龙格 -库塔方法进行迭代求解，能够得到较为精确的结果。
积分方程的数值解
积分方程是描述函数与积分之间的关 系的数学模型，常见于物理、工程等 领域。
VS
龙格-库塔方法也可以用于求解积分 方程的数值解，通过将积分方程转化 为常微分方程组，利用龙格-库塔方 法进行迭代求解，能够得到较为精确 的结果。
重要性及应用领域
龙格-库塔方法是数值分析中非常重要的内容， 它为解决常微分方程提供了一种有效的数值方 法。
在科学、工程和经济学等领域中，许多问题都 可以转化为求解常微分方程的问题，因此龙格库塔方法具有广泛的应用价值。
例如，在物理学、化学、生物学、金融学等领 域中，龙格-库塔方法被广泛应用于模拟和预测 各种动态系统的行为。
数值分析9-3：龙格-库塔方法
目录
• 引言 • 龙格-库塔方法概述 • 龙格-库塔方法在数值分析中的应用 • 龙格-库塔方法的实现与编程 • 龙格-库塔方法的改进与优化 • 结论与展望
01 引言
主题简介
龙格-库塔方法是一种用于求解常微 分方程的数值方法。
它通过构造一个离散化的时间序列来 逼近微分方程的解，并利用已知的离 散点来计算新的离散点，逐步逼近微 分方程的真实解。
02 龙格-库塔方法概述
定义与原理
93 龙格-库塔方法
r
K1 f ( xn , yn ),
i1 j 1
i 2,, r , （3.5）
K i f ( xn i h, yn h ij K j )
6
9.3.2
二阶显式R-K方法
对 r 2 的R-K方法，计算公式如下
yn 1 yn h(c1 K1 c2 K 2 ), K1 f ( xn , yn ), K f ( x h, y hK ). n 2 n 21 1 2
所耗费的计算量几乎相同.
龙格-库塔方法的推导基于泰勒展开方法，因而它要求 所求的解具有较好的光滑性质. 反之，如果解的光滑性差，那么，使用四阶龙格-库塔 方法求得的数值解，其精度可能反而不如改进的欧拉方法.
19
9.3.4
变步长的龙格-库塔方法
单从每一步看，步长越小，截断误差就越小，但随着 步长的缩小，在一定求解范围内所要完成的步数就增加了.
（3.6）
这里 c1 , c2 , 2 , 21均为待定常数. 希望适当选取这些系数，使公式阶数 p尽量高. 根据局部截断误差的定义，（3.6）的局部截断误差为
Tn1 y( xn1 ) y( xn )
（3.7）
7
h[c1 f ( xn , yn ) c2 f ( xn 2 h, yn 21hf n )],
这是8个未知数6个方程的方程组，解也不是惟一的. 所以 这是一簇公式.
满足条件（3.12）的公式（3.11）统称为三阶R-K公式.
一个常见的公式为
h yn 1 yn 6 ( K1 4 K 2 K 3 ), K1 f ( xn , yn ), h h K 2 f ( xn , yn K1 ), 2 2 K 3 f ( xn h, yn hK1 2hK 2 ).
龙格-库塔方法
h 按照数学学 习一般方法， yn 1 yn （k1 2k 2 2k3 +k 4） 6 继续上面的 做法，经过 k1 f ( xn , yn ) 较复杂的数 h 学推理和计 k （16） 2 f ( xn 1 , yn k1 ) 算，可以导 2 2 出四阶龙格 h －库塔格式， k f ( x 1 , yn k 2 ) 这里有一个 3 n 2 2 较有用的四 k f ( x , y hk ) 阶格式。 n 1 n 3 4
？
3、欧拉格式
yn 1 yn hf ( xn , yn ),
取点 xn 的斜率 k
*
(3)
n 0,1, 2,
f ( xn , yn )作为区间 [ xn , xn1 ] 上的平均斜率，精度低。
4、改进的欧拉格式
于是：
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn hf ( xn , yn ))] 2 h yn 1 yn [k1 k2 ] 2
龙格－库塔方法
考虑一维经典初值问题
dy f ( x , y ) , y( x0 ) y0 , x [a, b] dx
(1)
龙格－库塔方法
一、设计思想
加权平均斜率。
1、微分中值定理
如果 f (x) 在 [a, b]上连续，在 ( a, b) 内可导，则在( a, b) 内至 少存在一点 ，有 f (b) f (a) f ' ( ) ba
xnq xn qh, 0 p q 1
令
yn1 yn h[(1 )k1 k2 k3 ] (14)
k1 , k 2
三阶_龙格-库塔公式 详细推导过程
作者最近在学习《数值分析》这门课程时，发现多数教材只以二阶二段龙格库塔（R-K ）公式为例进行推导，对于三四阶R-K 公式只是简单一提，便给出结论。
作者查了一些资料，都未找到三阶R-K 公式的详细推导过程。
于是饶有兴致地尝试推导，功夫不负有心人，终于有所收获，特与大家分享。
作者以30岁的年龄、18岁的心态，终得此成果，很欣慰。
R-K 公式应用广泛，是一种高精度的单步法，读者在阅读本文之前，应该已经了解R-K 方法，并已熟悉二阶二段龙格库塔（R-K ）公式的推导过程。
主要参考资料：同济大学出版社：《高等数学》、《数值分析基础》。
三阶R-K 公式如下：k y iri in∑=++=11n yω 其中：f kn 1),(h hf y x nn ==),(1212k k b y a x nn hf ++= ｈ２ ),(23213133kk b k b y a x nn hf +++= ｈ将k 2按二元函数在处按Taylor 公式展开， （),(f n y x nn f =）)(221),(3222212122222212n121h ff hb ff b ha fh a ff hb fha fk b y a x O n n f yynxyxxy x nn ++++++=++）（ ｈ２然后k 3作类似的处理(注意：将其中的k 2用上式代替，并注意略去导致最终展开式中的O(h 4)项)：)()2)()222122),(322323122232231322331232232132y32313n232313)(h ff hb b ff h b b f b ha ff b ha fh a ff f h b b f f h a f h b f f hb fha fk b f h b y a x O n n n x n n n f yynyynxyxxy yyxnn ++++++++++++=++++（ ｈ将上述公式带入总公式即可，读者可自己尝试，此处不再赘述。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开(用二元函数泰勒展开)
y K 2 ] f ( x i , yi ) f ( xi ph, yi phK 1 )
d f ( x, y) dx 首先希望能确定系数 1、2、p，使得到的算法格式有 2阶 dy 精度，即在 yi y( xi ) 的前提假设下，使得 f x ( x, y) f y ( x, y) dx Ri y( xi 1 ) yi 1 O( h3 ) f x ( x, y) f y ( x, y) f ( x, y) y( x )
y( xi ) phy( xi ) O( h2 )
Step 2: 将 K2 代入第1式，得到
yi 1 y i h 1 y( x i ) 2[ y( x i ) phy( x i ) O ( h 2 )] y i (1 2 )h y( x i ) 2 ph2 y( x i ) O ( h 3 )
... ... K m f ( x i m h, y m1 hK 1 m 2 hK 2 ... m m 1 hK m 1 )
其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij ( i = 2, …, m; j = 1, …, i1 ) 均为待定系 数，确定这些系数的步骤与前面相似。
2
问题: 为获得更高的精度，应该如何进一步推广？
yi 1 y i h [ 1 K 1 2 K 2 ... m K m ] K 1 f ( x i , yi ) K 2 f ( x i 2 h, yi 21 hK 1 )
K 3 f ( x i 3 h, yi 31 hK 1 32 hK 2 )
考察改进的欧拉法，可以将其改写为： 斜率 一定取K1 K2 的平均值吗？
y i 1 K1 K2
1 1 yi h K 1 K 2 2 2 f ( xi , yi ) f ( x i h, y i hK 1 )
步长一定是一个h 吗？
将改进欧拉法推广为：
2
p
注：应用泰勒展式构造公式时，希望近 似值的泰勒展式与精确值的泰勒展式有尽可 能多的项相同，即：两者关于 h 的展式中的 系数尽可能多的项相同，这样局部截断误差 与 h 的更高次同阶，从而获得高精度的近似 公式。
二、Runge-Kutta 方法
建立高精度的单步递推格式。
单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜 率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所 能达到的最高精度为2阶。
由于龙格-库塔法的导出基于泰勒展开， 故精度主要受解函数的光滑性影响。对 于光滑性不太好的解，最好采用低阶算 法而将步长h 取小。
f ( xi h, yi hK 3 )
注： 龙格-库塔法的运算在于计算 Ki 的值，即 计算 f 的值。Butcher 于1965年给出了计 算量与可达到的最高精度阶数的关系：
每步须算Ki 的 2 3 4 5 6 7 n8 个数 可达到的最高 4 3 4 5 6 n2 O( h2 ) O( h ) O( h ) O( h ) O( h ) O( h ) O ( h ) 精度
Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较
yi 1 yi (1 2 )h y( xi ) 2 ph2 y( xi ) O( h3 )
h2 y( xi 1 ) y( xi ) hy( xi ) y( xi ) O( h3 ) 2
最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ：
y i 1 y i h ( K1 2K 2 2K 3 K 4 ) 6 K1 K2 K3 K4 f ( xi , yi )
h f ( xi h , y K1 ) i 2 2 h f ( xi h , y K2 ) i 2 2
3 要求 Ri y( xi 1 ) yi 1 O( h ) ，则必须有：
1 1 2 1 , 2 p 2
这里有 3 个未知 数， 2 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到，p 1, 1 2 1 就是改进的欧拉法。
第九章 常微分方程数值解 §2 龙格 - 库塔法 /* Runge-Kutta Method */
一、再看Taylor方法
二、Runge-Kutta 方法
dy f ( x , y ) x [a , b ] dx y ( a ) y0
复习
一元函数的Taylor展开
h y ( x h) y ( x) hy' ( x) y" ( x) 2 3 4 h h y ( x) y ( 4 ) ( x) 6 24
2
二元函数的Taylor展开
f ( x, y ) f ( x 0 , y 0 ) f x ( x0 , y 0 )( x x0 ) f y ( x0 , y 0 )( y y 0 ) f ( x , y )( x x ) 2 x 0 0 0 1 ( x0 , y 0 )( x x0 )( y y 0 ) 2 f xy 2 2 f y( x0 , y 0 )( y y 0 )
一、再看Taylor方法
h y ( x h) y ( x) hy' ( x) y" ( x) 2 h3 h 4 ( 4) y ( x) y ( x) 6 24 一般可取公式为如下形式
2
y n 1
h h ( p) y n hy' n y"n yn 2! p!