一、理论基础

四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“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: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为步数
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的常微分方程的求解matlab仿真相关推荐

基于matlab的ldpc编码的构造,基于LDPC编码的GMSK调制与解调及matlab仿真实现(含录像)...
基于LDPC编码的GMSK调制与解调及matlab仿真实现(含录像)(开题报告,论文10700字,程序代码,录像) 摘要随着无线通信技术的不断发展与进步,数字电视广播.移动视频点播等对数据吞吐量要 ...
基于二维切片图序列的三维立体建模MATLAB仿真
目录 1.算法概述 2.仿真效果预览 3.核心MATLAB程序 4.完整MATLAB程序 1.算法概述 isosurface 等值面函数调用格式: fv = isosurface(X,Y,Z,V,i ...
基于沙猫群优化算法的线性规划求解matlab程序
基于沙猫群优化算法的线性规划求解matlab程序 1 沙猫群优化算法沙猫的中文学名叫沙丘猫,俗名沙漠猫,与荒漠猫名字相似,但却是两种不同的猫科动物.沙猫生活在茫茫沙漠里,主要分布在分布于非洲北部,阿 ...
【LDPC-11】基于QC-LDPC的CDR系统LDPC编码实现与matlab仿真验证
目录 1.基于QC-LDPC的CDR系统LDPC编码理论概述 2.matlab编程实现 3.编程验证 1.基于QC-LDPC的CDR系统LDPC编码理论概述中国数字音频广播CDR是一种 ...
一个简单的基于形态学处理的报纸图像版面分割算法matlab仿真
目录一.理论基础二.部分MATLAB仿真三.仿真结论分析一.理论基础图像分割就是把图像分成若干个特定的.具有独特性质的区域并提出感兴趣目标的技术和过程.它是由图像处理到图像分析的关键步骤.现 ...
基于滑模控制的直接转矩控制的MATLAB仿真
模型是基于袁雷老师主编的<现代永磁同步电机控制原理机MATLAB仿真>一书为参考. 个人认为使用滑模控制的DTC控制器有几个原因:1 系统存在滑模态 2点击要求转矩与磁链小脉动3逆变器需要 ...
基于等波纹最佳逼近法的FIR数字滤波器实现matlab仿真
目录一.理论基础二.案例背景三.MATLAB核心代码四.仿真结论分析一.理论基础等波纹最佳逼近法,其本质是一种优化算法,该方法有效克服了基于窗函数的FIR滤波器设计方法以及基于频率抽样的F ...
【OFDM频域同步】基于OFDM数字电视地面广播系统中频域同步技术matlab仿真
1.软件版本 matlab2021a 2.本算法理论知识数字电视地面广播传输一直是各界研究的热点,采用OFDM多载波调制可以有效对抗地面传输中多径效应.本课题主要要求学生掌握数字电视地面广播传输系统 ...
【OMP信道估计】基于OMP压缩感知的信道估计算法的MATLAB仿真
1.软件版本 MATLAB2021a 2.本算法理论知识 3.核心代码 clc; clear; close all; warning off; addpath 'func\'CYC = 20; for ...

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

一、理论基础

二、核心程序

三、测试结果

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

最新文章

热门文章