MATLAB四阶Runge-Kutta代码无法显示Sigmoid曲线求助

在MATLAB中使用四阶Runge-Kutta（RK4）方法编写动力学模型代码，尝试将Sigmoid函数加入M1的更新逻辑，但绘制的M1曲线无法呈现Sigmoid形态。以下是主程序代码、各变量对应的函数文件及当前绘图情况：

主程序代码 #

clc;clear;
  
% 时间设置
t(1)=0;
dt=0.1; % 时间步长
t=0:dt:100; % 时间范围

% 初始化数组
T=zeros(length(t),1); % 时间数组
M1=zeros(length(t),1); % M1数组
M2=zeros(length(t),1); % M2数组
M3=zeros(length(t),1); % M3数组
O=zeros(length(t),1); % O数组
P=zeros(length(t),1); % P数组

% 初始条件（已注释）
    %M1 =10;
    %M2 = 0;
    %M3 = 0;
    %O =0;
    %P =0;

for j = 1:length((t))
    T(j+1)=T(j)+dt;
    M1(j+1)= M1(j)+1./(1+exp(-T(j))); % 尝试加入的Sigmoid更新
    
    % M2的RK4迭代
    k1M2 = dt*fRK4M2(M2(j),M1(j));
    k2M2 = dt*fRK4M2(M2(j)+k1M2/2,M1(j)+k1M2/2);
    k3M2 = dt*fRK4M2(M2(j)+k2M2/2,M1(j)+k2M2/2);
    k4M2 = dt*fRK4M2(M2(j)+k3M2,M1(j)+k3M2/2);
    M2(j+1) = M2(j)+1/6*(k1M2+2*k2M2+2*k3M2+k4M2);
    
    % M3的RK4迭代
    k1M3 = dt*fRK4M3(M1(j),M3(j));
    k2M3 = dt*fRK4M3(M3(j)+k1M3/2,M1(j)+k1M3/2);
    k3M3 = dt*fRK4M3(M3(j)+k2M3/2,M1(j)+k2M3/2);
    k4M3 = dt*fRK4M3(M3(j)+k3M3,M1(j)+k3M3);
    M3(j+1) = M3(j)+1/6*(k1M3+2*k2M3+2*k3M3+k4M3);
    
    % O的RK4迭代
    k1O = dt*fRK4O(O(j),M1(j));
    k2O = dt*fRK4O(O(j)+k1O/2,M1(j)+k1O/2);
    k3O = dt*fRK4O(O(j)+k2O/2,M1(j)+k2O/2);
    k4O = dt*fRK4O(O(j)+k3O,M1(j)+k3O);
    O(j+1) = O(j)+1/6*(k1O+2*k2O+2*k3O+k4O);
    
    % P的RK4迭代
    k1P = dt*fRK4P(P(j),M1(j));
    k2P = dt*fRK4P(P(j)+k1P/2,M1(j)+k1P/2);
    k3P = dt*fRK4P(P(j)+k2P/2,M1(j)+k2P/2);
    k4P = dt*fRK4P(P(j)+k3P,M1(j)+k3P/2);
    P(j+1) = P(j)+1/6*(k1P+2*k2P+2*k3P+k4P);
    
    % M1的RK4迭代（参数顺序错误）
    k1M1= dt*fRK4M1(M1(j),M2(j),M3(j),O(j),P(j));
    k2M1= dt*fRK4M1(M2(j)+k1M2/2,M3(j)+k1M3/2,O(j)+k1O/2,P(j)+k1P/2,M1(j)+k1M1/2);
    k3M1= dt*fRK4M1(M2(j)+k2M2/2,M3(j)+k2M3/2,O(j)+k2O/2,P(j)+k2P/2,M1(j)+k2M1/2);
    k4M1= dt*fRK4M1(M2(j)+k3M2/2,M3(j)+k3M3/2,O(j)+k3O/2,P(j)+k3P/2,M1(j)+k3M1/2);
    M1(j+1) = M1(j)+(1/6*(k1M1+(2*k2M1)+(2*k3M1)+k4M1));
end
  
% 绘图
figure;
plot (T,M1,'r','Linewidth',3)
xlabel('时间')
ylabel('M1')

figure 
plot (T,M2,'b','Linewidth',3)
xlabel('时间')
ylabel('M2')

figure 
plot (T,M3,'g','Linewidth',3)
xlabel('时间')
ylabel('M3')

figure 
plot (T,O,'b','Linewidth',5)
xlabel('时间')
ylabel('O')

%figure 
%plot (T,P,'r','Linewidth',5)
%xlabel('时间')
%ylabel('P')

函数文件（分文件存储） #

fRK4M1.m #

