我正在尝试用 Matlab 求解以下微分方程组。
我使用的参数在这里
.
我遇到的问题是，当我运行我的代码时，无论我为 d1 和 d2 设置的值如何，我都会得到相同的结果。准确地说，我得到的结果与系统的简化版本相同，其中 d1,d2 = 0。由于我试图复制其他人的结果，我知道情况不应该如此。有谁知道为什么会这样？
% Define u_1(t,j), u_2(t,i) as:
% u_1(t,1) = x(1)
% u_2(t,1) = x(2)
% u_1(t,2) = x(3)
% u_2(t,2) = x(4)
% So
% x(1)' = x(1)*(epsilon - (epsilon/k1)*x(1) - alpha*x(2)) + d1*(rho1(x(4))*x(3) - rho1(x(2))*x(1))
% x(2)' = x(2)*(gama + beta*x(1) - (gama/k2)*x(2)) + d2*(rho2(x(3))*x(4) - rho2(x(1))*x(2))
% x(3)' = x(3)*(epsilon - (epsilon/k1)*x(3) - alpha*x(4)) + d1*(rho1(x(2))*x(1) - rho1(x(4))*x(3))
% x(4)' = x(4)*(gama + beta*x(3) - (gama/k2)*x(4)) + d2*(rho2(x(1))*x(2) - rho2(x(3))*x(4))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsilon = 3;
alpha = 0.2;
gama = 0.4;
beta = 0.4;
k1 = 2;
k2 = 1;
m1 = 1;
m2 = 1;
d2 = 0.1;
d1 = 17;
% rho1(x(2)) = m1*exp(-x(2)/m1)
% rho1(x(4)) = m1*exp(-x(4)/m1)
% rho2(x(1)) = m2*exp(-x(1)/m2)
% rho2(x(3)) = m2*exp(-x(3)/m2)
g = @(t,x)[
x(1)*(epsilon - (epsilon/k1)*x(1) - alpha*x(2)) + d1*(m1*exp(-x(4)/m1)*x(3) - m1*exp(-x(2)/m1)*x(1));
x(2)*(gama + beta*x(1) - (gama/k2)*x(2)) + d2*(m2*exp(-x(3)/m2)*x(4) - m2*exp(-x(1)/m2)*x(2));
x(3)*(epsilon - (epsilon/k1)*x(3) - alpha*x(4)) + d1*(m1*exp(-x(2)/m1)*x(1) - m1*exp(-x(4)/m1)*x(3));
x(4)*(gama + beta*x(3) - (gama/k2)*x(4)) + d2*(m2*exp(-x(1)/m2)*x(2) - m2*exp(-x(3)/m2)*x(4))];
[t1,xa1] = ode45(@(t,x) g(t,x),[0 20],[4 4 4 4]);
[t2,xa2] = ode45(@(t,x) g(t,x),[0 20],[3 3 3 3]);
[t3,xa3] = ode45(@(t,x) g(t,x),[0 20],[2 2 2 2]);
[t4,xa4] = ode45(@(t,x) g(t,x),[0 20],[0.5 0.5 0.5 0.5]);
figure;
subplot(2,2,1)
plot(t1,xa1(:,1),"r",t1,xa1(:,2),"b")
title('y0 = 4')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','southwest')
subplot(2,2,2)
plot(t2,xa2(:,1),"r",t2,xa2(:,2),"b")
title('y0 = 3')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','southwest')
subplot(2,2,3)
plot(t3,xa3(:,1),"r",t3,xa3(:,2),"b")
title('y0 = 2')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','southwest')
subplot(2,2,4)
plot(t4,xa4(:,1),"r",t4,xa4(:,2),"b")
title('y0 = 0.5')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','southwest')
sgtitle('Labour Force and Free Jobs')
figure;
hold on
plot(t1,xa1(:,1),"r",t1,xa1(:,2),"b")
plot(t2,xa2(:,1),"r",t2,xa2(:,2),"b")
plot(t3,xa3(:,1),"r",t3,xa3(:,2),"b")
plot(t4,xa4(:,1),"r",t4,xa4(:,2),"b")
title('Labour Force and Free Jobs')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','northeast')
hold off