Ergebnis für URL: http://www.ummisco.ird.fr/perso/bacaer/chikungunya.sciclear;
//chikungunya
gamm=7/7; //incubation moustique
mu=7/30; //par semaine
delta=7/4; //incubation homme
alpha=7/7; //par semaine
bet=7/4; //par semaine
phi=2*%pi/12;
omega=2*%pi/52;
P=785000;
r=0.06;
epsi=(1-r)/(1+r);
pmax=1.2*P/bet; //par semaine
pmin=r*pmax;
p0=0.5*(pmin+pmax);
T=2*%pi/omega;
//order 2
//disp("R0 order 2 = "+string(a*b*(1-alpha*mu*e^2/2/(omega^2+(alpha+mu)^2))/alpha/mu));
//dichotomie
deff("[dqdt]=floquet(t,q)",["e=q(1)";
"i=q(2)";...
"E=q(3)";...
"I=q(4)";...
"dqdt(1)=-(gamm+mu)*e+(bet*p0/R0/P)*(1+epsi*cos(omega*t-phi))*I";...
"dqdt(2)=gamm*e-mu*i";...
"dqdt(3)=bet*i-delta*E";...
"dqdt(4)=delta*E-alpha*I"]);
R0min=0.5;
R0max=100;
while R0max-R0min>1e-7,
R0=(R0min+R0max)/2;
c1=ode([1;0;0;0],0,[0 T],floquet);
colonne1=c1(:,2);
c2=ode([0;1;0;0],0,[0 T],floquet);
colonne2=c2(:,2);
c3=ode([0;0;1;0],0,[0 T],floquet);
colonne3=c3(:,2);
c4=ode([0;0;0;1],0,[0 T],floquet);
colonne4=c4(:,2);
sp=spec([colonne1 colonne2 colonne3 colonne4]);
rho=max(real(sp));
if rho>1 then
R0min=R0;
else
R0max=R0;
end;
end;
disp("R0 dichotomy = "+string(R0));
//eigenvalue
if 0==0 then
N=100;
matrice=zeros(N,N);
dt=T/N;
deff("[y]=Ghat(x)","y=(bet^2*p0*gamm*delta/P)*(...
exp(-alpha*x)/(1-exp(-alpha*T))/(gamm+mu-alpha)/(mu-alpha)/(delta-alpha)...
+exp(-mu*x)/(1-exp(-mu*T))/(delta-mu)/(alpha-mu)/gamm...
+exp(-delta*x)/(1-exp(-delta*T))/(alpha-delta)/(gamm+mu-delta)/(mu-delta)...
+exp(-(mu+gamm)*x)/(1-exp(-(mu+gamm)*T))/(-gamm)/(delta-mu-gamm)/(alpha-mu-gamm))");
for i=1:N,
ti=(i-1)*dt;
f(i)=1+epsi*cos(omega*ti-phi);
end;
for i=1:N,
ti=(i-1)*dt;
for j=1:i-1,
tj=(j-1)*dt;
matrice(i,j)=f(i)*Ghat(ti-tj);
end;
for j=i:N,
tj=(j-1)*dt;
matrice(i,j)=f(i)*Ghat(ti-tj+T);
end;
end;
matrice=matrice*dt;
sp2=spec(matrice);
R0=max(real(sp2));
disp("R0 eigenvalue = "+string(R0));
end;
//Fourier
if 0==0 then
deff("[y]=G(n)","y=(bet^2*p0*gamm*delta/P)/(mu+n*%i*omega)/(alpha+n*%i*omega)/(delta+n*%i*
omega)/(gamm+mu+n*%i*omega)");
kappa=12;
Lambda=zeros(kappa+1,1);
c=zeros(kappa+2,kappa+1);
Lambda(1)=G(0);
c(1,1)=1;
for k=1:kappa,
Lambda(k+1)=real(G(1)*c(2,k));
for n=1:k,
c(n+1,k+1)=(0.5*G(n-1)*c(n,k)+0.5*G(n+1)*c(n+2,k)-c(n+1,k:-1:2)*Lambda(2:k))/(G(0)-G(n
));
end;
end;
R0=(epsi.^[0:kappa])*Lambda(1:kappa+1);
disp("R0 Fourier = "+string(R0));
end;
//Dye
if 0==0 then
N=4;
M=zeros(2*N+1,2*N+1);
deff("[y]=f(n)","y=1*(n==0)+(epsi/2)*(n==1)+(epsi/2)*(n==-1)");
deff("[y]=G(n)","y=(bet^2*p0*gamm*delta/P)/(mu+n*%i*omega)/(alpha+n*%i*omega)/(delta+n*%i*
omega)/(gamm+mu+n*%i*omega)");
fourier=zeros(2*N+1,2*N+1);
for n=-N:N,
fourier(N+n+1,N+n+1)=f(0)*G(n);
if n>-N then
fourier(N+n+1,N+n)=f(1)*G(n-1);
end;
if n
Usage: http://www.kk-software.de/kklynxview/get/URL
e.g. http://www.kk-software.de/kklynxview/get/http://www.kk-software.de
Errormessages are in German, sorry ;-)