在 MATLAB 中比较“ODE45”函数与 4 阶 Runge-Kutta 对流感疾病的结果

Question

我为一个非线性微分方程系统编写了一个代码，该系统对流感疾病进行建模。当我想用 4 阶 Runge-Kutta 和 ode45function 模拟这个系统时，我感到非常震惊。问题是这两种方法的结果是不一样的。事实上，当我为 4 阶 Runge-Kutta 算法编写代码时，状态会在每一步中变得越来越大，直到达到无穷大。但是当我用ode45function 模拟系统时，结果是合理合理的。我不认为我写的算法是不正确的。有人可以向我解释我做错了什么吗？如果代码中的某些内容无效？

我写的代码如下：

clc;clear;close all;

t = linspace(0,100,101)';
global alpha eta p f ep delta q beta kessi

kessi = 0.526;                % Transition rate for the exposed
alpha = 0.244;                % Recovery rate for the (symptomatic)
eta = alpha;                  % Infected
p = 0.667;                    % Recovery rate for the asymptomatic
f = 0.98;                     % 1-fatallity rate(one minus fatality rate)
ep = 0;                       % Infectivity reduction factor for the exposed
delta = 1;                    % Infectivity reduction factor for the asymptomatic
q = 0.5;                      % Contact reduction by isolation
beta = 1;                     % From CALCULATIONS

Influenza_Fun = @(t,x,u) [-beta*x(1)*(ep*x(2)+(1-q)*x(3)+delta*x(4))-u(1)*x(1)
                           beta*x(1)*(ep*x(2)+(1-q)*x(3)+delta*x(4))-kessi*x(2)
                            p*kessi*x(2)-alpha*x(3)-u(2)*x(3)
                            (1-p)*kessi*x(2)-eta*x(4)
                            f*alpha*x(3)+u(2)*x(3)+eta*x(4)+u(1)*x(1)];


% Initial Condtions
Ics = [15000 200 500 300 0]';                     % [S0 E0 I0 A0 R0]';
u = [0 0]';

[~,States] = ode45(@(t,x)Influenza_Fun(t,x,u),t,Ics);

for i=1:5
   plot(t,States(:,i),'LineWidth',2)
   hold on
   grid on

end
legend('S','E','I','A','R')
legend show

以上代码将基于ode45function 模拟系统。以下代码是我的 RK4 算法：

clc;clear;close all;

% { Control of Influenza disease based on Genetic Algorithm using desired Objective FUnction}

global alpha eta p f ep delta q beta kessi


kessi = 0.526;                % Transition rate for the exposed
alpha = 0.244;                % Recovery rate for the (symptomatic)
eta = alpha;                  % Infected
p = 0.667;                    % Recovery rate for the asymptomatic
f = 0.98;                     % 1-fatallity rate(one minus fatality rate)
ep = 0;                       % Infectivity reduction factor for the exposed
delta = 1;                    % Infectivity reduction factor for the asymptomatic
q = 0.5;                      % Contact reduction by isolation
beta = 1;                     % From CALCULATIONS

% Initial Condtions
Ics = [15000 200 500 300 0]';                     % [S0 E0 I0 A0 R0]';

max_days = 100;
y(:,1) = Ics;
t = zeros(1,max_days);

h = 1;
t(1) = 0;
LB = [0 0]';
UB = [1 1]';
g = @yMatrix; 
u = [0 0]';

Influenza_Fun = @(t,x,u) [-beta*x(1)*(ep*x(2)+(1-q)*x(3)+delta*x(4))-u(1)*x(1)
                           beta*x(1)*(ep*x(2)+(1-q)*x(3)+delta*x(4))-kessi*x(2)
                            p*kessi*x(2)-alpha*x(3)-u(2)*x(3)
                            (1-p)*kessi*x(2)-eta*x(4)
                            f*alpha*x(3)+u(2)*x(3)+eta*x(4)+u(1)*x(1)];


%% 2. mAiN Loop

for time = 1:max_days


    k1 = Influenza_Fun(t(time),y(:,time),u);
    k2 = Influenza_Fun(t(time)+(h/2),y(:,time)+k1,u);
    k3 = Influenza_Fun(t(time)+(h/2),y(:,time)+k2,u);
    k4 = Influenza_Fun(t(time)+h,y(:,time)+k3,u);

    % y(i+1) = y(i) + (h/6)*(k1+2k2+2k3+k4)
    y(:,time+1) = y(:,time)+(h/6)*(k1+2*k2+2*k3+k4);
    t(time+1) = t(time)+h;

end

y = y';

for i=1:size(y,2)

    plot(t,y(:,i),'LineWidth',2)
    hold on
    grid on


end

Answer 1

应该是

    k2 = Influenza_Fun(t(time)+(h/2),y(:,time)+(h/2)*k1,u);

等等

还可以使用某种方式来展示 ode45 实际使用的时间步长，例如使用 [ T,Y] = 输出格式，从 Refine 选项中扣除额外的分数。或者使用函数评估计数，至少文档中的最后一个调用变体应该显示这一点。在 RK4 中使用类似的步长或函数调用计数以获得可比较的结果。