0

我试图通过在 Matlab 中编写代码来解决微分方程组。我在这个论坛上发帖，希望有人能以某种方式帮助我。我有一个由 10 个耦合微分方程组成的系统。它是一种载体-宿主流行病模型，它捕捉了一种疾病在人群和昆虫种群之间的传播。由于它是一个简单的微分方程系统，因此我将求解器 ( ode45) 用于非刚性问题类型。

有 10 个微分方程，每个代表 10 个不同的状态变量。有两个函数具有相同的 10 个耦合 ODE 系统。一个被称为NoEffects_derivative_6_15_2012.m包含 ODE 的原始系统。调用另一个函数，OnlyLethal_derivative_6_15_2012.m它包含相同的 ODE 系统，从时间开始，提款率增加gamma=32 %days，提款率随时间呈指数衰减。

我使用ode45相同的初始条件来解决这两个系统。两个系统的时间向量也是相同的，从t0到tfinal。该向量tspan包含从t0到的时间值tfinal，每个增量为 0.25 天，总共有 157 个时间值。

解值存储在矩阵ye0和yeL中。这两个矩阵都包含 157 行和 10 列（对于 10 个状态变量值）。当我比较第 10 个状态变量的值时，对于time=tfinal矩阵中的ye0，并yeL通过绘制差异，我发现它在某些时间值内变为负数。（使用命令：）plot(te0,ye0(:,10)-yeL(:,10))。这不是预期的。t0对于从until开始的所有时间值tfinal，状态变量 10 的值应该更大，因为它是从 ODE 系统获得的解决方案，它没有应用增加的提款率。

我被告知我的 matlab 代码中有一个错误。我不确定如何找出那个错误。或者我正在使用的 matlab 中的求解器 ( ode45) 效率不高，确实会出现这种问题。任何人都可以帮忙。

我也尝试过ode23，ode113但遇到了同样的问题。图 (2) 显示了一条曲线，该曲线对于时间值 32 和 34 变为负数，这表明了一个未预期的结果。对于所有时间值，该曲线应始终具有正值。还有其他任何人可以建议的论坛吗？

这是主要的脚本文件：

clear memory; clear all;
global Nc capitalambda muh lambdah del1 del2 p eta alpha1 alpha2 muv lambdav global dims Q t0 tfinal gamma Ct0 b1 b2 Ct0r b3 H C m_tilda betaHV bitesPERlanding IC global tspan Hs Cs betaVH k landingARRAY muARRAY
Nhh=33898857; Nvv=2*Nhh; Nc=21571585; g=354; % number of public health centers in Bihar state %Fix human parameters capitalambda= 1547.02; muh=0.000046142; lambdah= 0.07; del1=0.001331871263014; del2=0.000288658; p=0.24; eta=0.0083; alpha1=0.044; alpha2=0.0217; %Fix vector parameters muv=0.071428; % UNIT:2.13 SANDFLIES DEAD/SAND FLY/MONTH, SOURCE: MUBAYI ET AL., 2010 lambdav=0.05; % UNIT:1.5 TRANSMISSIONS/MONTH, SOURCE: MUBAYI ET AL., 2010
Ct0=0.054;b1=0.0260;b2=0.0610; Ct0r=0.63;b3=0.0130;
dimsH=6; % AS THERE ARE FIVE HUMAN COMPARTMENTS dimsV=3; % AS THERE ARE TWO VECTOR COMPARTMENTS dims=dimsH+dimsV; % THE TOTAL NUMBER OF COMPARTMENTS OR DIFFERENTIAL EQUATIONS
gamma=32; % spraying is done of 1st feb of the year
Q=0.2554; H=7933615; C=5392890;
m_tilda=100000; % assumed value 6.5, later I will have to get it for sand flies or mosquitoes betaHV=66.67/1000000; % estimated value from the short technical report sent by Anuj bitesPERlanding=lambdah/(m_tilda*betaHV); betaVH=lambdav/bitesPERlanding; IC=zeros(dims+1,1); % CREATES A MATRIX WITH DIMS+1 ROWS AND 1 COLUMN WITH ALL ELEMENTS AS ZEROES
t0=1; tfinal=40; for j=t0:1:(tfinal*4-4) tspan(1)= t0; tspan(j+1)= tspan(j)+0.25; end clear j;
% INITIAL CONDITION OF HUMAN COMPARTMENTS q1=0.8; q2=0.02; q3=0.0005; q4=0.0015; IC(1,1) = q1*Nhh; IC(2,1) = q2*Nhh; IC(3,1) = q3*Nhh; IC(4,1) = q4*Nhh; IC(5,1) = (1-q1-q2-q3-q4)*Nhh; IC(6,1) = Nhh; % INTIAL CONDITIONS OF THE VECTOR COMPARTMENTS IC(7,1) = 0.95*Nvv; %80 PERCENT OF TOTAL ARE ASSUMED AS SUSCEPTIBLE VECTORS IC(8,1) = 0.05*Nvv; %20 PRECENT OF TOTAL ARE ASSUMED AS INFECTED VECTORS IC(9,1) = Nvv; IC(10,1)=0;
Hs=2000000; Cs=3000000; k=1; landingARRAY=zeros(tfinal*50,2); muARRAY=zeros(tfinal*50,2);

