常微分方程组-龙格库塔法

导论

        二阶龙格库塔法的思想及推导如下图所示，使用tn及tn+1两点处导数的平均值，且将tn+1处的导数值利用简单欧拉法代替，得到如下图中的公式。易知，龙格库塔方法的思想是通过采用不同点处的斜率加权平均作为递推斜率。

四阶龙格库塔方法

单步误差为，累积误差为

扩展到常微分方程组

clear
clc
T = 20; N = 10000; dt = T/N;
t = 0:dt:T;
y = zeros(3, length(t)); 
s = 10; b = 8/3; r = 28;
y(:, 1) = [-1.6; -0.5; 21];
for n = 1:length(t)-1
    k1 = fun(y(:,n),t(n));
    k2 = fun(y(:,n)+k1*dt/2,t(n)+dt/2);
    k3 = fun(y(:,n)+k2*dt/2,t(n)+dt/2);
    k4 = fun(y(:,n)+k3*dt,t(n)+dt);
    y(:,n+1) = y(:,n) + (k1+2*k2+2*k3+k4)/6*dt;
end
figure(1)
plot(t, y(1,:), t, y(2,:), t, y(3,:))
figure(2)
plot3(y(1,:), y(2,:), y(3,:))
function dydt = fun(y, t)
s = 10; b = 8/3; r = 28;
dydt = [s*(y(2)-y(1)); r*y(1)-y(2)-y(1)*y(3); -b*y(3)+y(1)*y(2)];
end