ode45를 사용하는 중에 예기치 않은 결과가 발생했습니다.

Question

Matlab에 코드를 작성하여 미분 방정식 시스템을 풀려고합니다. 누군가가 어떤 식 으로든 나를 도울 수 있기를 희망하면서이 포럼에 게시하고 있습니다. 10 개의 미분 방정식 시스템이 있습니다. 그것은 인류와 곤충 인구 사이에서 질병의 전파를 포착하는 벡터 숙주 전염병 모델입니다. 그것은 미분 방정식의 간단한 시스템이기 때문에, 나는 딱딱하지 않은 문제 유형에 대해 솔버 (ode45)를 사용하고 있습니다.

각각 10 개의 다른 상태 변수를 나타내는 10 개의 미분 방정식이 있습니다. 동일한 10 개의 ODE가 결합 된 두 가지 기능이 있습니다. 하나는 ODE의 원래 시스템을 포함하는 NoEffects_derivative_6_15_2012.m입니다. 다른 기능은 이며, ODE와 동일한 시스템을 포함하고 있으며, 시간당 증가하는 철수 율, gamma=32 %days 및 철수 율은 시간에 따라 기하 급수적으로 감소합니다.

동일한 초기 조건을 사용하여 두 시스템을 모두 해결하기 위해 ode45을 사용합니다. 시간 벡터는 두 시스템에서 동일하며, t0에서 tfinal으로 이동합니다. 벡터 tspan에는 t0에서 tfinal까지의 시간 값이 포함되며 각 시간 값은 0.25 일 증가하여 총 157 개의 시간 값을 생성합니다.

솔루션 값은 행렬 ye0 및 yeL에 저장됩니다. 이 행렬은 모두 157 개의 행과 10 개의 열 (10 개의 상태 변수 값)을 포함합니다. 차이를 플로팅하여 ye0 및 yeL 행의 time=tfinal에 대한 10 번째 상태 변수의 값을 비교할 때 일부 시간 값이 음수가된다는 것을 알았습니다. (명령 : plot(te0,ye0(:,10)-yeL(:,10)) 사용). 이것은 예상되지 않습니다. t0에서 tfinal까지의 모든 시간 값에 대해, 10 상태 변수의 값은 증가 된 철회 율을 가지지 않은 ODE 시스템에서 얻은 해이기 때문에 더 커야합니다.

내 MATLAB 코드에 버그가 있다고 들었습니다. 그 버그를 찾는 방법을 모르겠습니다. 아니면 (ode45) matlab에있는 솔버가 효율적이지 않고 이런 종류의 문제를 제공합니다. 누구든지 도와 줄 수 있어요.

나는 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

다음

주 스크립트 파일입니다 0

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 벡터로 검색/저장하는 것이 좋습니다. 또는 무슨 일이 일어나고 있는지 명확하게하는 인덱스를 사용하십시오.

예.

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 해석기와 관련하여 구현이 올바른지 확실하다면, 딱딱한 시스템이 없다고 확신하지 않는 한 딱딱한 문제에 대한 해결 방법을 시도해 볼 수도 있습니다.