[te0 ye0]=ode45(@NoEffects_derivative_6_15_2012,tspan,IC); [teL yeL]=ode45(@OnlyLethal_derivative_6_15_2012,tspan,IC);

figure(1) subplot(4,3,1); plot(te0,ye0(:,1),'b-',teL,yeL(:,1),'r-'); xlabel('time'); ylabel('S'); legend('susceptible humans'); subplot(4,3,2); plot(te0,ye0(:,2),'b-',teL,yeL(:,2),'r-'); xlabel('time'); ylabel('I'); legend('Infectious Cases'); subplot(4,3,3); plot(te0,ye0(:,3),'b-',teL,yeL(:,3),'r-'); xlabel('time'); ylabel('G'); legend('Cases in Govt. Clinics'); subplot(4,3,4); plot(te0,ye0(:,4),'b-',teL,yeL(:,4),'r-'); xlabel('time'); ylabel('T'); legend('Cases in Private Clinics'); subplot(4,3,5); plot(te0,ye0(:,5),'b-',teL,yeL(:,5),'r-'); xlabel('time'); ylabel('R'); legend('Recovered Cases');
subplot(4,3,6);plot(te0,ye0(:,6),'b-',teL,yeL(:,6),'r-'); hold on; plot(teL,capitalambda/muh); xlabel('time'); ylabel('Nh'); legend('Nh versus time');hold off;
subplot(4,3,7); plot(te0,ye0(:,7),'b-',teL,yeL(:,7),'r-'); xlabel('time'); ylabel('X'); legend('Susceptible Vectors');
subplot(4,3,8); plot(te0,ye0(:,8),'b-',teL,yeL(:,8),'r-'); xlabel('time'); ylabel('Z'); legend('Infected Vectors');
subplot(4,3,9); plot(te0,ye0(:,9),'b-',teL,yeL(:,9),'r-'); xlabel('time'); ylabel('Nv'); legend('Nv versus time');
subplot(4,3,10);plot(te0,ye0(:,10),'b-',teL,yeL(:,10),'r-'); xlabel('time'); ylabel('FS'); legend('Total number of human infections');
figure(2) plot(te0,ye0(:,10)-yeL(:,10)); xlabel('time'); ylabel('FS(without intervention)-FS(with lethal effect)'); legend('Diff. bet. VL cases with and w/o intervention:ode45');

函数文件：NoEffects_derivative_6_15_2012

