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

目录

一、理论基础

二、核心程序

三、测试结果

---

一、理论基础

       常微分方程的求解是数学中的重要问题，而龙格-库塔法是一种常用的数值方法，它可以用于求解一阶常微分方程。下面我们将介绍龙格-库塔法的基本原理和数学公式。

一、龙格-库塔法的基本原理
        龙格-库塔法是一种基于数值逼近的方法，它通过构造一组离散点来逼近函数的曲线，从而得到函数值的近似解。具体来说，它通过构造一组泰勒级数的和来逼近常微分方程的解。
考虑一阶常微分方程：
y'=f(x,y)
其中f是一个给定的函数，y是未知函数，x是自变量。
龙格-库塔法的基本思想是构造一组离散点(x_i, y_i)，使得这些点逼近方程的解曲线。为了做到这一点，我们选择一个步长h，并计算一系列点(x_i, y_i)，其中i=0,1,2,...,n，使得每个点满足方程：
y_i+1=y_i+hf(x_i,y_i)
其中y_i+1表示下一个离散点(x_i+h, y_i+h)的y坐标，x_i表示当前离散点的x坐标。
这样，通过逐步迭代，我们可以得到一组离散点，这些点逼近方程的解曲线。通过适当的调整步长h，我们可以控制解的精度。

二、龙格-库塔法的数学公式
龙格-库塔法有很多种实现方式，其中最常用的是四阶龙格-库塔法。它使用了四个参数a, b, c, d来构造泰勒级数的和，具体公式如下：
y_0=f(x0)
y_1=y0+haf(x0,y0)
y_2=y1+hbf(x1,y1)
y_3=y2+hcf(x2,y2)
y_4=y3+hd*f(x3,y3)
其中，x0, x1, x2, x3分别表示初始离散点以及后续三个离散点的x坐标，y0, y1, y2, y3分别表示相应离散点的y坐标，h表示步长。
       通过调整参数a, b, c, d的值，我们可以得到不同精度的解。此外，如果步长h选择合适，四阶龙格-库塔法可以提供O(h^4)阶的精度。
       需要注意的是，龙格-库塔法是一种数值方法，它的解是近似解。因此，在求解常微分方程时，我们需要选择合适的步长h和迭代次数n，以确保解的精度满足要求。同时，还需要进行适当的误差分析和收敛性判断。

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

龙格库塔法是自洽的，如果

如果要求方法有精度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