三、四阶Runge-Kutta法求解常微分方程
一、龙格库塔法的思想
根据第九章的知识可知道，Euler方法的局部截断误差是
，而当用Euler方法估计出
再用梯形公式
进行校正，即采用改进Euler方法得出数值解的截断误差为
。
由Lagrange微分中值定理
记
，得到
这样只要给出一种计算
的算法，就能得到相应的
计算公式
六西格玛计算公式下载结构力学静力计算公式下载重复性计算公式下载六西格玛计算公式下载年假计算公式
。
用这种观点的来分析Euler方法和改进Euler方法，Euler方法的迭代公式可改写为
改进Euler方法的预报-校正公式可改写为
Euler方法实际上是用一个点处的值
近似
，而改进Euler方法是用两个点处的值
，和
，做算术平均值近似
自然改进Euler方法要优于Euler方法。
因此，可以想到假如在
内多预报几个点值
，并用他们的加权平均值作为
的近似值，则有可能构造出具有更高精度的计算公式，这就是Runge-Kutta法的基本思想。
二、四阶龙格库塔法
由Runge-Kutta的基本思想，构造四阶Runge-Kutta法是利用
的加权平均值来近似
，因此令
使得
即其总体截断误差为
。
采用泰勒公式展开，经过复杂的推导，得到一个具有13个参数，11个方程的线性方程组。由于方程的个数少于未知量的个数，因此方程有无穷多个解。可以根据情况得到几种常用的解，即得到相应的四阶公式。最常见的四阶公式如式（6）：
，
也称为标准四阶Runge-Kutta法。
三、四阶龙格库塔法程序说明及应用
3.1龙格库塔的计算程序
function [x,y] =Runge(ydot_fun,x0,y0,h,N )
x=zeros(1,N+1);y=zeros(length(y0),N+1);
x(1)=x0;y(:,1)=y0;
for n=1:N
x(n+1)=x(n)+h;
k1=h*feval(ydot_fun,x(n),y(:,n));
k2=h*feval(ydot_fun,x(n)+1/2*h,y(:,n)+1/2*k1);
k3=h*feval(ydot_fun,x(n)+1/2*h,y(:,n)+1/2*k2);
k4=h*feval(ydot_fun,x(n)+h,y(:,n)+k3);
y(:,n+1)=y(:,n)+1/6*(k1+2*k2+2*k3+k4);
end
3.2 程序解释及使用
该算法可以对一阶微分方程，一阶微分方程组进行有效的求解。ydot_fun为一阶微分方程的函数，x0为初始点，y0为初始向量，h为步长，N为区间的等分数，x为Xn构成的向量，y为Yn构成的矩阵。
程序调用方法：
1，先编写要求解的一阶微分方程或方程组的函数文件ydot_fun.m文件，将该文件和Runge文件放到同一个目录下。
2. 调用求解程序，[x,y]=Runge(@dot_fun,x0,y0,h,N)，运行后即可得出结果。
或者用内部函数调用：
输入：ydot_fun=(x,y)[]
[x,y]= Runge(ydot_fun,x0,y0,h,N)
3.3实例求解
课本304页9.2题目：
用标准4级4阶R-K法求解
，
，
取步长h=0.1，计算
的近似值，并与解析解
作比较。
解：首先将三阶方程改写成微分方程组的形式：
令
得如下微分方程组
在ydot_fun.m文件中编写待求解微分方程组，调用计算程序，保留5位小数得：
表3-1 三阶微分方程求解结果
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
-1.00000
-0.68948
-0.35572
0.00496
0.39673
0.82436
1.29327
1.80962
2.38042
3.01363
3.71827
3.00000
3.21569
3.46568
3.75481
4.08855
4.47308
4.91538
5.42337
6.00596
6.67323
7.43655
2.00000
2.32086
2.68708
3.10467
3.58038
4.12180
4.73750
5.43712
6.23150
7.13284
8.15483
表3-1中第二行
为原三阶微分方程对应的数值解，第二行
为其一阶导数值，第三行
为其二阶导数值。
由结果可知，数值解y(1)=3.7183，其对应的精确解析解
3.71828188的相对误差为4.8947e-6，可知四阶龙格库塔法具有很高的代数精度。
浙公网安备 33021202002483