function dx = NoEffects_derivative_6_15_2012( t , x )
global Nc capitalambda muh del1 del2 p eta alpha1 alpha2 muv global dims m_tilda betaHV bitesPERlanding betaVH
dx       = zeros(dims+1,1); % t % dx 
dx(1,1)  = capitalambda-(m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/(x(7,1)+x(8,1))-muh*x(1,1); 
dx(2,1)  = (m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/(x(7,1)+x(8,1))-(del1+eta+muh)*x(2,1); 
dx(3,1)  = p*eta*x(2,1)-(del2+alpha1+muh)*x(3,1); 
dx(4,1)  = (1-p)*eta*x(2,1)-(del2+alpha2+muh)*x(4,1);
dx(5,1)  = alpha1*x(3,1)+alpha2*x(4,1)-muh*x(5,1); 
dx(6,1)  = capitalambda -del1*x(2,1)-del2*x(3,1)-del2*x(4,1)-muh*x(6,1); 
dx(7,1)  = muv*(x(7,1)+x(8,1))-bitesPERlanding*betaVH*x(7,1)*x(2,1)/(x(6,1)+Nc)-muv*x(7,1);
%dx(8,1) = lambdav*x(7,1)*x(2,1)/(x(6,1)+Nc)-muvIOFt(t)*x(8,1); 
dx(8,1)  = bitesPERlanding*betaVH*x(7,1)*x(2,1)/(x(6,1)+Nc)-muv*x(8,1); 
dx(9,1)  = (muv-muv)*x(9,1); 
dx(10,1) = (m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/x(9,1);

函数文件：OnlyLethal_derivative_6_15_2012 

function dx=OnlyLethal_derivative_6_15_2012(t,x) 
global Nc capitalambda muh del1 del2 p eta alpha1 alpha2 muv global dims m_tilda betaHV bitesPERlanding betaVH k muARRAY
dx=zeros(dims+1,1);
% the below code saves some values into the second column of the two arrays % t muARRAY(k,1)=t; muARRAY(k,2)=artificialdeathrate1(t); k=k+1; 
dx(1,1)= capitalambda-(m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/(x(7,1)+x(8,1))-muh*x(1,1);    
dx(2,1)= (m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/(x(7,1)+x(8,1))-(del1+eta+muh)*x(2,1);
dx(3,1)=p*eta*x(2,1)-(del2+alpha1+muh)*x(3,1);
dx(4,1)=(1-p)*eta*x(2,1)-(del2+alpha2+muh)*x(4,1);
dx(5,1)=alpha1*x(3,1)+alpha2*x(4,1)-muh*x(5,1);
dx(6,1)=capitalambda -del1*x(2,1)-del2*( x(3,1)+x(4,1) ) - muh*x(6,1);
dx(7,1)=muv*( x(7,1)+x(8,1) )- bitesPERlanding*betaVH*x(7,1)*x(2,1)/(x(6,1)+Nc) - (artificialdeathrate1(t) + muv)*x(7,1);
dx(8,1)= bitesPERlanding*betaVH*x(7,1)*x(2,1)/(x(6,1)+Nc)-(artificialdeathrate1(t) + muv)*x(8,1);
dx(9,1)= -artificialdeathrate1(t) * x(9,1);
dx(10,1)= (m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/x(9,1);

函数文件：artificialdeathrate1

function art1=artificialdeathrate1(t) 
global Q Hs H Cs C   
art1= Q*Hs*iOFt(t)/H + (1-Q)*Cs*oOFt(t)/C ;

函数文件：iOFt 

function i = iOFt(t) 
global gamma tfinal Ct0 b1
if t>=gamma && t<=tfinal 
    i = Ct0*exp(-b1*(t-gamma)); 
else 
    i =0; 
end

函数文件：oOFt

function o = oOFt(t)
global gamma Ct0 b2 tfinal    
if (t>=gamma && t<=tfinal) 
    o = Ct0*exp(-b2*(t-gamma)); 
else 
    o = 0; 
end
4

1 回答 1

2

如果您的工作代码甚至与您发布的代码一样混乱，那么恕我直言，您应该首先解决的问题。

我为你清理了一点iOFt，oOFt因为这些很容易处理。我尽力了NoEffects_derivative_6_15_2012。我个人对您的代码所做的更改是使用体面的索引。你有 10 个变量，如果你让你的代码休息几个星期或几个月，你就无法记住例如状态 7 是什么。因此(7,1)，您可能不想使用 ，而是使用详细名称重写 ODE，然后将它们检索/存储在x和dx向量中。或者使用索引来明确正在发生的事情。

例如

function ODE(t,x)
insectsInfected = x(1);
humansInfected  = x(2);
%etc

dInsectsInfected = %some function of the rest
dHumansInfected  = %some function of the rest
% etc

dx = [dInsectsInfected; dHumansInfected; ...];

或者

function ODE(t,x)
iInsectsInfected = 1;
iHumansInfected  = 2;
%etc

dx(iInsectsInfected) = %some function of x(i...)
dx(iHumansInfected)  = %some function of x(i...)
%etc

当你不做这些事情时，你最终可能会在某些公式中使用x(6,1)而不是例如x(3,1)，并且可能需要几个小时才能发现这样的事情。如果您使用冗长的名称，则键入时间会更长，但它会使调试变得容易得多，并且如果您理解您的方程式，那么当发生此类错误时应该会更加明显。

此外，不要犹豫，在公式中添加空格，这会使阅读更容易。如果您有一些有意义的子表达式（例如，如果(1-p)*eta*x(2,1)是死于疾病的昆虫的数量，只需将其放入变量中dyingInsects并在它出现的任何地方使用）。如果你对齐你的作业（就像我在上面所做的那样），这可能会增加更容易阅读和理解的代码。

关于 ODE 求解器，如果你确定你的实现是正确的，我也会尝试一个求解器来解决僵硬问题（除非你绝对确定你没有僵硬的系统）。

于 2012-06-24T07:27:07.983 回答