function M1 =fRK4M1(M1,M2,M3,O,P) 
  delta=50;
  K1= 10^-4;
  Ko=0.1;
  n=3;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  K2=5*10^-4;
  K3=10^-3;
  gamma=75;
 
  M1=(delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2.*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1)); 
end

fRK4M2.m #

function M2 =fRK4M2(M2,M1)
M1=10; % 硬编码覆盖传入参数
mu_2=10^-3;
K1= 10^-4;
K2=5*10^-4;     

M2 = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);
end

fRK4M3.m #

function M3 =fRK4M3(M3,M1)
M2=0; % 硬编码覆盖传入参数
mu_3=10^-3;
K2=5*10^-4;
K3=10^-3;

M3=(K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);
end

fRK4O.m #

function O =fRK4O(O,M1)
M3=0; % 硬编码覆盖传入参数
Ko=0.1;
mu_o=10^-4;
K3=10^-3;

O = (K3*M1*M3)-(Ko*M1*O)-(mu_o*O);
end

fRK4P.m #

function P =fRK4P(P,M1)
O=0; % 硬编码覆盖传入参数
Ko=0.1;
mu_p= 10^-5;

P  = (Ko*M1*O)-(mu_p*P);
end

当前绘图情况 #

• M1、M2、M3曲线：M1无Sigmoid形态，走势异常；M2、M3曲线无预期动力学变化
• O、P曲线：几乎维持在0值，无明显波动

以下是针对性的代码修改步骤：

1. 移除M1的重复更新，整合Sigmoid逻辑 #

主程序中先执行了Sigmoid更新，随后又用RK4重新计算M1(j+1)，导致Sigmoid效果被完全覆盖。需删除该行，并将Sigmoid项加入M1的微分方程：
修改fRK4M1.m：

function dM1 = fRK4M1(M1,M2,M3,O,P) 
  delta=50;
  K1= 10^-4;
  Ko=0.1;
  n=3;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  K2=5*10^-4;
  K3=10^-3;
  gamma=75;
  % 加入Sigmoid项，可调整缩放因子控制拐点位置和斜率
  sigmoid_term = 1/(1+exp(-(M1-30)/5)); 
 
  dM1 = delta*M1*(1-(M1/gamma)) - 2*K1*M1^2 - K2*M1*M2 - (Oa-n)*K3*M1*M3 - (Pa-Oa)*Ko*M1*O - mu_1*M1 + sigmoid_term; 
end

2. 修正RK4调用的参数顺序错误 #

fRK4M1的参数顺序为(M1,M2,M3,O,P)，但主程序中调用时参数顺序完全混乱，修改如下：

% 修正M1的RK4迭代参数顺序
k1M1= dt*fRK4M1(M1(j),M2(j),M3(j),O(j),P(j));
k2M1= dt*fRK4M1(M1(j)+k1M1/2, M2(j)+k1M2/2, M3(j)+k1M3/2, O(j)+k1O/2, P(j)+k1P/2);
k3M1= dt*fRK4M1(M1(j)+k2M1/2, M2(j)+k2M2/2, M3(j)+k2M3/2, O(j)+k2O/2, P(j)+k2P/2);
k4M1= dt*fRK4M1(M1(j)+k3M1/2, M2(j)+k3M2/2, M3(j)+k3M3/2, O(j)+k3O/2, P(j)+k3P/2);
M1(j+1) = M1(j)+(1/6*(k1M1+2*k2M1+2*k3M1+k4M1));

3. 删除函数内部的硬编码变量 #

所有函数中硬编码的M1=10;、M2=0;等语句会覆盖传入的动态参数，必须删除：

• fRK4M2.m：删除M1=10;
• fRK4M3.m：删除M2=0;
• fRK4O.m：删除M3=0;
• fRK4P.m：删除O=0;

4. 启用初始条件并修正时间数组 #

取消初始条件注释，同时删除冗余的T数组构建，直接使用已定义的t数组：

% 启用初始条件
M1(1) =10;
M2(1) = 0;
M3(1) = 0;
O(1) =0;
P(1) =0;

% 调整循环范围避免数组越界
for j = 1:length(t)-1
    % 删除T数组构建代码
    % T(j+1)=T(j)+dt;
    
    % ... 其余RK4计算代码不变 ...
end

% 绘图时使用t数组
figure;
plot(t, M1,'r','Linewidth',3)
xlabel('时间')
ylabel('M1')
% 其余绘图同理替换T为t

5. 调整Sigmoid项的输入与系数 #

若M1仍未呈现预期形态，可微调fRK4M1中的Sigmoid参数：

• 修改(M1-30)/5中的数值：30控制拐点位置，5控制斜率
• 给Sigmoid项乘以系数（如5*sigmoid_term）增强其影响

内容的提问来源于stack exchange，提问作者Cindy