1、龙格-库塔方法及其matlab实现摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具达到求解目的。龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。MatLab软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。MatLab是英文MATrix LABoratory（矩阵实验室）的缩写。在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。关键词：龙格-库塔 matlab 微分方程1. 前言1.1：知识背景龙格库塔法（Runge-Ku
2、tta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔龙格和马丁威尔海姆库塔于1900年左右发明。通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库程序的“通俗易用”的接口，此即用F
3、ORTRAN编写的萌芽状态的MATLAB。 经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。从这时起，MATLAB的内核采用C语言编写，而且除原有的数值计算能力外，还新增了数据图视功能。 MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行的可靠性，使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND和SIMNON，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。在时间进入20世纪九十年代的时候，MATLAB已经成为国际控制界公
4、认的标准计算软件。到九十年代初期，在国际上30几个数学类科技应用软件中，MATLAB在数值计算方面独占鳌头，而Mathematica和Maple则分居符号计算软件的前两名。Mathcad因其提供计算、图形、文字处理的统一环境而深受中学生欢迎。1.2研究的意义精确求解数值微分方程，对龙格库塔的深入了解与正确运用，主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。利用matlab强大的数值计算功能，省去认为计算的过程，达到快速精确求解数值微分方程。在实际生活中可以利用龙格库塔方法和matlab的完美配合解决问题。1.3研究的方法对实例的研究对比，实现精度的要求，龙格
5、库塔是并不是一个固定的公式，所以只是对典型进行分析2. 龙格-库塔方法2.1龙格-库塔公式在一阶精度的的拉格朗日中值定理有：对于函数y=f（x，y） y=f(x,y)y(n+1)=y(n)+h*K1K1=f(xn, yn)这就是一阶龙格-库塔方法形如 y(n+1)=y(n)+h*i=1rcikik1 =f（xn，yn） ki=f（xn+hai ，yn+h*j=1i-1bijki） i=2r故二阶龙格-库塔公式y(n+1)=y(n)+h（c1k1+c2k2）k1= f（xn，yn） (2)k2= f（xn+ha2 ，yn+ha2 k1）将y（x）在xn处展成幂级数y（xn+1）=y(xn)+hy
6、(xn)+h22y (xn)+o（h3）y（x）= f（x，y（x）yx= fx（x，y（x）+ fy（x，y（x）f(x，y（x）)y（xn+1）=y(xn)+hf+h22(fx+fyf)+ o（h3） (3)将（2）式中的k2在（xn，yn）点展成幂级数k2= f（xn+ha2 ，yn+ha2 k1） =f+ha2fx+ ha2fyf+ o（h2）将k1，k2代入（2）式，得yn+1=yn +h（c1+c2）f+ha2c2（fx+fyf）+ o（h3） （4） 对比（3）（4），当y(xn)= yn时只有c1+c2=1，a2c2=12 （5）形如(2)存在常数满足（5）式，局部截断误差为o
7、（h3）的求解方法称为二阶龙格-库塔法。满足（5）式，若取c1=12，则得到c2=12，a2=1，则公式则恰为预估-校正法公式若取c1=0，则c2=1，a2=12， yn+1=yn+hk2 k1= f（xn，yn） (6) k2= f（xn+h2，yn+h2k） n=0,1N-1由（5）式，可知龙格-库塔法不是唯的三阶龙格-库塔法 yn+1=yn+h（c1k1+c2k2+ c3k3） k1= f（xn，yn） k2= f（xn+ha2，yn+ha2k1） （7） k3= f（xn+ha3，yn+hb31k1+hb32k2）若c1，c2， c3，a2，a3，b31, b32且满足b31+ b32
8、=a3，并使得局部截断误差为o（h4）。类似二阶龙格-库塔法推导的 c1+c2+ c3=1 a2c2+a3c3=12 a2b32c3=16 （8） a22c2+a32c3=13形如（7），常数满足（8），局部截断误差为o（h4）的求解方法称为三阶龙格-库塔法在（8）式中若取c1=16，c3=16，则得c2=23，a2=12，a3=1，b31=-1，b32=2代入（7）中得三阶龙格-库塔法公式 yn+1=yn+h6（k1+4k2+ k3） k1= f（xn，yn） k2= f（xn+h2，yn+h2k1） （9） k3= f（xn+ha3，yn-hk1+2hk2）四阶龙格库塔法的推导类似于三阶龙
9、格-库塔法，但相对复杂这里不再进行推导，公式如下 yn+1=yn+h6（k1+2k2+2 k3+k4） k1= f（xn，yn） k2= f（xn+h2，yn+h2k1） （10） k3= f（xn+h2，yn+h2k2） k4= f（xn+h，yn+h k3） n=0,1N-1这就是标准四阶龙格库塔公式2.1 对实例的研究 利用龙格-库塔法求解方程 y=8-3yy0=2的数值，其中h=0.2，计算y(0.4)的近似值。至少保留四位小数。解：f(x，y)83y，利用四阶龙格-库塔公式有 yn+1=yn+h6（k1+2k2+2 k3+k4） k1= f（xn，yn）=8-3yn k2= f（xn+h2，yn+h2k1）=5.6-2.1yn k3= f（xn+h2，yn+h2k2）=6.32-2.37yn k4= f（xn+h，yn+h k3）=4.2081.578yn n=0,1N-1 yn+1=yn1.20160.5494yn当x00，y02，y(0.2)y11.20160.5494y01.20160.549422.3004y(0.4)y21.20160.5494y11.20160.54942.30042.4654
