二阶龙格库塔方法
python龙格库塔法求解二阶方程
Python是一种高级编程语言，广泛应用于科学计算和工程领域。
而在科学计算中，求解二阶常微分方程是一个常见的问题。
龙格库塔法（Runge-Kutta method）是一种常用的数值求解方法，可以用来求解常微分方程的初值问题。
在Python中，我们可以利用其强大的数值计算库来实现龙格库塔法求解二阶常微分方程。
在本文中，我们将介绍如何使用Python中的龙格库塔法来求解二阶常微分方程。
文章将从以下几个方面展开讲解：1. 二阶常微分方程的基本概念2. 龙格库塔法的原理与公式推导3. Python中如何实现龙格库塔法求解二阶常微分方程4. 一个具体的求解例子及代码实现5. 总结与展望一、二阶常微分方程的基本概念二阶常微分方程是指具有形如y''(t) = f(t, y, y')(t)的形式，其中t是自变量，y是未知函数，f是已知函数。
求解二阶常微分方程的目标是找到一个满足方程的函数y(t)。
通常情况下，需要给出该方程的初值条件，即y(0)和y'(0)，以求得方程的具体解。
二、龙格库塔法的原理与公式推导龙格库塔法是一种数值求解常微分方程初值问题的方法，通过迭代计算来逼近方程的解。
它的基本思想是将微分方程的解按照一定的步长进行逼近，然后利用逼近值来计算下一个逼近值，直到达到所需的精度。
龙格库塔法的一般形式为：k1 = h * f(tn, yn)k2 = h * f(tn + 1/2 * h, yn + 1/2 * k1)k3 = h * f(tn + 1/2 * h, yn + 1/2 * k2)k4 = h * f(tn + h, yn + k3)yn+1 = yn + 1/6 * (k1 + 2*k2 + 2*k3 + k4)其中，h是步长，t是自变量，y是未知函数，f是已知函数，k1、k2、k3、k4是中间变量。
三、Python中如何实现龙格库塔法求解二阶常微分方程在Python中，可以利用其强大的数值计算库来实现龙格库塔法求解二阶常微分方程。
matlab龙格库塔方法求解二元二阶常微分方程组
matlab龙格库塔方法求解二元二阶常微分方程组文章标题：深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用摘要：在科学与工程领域，常常需要求解复杂的微分方程组，而matlab作为一种强大的数学工具，提供了许多求解微分方程组的方法。
本文将深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用，以便读者全面理解该方法并能灵活应用于实际问题中。
正文：一、介绍龙格库塔方法龙格-库塔法（Runge-Kutta methods）是一种数值求解常微分方程的方法，通过将微分方程的解进行离散化，将微分方程转化为差分方程，从而进行数值求解。
龙格库塔方法通过迭代计算，能够得到微分方程的数值解，广泛应用于科学计算和工程技术领域。
二、matlab中的龙格库塔方法在matlab中，龙格库塔方法通过ode45函数实现，该函数能够对一阶或高阶常微分方程进行数值求解。
用户可以通过设定初始条件、微分方程表达式，以及积分区间等参数，快速得到微分方程的数值解。
ode45函数采用自适应步长的方式进行求解，能够有效解决微分方程解的数值稳定性和精确度问题。
三、龙格库塔方法在求解二元二阶常微分方程组中的应用考虑如下形式的二元二阶常微分方程组：$$\begin{cases}y_1' = f_1(t, y_1, y_2) \\y_2' = f_2(t, y_1, y_2)\end{cases}$$其中，$y_1(t)$和$y_2(t)$是未知函数，$f_1(t, y_1, y_2)$和$f_2(t,y_1, y_2)$分别表示其对应的函数表达式。
通过matlab中的ode45函数，可以将该二元二阶常微分方程组转化为一阶常微分方程组的形式，然后利用龙格库塔方法进行数值求解。
设定初始条件$y_1(0) = y1_0, y_2(0) = y2_0$，对应的一阶方程组为：$$\begin{cases}u_1' = u_3 \\u_2' = u_4 \\u_3' = f_1(t, u_1, u_2) \\u_4' = f_2(t, u_1, u_2)\end{cases}$$其中，$u_1(t) = y_1(t), u_2(t) = y_2(t), u_3(t) = y_1'(t), u_4(t) =y_2'(t)$，通过ode45函数求解该一阶常微分方程组即可得到原二元二阶常微分方程组的数值解。
龙格库塔求解二阶常微分方程
龙格库塔求解二阶常微分方程一、前言在数学领域中，常微分方程是一种非常重要的数学工具。
它可以用来描述许多自然现象，如物理学、化学、生物学等。
而龙格库塔法是一种求解常微分方程的方法之一。
本文将介绍如何使用龙格库塔法求解二阶常微分方程。
二、什么是二阶常微分方程二阶常微分方程是指形如y''+p(t)y'+q(t)y=f(t)的方程，其中y表示未知函数，p(t)、q(t)和f(t)都是已知函数。
通常情况下，我们需要找到一个特定的y函数来满足这个方程。
三、什么是龙格库塔法龙格库塔法是一种数值求解常微分方程的方法。
它基于欧拉方法，但更准确和稳定。
欧拉方法使用线性逼近来估计未知函数值，而龙格库塔法使用高阶多项式逼近来更准确地估计未知函数值。
四、龙格库塔法的基本思想1. 将时间区间[0,T]平均分成N个小区间；2. 在每个小区间内进行递推计算，直到得到所有时间点上的y值；3. 每次递推计算都使用龙格库塔公式来估算y的值。
五、龙格库塔法的公式对于二阶常微分方程y''+p(t)y'+q(t)y=f(t)，我们可以使用以下公式来递推计算：1. k1=h*y'(t)l1=h*f(t)-p(t)*k1/2-q(t)*y(t);2. k2=h*(y'(t)+l1/2)l2=h*f(t+h/2)-p(t+h/2)*k2/2-q(t+h/2)*(y(t)+k1/2);3. k3=h*(y'(t)+l2/2)l3=h*f(t+h/2)-p(t+h/2)*k3/2-q(t+h/2)*(y(t)+k2/2);4. k4=h*(y'(t)+l3)l4=h*f(t+h)-p(t+h)*k4-q(t+h)*(y(t)+k3);其中，h表示时间步长，t表示当前时间点，f表示已知函数，p和q都是已知常数。
通过这些公式，我们可以得到下一个时间点上的y值。
六、龙格库塔法的实现为了更好地理解龙格库塔法，我们可以看一下具体的实现过程。
二阶微分方程龙格库塔法
二阶微分方程龙格库塔法
1.什么是二阶微分方程？
二阶微分方程是指二阶导数出现的方程，其通用形式为
y''+P(x)y'+Q(x)y=F(x)，其中P(x)、Q(x)、F(x)均为已知函数，y是所求函数。
2.为什么要用龙格库塔法？
求解二阶微分方程是数学、物理等领域中常见的问题，然而大多数情况下无法直接解题，所以需要使用数值方法。
龙格库塔法是一种数值方法，被广泛应用于求解微分方程，其优点是精度高、计算速度快。
3.龙格库塔法的原理
龙格库塔法是一种迭代方法，将微分方程看作初值问题，从初始点出发，采用一定的步长h，在每个点上用插值公式逼近y(x+h)，以此求得y(x+h)的近似值，一步步逼近所要求的精度。
4.龙格库塔法的步骤
（1）确定步长h和积分区间[a,b]。
（2）用初值y(x0)=y0，y'(x0)=y'0求出y(x0+h)的近似值。
（3）根据龙格库塔公式求得y(x0+2h)的近似值。
（4）对于连续求解的情况，重复以上步骤，直到求得所需的精度或者达到指定的终点。
5.龙格库塔公式
龙格库塔法的精度与所采用的公式有关，一般采用二阶或四阶的龙格库塔公式。
二阶龙格库塔公式为：
y0=y(x0)
k1=h*f(x0,y0)
k2=h*f(x0+h,y0+k1)
y(x0+2h)=y0+1/2(k1+k2)
其中，f(x,y)是待求函数。
6.总结
龙格库塔法是求解微分方程的一种常用数值方法，可以高精度、高效地解决二阶微分方程的问题。
该方法所需的计算量较小，易于编写程序实现，在实际应用中具有广泛的用途。
二阶龙格库塔方法
2012-2013（1）专业课程实践论文二阶Runge-Kutta方法董文峰，0818180123，R数学08-1班由改进的Euler 方法得到：()),(),(21121211K y h x hf K y x hf K K K y y i i i i i i ++==++=⎪⎩⎪⎨⎧+凡满足条件式有一簇形如上式的计算格式，这些格式统称为二阶龙格—库塔格式。
因此改进的欧拉格式是众多的二阶龙格—库塔法中的一种特殊格式。
若取1,0,2121212====c c b a ，就是另一种形式的二阶龙格 - 库塔公式。
⎪⎪⎩⎪⎪⎨⎧++==+=+)2,2(),(12121K hy h x f K y x f K hK y y n n n n n n （1）此计算公式称为变形的二阶龙格—库塔法。
二级龙格-库塔方法是显式单步式,每前进一步需要计算两个函数值。
由上面的讨论可知,适当选择四个参数y0,a,b,n ，可使每步计算两次函数值的二阶龙格-库塔方法达到二阶精度。
下面以式子（1）为依据利用VC++6.0编译程序进行问题的求解。
#include<stdlib.h>#include<stdio.h>/*n表示几等分，n+1表示他输出的个数*/int RungeKutta(double y0,double a,double b,int n,double *x,double *y,double (*function)(double,double)){double h=(b-a)/n,k1,k2;int i;x[0]=a;y[0]=y0;for(i=0;i<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);y[i+1]=y[i]+h*k2;}return 1;}double function(double x,double y){return y-2*x/y;}int main(){ int i;double x[6],y[6];printf("用二阶龙格-库塔方法\n");RungeKutta(1,0,1,5,x,y,function);for( i=0;i<6;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);return 1;}例1. 取h=0.2,用二阶Runge-Kutta 公式求解⎩⎨⎧=<<+='1)0(10,y x y x y解：1.将程序中的return y-2*x/y 改成 return x+y2.输入y0,a,b,n 分别为 1，0，1，53.结果为例2. 取h=0.2，用二阶Runge-Kutta 公式求解⎩⎨⎧=<<+='1)0(10),1/(3y x x y y解：1. 将程序中的return y-2*x/y 改成 return 3*y/(1+x)2.输入y0,a,b,n 分别为1,0,1,53.结果为。
龙格库塔法介绍
yn
hf
(xn, yn ))]，
(x, y,h) 1[ f (x, y) f (x h, y hf (x, y))],
2
|
( x,
y1,
h)
(x,
y2 ,
h)
|
[L
2
L 2
(1
hL)]
|
y1
y2
|,
L
L(1
h0L)，h 2
h0.
类似地,不难验证其他龙格 库塔方法的收敛性.
这里c1,c2,c3,2,3, 21, 31, 32均为待定参数.
Tn1 y(xn1) yn1 O(h4 )
(3.11)
c1 c2 c3 1
2
21
3 31 32
c22
c33
1 2
cc232223c2332
将步长折半，从xn用两步求xn1处的近似值，则有
y(xn1)
h
yn21
2c
h 2
5
.
从而
h
y ( xn 1) y ( xn 1)
yn21 ynh1
1， 16
得到事后估计式：
y ( xn 1)
h
yn21
1 15
(
h
yn21
ynh1).
通过检查步长折半前后计算结果的偏差，
y(x) (x, y(x),0) 0 p 1 单步法(4.1)收敛. 定义4 若单步法(4.1)增量函数(x, y,h)是否满足
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('自编的龙格库塔函数')运行结果如下：。
龙格-库塔方法-文档资料
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 阶精度，只需
第二节 龙格-库塔方法
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种三阶方法： 常见的 种三阶方法： 种三阶方法
matlab ppt_13_龙格-库塔法
这就是改进的欧拉公式，通常写成 h y = y + ( K1 + K2 ) i i+1 2 K1 = f ( x i , y i ) K2 = f (xi + h, yi + hK1) 由此可见，为了提高精度，可以多取几点的斜率值作加权平均当作平均斜 率kave，这就是龙格库塔法基本思想. 一般的显式龙格库塔方法可以写成
2. 二阶龙格――库塔法
在改进的欧拉公式中将xi+1替换为区间中的任意一点 xi+p = xi + ph (0 < p ≤ 1)
将它的斜率记为K2。利用欧拉公式得y (xi+p)的预报值 y i+p = yi + phK1 求得 K2 = f (xi+p, y i+p) 平均斜率则取 Keve ≈ λ1K1 + λ2K2 其中λ1, λ2是待定的系数。
(h)
y (xi+1) − yi+1 由此可得下列事后估计式
2 y (xi+1) − yi+1 ≈
(h)
≈
1 16
(h)
1 (h (h) 2) (yi+1 − yi+1) 15
这样，可以通过检查步长折半前后两次计算结果的偏差 ∆=
(h 2) |yi+1
− yi+1|
(j )
来判断所选取的步长是否合适。具体可分以下两种情况： 1. 对于给定的精度 ,如果∆ > ，就反复将步长折半计算，直至∆ < 为 (h 2) 止，这时取步长折半后的“新值”yi+1 作为结果。 2. 如果∆ < ，就反复将步长加倍，直至∆ > 为止，这时取步长加倍前 的“旧值”作为结果。 这种通过步长折半或加倍的计算的方法就叫变步长方法。表面上看，为了 选择步长，每一步的计算量似乎增加了，但从总体上考虑往往是合算的。
二阶runge-kutta法介数证明
二阶runge-kutta法介数证明
二阶Runge-Kutta法是一种常用的数值解常微分方程的方法。
该方法基于以下假设：
1. 求解的常微分方程为一阶方程：dy/dt = f(t, y)。
2. 采用离散化的方式对时间轴进行分割，将时间区间[tn, tn+1]划分为n个子区间，每个子区间的长度为h = (tn+1 - tn)/n。
3. 求解的方程在每个子区间上取常值，即在子区间[tn, tn+1]上可以近似为：y(t) ≈ ytn。
4. 采用迭代的方式，首先将当前y的值初始化为y0，然后通过迭代得到下个时间点的近似值。
5. 第二次迭代的近似值将作为下个时间点的初始值。
基于以上假设，二阶Runge-Kutta法可以描述为以下步骤：
Step 1: 给定初始条件 y0 和迭代时间点 t0，设定步长 h。
Step 2: 计算 k1 = hf(t0, y0)。
Step 3: 计算 k2 = hf(t0 + h/2, y0 + k1/2)。
Step 4: 计算下个时间点的近似值 y1 = y0 + k2。
Step 5: 更新时间点和迭代值，即 t1 = t0 + h，并取 y0 = y1。
Step 6: 重复步骤2到步骤5，直到达到所需的时间点。
通过以上的步骤，二阶Runge-Kutta法可以较好地逼近求解的常微分方程。
其中，k1为在当前时间点t0处的斜率值，k2则是在下个时间点t0 + h/2处的斜率值。
该方法利用了斜率在当前时间点和下个时间点之间的平均变化率来逼近方程的解。
二阶常系数微分方程的龙格库塔法MATLAB实现
二阶常系数微分方程的龙格库塔发MATLAB实现
作者：头铁的小甘
上一篇我们讲到一阶常系数微分方程利用MATLAB写四阶龙格库塔法进行求解。
无论如何，对于线性系统，系统都可以表示为输出的n阶导数和等于输入的m阶导数和进行表示，因此两边都是高阶微分，那么怎样把高阶微分方程变成一阶微分方程组合，再使用四阶龙格库塔法进行求解成为我们需要解决的问题。
因此其中一个最常用的思路是将高阶导数变成一阶微分的线性方程组表示。
下面我们解决一个二阶微分方程：
ð2U c ðt =−
RC
L
ðU c
ðt
−
1
L
U c+U i
该方程式对应物理电路的RLC串联，或者受迫振动的弹簧振子系统。
现在我们用线性方程组进行降阶。
令
x1=ðU c ðt
x2=U c
利用线性代数对变化后方程进行表示
(ẋ1
ẋ2
)=(
−
RC
L，−
1
L
1，0
)(
x1
x2
)+(
1
L
)U i
到该式已经完成利用一阶微分方程表示二阶微分方程，更高阶依次类推。
MATLAB代码实现
图1
图2
图3
图1表示时间轴定义的主程序，图2 为四阶龙格库塔算法，图3为二阶微分方程用一阶微
分方程组进行表示、
结果表示
下一篇：一阶空间域、时域常系数微分方程四阶龙格库塔法MATLAB求解。
(完整版)第二节龙格-库塔方法
因为单步法是 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
数值分析72龙格—库塔方法剖析
半计算,直至Δ<ε为止,这时取最终得到的 作为结yn果h21；
2.如果Δ<ε，我们将反复将步长作加倍计算，直 至Δ>ε为止，这时再将步长折半计算一次，就得到所
要的结果. 这种通过加倍或折半处理步长的方法称为变步长
方法.表面上看，为了选择步长，每一步的计算量增 加了，但总体考虑往往是合算的.
其中
y( xn1)
yn
hyn
h2 2
yn
h3 3!
ynO(h4 ),来自ynf ( xn, yn )
fn,
yn
d dx
f ( xn , y( xn ))
f x( xn , yn )
fn f y( xn , yn ),
yn
f xx 2 fn f xy
f
2 n
f
yy
f y[ f x
fn f y].
在选择步长时，需要考虑两个问题：
1. 怎样衡量和检验计算结果的精度？
2. 如何依据所获得的精度处理步长？
我们考察四阶R-K公式(3.13) ，从节点xn出发，
先以h为步长求出一个近似值，记为 yn(h)1，由于公式
的局部截断误差为O(h5)，故
y( xn1 )
y(h) n1
ch5 ,
(3.14)
然后将步长折半，即取为步长
yn1
yn
h 2 (K1
K2 ),
K1 f ( xn , yn ),
K
2
f ( xn1 , yn hK1 ).
如取a=1，则c1=0, c2=1, λ2=μ21=1/2. 得计算公式
yn1
yn
龙格库塔法
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).
三、三、四阶龙格—库塔法
三阶龙格—库塔法
数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程
Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。
由角动量定理我们知道εJ M =化简得到0sin 22=+θθl g dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为，这样比较容易解。
实际上这是一个解二阶常微分方程的问题。
θ在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上，使用能量法建立方程W T =h mg ∆=2J 21ω化简得到四阶龙格库塔方法使用四级四阶经典显式Rungkutta公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。
所以比欧拉稳定。
运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。
通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。
三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期legend('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%set(0,'RecursionLimit',500)h=0.01;a=0;b=25;x=a:h:b;RK_y(1)=0; %初值%RK_z(1)=0;初值for i=1:length(x)-1K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1 and L1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);K4=RK_z(i)+h*L3;% K4L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b+');xlabel('Variable x');ylabel('Variable y');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')%周期fprintf('周期等于%d',T*h)两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;%欧拉的初值Euler_z(1)=0;RK_y(1)=0;%龙格库塔初值RK_z(1)=0;for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 andL1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);% K2 and L2K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);% K3 and L3K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); K4 and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);%角度的峰值也就是πfprintf('角度峰值等于%d',y)fprintf('\n')%周期fprintf('周期等于%d',T*h)xlabel('时间');ylabel('角度');legend('欧拉','RK4');。
runge kutta 二阶例题
一、介绍Runge-Kutta 方法是求解常微分方程数值解的一种常用方法，广泛应用于工程、物理、经济等领域。
其中二阶 Runge-Kutta 方法是一种比较简单但准确度较高的求解方法。
本文将通过一个具体的例题来介绍二阶 Runge-Kutta 方法的应用过程，以帮助读者更好地理解和掌握这一数值计算方法。
二、例题描述考虑如下的初值问题：\(\frac{dy}{dt}=t-y, \quad y(0)=1\)我们希望利用二阶 Runge-Kutta 方法求解该初值问题的数值解。
三、方法步骤1. 离散化区间首先我们需要确定求解的区间范围，以及离散化的步长。
在本例中，我们以区间[0, 1]为例，取步长h=0.1进行离散化。
2. 迭代计算根据二阶 Runge-Kutta 方法的迭代计算公式，我们可以按照以下步骤进行数值解的计算:\(k_1=h(t_n, y_n)\)\(k_2=h(t_n+\frac{h}{2}, y_n+\frac{k_1}{2})\)\(y_{n+1}=y_n+k_2\)四、数值计算根据以上步骤，我们可以开始进行数值计算。
在本例中，我们可以依次按照步长h=0.1进行计算，并利用迭代公式计算出相应的数值解。
五、结果与分析通过计算得到的数值解可以与真实解进行比较，分析其误差大小。
我们也可以将计算得到的数值解进行绘图，以直观地观察数值解的走势。
通过对结果的分析，我们可以评价二阶 Runge-Kutta 方法在该初值问题下的适用性和准确性。
六、总结通过本文的介绍和实例计算，读者可以更加深入地理解二阶 Runge-Kutta 方法的应用过程，并掌握相应的计算步骤和技巧。
在实际应用中，可以根据具体的问题需求选择合适的数值计算方法并进行相应的调整和优化，以获得更加准确和可靠的计算结果。
希望本文能够为读者在理论和实践中对二阶 Runge-Kutta 方法的应用提供一定的帮助和指导。
二阶 Runge-Kutta 方法是常微分方程数值解的一种重要方法，在工程、物理、经济等领域有着广泛的应用。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2012-2013（1）专业课程实践论文二阶Runge-Kutta方法
董文峰，0818180123，R数学08-1班
一、算法理论
由改进的Euler 方法得到：
()),(),(2
1
12
1211K y h x hf K y x hf K K K y y i i i i i i ++==++
=
⎪⎩⎪⎨⎧+
凡满足条件式有一簇形如上式的计算格式，这些格式统称为二阶龙格—库塔格式。
因此改进的欧拉格式是众多的二阶龙格—库塔法中的一种特殊格式。
若取1,0,2
1
21212====c c b a ，就是另一种形式的二阶龙格 - 库塔公式。
⎪⎪⎩⎪
⎪⎨⎧
++==+=+)
2,2()
,(12121K h
y h x f K y x f K hK y y n n n n n n （1）
此计算公式称为变形的二阶龙格—库塔法。
二级龙格-库塔方法是显式单步式,每前进一步需要计算两个函数值。
由上面的讨论可知,适当选择四个参数y0,a,b,n ，可使每步计算两次函数值的二阶龙格-库塔方法达到二阶精度。
下面以式子（1）为依据利用VC++6.0编译程序进行问题的求解。
二、算法框图
三、算法程序
#include<stdlib.h>
#include<stdio.h>
/*n表示几等分，n+1表示他输出的个数*/
int RungeKutta(double y0,double a,double b,int n,double *x,double *y,double (*function)(double,double))
{
double h=(b-a)/n,k1,k2;
int i;
x[0]=a;
y[0]=y0;
for(i=0;i<n;i++)
{
x[i+1]=x[i]+h;
k1=function(x[i],y[i]);
k2=function(x[i]+h/2,y[i]+h*k1/2);
y[i+1]=y[i]+h*k2;
}
return 1;
}
double function(double x,double y)
{
return y-2*x/y;
}
int main()
{ int i;
double x[6],y[6];
printf("用二阶龙格-库塔方法\n"); RungeKutta(1,0,1,5,x,y,function);
for( i=0;i<6;i++)
printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); return 1;
}
四、算法实现
例1. 取h=0.2,用二阶Runge-Kutta 公式求解⎩⎨⎧=<<+='1)0(1
0,y x y x y
解：1.将程序中的return y-2*x/y 改成 return x+y
2.输入y0,a,b,n 分别为 1，0，1，5
3.结果为
例2. 取h=0.2，用二阶Runge-Kutta 公式求解⎩⎨⎧=<<+='1)0(1
0),1/(3y x x y y
解：1. 将程序中的return y-2*x/y 改成 return 3*y/(1+x)
2.输入y0,a,b,n 分别为1,0,1,5
3.结果为。