使用ode45獲取意外結果

Question

我想通過在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個狀態變量值）。當我通過繪製差異來比較time=tfinal，矩陣ye0和yeL中的第10個狀態變量的值時，我發現它對某些時間值變爲負值。 （使用命令：plot(te0,ye0(:,10)-yeL(:,10))）。這不是預期的。對於從t0到tfinal的所有時間值，10狀態變量的值應該更大，因爲它是從ODE系統獲得的解決方案，該系統沒有對其應用增加的提取率。

我被告知我的matlab代碼中有一個錯誤。我不知道如何找出那個bug。或者，也許我正在使用的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 

Answer 1

如果你的工作代碼甚至是遠程爲您發佈的代碼凌亂，那麼就應該恕我直言你應該首先解決的問題。

我清理了iOFt,oOFt有點適合你，因爲這些都很容易處理。我盡我所能在NoEffects_derivative_6_15_2012。我親自改變你的代碼是使用像樣的索引。你有10個變量，如果你讓你的代碼休息幾周或幾個月，你就不會記住，例如，你會記得狀態7。因此，您可能不想使用(7,1)，而是使用詳細名稱重寫ODE，然後將它們檢索/存儲在x和dx向量中。或者使用清楚發生什麼的索引。

E.g.

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求解器，如果你確定你的實現是正確的，我也會嘗試解決僵硬問題（除非你確定你沒有僵硬的系統）。