变步长四阶龙格库塔法原理
龙格库塔方法
• 可以通过检查步长折半前后两次计
算结果的偏差
=
(h)
y2 n1
y(h) n1
变步长方
来判断所选取的步长是否合法适
（1）对于给定的精度，如果 ,则反复将步长折半 进行计算直至 为止，这时取折半以后的“ 新值” 作为结果；
（2）如果 ，则反复将步长加倍，直至 为止， 这时取步长加倍前的“ 老值” 作为结果
四阶RungeKutta方法
这样，下一个值(yn+1)由现在的值(yn)加上时间间隔 (h)和一个估算的斜率的乘积 决定。该斜率是以下斜率的加权平均：
k1是时间段开始时的斜率； k2是时间段中点的斜率，通过欧拉法采用斜率 k1 来 决定 y在点 tn + h/2的值； k3也是中点的斜率，但是这次采用斜率 k2决定 y值； k4是时间段终点的斜率，其 y值用 k3 决定。 当四个斜率取平均时，中点的斜率有更大的权值：
y=y0;z=zeros(c,6); %z生成c行，6列的零矩阵存放结 果；
不断迭代运算： for x=a:h:b
if i1<=c k1=feval(yy,x,y); k2=feval(yy,x+h/2,y+(h*k1)/2); k3=feval(yy,x+h/2,y+(h*k2)/2); k4=feval(yy,x+h,y+h*k3); y=y+(h/6)*(k1+2*k2+2*k3+k4);
五、变步长Runge-Kutta方
法• 从每一步看，步长越小，截断误差
越小；但随着步长的缩小，在一定 求解范围内所要完成的步数就会增 加，步数的增加不但引起计算量的 增大，而且可能导致舍入误差的严 重积累，因此需要选择步长
四阶龙格—库塔法的原理及其应用
《四阶龙格—库塔法的原理及其应用》
龙格—库塔法(又称龙格库塔法)是由一系列有限的、独立的可能解组成的无穷序列,这些解中每个都与原来的数列相差一个常数。
它是20世纪30年代由匈牙利著名数学家龙格和库塔提出的,故得此名。
1.它的基本思想是:在n 阶方阵M 上定义一个函数,使得当n 趋于无穷时,它在m 中所表示的数值为M 的某种特征值,从而构造出一族具有某种特性的可计算函数f (x)= Mx+ C (其中C 为任意正整数)。
例如,若f (x)=(a-1) x+ C,则称之为(a-1) x 的龙格—库塔法。
2.它的应用很广泛,可以求解各类问题,且能将大量的未知数变换成少数几个已知数,因此它是近似计算的一种重要工具。
3.
它的优点主要有:(1)可以将多项式或不等式化成比较简单的形式;(2)对于同一问题可以用不同的方法来解决,并取得同样的结果;(3)适合处理高次多项式或者不等式,尤其适合处理多元函数的二次型。
4阶经典龙格库塔公式求解微分方程
4阶经典龙格库塔公式求解微分方程4阶龙格库塔法（也称为四阶Runge-Kutta方法）是一种常用于求解常微分方程的数值方法。
它是由德国数学家卡尔·龙格以及他的学生马丁·威尔海姆·库塔于1901年独立提出的。
以下是详细的介绍。
1.问题描述我们考虑一个典型的常微分方程初值问题：dy/dx = f(x, y)，y(x0) = y0其中，x0是给定的初始点，y(x)是我们要求解的函数，f(x,y)是已知的函数。
2.方法原理四阶龙格库塔方法的基本思想是通过使用全区间的函数信息来估计下一步的函数值。
在每个步骤中，我们计算四个不同的估计量，然后将它们组合起来以获取最终的解。
具体来说，我们首先计算在初始点x0上的斜率k1=f(x0,y0)。
然后我们计算在x0+h/2处的斜率k2=f(x0+h/2,y0+h/2*k1)，其中h是步长。
以此类推，我们计算k3和k4分别在x0+h/2和x0+h处的斜率。
然后，我们通过加权组合这些斜率来计算下一个点的函数值y1：y1=y0+(h/6)*(k1+2*k2+2*k3+k4)。
3.算法步骤以下是使用四阶龙格库塔法求解常微分方程的详细步骤：步骤1：给定初始条件 y(x0) = y0，选择步长 h 和积分终点 x_end。
步骤2：计算积分步数n：n = (x_end - x0)/h。
步骤3：设置变量x=x0和y=y0，并将步长设置为h。
步骤4：对每个步数i=1,2,...,n，执行以下计算：4.1：计算斜率k1=f(x,y)。
4.2：计算斜率k2=f(x+h/2,y+h/2*k1)。
4.3：计算斜率k3=f(x+h/2,y+h/2*k2)。
4.4：计算斜率k4=f(x+h,y+h*k3)。
4.5：计算下一个点的函数值y1=y+(h/6)*(k1+2*k2+2*k3+k4)。
4.6：将x更新为x+h，y更新为y1步骤5：重复步骤4，直到达到积分终点 x_end。
龙格-库塔方法基本原理
c1 c2 c3 1
a2c2
a3c3
1 2
,
a22c2 a32c31 3,b221c2
(b31
b32 )2 c3
1 3
a 2 b32 c 2
1 6
,
b21c2
(b31
b32 )c3
1 2
a2b21c2
a3c3 (b31
b32 )
1 3
1 b21b32c3 6
2021/5/25
6
同理，改进Euler公式可改写成
y
i
1
yi
1 2
K1
1 2
K2
K
1
hf
( xi ,
yi )
K
2
hf
( xi
h, yi
K1)
局部截断误差为O(h3)
上述两组公式在形式上共同点:都是用f(x,y)在某些 点上值的线性组合得出y(xi+1)的近似值yi+1, 且增加计 算的次数f(x,y)的次数,可提高截断误差的阶。如欧拉 法:每步计算一次f(x,y)的值,为一阶方法。改进欧拉法 需计算两次f(x,y)的值，为二阶方法。
17
若取 a2b2 塔公式。
11 2,
c1，就0,是c2另一1种形式的二阶龙格
-
库
yi1 yi K2
K1 hf(xi, yi)
K2
hf(xi
1h, 2
yi
12K1)
此计算公式称为变形的二阶龙格—库塔法。式中
xi
1为h区间
2
x的i , x中i1点。也称中点公式。
Q:为获得更高的精度，应该如何进一步推广？
matlab利用四阶runge-kutta算法求解原理
matlab利用四阶runge-kutta算法求解原理四阶Runge-Kutta（RK4）方法是一种常用的数值求解常微分方程（ODEs）的方法。
下面是RK4方法的基本原理，以及如何在MATLAB中实现：###基本原理：1.ODE表示：我们考虑形如dy/dx=f(x,y)的常微分方程，其中y是未知函数，f是给定的函数。
2.离散化：我们将x轴上的区间分成若干小步长h。
我们的目标是找到每一步上的y值。
3.四阶Runge-Kutta公式：这个方法的核心是通过四个中间步骤来逼近每一步的斜率，然后计算新的y值。
具体的步骤如下：-k1=h*f(x_n,y_n)-k2=h*f(x_n+h/2,y_n+k1/2)-k3=h*f(x_n+h/2,y_n+k2/2)-k4=h*f(x_n+h,y_n+k3)其中，x_n和y_n是当前步的x和y值，h是步长。
新的y值计算为：y_{n+1}=y_n+(k1+2*k2+2*k3+k4)/6###在MATLAB中的实现：在MATLAB中，你可以使用以下的代码来实现四阶Runge-Kutta算法：```matlabfunction[x,y]=runge_kutta_4th_order(f,x0,y0,h,x_end)x=x0:h:x_end;y=zeros(size(x));y(1)=y0;for i=1:(length(x)-1)k1=h*f(x(i),y(i));k2=h*f(x(i)+h/2,y(i)+k1/2);k3=h*f(x(i)+h/2,y(i)+k2/2);k4=h*f(x(i)+h,y(i)+k3);y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6;endend```这个函数接受一个ODE的右侧函数f，初始值x0和y0，步长h，以及求解的终点x_end。
返回的x和y包含了在给定区间内的解。
你可以调用这个函数并提供你自己的ODE右侧函数f。
第四讲龙格-库塔方法
龙格-库塔方法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 。
4阶runge-kutta原理
4阶Runge-Kutta方法是一种数值求解常微分方程的方法，它通过迭代的方式逐步逼近微分方程的解。
本文将从原理、推导以及应用等方面对4阶Runge-Kutta方法进行详细解读。
1. 原理4阶Runge-Kutta方法是数值分析中常用的数值解常微分方程的方法之一。
它的核心思想是利用哈密顿显式中点法求解微分方程。
该方法通过将微分方程的解离散化，然后通过计算每一步的斜率来逐步逼近方程的解，最终得到数值解。
2. 推导假设我们要求解如下的一阶常微分方程初值问题：$\frac{dy}{dx} = f(x, y)$$y(x_0) = y_0$其中$f(x, y)$是关于$x$和$y$的函数，$y_0$是初值，$x_0$是初始点。
现在我们希望通过4阶Runge-Kutta方法来求解上述方程。
我们将自变量$x$进行离散化，即将其分成$n$个小区间，每个小区间长度为$h$，即$x_i = x_0 + ih$，$i=0,1,2,...,n$。
然后我们利用下面的迭代公式来计算每一步的$y$的近似值：$k_1 = h f(x_i, y_i)$$k_2 = h f(x_i + \frac{h}{2}, y_i + \frac{k_1}{2})$$k_3 = h f(x_i + \frac{h}{2}, y_i + \frac{k_2}{2})$$k_4 = h f(x_i + h, y_i + k_3)$$y_{i+1} = y_i + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$式中，$k_1$、$k_2$、$k_3$、$k_4$分别表示斜率的四个近似值，$y_{i+1}$表示下一个点的近似值。
3. 应用4阶Runge-Kutta方法在实际工程问题中有着广泛的应用。
它不仅可以用来解决一阶常微分方程，还可以推广到高阶微分方程、常微分方程组以及偏微分方程等更复杂的问题。
由于该方法的高精度和稳定性，它也被广泛应用于科学计算领域，例如物理学、工程学、生物学和经济学等各个领域。
龙格-库塔方法-文档资料
c3a 2b32
c3a3 1
6
; 2
O (h4)
常见的2种三阶方法：
库塔三阶方法
h yn1yn6(k14k2k3)
k1
f(xn,yn);
k2
hh f(xn2,yn2k1)
k 3 f(x n h ,y n h k 1 2 h k 2 ) •5
四级方法：N = 4
局部截断误差 O ( h 5 )
可见误差随着 x n 的增加呈指数函数递减
当 f y 0 时，微分方程是不稳定的； 而 f y 0 时，微分方程是稳定的。
上面讨论的稳定性，与数值方法和方程中 f 有关
•21
实验方程： y y C ,R e () 0
D e f 3 对单步法 yn 1ynh(xn,yn,h )应用实验方程，
e n 1 e n h [ ( x n ,y ( x n ) , h ) ( x n ,y n , h ) ] T n 1
•15
因为单步法是 p 阶的：h0,0hh0满足|Tn1|Chp1
|e n 1| |e n| h L |e n| C h p 1|en |
其中 1hL,C hp1
•18
三、绝对稳定性 /*Absolute Stibility*/ 计算过程中产生的舍入误差对计算结果的影响
首先以Euler公式为例，来讨论一下舍入误差的传播:
yn1ynhf(xn,yn)
设实际计算得到的点 x n 的近似函数值为 yn yn n，
其中 y n 为精确值， n 为误差
yn1ynhf(xn,yn)
通过适当选取参数 1,2和 p 的值，使得公式具有 2阶精度！！
•3
由泰勒公式展开，要使公式具有 2 阶精度，只需
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
变步长的龙格变步长的龙格-库塔方法
单从每一步看，步长越小，截断误差就越小，但随着 步长的缩小，在一定求解范围内所要完成的步数就增加了. 步数的增加不但引起计算量的增大，而且可能导致舍 入误差的严重积累. 因此同积分的数值计算一样，微分方程的数值解法也 有个选择步长的问题. 在选择步长时，需要考虑两个问题： 1° 怎样衡量和检验计算结果的精度？
y(xn+1) y
h ( ) 2 n+1
h ≈ 2c , 2
5
（3.15）
比较（3.14）式和（3.15）式我们看到，步长折半后， 1 误差大约减少到 ， 16
3
即有
y(xn+1) y y(xn+1) y
h ( ) 2 n+1 (h) n+1
1 ≈ . 16
由此易得下列事 h为步长求出一个近似值 ynh) ， +1
2
由于公式的局部截断误差为 O(h5 )，故有
( y(xn+1) ynh) ≈ ch5 , +1
（3.14）
然后将步长折半，即取 h 为步长从 xn 跨两步到 xn+1， 2 5 h 再求得一个近似值 y( 2) ，每跨一步的截断误差是 c h ， n+1 2 因此有
1
2° 如何依据所获得的精度处理步长？ 考察经典的四阶龙格-库塔公式
h yn+1 = yn + 6 (K1 + 2K2 + 2K3 + K4 ), K1 = f (xn , yn ), h h K2 = f (xn + , yn + K1), （3.13） 2 2 K = f (x + h , y + h K ), n n 2 3 2 2 K4 = f (xn + h, yn + hK3 ).
h ( ) 2 n+1
作为结果；
2. 如果 < ε ，反复将步长加倍，直到 > ε为止， 这时再将步长折半一次，就得到所要的结果. 这种通过加倍或折半处理步长的方法称为变步长方法 变步长方法. 变步长方法 表面上看，为了选择步长，每一步的计算量增加了， 但总体考虑往往是合算的.
5
h ( ) 2 n+1 h ( ) 1 ( 2 ≈ [ yn+1 ynh) ]. +1 15
这样，可以通过检查步长，折半前后两次计算结果的偏差
= y
h ( ) 2 n+1 ( ynh) +1
来判定所选的步长是否合适. 具体地说，将区分以下两种情况处理：
4
1. 对于给定的精度 ε，如果 > ε ，反复将步长折半 进行计算，直至 < ε 为止. 这时取最终得到的 y