用4阶龙格库塔公式求解微分方程组

公式真的很难编辑，还是直接放源码。

龙格库塔法:（runge-kutta)是一种高阶显示一步法，而且不需要计算导数。

间接使用泰勒展开式的方法：用f(x,y)在一些点上函数值的线性组合替代y(x)的各阶导数，构造yn+1的表达式；

比较这个表达式与y(xn+1)在x=xn处泰勒展开式，确定yn+1的表达式中的系数，使yn+1的表达式与y(xn+1)泰勒展开式前面的若干项相等，从而具有对应的泰勒展开式的精度阶次。

1.四阶龙格库塔公式函数文件

function varargout = RungeKutta_odes(odefun,xspan,y0,h)
% RungeKutta_odes求解显示微分方程组，用4阶龙格库塔公式
% 调用格式[t,y] = RungeKutta_odes(odefun,xspan,y0,h)
% 调用格式sol = RungeKutta_odes(odefun,xspan,y0,h)
% 输入参数，odefun定义微分方程组的M文件或者匿名形式
% xspan:求解区间，y0：初值，h:步长，缺省0.01
%varargout是可变参数，输出，可以是一个参数，两个参数，三个参数
%varargout是根据调用的参数确定的
%for example 
%[t,y] = RungeKutta_odes(odefun,xspan,y0,h)
%[t,y,m] = RungeKutta_odes(odefun,xspan,y0,h)

%%
%输入参数控制
if nargin == 3
    h = 0.01; % 缺省，步长精度
elseif nargin < 3
    error('此方法是求解带有初值问题的微分方程...')
elseif nargin > 4
    error('输入参数个数too many...')
end

%%
%变量的初始化
t = xspan(1):h:xspan(2); % 根据步长确定子区间，
n = length(t); % 统计变量t的数量，用于循环求解x的值
y = zeros(length(y0),n); % 根据初值判断微分方程组的个数，一行代表一个微分方程解
y(:,1) = y0(:); % 第一列存初值，

%%
%循环求解，使用4阶龙格库塔公式,核心代码
for i = 2:n
    k1 = odefun(t(i-1), y(:,i-1)); %y是所有rows,columns：i-1,公式：K1=f(tn,yn,zn)
    k2 = odefun(t(i-1)+h/2, y(:,i-1)+(h/2)*k1); % K2=f(tn+h/2,yn+h/2*K1,zn+h/2*L1)
    k3 = odefun(t(i-1)+h/2, y(:,i-1)+(h/2)*k2);
    k4 = odefun(t(i-1)+h, y(:,i-1)+h*k3);
    y(:,i) = y(:,i-1)+(h/6)*(k1+2*k2+2*k3+k4);
end

%%
%可变输出参数的判断
if nargout == 2
    varargout{1} = t; % varargout是一个数组，如果输出变量是两个，
    varargout{2} = y;
elseif nargout == 1
    varargout{1}.t = t;
    varargout{1}.y = y;
elseif varargout > 2
    error('输出参数过多too many');
end

%%
%数值可视化
leg = {};
for i = 1:length(y0)
    plot(t,y(i,:))%,'Marker','-','MarkerEdgeColor','green')
    hold on
    leg{i} = strcat('y_',num2str(i),'(t)');
end
hold off
grid on
xlabel('t');ylabel('y');

title('微分方程（组）的数值解曲线')
legend(leg)
legend('boxoff')


end

2.所求函数文件

function dy = dyfun2(t,y)
    dy = [y(2); % y(1)' = y(2)
        y(3); % y(2)' = y(3)
        y(3)/t-3*y(2)/(t^2)+2*y(1)/(t^3)+9*t^3*sin(t)]; % y(3)' 
end

3.主文件，调用两个函数

clc;clear;close all


xspan = [0.1,60];
%xspan = [0,20];
y0 = [1;1;1];
h = 0.02;
% 调用函数
sol = RungeKutta_odes(@dyfun2,xspan,y0,h)