四阶龙格库塔（Runge-Kutta）求解微分方程-多种编程语言
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了四阶龙格库塔（Runge-Kutta）求解微分方程-多种编程语言相关的知识，希望对你有一定的参考价值。
✅作者简介：热爱科研的算法开发者，Python、Matlab项目可交流、沟通、学习。
基于龙格-库塔法Runge-Kutta的常微分方程的求解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
以上是关于四阶龙格库塔（Runge-Kutta）求解微分方程-多种编程语言的主要内容，如果未能解决你的问题，请参考以下文章
