手搓SIP模型

一、引言

在众多的动力学模型中，SIP 模型是一种用于描述系统动态变化的重要模型。它可以用来研究各种复杂系统中不同状态变量之间的相互作用和演化过程。本博客将围绕一个使用 MATLAB 实现的 SIP 模型展开，详细介绍模型的参数设置、平衡点计算、轨迹绘制以及全局相图分析。

二、模型介绍

SIP 模型主要涉及三个状态变量：S、I 和 P。在这个模型中，这些变量随时间的变化遵循一定的动力学规则。模型的参数设置如下：

r     = 2;
alpha = 1;
beta  = 1;
c     = 1;
d     = 0.25;
mu    = 0.4;    %  [0.4,0.8] 之间尝试不同值

• r：代表系统的某个基础增长率。
• alpha 和 beta：用于描述不同状态变量之间的相互作用系数。
• c 和 d：与状态变量 P 的变化相关。
• mu：一个关键参数，可在 [0.4, 0.8] 范围内取值，它会影响系统的动态行为。

三、平衡点计算

平衡点是系统在长时间演化后可能达到的稳定状态。通过对模型的分析，我们可以计算出三个平衡点：

S1 = r/alpha;
I1 = 0;
P1 = 0;

S2 = mu/beta;
I2 = - (alpha*mu - beta*r)/(beta*mu);
P2 = 0;

S3 = r*(c - d) / (c*alpha - alpha*d + beta*d);
I3 = d/(c - d);
P3 = - c*(c*alpha*mu - c*beta*r - alpha*d*mu + beta*d*mu + beta*d*r) ...
     / ( (c*alpha - alpha*d + beta*d)*(c - d) );

eqs = [S1,I1,P1; S2,I2,P2; S3,I3,P3];

这三个平衡点分别存储在 eqs 矩阵中，它们代表了系统可能的三种稳定状态。 

四、映射函数

映射函数 F 描述了状态变量 S、I 和 P 在每个时间步的更新规则：

F = @(S,I,P) deal( ...
    S + r - beta*S*I - alpha*S, ...
    I + beta*S*I - (P./(1+I)).*I - mu*I, ...
    P + P.*(-d + c*I./(1+I)) ...
);

通过这个函数，我们可以根据当前的状态变量值计算出下一个时间步的状态变量值。

五、平衡点附近的小范围轨迹

为了研究系统在平衡点附近的行为，我们在每个平衡点附近添加微小扰动，并观察系统的演化轨迹：

delta = 0.02;    % 微小扰动幅度
Nstep = 30;      % 步数

figure('Name','平衡点附近轨迹','Position',[100 100 900 300]);
for k = 1:3
    subplot(1,3,k);
    hold on; grid on;
    title(sprintf('Equilibrium #%d',k));
    xlabel('S'); ylabel('I'); zlabel('P');
    
    % 三条从三个微扰初始点开始的轨迹
    for sign = [-1,1]
        init = eqs(k,:) + sign*delta;
        traj = zeros(Nstep+1,3);
        traj(1,:) = init;
        S=init(1); I=init(2); P=init(3);
        for n=1:Nstep
            [Snew,Inew,Pnew] = F(S,I,P);
            traj(n+1,:) = [Snew,Inew,Pnew];
            S=Snew; I=Inew; P=Pnew;
        end
        plot3(traj(:,1),traj(:,2),traj(:,3),'-o','LineWidth',1);
    end
    plot3(eqs(k,1),eqs(k,2),eqs(k,3),'kp','MarkerSize',12,'MarkerFaceColor','y');
    view(3);
    axis tight;
end

通过绘制这些轨迹，我们可以直观地看到系统在平衡点附近的稳定性。如果轨迹逐渐靠近平衡点，说明该平衡点是稳定的；反之，如果轨迹远离平衡点，则说明该平衡点是不稳定的。

六、全局相图

全局相图展示了系统从给定初始值开始的整体演化过程：

S0 = 0.3; I0 = 0.2; P0 = 0.5;
T  = 100;
[y, t] = ode45(F, 1: T, [S0, I0, P0]);
traj = zeros(T+1,3);
traj(1,:) = [S0,I0,P0];
S=S0; I=I0; P=P0;
for n=1:T
    [Snew,Inew,Pnew] = F(S,I,P);
    traj(n+1,:) = [Snew,Inew,Pnew];
    S=Snew; I=Inew; P=Pnew;
end

figure('Name','全局相图','Position',[200 200 600 500]);
plot3(traj(:,1),traj(:,2),traj(:,3),'-', 'LineWidth',1.5);
hold on; grid on;
plot3(eqs(:,1),eqs(:,2),eqs(:,3),'kp','MarkerSize',10,'MarkerFaceColor','r');
legend('从 (0.3,0.2,0.5) 出发','平衡点','Location','best');
xlabel('S'); ylabel('I'); zlabel('P');
title('SIP 模型三维相图');
view(45,25);
axis tight;

全局相图可以帮助我们了解系统在不同初始条件下的长期行为，以及平衡点在整个系统中的作用。

七、结果解读

• 平衡点附近轨迹：通过观察平衡点附近的轨迹，我们可以判断每个平衡点的稳定性。如果轨迹最终收敛到平衡点，说明该平衡点是稳定的；如果轨迹发散，则说明该平衡点是不稳定的。这有助于我们理解系统在小扰动下的行为。
• 全局相图：全局相图展示了系统从给定初始值开始的整体演化过程。我们可以看到系统的轨迹如何在三维空间中移动，以及它是否会趋近于某个平衡点。这对于预测系统的长期行为非常有帮助。

通过对 SIP 模型的分析，我们可以深入了解系统的动态特性，为进一步的研究和应用提供理论基础。同时，通过调整模型参数（如 mu），我们可以研究不同参数对系统行为的影响，从而更好地控制和优化系统。

看下效果

 Over！