用 MatLab 求解这个微分方程组

Question

我正在尝试用 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

Answer 1

这里的问题是，您仅在所有 x 值都相等时才进行测试。查看g 中包含d1 的第一个词

d1*(m1*exp(-x(4)/m1)*x(3) - m1*exp(-x(2)/m1)*x(1))

如果所有x 值都相等，则该等式简化为d1*0。对于函数g 中带有d1 或d2 的所有其他术语也是如此。因此，更改这些值对结果没有影响。如果您的 x 值不同，则 d1 和 d2 值将影响结果。