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

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

时间：2022-10-08 19:27:13

一、理论基础

二、核心程序

三、测试结果

一、理论基础

四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“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:Stepu(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);endfigure;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为步数clcformat long%算步长hh = (b - a)/m; %建立1*m+1矩阵，并赋给T,YT = 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%将计算结果赋给YY1=Y';