二阶龙格-库塔源程序1
龙格-库塔方法(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值(即由前面(i-1)个已算出的表示，故公式是显式的.例)，而ki如当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 。
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函数求解该一阶常微分方程组即可得到原二元二阶常微分方程组的数值解。
龙格库塔方程
龙格库塔方程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 方法可以获得任意高的精度。
这个特性非常适用于涉及长时间尺度和小尺度特征的问题，例如天气预报，需要同时精确地处理地球的自转和大气的扰动。
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('自编的龙格库塔函数')运行结果如下：。
计算方法 15 龙格库塔-常微分方程
以上公式就是三阶龙格-库塔公式
计算方法（2016/2017 第一学期）
西南科技大学
制造科学与工程学院
11
四阶龙格-库塔公式
四阶龙格-库塔公式：
h yn1 yn 6 ( k1 2k2 2k3 k4 ) k1 f ( xn , yn ) h hk1 ) k2 f ( xn , yn 2 2 h hk2 k3 f ( xn 2 , yn 2 ) k4 f ( xn h, yn hk3 )
以上公式就是四阶龙格-库塔公式
也称作经典龙格-库塔公式
计算方法（2016/2017 第一学期） 西南科技大学 制造科学与工程学院
12
高阶龙格-库塔公式
高阶龙格-库塔公式：
yn1 yn h( λ1k1 λ2 k2 λ3 k3 λm km ) k1 f ( xn , yn ) k2 f ( xn 2 h, yn 21hk1 ) k3 f ( xn 3 h, yn 31hk1 32 hk2 ) km f ( xn m h, yn m1hk1 m ( m 1) hkm 1 ) i (i 2,, m) 和 ij (i 2,, m; 其中 i (i 1,, m)，
14
龙格-库塔公式
由于龙格-库塔法的导出基于泰勒展开，故精度
主要受解函数的光滑性影响。
对于光滑性不好的解，最好采用低阶算法来求解，
而将步长 h 取小。
计算方法（2016/2017 第一学期）
西南科技大学
制造科学与工程学院
15
例题
例：利用四阶龙格-库塔方法求解微分方程的
C语言 微分方程组龙格-库塔法
/* By wende Isaac微分方程组龙格－库塔法*/#include<stdio.h>#include<math.h>double Ft(int j, double x ,double y[]);void rkutta4g(double x0, double y0[], double ym[], int M,double h);int main(){int N=5000,i,j,ipm=0,ip=100,M=50;double H=0.05;double x[100],y[100][10];double x0,y0[M],ym[M];x[0]=0.0;y[0][0]=1.0;for(i=1;i<M;i++){y[0][i]=0;}for(i=1;i<=N;i++){x[i]=x[i-1]+H;ipm++;x0=x[i-1];for(j=0;j<M;j++){y0[j]=y[i-1][j];}double ym[M];rkutta4g(x0,y0,ym,M,H);for(j=0;j<M;j++){y[i][j]=ym[j];// printf("ym[%d]= %f \n",j,ym[j]);}if(ipm>=ip){printf(" T: %.f(sec) ",x[i]);for (j=0;j<M-1;j++){printf(" C%d= %0.5f(m) ",j+1,y[i][j]);}printf(" C%d= %f(m) \n",M,y[i][j]);ipm=0;}}return 0;}void rkutta4g(double x0, double y0[], double ym[], int M,double h) {double s[4]={0,0.5,0.5,1},k[M][5],c[4]={1.0/6,1.0/3,1.0/3,1.0/6};int i,j,l;for(i=0;i<M;i++){k[i][0]=0;}for(j=0;j<M;j++){ym[j]=y0[j];// printf("ym[%d]=%f\n",j,ym[j]);}double xt,yt[5];for(j=0;j<4;j++){xt=x0+s[j]*h; //printf("xt=%f \n",xt);for(i=0;i<M;i++){yt[i]=y0[i]+s[j]*k[i][j];// printf("yt[i]=%f \n",yt[i]);}for(i=0;i<M;i++){k[i][j+1]=h*Ft(i,xt,yt);// printf("k[%d][%d] =%f \n",i,j+1,k[i][j+1]);}for(i=0;i<M;i++){ym[i]=ym[i]+c[j]*k[i][j+1];}}}double Ft(int j, double x ,double y[]){double k0,k1,k2,k3,ft;switch(j){case 0:ft=导数方程1;break;case 1:ft=导数方程2;break;case 2:ft=导数方程13;break;case 3:ft=导数方程4;break;case 4:ft=导数方程5;break;default :printf("error operation \n");}// printf("ft= %f \n", ft);return ft;}。
龙格-库塔方法-文档资料
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 阶精度，只需
微分方程常用的两种数值解法：欧拉方法与龙格—库塔法
资料范本本资料为word版本，可以直接编辑和打印，感谢您的下载微分方程常用的两种数值解法：欧拉方法与龙格—库塔法地点：__________________时间：__________________说明：本资料适用于约定双方经过谈判，协商而共同承认，共同遵守的责任与义务，仅供参考，文档可直接下载或修改，不需要的部分可直接删除，使用时请详细阅读内容四川师范大学本科毕业论文四川师范大学教务处二○一○年五月微分方程常用的两种数值解法：欧拉方法与龙格—库塔法学生姓名：xxx 指导教师：xx【内容摘要】微分方程是最有生命力的数学分支，在自然科学的许多领域中，都会遇到常微分方程的求解问题。
当前计算机的发展为常微分方程的应用及理论研究提供了非常有力的工具，利用计算机解微分方程主要使用数值方法，欧拉方法和龙格——库塔方法是求解微分方程最典型常用的数值方法。
本文详细研究了这两类数值计算方法的构造过程，分析了它们的优缺点，以及它们的收敛性，相容性，及稳定性。
讨论了步长的变化对数值方法的影响和系数不同的同阶龙格—库塔方法的差别。
通过编制C程序在计算机上实现这两类方法及对一些典型算例的结果分析比较，能更深切体会它们的功能，优缺点及适用场合，从而在实际应用中能对不同类型和不同要求的常微分方程会选取适当的求解方法。
关键词：显式单步法欧拉（Euler）方法龙格—库塔（Runge—Kutta）方法截断误差收敛性Two commonly used numerical solution of differential equations：Euler method and Runge - Kutta methodStudent Name: Xiong Shiying Tutor：Zhang Li【Abstract】The differential equation is the most vitality branch in mathematics. In many domains of natural science, we can meet the ordinary differential equation solution question. Currently, the development of computer has provided the extremely powerful tool for the ordinary differential equation application and the fundamental research, the computer solving differential equation mainly uses value method. The Euler method and the Runge—Kutta method are themost typical commonly value method to solve the differential equation. This article dissects the structure process of these two kinds of values commonly value method to solve the analyses their good and bad points, to their astringency, the compatibility, and the stabilityhas made the proof. At the same time, the article discuss the lengthof stride to the numerical method changing influence and thedifference of the coefficient different same step Runge—kutta method. Through establishing C program on the computer can realize these two kind of methods, Anglicizing some models of calculate example result can sincerely realize their function, the advantage and disadvantage points and the suitable situation, thus the suitable solution method can be selected to solve the different type and the different request ordinary differential equation in the practical application .Keywords: Explicit single-step process Euler method Runge—Kutta method truncation error convergence目录微分方程常用的两种数值解法：欧拉方法与龙格—库塔法前言常微分方程的形成与发展是和力学、天文学、物理学以及其他科学技术的发展密切相关的。