目录
一、理论基础
二、核心程序
三、测试结果
一、理论基础
四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。令初值问题表述如下。
这样，下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均：
k1是时间段开始时的斜率；
k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值；
k3也是中点的斜率，但是这次采用斜率k2决定y值；
k4是时间段终点的斜率，其y值用k3决定。
当四个斜率取平均时，中点的斜率有更大的权值：
RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。
注意上述公式对于标量或者向量函数(y可以是向量)都适用。
显式龙格库塔法
显示龙格－库塔法是上述RK4法的一个推广。它由下式给出
(注意：上述方程在不同著述中由不同但却等价的定义)。
要给定一个特定的方法，必须提供整数s (阶段数)，以及系数 aij (对于1 ≤ j < i ≤ s), bi (对于i = 1, 2, ..., s)和ci (对于i = 2, 3, ..., s)。这些数据通常排列在一个助记工具中，称为龙格库塔表：
| c2 | a21 | |||
| c3 | a31 | a32 | ||
| | ||||
| cs | as1 | as2 | | as,s ? 1 |
| b1 | b2 | | bs ? 1 | bs |
龙格库塔法是自洽的，如果
如果要求方法有精度p则还有相应的条件，也就是要求舍入误差为O(hp+1)时的条件。这些可以从舍入误差本身的定义中导出。
二、核心程序
clc;
clear;
close all;
warning off;
pack;
u = zeros(9,1);
Step = 3000;
R1 = zeros(Step,1);
R2 = zeros(Step,1);
R3 = zeros(Step,1);
y = zeros(3,Step);
y(:,1)= [450;541;600];
R1(1) = y(1,1);
R2(1) = y(2,1);
R3(1) = y(3,1);
h = 0.4;
for j = 2:Step
u(1) = y(1,j-1);
u(2) =y(2,j-1);
u(3) = y(3,j-1);
u(4) = 574;
u(5) = 470;
u(6) = 27.5;
u(7) = 283.4;
u(8) = 731.9;
u(9) = 950;
T_s_jw_out = u(1);
T_s_out = u(2);
T_j = u(3);
D_s_jw_in = u(4);
T_s_jw_in = u(5);
D_w = u(6);
T_w = u(7);
D_y_in = u(8);
T_y_in = u(9);
I_jw = 5000;
I_s = 5000;
c_j = 435;
c_p_y = 10;
k1 = 60;
A1 = 5.9;
k2 = 2000;
A2 = 5.9;
%*************************************************************************************************************************************
Ky1_1 = [(-(D_s_jw_in+D_w)*h_s(T_s_jw_out)+D_s_jw_in*h_s(T_s_jw_in)+D_w*h_w(T_w))/I_jw];
Ky1_2 = [(-(D_s_jw_in+D_w)*h_s(T_s_jw_out+h/2)+D_s_jw_in*h_s(T_s_jw_in)+D_w*h_w(T_w))/I_jw];
Ky1_3 = [(-(D_s_jw_in+D_w)*h_s(T_s_jw_out+h)+D_s_jw_in*h_s(T_s_jw_in)+D_w*h_w(T_w))/I_jw];
Ky1_4 = [(-(D_s_jw_in+D_w)*h_s(T_s_jw_out+3*h/2)+D_s_jw_in*h_s(T_s_jw_in)+D_w*h_w(T_w))/I_jw];
y(1,j)= y(1,j-1) + h*(Ky1_1+Ky1_2+Ky1_2+Ky1_3+Ky1_3+Ky1_4)/6;
%*************************************************************************************************************************************
%*************************************************************************************************************************************
Ky2_1 = [((D_s_jw_in+D_w)*h_s(T_s_jw_out)-(D_s_jw_in+D_w)*h_s(T_s_out)-k2*A2*T_s_out+k2*A2*T_j)/I_s];
Ky2_2 = [((D_s_jw_in+D_w)*h_s(T_s_jw_out+h/2)-(D_s_jw_in+D_w)*h_s(T_s_out+h/2)-k2*A2*(T_s_out+h/2)+k2*A2*(T_j+h/2))/I_s];
Ky2_3 = [((D_s_jw_in+D_w)*h_s(T_s_jw_out+h)-(D_s_jw_in+D_w)*h_s(T_s_out+h)-k2*A2*(T_s_out+h)+k2*A2*(T_j+h))/I_s];
Ky2_4 = [((D_s_jw_in+D_w)*h_s(T_s_jw_out+3*h/2)-(D_s_jw_in+D_w)*h_s(T_s_out+3*h/2)-k2*A2*(T_s_out+3*h/2)+k2*A2*(T_j+3*h/2))/I_s];
y(2,j)= y(2,j-1) + h*(Ky2_1+Ky2_2+Ky2_2+Ky2_3+Ky2_3+Ky2_4)/6;
%*************************************************************************************************************************************
%*************************************************************************************************************************************
Ky3_1 = [(k2*A2*T_s_out-(((k1*A1+k2*A2)*D_y_in*c_p_y+k1*A1*k2*A2)*T_j-k1*A1*D_y_in*c_p_y*T_y_in)/(k1*A1+D_y_in*c_p_y))/(c_j*A1*78.5)];
Ky3_2 = [(k2*A2*(T_s_out+h/2)-(((k1*A1+k2*A2)*D_y_in*c_p_y+k1*A1*k2*A2)*T_j-k1*A1*D_y_in*c_p_y*T_y_in)/(k1*A1+D_y_in*c_p_y))/(c_j*A1*78.5)];
Ky3_3 = [(k2*A2*(T_s_out+h)-(((k1*A1+k2*A2)*D_y_in*c_p_y+k1*A1*k2*A2)*T_j-k1*A1*D_y_in*c_p_y*T_y_in)/(k1*A1+D_y_in*c_p_y))/(c_j*A1*78.5)];
Ky3_4 = [(k2*A2*(T_s_out+3*h/2)-(((k1*A1+k2*A2)*D_y_in*c_p_y+k1*A1*k2*A2)*T_j-k1*A1*D_y_in*c_p_y*T_y_in)/(k1*A1+D_y_in*c_p_y))/(c_j*A1*78.5)];
y(3,j)= y(3,j-1) + h*(Ky3_1+Ky3_2+Ky3_2+Ky3_3+Ky3_3+Ky3_4)/6;
%*************************************************************************************************************************************
R1(j)= y(1,j);
R2(j)= y(2,j);
R3(j)= y(3,j);
end
figure;
subplot(311);
plot(R1(1:200));
xlabel('仿真时间');
ylabel('T-s-jw-out');
legend('龙格库塔仿真结果');
subplot(312);
plot(R2);
xlabel('仿真时间');
ylabel('T-s-out');
legend('龙格库塔仿真结果');
subplot(313);
plot(R3);
xlabel('仿真时间');
ylabel('T-j');
legend('龙格库塔仿真结果'); %自编四阶龙格库塔法解微分方程，并与ode45的计算结果比较
function Y1 = func_4RGKT(f,a,b,ya,m)
%f为要求的函数，a,b分别为上下限，ya为y的初值，m为步数
clc
format long
%算步长h
h = (b - a)/m;
%建立1*m+1矩阵，并赋给T,Y
T = zeros(1,m+1);
Y = zeros(1,m+1);
%给初值赋值
T(1) = a;
Y(1) = ya;
for j=1:m
tj = T(j);
yj = Y(j);
k1 = h*feval（f,tj,yj);
k2 = h*feval（f,tj+h/2,yj+k1/2);
k3 = h*feval（f,tj+h/2,yj+k2/2);
k4 = h*feval（f,tj+h,yj+k3);
Y(j+1) = yj + (k1 + 2*k2 + 2*k3 + k4)/6;
T(j+1) = a + h*j;
end
%将计算结果赋给Y
Y1=Y'; 三、测试结果
从上面的仿真结果为当迭代次数大于40的时候，采用龙格库塔算法的精度非常接近真实的值，因此，在实际仿真过程中，我们一般将迭代次数设置为至少40。