龙格库塔方法及其matlab实现

资料范本
本资料为word版本，可以直接编辑和打印，感谢您的下载
龙格库塔方法及其matlab实现
地点：__________________
时间：__________________
说明：本资料适用于约定双方经过谈判，协商而共同承认，共同遵守的责任与义务，仅供参考，文档可直接下载或修改，不需要的部分可直接删除，使用时请详细阅读内容
龙格-库塔方法及其matlab实现摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具达到求解目的。

龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。

MatLab软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。

MatLab 是英文MATrix LABoratory（矩阵实验室）的缩写。

在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。

关键词：龙格-库塔 matlab 微分方程
前言
1.1：知识背景
龙格－库塔法（Runge-Kutta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。

这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。

通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。

该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。

Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。

经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。

从这时起，MATLAB的内核采用C语言编写，而且除原有的数值计算能力外，还新增了数据图视功能。

MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行的可靠性，使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND和SIMNON，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。

在时间进
入20世纪九十年代的时候，MATLAB已经成为国际控制界公认的标准计算软件。

到九十年代初期，在国际上30几个数学类科技应用软件中，MATLAB在数值计算方面 HYPERLINK
"/s?wd=%E7%8B%AC%E5%8D%A0%E9%B3%8C%E5%A4%B4&hl_ta g=textlink&tn=SE_hldp01350_v6v6zkg6" \t "_blank" 独占鳌头，而Mathematica和Maple则分居符号计算软件的前两名。

Mathcad因其提供计算、图形、文字处理的统一环境而深受中学生欢迎。

1.2研究的意义
精确求解数值微分方程，对龙格库塔的深入了解与正确运用，主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。

利用matlab强大的数值计算功能，省去认为计算的过程，达到快速精确求解数值微分方程。

在实际生活中可以利用龙格库塔方法和matlab的完美配合解决问题。

1.3研究的方法
对实例的研究对比，实现精度的要求，龙格库塔是并不是一个固定的公式，所以只是对典型进行分析
龙格-库塔方法
2.1龙格-库塔公式
在一阶精度的的拉格朗日中值定理有：
对于函数y=f（x，y）
y'=f(x,y)
y(n+1)=y(n)+h*K1
K1=f(xn, yn)
这就是一阶龙格-库塔方法
形如 y(n+1)=y(n)+h*i=1rciki
k1 =f（xn，yn）
ki=f（xn+hai ，yn+h*j=1i-1bijki）
i=2…r
故二阶龙格-库塔公式
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'(xn)+h22y'’ (xn)+o（h3）
y'（x）= f（x，y（x））
y'’x= fx‘（x，y（x））+ fy‘（x，y（x））·f(x，y（x）)
y（xn+1）=y(xn)+hf+h22(fx‘+fy‘f)+ o（h3） (3)
将（2）式中的k2在（xn，yn）点展成幂级数
k2= f（xn+ha2 ，yn+ha2 k1）
=f+ha2fx‘+ ha2fy‘f+ o（h2）
将k1，k2代入（2）式，得
yn+1=yn +h（c1+c2）f
+ha2c2（fx‘+fy‘f）+ o（h3）（4）
对比（3）（4），当y(xn)= yn时
只有c1+c2=1，a2c2=12 （5）形如(2)存在常数满足（5）式，局部截断误差为o（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,1…N-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=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）
四阶龙格库塔法的推导类似于三阶龙格-库塔法，但相对复杂这里不再进行推导，公式如下
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,1…N-1
这就是标准四阶龙格库塔公式
2.1 对实例的研究
利用龙格-库塔法求解方程
y'=8-3yy0=2的数值，其中h=0.2，计算y(0.4)的近似值。

至少保留四位小数。

解：f (x，y)＝8－3y，利用四阶龙格-库塔公式有
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.208－1.578yn
n=0,1…N-1
yn+1=yn1.2016＋0.5494yn
当x0＝0，y0＝2，
y(0.2)≈y1＝1.2016＋0.5494y0＝1.2016＋0.5494×2＝2.3004
y(0.4)≈y2＝1.2016＋0.5494y1＝1.2016＋0.5494×2.3004＝2.4654。