% Corrigé de la séance 11 -  Décembre 2011
%
%      METHODE D'EULER EXPLICITE
%
% Pour q=10, le graphe de la figure 2 de la "note" montre tros droites
% rouges qui correspondent à trois comportements différents possibles de
% l'évolution de la population (réduite !) de chenilles
% Pour r=0.3 il n'y a qu'une population d'équilibre, proche de u=0.32, elle
% correspond au "confinement". Quel que soit l'effectif de la population à
% l'instant initial, elle tendra inéluctablement vers le confinement.   
% Pour r=0.6 il n'y a qu'une population d'équilibre, proche de u=8, elle
% correspond au phénomène de  "pullulation". Quel que soit l'effectif de 
% la population à l'instant initial, elle tendra inéluctablement vers la
% pullulation.
% Pour r=0.45 il y a trois populations d'équilibre celle qui est proche de
% u=0.5 et celle qui est proche de u=6.8 sont stables, tandis que la
% population d'équilibre intermédiaire proche de u=2.7 est instable.
% Si l'effectif initial est inférieur à la valeur d'équilibre intermédiaire
% la population évolue vers le confinement, s'il est supérieur elle évolue
% vers la pullulation.
% On se limte aux trois valeurs de r : 0.3, 0.45 et 0.6
% Pour chaque valeur de r on simule l'évolution de la population pour
% quatre valeurs initiales, situées de part et d'autre de la position 
% d'équilibre intermédiaire correspondant à r=0.45
clc;clear
f   = @(u,r,q) (r*(1-u/q)-u/(1+u*u))*u;     % Second membre de l'EDO
q = 10;
to=0; tf=200;
N=500; 
t=linspace(to,tf,N+1);
u=zeros(length(t));
h=(tf-to)/N;
figure(1)
hold on
axis ([0 200 0 11])
uo=[0.2 1.5 2.8 10];
R=[0.3 0.45 0.6];
col=['b' 'r' 'g'];
for j=1:length(R)
    r=R(j);
    for i=1:length(uo)
        u(1)=uo(i);
        for n=1:N
            u(n+1)=u(n)+h*f(u(n),r,q);
        end
        plot(t,u,col(j))
        pause(0.2)
    end
    pause(1)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Corrigé de la séance 11 -  Décembre 2011
%
%      METHODE D'EULER IMPLICITE
%
% Pour q=10, le graphe de la figure 2 de la "note" montre tros droites
% rouges qui correspondent à trois comportements différents possibles de
% l'évolution de la population (réduite !) de chenilles
% Pour r=0.3 il n'y a qu'une population d'équilibre, proche de u=0.32, elle
% correspond au "confinement". Quel que soit l'effectif de la population à
% l'instant initial, elle tendra inéluctablement vers le confinement.   
% Pour r=0.6 il n'y a qu'une population d'équilibre, proche de u=8, elle
% correspond au phénomène de  "pullulation". Quel que soit l'effectif de 
% la population à l'instant initial, elle tendra inéluctablement vers la
% pullulation.
% Pour r=0.45 il y a trois populations d'équilibre celle qui est proche de
% u=0.5 et celle qui est proche de u=6.8 sont stables, tandis que la
% population d'équilibre intermédiaire proche de u=2.7 est instable.
% Si l'effectif initial est inférieur à la valeur d'équilibre intermédiaire
% la population évolue vers le confinement, s'il est supérieur elle évolue
% vers la pullulation.
% On se limte aux trois valeurs de r : 0.3, 0.45 et 0.6
% Pour chaque valeur de r on simule l'évolution de la population pour
% quatre valeurs initiales, situées de part et d'autre de la position 
% d'équilibre intermédiaire correspondant à r=0.45
clc;clear
f   = @(u,r,q) (r*(1-u/q)-u/(1+u*u))*u;     % Second membre de l'EDO
fpu = @(u,r,q) r*(1-2*u/q)-2*u/((1+u*u)^2); % Dérivée de f par rapport à u
                                            % pour la méthode de Newton
