龙格—库塔四阶五级方程MATLAB求解方程

龙格-库塔一般形式 

 四级五阶龙格库塔

 代码

脚本文件

%m文件
%y'=3*y+x,y(0)=1,x从0~2
h=0.01;
X=0:h:2;
y0=1;    
n=length(X);
[x,y]= r_k_xy(X,y0,h,n);
for j=1:n
dy=r_k_f(x,y);
end
%eval函数：和文本变量一起使用，为有利的文本宏工具
String1='[x;y]';
String2='[x;dy]';
output1=eval(String1)
output2=eval(String2)
%绘图
figure;
plot(x,y,'r^','MarkerSize',3);xlim([0 2]);hold on;
plot(x,dy,'b-','LineWidth',2);
legend('y''=3*y+x','y=y(x)');
xlabel('x');
title('y-x-y''');

 函数文件

%方程
%y'=3*y+x,y(0)=1,x从0~2
%即dy=3*y+x
function [x,y]= r_k_xy(X,y0,h,n)
x=0:h:2;
y=zeros(1,n);
y(1)=y0
for j=1:n-1
    k1=r_k_f(x(j),y(j));
    k2=r_k_f(x(j)+h/2,y(j)+h*k1/2);
    k3=r_k_f(x(j)+h/2,y(j)+h*(k1+k2)/4);
    k4=r_k_f(x(j)+h/2,y(j)+h*(k1+k2+k3)/6);
    k5=r_k_f(x(j)+1/2,y(j)+h*(k1+k2+k3+k4)/8);
    y(j+1)=y(j)+h*(k1+k2+k3+k4+k5)/5;
end
end

%% 斜率
function dy=r_k_f(x,y)
dy=3*y+x;
end

 结果

 2022.4.26