基于龙格-库塔法Runge-Kutta的常微分方程的求解matlab仿真

Posted 2022-11-26 fpga和matlab

目录

一、理论基础

二、核心程序

三、测试结果

---

一、理论基础

       四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“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)。这些数据通常排列在一个助记工具中，称为龙格库塔表：

	0
	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。

 A16-21