q = 10;
to=0; tf=200;
N=10000; 
t=linspace(to,tf,N+1);
u=zeros(length(t),1);
h=(tf-to)/N;
figure(1)
hold on
axis ([0 200 0 11])
uo=[0.22 1.5 2.8 10];
R=[0.3 0.45 0.6];
col=['b' 'r' 'g'];
for j=1:length(R)
    r=R(j);
    for i=1:length(uo)
        u(1)=uo(i);
% Valeur de démarrage calculée par Euler explicite
        u(2)=u(1)+h*f(u(1),r,q);
for n=2:N
    z=u(n);
    for k=1:4
        z=z-(z-u(n-1)-h/3*(f(u(n-1),r,q)+4*f(u(n),r,q)+f(z,r,q)))/...
            (1-h/3*fpu(z,r,q));
    end
    u(n+1)=z;
end
        plot(t,u,col(j))
        pause(0.2)
    end
    pause(1)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
%      METHODE DE HEUN (RUNGE-KUTTA A DEUX PARAMETRES)
%

clc;clear
f   = @(u,r,q) (r*(1-u/q)-u/(1+u*u))*u;     % Second membre de l'EDO
q = 10;
to=0; tf=200;
N=500; 
t=linspace(to,tf,N+1);
u=zeros(length(t));
h=(tf-to)/N;
figure(1)
hold on
axis ([0 200 0 11])
uo=[0.2 1.5 2.8 10];
R=[0.3 0.45 0.6];
col=['b' 'r' 'g'];
for j=1:length(R)
    r=R(j);
    for i=1:length(uo)
        u(1)=uo(i);
        for n=1:N
            k1=f(u(n),r,q);
            k2=f(u(n)+h*k1,r,q);
            u(n+1)=u(n)+h*(k1+k2)/2;
        end
        plot(t,u,col(j))
        pause(0.2)
    end
    pause(1)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        

%
%      METHODE DE RUNGE-KUTTA A 4 PARAMETRES
%

clc;clear
f   = @(u,r,q) (r*(1-u/q)-u/(1+u*u))*u;     % Second membre de l'EDO
q = 10;
to=0; tf=200;
N=500; 
t=linspace(to,tf,N+1);
u=zeros(length(t));
h=(tf-to)/N;
figure(1)
hold on
axis ([0 200 0 11])
uo=[0.2 1.5 2.8 10];
R=[0.3 0.45 0.6];
col=['b' 'r' 'g'];
for j=1:length(R)
    r=R(j);
    for i=1:length(uo)
        u(1)=uo(i);
        for n=1:N
           k1 = f(u(n),r,q);
           k2 = f(u(n)+h/2*k1,r,q); 
           k3 = f(u(n)+h/2*k2,r,q);
           k4 = f(u(n)+h*k3,r,q);
           u(n+1) = u(n) + h/6*(k1 + 2*(k2+k3)+k4);   
        end                          
        plot(t,u,col(j))
        pause(0.2)
    end
    pause(1)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
%
%      METHODE DE RUNGE-KUTTA A 4 PARAMETRES
%
%      RHO VARIE EN FONCTION DU TEMPS (REDUIT)
clc;clear
f   = @(t,u,q) ((0.45+0.20*sin(t/15))*(1-u/q)-u/(1+u*u))*u; 
q = 10;
to=0; tf=200;
N=1000; 
t=linspace(to,tf,N+1);
u=zeros(length(t));
h=(tf-to)/N;
figure(1)
hold on
axis ([0 200 0 11])
uo=[0.2 1.5 2.8 10];


    for i=1:length(uo)
        u(1)=uo(i);
        for n=1:N
           k1 = f(t(n),u(n),q);
           k2 = f(t(n)+h/2,u(n)+h/2*k1,q); 
           k3 = f(t(n)+h/2,u(n)+h/2*k2,q);
           k4 = f(t(n+1),u(n)+h*k3,q);
           u(n+1) = u(n) + h/6*(k1 + 2*(k2+k3)+k4);   
        end                          
        plot(t,u,'r')
        pause(0.2)
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        

