龙格库塔方法及其matlab实现

龙格-库塔方法及其matlab实现

摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具

达到求解目的。龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。MatLab软件是由美国Mathworks公司推出的用于数值计算和图

形处理的科学计算系统环境。MatLab是英文MATrix LABoratory（矩阵实验室）的缩写。

在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文

件管理等各项操作。

关键词：龙格-库塔 matlab 微分方程

1.前言

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在数值计算方面独

占鳌头，而Mathematica和Maple则分居符号计算软件的前两名。Mathcad因其提供计算、图形、文字处理的统一环境而深受中学生欢迎。

1.2研究的意义

精确求解数值微分方程，对龙格库塔的深入了解与正确运用，主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。利用matlab 强大的数值计算功能，省去认为计算的过程，达到快速精确求解数值微分方程。在实际生活中可以利用龙格库塔方法和matlab 的完美配合解决问题。

1.3研究的方法

对实例的研究对比，实现精度的要求，龙格库塔是并不是一个固定的公式，所以只是对典型进行分析

2. 龙格-库塔方法

2.1龙格-库塔公式

在一阶精度的的拉格朗日中值定理有：

对于函数y=f （x ，y ）

y'=f(x,y)

y(n+1)=y(n)+h*K1

K1=f(x n , y n )

这就是一阶龙格-库塔方法

形如 y(n+1)=y(n)+h*∑c i r i=1k i

k 1 =f （x n ，y n ）

k i =f （x n +ha i ，y n +h*∑b ij i−1j=1k i ）

i=2…r

故二阶龙格-库塔公式

y(n+1)=y(n)+h （c 1k 1+c 2k 2）

k 1= f （x n ，y n ） (2)

k 2= f （x n +ha 2 ，y n +h a 2 k 1）

将y （x ）在x n 处展成幂级数

y （x n+1）=y(x n )+h y ′(x n )+h 22y ′’ (x n )+o （h 3）

y ′（x ）= f （x ，y （x ））

y ′’(x )= f x ‘（x ，y （x ））+ f y ‘（x ，y （x ））·f(x ，y （x ）)

y （x n+1）=y(x n )+hf+h 22(f x ‘+f y ‘f )+ o （h 3） (3) 将（2）式中的k 2在（x n ，y n ）点展成幂级数

k 2= f （x n +ha 2 ，y n +h a 2 k 1）

=f+h a 2f x ‘+ h a 2f y ‘f + o （h 2）

将k 1，k 2代入（2）式，得

y n+1=y n +h （c 1+c 2）f

+h a 2c 2（f x ‘+f y ‘f ）+ o （h 3） （4）

对比（3）（4），当y(x n )= y n 时

只有c 1+c 2=1，a 2c 2=12 （5） 形如(2)存在常数满足（5）式，局部截断误差为o （h 3）的求解方法称为二阶龙格-库塔法。

满足（5）式，若取c 1=12，则得到c 2=12，a 2=1，则公式则恰为预估-校正法公式 若取c 1=0，则c 2=1，a 2=12，

y n+1=y n +hk 2

k 1= f （x n ，y n ） (6)

k 2= f （x n +h 2，y n +h 2k ）

n=0,1…N-1

由（5）式，可知龙格-库塔法不是唯的

三阶龙格-库塔法

y n+1=y n +h （c 1k 1+c 2k 2+ c 3k 3）

k 1= f （x n ，y n ）

k 2= f （x n +h a 2，y n +ha 2k 1） （7）

k 3= f （x n +h a 3，y n +hb 31k 1+hb 32k 2）

若c 1，c 2， c 3，a 2，a 3，b 31, b 32且满足b 31+ b 32=a 3，，并使得局部截断误差为o

（h 4）。类似二阶龙格-库塔法推导的

c 1+c 2+ c 3=1