0

Matlab でコードを記述して微分方程式系を解こうとしています。誰かが私を助けてくれることを願って、私はこのフォーラムに投稿しています。10 個の結合微分方程式のシステムがあります。これは、人間の集団と昆虫の集団の間での病気の伝染を捉えたベクター宿主流行モデルです。これは微分方程式の単純なシステムであるode45ため、非スティッフな問題タイプにはソルバー ( ) を使用しています。

10 個の微分方程式があり、それぞれが 10 個の異なる状態変数を表しています。10 連成 ODE の同じシステムを持つ関数が 2 つあります。NoEffects_derivative_6_15_2012.mODE の元のシステムを含む1 つが呼び出されます。もう 1 つの関数が呼び出さOnlyLethal_derivative_6_15_2012.mれます。この関数には、同じ ODE のシステムが含まれており、時間から始まる引き出し率が増加しgamma=32 %days、その引き出し率は時間とともに指数関数的に減衰します。

ode45同じ初期条件を使用して、両方のシステムを解決するために使用します。時間ベクトルも両方のシステムで同じで、 から に進みt0ますtfinal。このベクトルには からまでtspanの時間値が含まれており、それぞれが 0.25 日の増分で、合計 157 の時間値になります。t0tfinal

解の値は、行列ye0およびに格納されyeLます。これらの行列はどちらも、157 行と 10 列 (10 個の状態変数値) を含みます。10 番目の状態変数の値をtime=tfinal行列ye0で比較しyeL、差をプロットすると、いくつかの時間値で負になっていることがわかります。(コマンドを使用: plot(te0,ye0(:,10)-yeL(:,10)))。これは想定外です。t0からまでのすべての時間値について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 に答える