%
%      METHODE DE RUNGE-KUTTA A 4 PARAMETRES
%
%      RHO ET Q VARIENT EN FONCTION DU TEMPS (REDUIT)
clc;clear
f   = @(t,u) ((0.45+0.20*sin(t/15))*(1-u/(10-2*cos(t/12)))-u/(1+u*u))*u; 
q = 10;
to=0; tf=500;
N=2000; 
t=linspace(to,tf,N+1);
u=zeros(length(t));
h=(tf-to)/N;
figure(1)
hold on
axis ([0 500 0 11])
uo=[0.2 1.5 2.8 10];


    for i=1:length(uo)
        u(1)=uo(i);
        for n=1:N
           k1 = f(t(n),u(n));
           k2 = f(t(n)+h/2,u(n)+h/2*k1); 
           k3 = f(t(n)+h/2,u(n)+h/2*k2);
           k4 = f(t(n+1),u(n)+h*k3);
           u(n+1) = u(n) + h/6*(k1 + 2*(k2+k3)+k4);   
        end                          
        plot(t,u,'r')
        pause(0.2)
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
 
        
          
        
                
                        
                
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation du modèle SIR de Kermack et McKendrick (1927).
% cf.  Mathematical biology - 
% Ed. Springer Verlag (1990) par J.D.Murray (page 613).
%
% Intégration par la méthode d'Euler explicite.
%
% Le taux de propagation de l'épidémie est appelé beta. 
% Le taux de guérison est noté gamma et 
% 1/gamma est la durée de l'infection.
%
% On peut discuter l'influence de beta et gamma.
%
% S : nombre d'individus sains (susceptibles);
% I : nombre de malades (infectives);
% R : nombre de guéris ou décédés (recovered).
%
% NB : gamma peut aussi inclure le taux de décès, il s'agit en fait du taux
% de passage de la catégorie I à la catégorie R (mort et/ou guéri)
%
clf;
figure(1)
hold on
axis([0 100 0 10000]); title('Modèle de Kermack et McKendrick',...
    'FontSize',14)
xlabel('Temps (jours)','FontSize',14);...
    ylabel('S I R  (nombre d''individus)','FontSize',14)
dt = 0.05;
Tmax=100;
beta=1/20000;
gamma=1/3;
I=1;                    % Un seul malade est introduit dans la population
R=0;                    % Au départ, personne n'est immunisé
S=10000-I;              % La population initiale est de 10000 individus
for t=0:dt:Tmax
% On définit deux variables auxiliaires    
    T1 = beta*S*I;
    T2 = gamma*I;
% Heun constantes k1
    k1S = -T1;
    k1I =  T1 - T2;
    k1R =  T2;
% Heun constantes k2
    k2S = -beta*(I+dt*k1I)*(S+dt*k1S);
    k2I = -k2S - gamma*(I+dt*k1I);
    k2R = gamma*(I+dt*k1I);
	S = S + dt/2*(k1S+k2S);
	I = I + dt/2*(k1I+k2I);
	R = R + dt/2*(k1R+k2R);
%
% Petite astuce, on utilise des complexes pour représenter les résultats
%
pointS = t + i*S;
pointR = t + i*R;
pointI = t + i*I;
plot(pointS,'g','LineWidth',1.5);  % Individus sains (susceptibles) vert
plot(pointI,'r','LineWidth',1.5);  % Individus malades (infectives) rouge
plot(pointR,'b','LineWidth',1.5);  % Individus guéris ou décédés (Removed) 
end
grid on
gtext('Sains','FontSize',14)
gtext('Infectés (malades contagieux)','FontSize',14)
gtext('Guéris','FontSize',14)