计算方法上机作业——龙格库塔法
龙格-库塔方法（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是微分方程的表达式。
通过多次迭代，可以逐渐逼近微分方程的解。
龙格库塔算法的优点在于其精确度较高，可以通过调整步长来控制计算的精度和效率。
此外，该算法具有较好的数值稳定性，可以有效处理非线性、刚性或高阶微分方程等复杂问题。
因此，在科学和工程计算中，龙格库塔算法被广泛应用于各种数值模拟和求解问题。
需要注意的是，龙格库塔算法并非万能的，对于一些特殊的问题，可能存在数值不稳定性或计算精度不够的情况。
此外，算法的步长选择也需要根据具体问题进行调整，过小的步长会增加计算量，而过大的步长可能导致精度下降。
因此，在使用龙格库塔算法时，需要根据具体问题的特点和要求来选择合适的步长和算法参数，以获得满意的计算结果。
总结起来，龙格库塔算法是一种常用的数值解微分方程的方法，具有较高的精度和稳定性。
通过离散化和迭代的方式，可以逐步逼近微分方程的解。
龙格-库塔(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 ))处的函数值.
数值 四阶龙格库塔
数值计算方法第五次上机实习报告班级: 姓名: 学号一. 实习内容：常微分方程的数值解法二．主要算法：四阶龙格—库塔公式为：[]()()()()112341213243/622,/2,/2*,/2,/2*,,n n n n n n n n n n y y h k k k k k f x y k f x h y h k k f x h y h k k f x h y hk +⎧⎫=++++⎪⎪=⎪⎪⎪⎪=++⎨⎬⎪⎪=++⎪⎪⎪⎪=++⎩⎭计算过程：先分别求出k1,k2,k3,k4，,然后再将其代入到1n y =中即可。
实例：用经典四阶龙格——库塔方法求初值解问题：()'2,01x y y y y ⎧⎫=-⎪⎪⎨⎬⎪⎪=⎩⎭在[0,1]上的数值解（取h=0.2） 三.上机过程：1.程序代码：public static void ModEler(double x0, double y0,double xn, int n) {double yp , yc , x = x0, y = y0, h = (xn - x0) / n;System.out.println("x[0] = "+QuSiWei(x)+", y[0] ="+QuSiWei(y));for(int i=1;i<=n;i++) {yp = y + h * f1(x,y);x = x0 + i * h;yc = y + h * f1(x,yp);y = (yp + yc) / 2.0;System.out.println("x["+i+"] = "+QuSiWei(x)+", y["+i+"] = "+QuSiWei(y));}}public static double f1(double x, double y) {return y - 2 * x / y;}2.输出结果：3.编译中出现的问题：开始时四阶龙格——库塔算法理解起来有些不深，算法出了些问题。
第三部分龙格-库塔方法
内江师范学院数学与信息科学学院 吴开腾 制作
于是有
其中
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 )
计算方法上机实习六
计算方法上机实习六
实习内容六 龙格-库塔方法的使用
1、方法介绍 给定一阶常微分方程的初值问题 y´=f (x,y)， y(0)=y0 四阶龙格-库塔方法为： yn+1=yn +(K1+2K2+2K3+K4) /6 其中： K1=hf (xn,yn) K2=hf (xn+h/2, yn + K1/2) K3=hf (xn+h/2, yn + K2/2) K4=hf (xn+h, yn + K3) 2、计算步骤 ① 输入 x0, x1,h,N； ② 使用经典龙格－库塔公式计算出 y1； ③ 输出 x1, y1，并使 x1 x0 , y1 y0 转到 ② 直至 n > N 结束。 3、流程
开始 输入 x0, y0,h, N 1n
x0 +h ;hk1/2) k2 f(x0+h/2,y0 +hk2/2) k3, f(x1,y0+hk3) k4 y0+h(k1+2k2++2k3+k4)/6 y1
输出 x1, y1
n+1 n
n=N ?
n
y 结束
x1 x0 y 1 y0
4、实例 取步长 h=0.2，用经典格式求解初值问题
y 2xy y(0) 1
5 实习要求及实习报告
0 x 1
要求按以上过程完成实习内容，将结果输出到文件中，完成实习报告。实习报告 包括： 实习内容，要求，基本原理，运行结果，结果分析和总结。
龙格-库塔方法-文档资料
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)。
计算方法上机作业——龙格库塔法
龙格库塔法（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₄)/6
y₁=ℎ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.返回结果：最终求得的函数值。
需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重
要影响。
一般来说，步长越小，计算结果越精确，但计算时间也越长。
龙格库塔法是一种常用的数值计算方法，能够有效地求解常微分方程，并且较为精确。
在应用领域中，龙格库塔法被广泛用于物理学、工程学、
计算机科学等领域，如控制系统的设计、电路仿真、气象预报等。
总结起来，龙格库塔法是一种求解常微分方程的数值方法，通过数值逼近的方式来计算微分方程的近似解。
它具有精度高、稳定性好等优点，在实际应用中有着广泛的应用。