前言
数值分析中,龙格-库塔法(Runge-Kutta methods)是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法,其中包括著名的欧拉法,用于数值求解微分方程。由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。
在各种龙格-库塔法当中有一个方法十分常用,以至于经常被称为“RK4”或者就是“龙格-库塔法”。该方法主要是在已知方程导数和初值信息,利用计算机仿真时应用,省去求解微分方程的复杂过程。
令初值问题表述如下。
则,对于该问题的RK4由如下方程给出:
其中
这样,下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积所决定。该斜率是以下斜率的加权平均:
k1是时间段开始时的斜率;
k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn+h/2的值;
k3也是中点的斜率,但是这次采用斜率k2决定y值;
k4是时间段终点的斜率,其y值用k3决定。
当四个斜率取平均时,中点的斜率有更大的权值:
RK4法是四阶方法,也就是说每步的误差是h阶,而总积累误差为h阶。
注意上述公式对于标量或者向量函数(y可以是向量)都适用。
matlab代码实现
问题:dy/dt=y-t^2+1 ; 0<=t<=2 ; y(0)=0.5;
clear
clc
close all
f = @(t,y) (y-t^2+1);
a = input('开始时间, a: ');
b = input('结束时间, b: ');
n = input('步数, n: '); % n=(b-a)/h
alpha = input('初始值, alpha: ');
h = (b-a)/n;
t = a;
w = alpha;
fprintf(' t w\n');
fprintf('%5.3f %11.7f\n', t, w);
for i = 1:n
k1 = h*f(t,w);
k2 = h*f(t+h/2.0, w+k1/2.0);
k3 = h*f(t+h/2.0, w+k2/2.0);
k4 = h*f(t+h,w+k3);
w = w+(k1+2.0*(k2+k3)+k4)/6.0;
t = a+i*h;
fprintf('%5.3f %11.7f\n', t, w);
plot(t,w,'b*'); grid on;
xlabel('t values'); ylabel('w values');
hold on;
end
开始时间, a: 0
结束时间, b: 2
步数, n: 10
初始值, alpha: 0.5
t w
0.000 0.5000000
0.200 0.8292933
0.400 1.2140762
0.600 1.6489220
0.800 2.1272027
1.000 2.6408227
1.200 3.1798942
1.400 3.7323401
1.600 4.2834095
1.800 4.8150857
2.000 5.3053630