龙格库塔法
龙格库塔法解微分方程组
龙格库塔法解微分方程组引言微分方程组是数学中经常遇到的问题，在物理、工程和自然科学中都有广泛的应用。
为了求解微分方程组，我们需要利用数值方法来逼近解析解。
本文将介绍一种常用的数值方法——龙格库塔法（Runge-Kutta method），并探讨如何利用该方法来解微分方程组。
龙格库塔法概述龙格库塔法是一种迭代数值方法，用于求解常微分方程的初值问题。
它的主要思想是将微分方程的解进行离散化，将其转化为一系列的逼近值。
龙格库塔法的基本步骤如下：1.确定步长h和迭代次数n。
2.初始化初始条件，并假设第一个逼近值为y(xi)。
3.依次计算每个逼近值，直到得到y(xi+n*h)为止。
在每个迭代步骤中，龙格库塔法根据前一步的逼近值来计算下一步的逼近值。
该方法具有高阶精度和较好的稳定性，在实际应用中广泛使用。
单一微分方程的龙格库塔法首先，我们来看如何利用龙格库塔法来解一阶常微分方程。
以方程dy/dx = f(x, y)为例，其中f(x, y)为给定的函数。
步骤一：确定步长和迭代次数选择合适的步长h和迭代次数n来进行数值计算。
步长h决定了离散化的精度，而迭代次数n决定了逼近解的数目。
步骤二：初始化条件并计算逼近值设初始条件为y(x0) = y0，其中x0为起始点，y0为起始点处的函数值。
我们先通过欧拉法计算出y(x0 + h)的逼近值，然后再通过该逼近值来计算下一个逼近值。
逼近值的计算公式如下：k1 = h * f(x0, y0)k2 = h * f(x0 + h/2, y0 + k1/2)k3 = h * f(x0 + h/2, y0 + k2/2)k4 = h * f(x0 + h, y0 + k3)y(x0 + h) = y0 + 1/6 * (k1 + 2k2 + 2k3 + k4)步骤三：重复步骤二直到得到y(xi+n*h)依次利用上一步计算出的逼近值来计算下一个逼近值，直到得到y(xi+n*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)法
龙格-库塔（Runge-Kutta）法 1.1 龙格-库塔(Runge-Kutta)法的基本思想
Euler公式可改写成
yi1 yi hK1 K1 f ( xi , yi )
则yi+1的表达式y(xi+1)与的Taylor展开式的前两项 完全相同,即局部截断误差为 O(h 2 ) 。
为了进一步提高精度，设除 xi p 外再增加一点
xiq xi qh ( p q 1)
并用三个点 xi ,xi p , xiq 的斜率k1,k2,k3加权平均
得出平均斜率k*的近似值，这时计算格式具有形式:
yi1 yi h(1 )k1 k2 k3
k1 f (xi , yi ) k2 f (xi ph, yi phk1 )
格式。
若取 1 0 ,则 2 法的计算公式为
1,
p
1 2
，此时二阶龙格-库塔
ky1i
1
f
yi hk2 ( xi , yi )
k
2
h
f
(
x
i
1
,
yi
2
2 k1 )
i 0,1,2, n 1
此计算公式称为变形的二阶龙格—库塔法。式中
x 1 i 2
为区间
xi , xi1
的中点。
1.3 三阶龙格-库塔法
拉法,将 xi p 视为 xi1，即可得
k2 f (xi ph, yi phk1 ) 对常微分方程初值问题(7.1)式的解 y=y(x),根据微 分中值定理，存在点 (xi , xi1 ) ，使得
也即
y(xi1 ) y(xi ) y( )( xi1 xi )
y( xi1 ) y( xi ) hK
龙格库塔法
一、高阶泰勒法
假设初值问题
龙格—库塔法 龙格 库塔法
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）是一种常用于求解常微分方程（Ordinary Differential Equation，ODE）的数值解法。
它是由德国数学家卡尔·龙格（Carl Runge）和马丁·威尔海姆·库塔（Martin Wilhelm Kutta）在20世纪初提出的。
龙格库塔法的基本思想是通过数值逼近来计算微分方程的近似解。
在讲解龙格库塔法之前，我们先来简单回顾一下ODE的一阶常微分方程的基本形式：y′(y)=y(y,y)其中，y(y,y)是已知函数。
龙格库塔法的核心是使用差分逼近计算函数的斜率。
假设我们要求解的方程为：y′(y)=y(y,y),y(y)=y₀所需计算的点为y₀,y₁,...,yy，对应的函数值为y₀,y₁,...,yy，其中y是步长的个数。
龙格库塔法通过递推关系式来计算估计值，并不断更新当前点的函数值。
接下来以龙格库塔法的经典四阶形式为例进行说明。
该方法的基本方程如下：yy+1=yy+(y₁+2y₂+2y₃+y₄)/6y₁=ℎy(yy,yy)y₂=ℎy(yy+ℎ/2,yy+y₁/2)y₃=ℎy(yy+ℎ/2,yy+y₂/2)y₄=ℎy(yy+ℎ,yy+y₃)其中y表示当前步骤，ℎ表示步长，yy表示当前点的函数值，y₁,y₂,y₃和y₄则表示对应的斜率。
使用龙格库塔法，我们可以通过不断递归计算来求得指定区间（例如[y,y]）上的函数值。
具体步骤如下：1.确定求解区间[y,y]和初始点(y₀,y₀)以及步长ℎ。
2.初始化：设置yy=y₀，yy=y₀。
3.对所有y=0,...,y−1：计算y₁,y₂,y₃和y₄，根据上述递推关系式。
根据递推关系式计算yy+1更新当前点的函数值，即yy+1=y(yy+1)。
更新当前点的y值，即yy+1=yy+ℎ。
4.返回结果：最终求得的函数值。
需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重要影响。
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方法不是通过求导数的方法构造近似公式, 而是通过计算不同点上的函数值, 并对这些函数值作 线性组合, 构造近似公式, 再把近似公式与解的泰勒 展开式进行比较, 使前面的若干项相同 , 从而使近似 公式达到一定的阶数.
龙格库塔方程
龙格库塔方程1.介绍龙格-库塔（RK）方法是求解常微分方程（ODE）最常见的数值方法之一。
对于大多数非线性ODE问题，解析解并不存在或难以获得，因此需要使用数值方法来近似计算解。
RK方法通过迭代逼近ODE的解来得到精确性可控、收敛性好、易实现的数值解。
RK方法的基本思想是将ODE中的一阶导数转化为一组计算步骤，以得到相邻时间点之间的函数值和一阶导数的近似值，然后将其结合起来得到一个更精确的解。
2.RK方法的推导RK方法的推导过程是基于欧拉方法的，欧拉方法是RK方法的一阶近似。
假设有ODE$\frac{dx}{dt}=f(x,t)$，欧拉方法的迭代公式为$$x_{n+1}=x_n+hf(x_n,t_n)$$其中$h$是时间步长，$t_n=n*h$。
这个公式的意思是，从$x_n$开始，用一阶导数$f(x_n,t_n)$来列出切线，然后沿着切线向前移动$h$个单位，得到$x_{n+1}$。
更高阶的RK方法则基于更精细的近似。
例如，经典的四阶RK方法（RK4）迭代公式为：\begin{align*}k_1&=f(x_n,t_n)\\k_2&=f(x_n+\frac{h}{2}k_1,t_n+\frac{h}{2})\\k_3&=f(x_n+\frac{h}{2}k_2,t_n+\frac{h}{2})\\k_4&=f(x_n+h k_3,t_n+h)\\x_{n+1}&=x_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\end{align*}其中，$k_1$是欧拉方法的一阶导数解，依次计算得到更高阶的导数近似值$k_2-k_4$。
3.RK方法的优势RK方法与其他数值方法相比具有众多优点。
首先，RK方法的精度可控。
通过增加迭代次数或者近似阶次，RK 方法可以获得任意高的精度。
这个特性非常适用于涉及长时间尺度和小尺度特征的问题，例如天气预报，需要同时精确地处理地球的自转和大气的扰动。
龙格库塔 高阶常微分方程组
龙格库塔高阶常微分方程组龙格库塔法（Runge-Kutta method）是一种数值解常微分方程的方法，它由卡尔·龙格（Carl Runge）和马丁·康托尔·库塔（Martin Kutta）分别独立发现。
这种方法常常被用来解决高阶常微分方程组，其优点在于精度高且适用范围广。
1. 龙格库塔法的基本原理我们来简要介绍一下龙格库塔法的基本原理。
对于一个一阶常微分方程y′=f(x,y)，我们可以通过欧拉法进行数值解，其迭代公式为y_(n+1) = y_n + hf(x_n, y_n)。
而龙格库塔法则通过多个步骤的迭代来提高精度。
常见的四阶龙格库塔法的迭代公式如下：k1 = hf(x_n, y_n)k2 = hf(x_n+1/2h, y_n+1/2k1)k3 = hf(x_n+1/2h, y_n+1/2k2)k4 = hf(x_n+h, y_n+k3)y_(n+1) = y_n + 1/6(k1 + 2k2 + 2k3 + k4)2. 高阶常微分方程组的数值解对于高阶常微分方程组，我们可以通过变量替换的方式将其转化为一组一阶微分方程，然后利用龙格库塔法进行数值解。
以二阶常微分方程为例，我们可以令y1 = y，y2 = y'，然后构造一组一阶微分方程：y1' = y2y2' = f(x,y1,y2)这样，我们就可以利用龙格库塔法对其进行数值解了。
3. 个人观点和理解在我看来，龙格库塔法是一种非常有效的数值解法，尤其适用于高阶常微分方程组。
通过多步迭代的方式，可以大大提高数值解的精度，从而在实际问题中得到更为可靠的结果。
当然，龙格库塔法也存在一定的局限性，比如对于一些特殊形式的方程可能不够有效，需要结合其他方法进行求解。
总结回顾通过本文的介绍，我们对龙格库塔法有了更深入的了解。
从基本原理到在高阶常微分方程组中的应用，我们了解了其在数值计算中的重要性。
我也分享了自己对这一方法的个人观点和理解。
龙格库塔法介绍
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
常用的经典四阶龙格 库塔方法:
第四讲龙格-库塔方法
第四讲龙格-库塔⽅法龙格-库塔⽅法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个，所应满⾜的⽅程为(3.2.10)这是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）就是一种离散化的数值方法，其原理可以概括为：通过相应的递推公式，将微分方程在离散时间点上进行逼近，从而得到近似的解。
具体来说，假设要求解如下形式的一阶常微分方程：$$ y'=f(t,y) $$其中，$f(t,y)$是已知的函数，$y(t)$是未知函数，并且已知初值$y(t_0)=y_0$。
为了离散化这个方程，我们可以采用以下的递推公式：$$ \begin{aligned} y_1 &=y_0 + h\varphi_1 \\ y_2 &=y_0 +h\varphi_2 \\ \cdots &=\cdots \\ y_n &=y_0 + h\varphi_n \end{aligned} $$其中，$h$是离散时间点的时间步长，$t_n=t_0+nh$，$\varphi_i$是与$t_i$有关的递推公式。
根据龙格库塔方法的不同级别，$\varphi_i$也有不同的形式。
二、龙格库塔方法的应用由于龙格库塔方法的较高精度和鲁棒性，以及易于实现等特点，它在各个领域都有着广泛的应用。
1. 数学领域在数学领域，龙格库塔方法可以用于求解常微分方程、偏微分方程、常微分方程组等等，特别是对于复杂的高阶微分方程，龙格库塔方法更是可以发挥其优势。
2. 物理学领域在物理学领域，各种微分方程是研究物理过程的基础。
龙格库塔方法在求解各种物理问题时也得到了广泛的应用，如天体力学、流体力学、电磁场问题等等。
3. 经济学领域在经济学领域，许多经济问题可以通过微分方程的形式进行建模，并采用龙格库塔方法进行数值求解。
龙格库塔法
2020/4/25
10
令 y(xi1) yi1 对应项的系数相等，得到
c1 c2 1 ,
a2c2
1 2
,
b21c2
1 2
这里有 4 个未知 数，3 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库塔格式。
2020/4/25
• 1.在运动流体的整个空间，可绘出一系列的流线，称为流 线簇。流线簇的疏密程度反映了该时刻流场中速度的不同。 2.当为非定常流时，流线的形状随时间改变：对于定常流， 流线的形状和位置不随时间而变化。 3.定常流的流线和迹线重合。 4.一般情况下，流线不能相交，不能折转，只能是一条光 滑曲线。
龙格库塔法
2020/4/25
5
引入记号
y(xi1) y(xi ) K
K hy(i) hf i, y(i)
yi1 yi K
K可以认为是y y(x)在区间[xi , xi1]上的平均斜率
y
只要使用适当的方法求 出y(x)在区
y y(x)
间[xi , xi1]上平均斜率的近似值 K
K
就可得到相应的Runge-Kutta方法
2!
3!
4!
2020/4/25
20
龙格—库塔方法的推导基于Taylor展开方法，因而它要求所求的解具有
较好的光滑性。如果解的光滑性差，那么，使用四阶龙格—库塔方法求得的数 值解，其精度可能反而不如改进的欧拉方法。在实际计算时，应当针对问题的 具体特点选择合适的算法。对于光滑性不太好的解，最好采用低阶算法而将步 长h 取小。
23 4 5 6
第7-3龙格-库塔方法
一、Taylor展开法
设 y f ( x, y)Βιβλιοθήκη y( x0 ) y0
(1)
在[a,b]上有解 y( x),将y( xn1 )在xn处泰勒展开
y( xn1 )
y( xn )
hy( xn )
h2 2!
y( xn )
h3 3!
y( xn )
截取有限项作为 y( xn1 ) 的近似值,有
四、三阶及四阶龙格－库塔公式
三阶龙格—库塔公式有：
yn1 k1 k2
f f
h yn 6 (k1 ( xn , yn )
1 ( xn 2 h,
4k2 k3 h
yn 2 k1 )
)
k3 f ( xn h, yn hk1 2hk2 )
yn1 k1 k2
f f
1 (h n! x
k )n y
f ( x0, y0 )
(n
1 (h 1)! x
k
)n1 y
f
( x0
h,
y0
k)
(0 1)
其中记号
(h
x
k
y
)
f
(
x0
,
y0
)＝hf
x
(
x0
,
y0
)
kf
y
(
x0
,
y0
)
(h
x
k
)2 y
f
(
x0 ,
y0
)
h2
f
xx (
x0 ,
y0
)
2hkf
f ( xn , yn ) h(c2 f x ( xn , yn ) a21k1 f y ( xn , yn )) O(h2 ) f ( xn , yn ) h(c2 f x ( xn , yn ) a21 f y ( xn , yn ) f ( xn , yn )) O(h2 )
龙格库塔法
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 Method）是一种常用的数值求解微分方程的方法，适用于解决轨道参数的计算问题。
轨道参数如卫星的位置、速度、加速度等，可以通过求解对应的微分方程来确定。
在实际的轨道计算中，常常需要求解二阶或更高阶的微分方程，而龙格库塔法是一种高阶的数值求解方法，具有较高的精度和稳定性。
它的基本思想是将微分方程在一段很小的时间间隔内进行逐步的近似计算，从而得到轨道参数的数值解。
龙格库塔法的一般步骤如下：1. 给定初始条件，即初始位置和速度。
2. 将时间间隔划分为若干个小的步长，通常称为时间步。
3. 对于每一个时间步，依次进行以下计算：- 计算当前时刻的加速度，通过已知的物理力学方程或其他相关公式得到。
- 根据当前时刻的位置、速度和加速度，利用龙格库塔法的公式进行数值计算，得到下一个时刻的位置和速度。
- 更新时间，进入下一个时间步。
4. 重复步骤3，直到达到所需的计算精度或时间点。
需要注意的是，龙格库塔法的精度与时间步的选择有关，一般来说，时间步越小，计算结果越精确，但计算量也会随之增加。
因此，在实际应用中，需要根据具体情况选择合适的时间步和计算精度。
龙格库塔法是一种广泛应用于轨道计算、天体力学等领域的数值方法，它在解决轨道参数问题上具有一定的优势。
但需要注意的是，龙格库塔法是一种近似计算方法，其结果仍然存在误差，因此在使用时需要进行适当的误差估计和控制。
总之，通过龙格库塔法，我们可以较为准确地求解轨道参数的微分方程，得到物体在不同时间点的位置和速度等信息。
这对于轨道设计、卫星运行轨迹的预测等应用具有重要意义。
第3讲(龙格-库塔方法)
易见，它与二阶泰勒级数方法仅相差 O( h3 )!
这一分析给我们提供了一个重要信息，那就是 我们所遇到的泰勒级数方法中求导数的困难是可以 克服的，改进的欧拉方法就没有用到导数，而是借 助于函数在某些点处的值 (复合函数的思想）。
又 y( x ) df ( x , y( x )) f x f y y f x f y f dx
故二阶泰勒级数方法为 h2 yi 1 yi h f ( xi , yi ) ( f x ( xi , yi ) f ( xi , yi ) f y ( xi , yi )) 2! 更高阶方法更复杂，主要是求导复杂！
yi 1
h2 hk ( k ) yi h y y yi i i 2! k!
这样的数值方法称为k 阶泰勒级数方法。
yi 1
h2 hk ( k ) yi h y y yi i i 2! k!
泰勒级数方法也是单步法，且其局部截断误差为
h2 hk ( k ) LTE y( xi 1 ) y( xi ) hy( xi ) y( xi ) y ( x i ) 2! k!
第二节 龙格-库塔方法
（Runge-Kutta）
根据局部截断误差与整体误差的关系可知， 局部截断误差的阶是衡量一个方法优劣的重要依据。 考虑用提高局部截断误差的阶来提高数值方法的 精度。 泰勒级数法 龙格―库塔方法
一、泰勒级数方法
d y f ( x, y ), x I 如果初值问题 d x 的精确解 y(x) 在 I y( x ) y 0 0
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
北京航空航天大学自动控制原理实验报告
目录
实验四控制系统数字仿真 (2)
一、实验目的 (2)
二、实验方法 (2)
1、四阶龙格-库塔法 (2)
2、控制系统数字仿真 (3)
三、实验内容 (3)
四、实验步骤 (4)
1、确定K值 (4)
2、仿真过程 (4)
3、用Simulink工具箱进行仿真 (7)
五、总结与分析 (9)
实验四 控制系统数字仿真
一、 实验目的
通过实验掌握四阶龙格库塔算法，并进行数字系统仿真实验，观察分析系统参数变换对实验结果及系统性能的影响。
二、 实验方法 1、四阶龙格-库塔法
若一阶微分方程如下：{ẏ(t )=f(t,y(t))
y (t )=y 0 则在t n+1(t n+1>t 0)处，
y(t n+1)的近似值为：
y n+1
=y n +h
6
(k 1+2k 2+2k 3+k 4)
式中：h=t n+1−t n k 1= f(t n ,y n )
k 2= f(t n +1
2h,y n +1
2hk 2)
k 4= f(t n +12
h,y n +12
hk 1) k 3= f((t n +1
2
h,y n +1
2
hk 3)
n=0,1,2,3......
如果微分方程是如下形式的向量微分方程：{
X (t )=F(t,x (t ),u(t))
X (t )=X 0 其中，X (t )为m 维向量，t,u(t)均为标量，则在t n+1处（t n+1>t n ）, X (t )的近似值为：
X n+1
=X n +h
6
(K 1+2K 2+2K 3+K 4)
式中：h=t n+1−t n K 1= F(t n ,X n ,u(t n ))
K 2= F(t n +1
2h,X n +1
2hK 1,u(t n ))
K 4= F(t n +12
h,X n +12
hK 2,u(t n )) K 3= F((t n +1
2
h,X +1
2
hK 3,u(t n ))
n=0,1,2,3,……
2、控制系统数字仿真
设系统的闭环传递函数为：φ(s )=
y(s)u(s)=c 1s n+1+c 2s n+2+⋯+c n+1+c n
s n +as n−1+⋯+a n−1+a n
引入中间变量V(s)则上式可化为：y(s)u(s)
=
y(s)v(s)
×
v(s)u(s)
令
v(s)u(s)
=
1
s n +as n−1+⋯+a
n−1s+a n
，
y(s)v(s)
=c 1s n+1+c 2s n+2+⋯+
c n+1s +c n
经过运算、整理可得如下向量形式：
{
ẋ(t )=AX (t )+bu(t)
y (t )=cX(t)x (0)=0
所以，可得F(t,X (t ),u (t ))=AX (t )+bu(t)
三、实验内容
已知系统结构如图1所示：
图1 系统结构图
若输入为单位阶跃函数，计算当超调量分别为非作5%，25%和50%时K的取值，（用主导极点方法估算），并根据确定的K值在计算机上进行数字仿真。
四、实验步骤
1、确定K值
用主导极点法估算K值，将三阶系统近似为二阶系统处理，已知二阶系统超调量σ%与阻尼比ζ的关系如下：
σ%=e−πζ/√1−ζ2×100%
故由此可求得超调量σ%与阻尼比ζ的对应值，列表如下：
由二阶系统计算公式可设主导极点S1,2=−ζωn±j√1−ζ2ωn，将主导极点S1,2代入高阶系统特征方程D(s)=S3+10S2+25S+k= 0中，再由幅值条件和相角条件可解得对应极点和K值，列表如下：
2、仿真过程
用Matlab编写龙格库塔程序，将计算得到的K值带入程序，得到超调量，同时将阶跃响应曲线输出。
程序如下：（替换K值为上述计算值即可）
A=[0 1 0;0 0 1;-K -25 -10];
b=[0 0 1];
c=[K 0 0];
X=zeros(3,1);
t=0:0.01:10;
n=length(t);
h=0.01;
for i=1:n
K1=A*X+b;
K2=A*(X+(h/2)*K1)+b;
K3=A*(X+(h/2)*K2)+b;
K4=A*(X+h*K3)+b;
X=X+(h/6)*(K1+2*K2+2*K3+K4);
y(i)=c*X;
end
plot(y);
s=1001;while y(s)>0.95&y(s)<1.05;s=s-1;end;
t=(s-1)*0.005;
max(y)-1
仿真波形如下：（依次为K=31.1、K=59.5、K=103）
龙格库塔仿真结果如下：
3、用Simulink工具箱进行仿真
将龙格库塔仿真波形与Simulink仿真波形进行比较，验证龙格库塔程序的正确性。
Simulink仿真模型如图：
在Fransfer Fcn模块中改变分子值即可实现不同K值下的仿真。
经两种波形对比，验证了龙格库塔仿真程序的正确性，下面以K=59.5为例展示仿真结果。
K=59.5时对比图：（上图为龙格库塔仿真结果，下图为Simulink 仿真结果）
五、总结与分析
由龙格库塔仿真波形与Simulink仿真波形的比较可知龙格库塔算法本身误差较小，龙格库塔仿真的误差主要由主导极点近似处理高阶系统带来的。
对比龙格库塔仿真结果与理论值，可知随系统增益K增加，主导极点近似法带来的误差增加。
本次实验难度较大，主要原因是对Matlab掌握不足，但在实验后见识到Matlab强大的矩阵处理能力和Simulink仿真工具仿真自控系统的便捷。
如果是用C语言编写龙格库塔仿真程序将十分繁琐，增加仿真难度。