// Etude du modèle SIR et de méthodes numériques pour les équations différerentielles

// Paramètres scientifiques

R0 = 3          // Taux de reproduction de base du virus
Tr = 10         // Temps de rétablissement (en jours)
N = 100000      // Nombre d'individus
t_immu = 0.1    // Taux d'immunité initiale: Proportion d'individus initialement immunisées au début de la simulation
I0 = 3          // Nombre d'individus malades au début de la simulation

// Paramètres mathématiques

T = 150         // Durée de la simulation (en jours)
h = 5           // Pas de temps pour l'intégration numérique


// Condition initiale

X0  = [N-I0-t_immu*N;I0;t_immu*N]


// Système différentiel

function Xdot = F(t,X)
    // Fonction donnant le champ de vecteurs associé à l'équation différentielle  X' = F(t,X)
    // à l'aide du modèle SIR
    // Paramètres d'entrée:
    //  - t: Temps
    //  - X: Variable d'espace (vecteur de format 3 x 1): X = [S;I;R], où:
    //      > S: Nombre d'individus sains
    //      > I: Nombre d'individus infectés
    //      > R: Nombre d'individus infectés
    S = X(1,1)
    I = X(2,1)
    R = X(3,1)
    Xdot = [-(R0/(Tr*N))*S*I;(1/Tr)*I*((R0/N)*S-1);(1/Tr)*I]
    
endfunction



// Intégration du système différentiel

function X = Exact(X0,T,h)
    // Fonction donnant une bonne approximation de la solution exacte à l'aide de l'intégrateur de Scilab
    // Paramètres d'entrée:
    //  - X0: Condition initiale (vecteur de format 3 x 1)
    //  - T: Durée de la simulation
    //  - h: Pas de simulation
    X = ode(X0,0,0:h:T,F)
endfunction


function X = Euler(X0,T,h)
    // Fonction donnant une approximation de la solution à l'aide de la méthode d'Euler Explicite
    // Paramètres d'entrée:
    //  - X0: Condition initiale (vecteur de format 3 x 1)
    //  - T: Durée de la simulation
    //  - h: Pas de simulation
    X = X0
    x = X0
    for t = h:h:T
        x = x + h*F(t,x);
        X = [X,x];
    end
endfunction


function X = RK2(X0,T,h)
    // Fonction donnant une approximation de la solution à l'aide de la méthode de Runge-Kutta 2
    // Paramètres d'entrée:
    //  - X0: Condition initiale (vecteur de format 3 x 1)
    //  - T: Durée de la simulation
    //  - h: Pas de simulation
    X = X0
    x = X0
    for t = h:h:T
        k = x + (h/2)*F(t,x)
        x = x + h*F(t,k);
        X = [X,x];
    end
endfunction


// Tracé des solutions

X_Exact = Exact(X0,T,h)
X_Euler = Euler(X0,T,h)
X_RK2 = RK2(X0,T,h)
TT = 0:h:T

scf();
plot(0:1:T,Exact(X0,T,1)(1,:),'color','blue','thickness',2)
plot(0:1:T,Exact(X0,T,1)(2,:),'color','red','thickness',2)
plot(0:1:T,Exact(X0,T,1)(3,:),'color',color(0,210,0),'thickness',2)
legend('S','I','R')
title('Solution du système différentiel SIR')
xlabel("Temps (jours)")
ylabel("Population")
gca().grid=[1 1 1]


scf();
plot(TT,X_Exact(1,:),'color',color(0,0,250),'thickness',2)
plot(TT,X_Exact(2,:),'color',color(250,0,0),'thickness',2)
plot(TT,X_Exact(3,:),'color',color(0,200,0),'thickness',2)
plot(TT,X_RK2(1,:),'color',color(0,150,250),'thickness',2)
plot(TT,X_RK2(2,:),'color',color(250,150,0),'thickness',2)
plot(TT,X_RK2(3,:),'color',color(0,235,0),'thickness',2)
plot(TT,X_Euler(1,:),'color',color(0,250,250),'thickness',2)
plot(TT,X_Euler(2,:),'color',color(250,250,0),'thickness',2)
plot(TT,X_Euler(3,:),'color',color(0,250,0),'thickness',2)
legend('S (ode)','I (ode)','R (ode)','S (RK2)','I (RK2)','R (RK2)','S (Euler)','I (Euler)','R (Euler)')
title('Solutions du système différentiel SIR selon les différentes méthodes numériques: Intégrateur Scilab (ode), Euler et RK2')
xlabel("Temps (jours)")
ylabel("Population")
gca().grid=[1 1 1]


scf();
deff("[xdot] = fSIR(t,x)",["xd1 = -(R0/(Tr*N))*x(1)*x(2)";"xd2 = (1/Tr)*x(2)*((R0/N)*x(1)-1)";"xdot = [ xd1 ; xd2 ]"])
MinMax = [min(Exact(X0,T,1)(1,:)),max(Exact(X0,T,1)(1,:)),min(Exact(X0,T,1)(2,:)),max(Exact(X0,T,1)(2,:))]
minS = MinMax(1)
maxS = MinMax(2)
minI = MinMax(3)
maxI = MinMax(4)
xf= minS:(maxS-minS)/20:maxS;
yf= minI:(maxI-minI)/20:maxI;
fchamp(fSIR,0,xf,yf)
plot(Exact(X0,T,1)(1,:),Exact(X0,T,1)(2,:),'color','red','thickness',2)
title('Portrait de phase (S,I) de la solution')
xlabel("S")
ylabel("I")
gca().grid=[1 1 1]




//Ordre de convergence pour les méthodes numériuques

function e = Erreur_Euler(X0,T,h)
    // Fonction donnant l'erreur globale entre la vraie solution (approchée par ode) et la solution
    // calculée numériquement avec la méthde d'Euler Explicite
    //  - X0: Condition initiale (vecteur de format 3 x 1)
    //  - T: Durée de la simulation
    //  - h: Pas de simulation
    X_Exact = Exact(X0,T,h)
    X_Euler = Euler(X0,T,h)
    err = X_Exact - X_Euler
    err_loc = sqrt(err(1,:)**2 + err(2,:)**2 + err(3,:)**2)
    e = norm(err_loc,'inf')
endfunction


function e = Erreur_RK2(X0,T,h)
    // Fonction donnant l'erreur globale entre la vraie solution (approchée par ode) et la solution
    // calculée numériquement avec la méthde de Runge-Kutta 2
    // Paramètres d'entrée:
    //  - X0: Condition initiale (vecteur de format 3 x 1)
    //  - T: Durée de la simulation
    //  - h: Pas de simulation
    X_Exact = Exact(X0,T,h)
    X_Euler = RK2(X0,T,h)
    err = X_Exact - X_Euler
    err_loc = sqrt(err(1,:)**2 + err(2,:)**2 + err(3,:)**2)
    e = norm(err_loc,'inf')
endfunction


function [ord,ord2,HH,Err] = Ord_Conv_Euler(X0,T,h,n)
    // Fonction donnant l'ordre de convergence de la méthode d'Euler Explicite, ainsi que la liste des pas de temps choisis et les erreurs correspondantes (globales)
    // Paramètres d'entrée:
    //  - X0: Condition initiale (vecteur de format 3 x 1)
    //  - T: Durée de la simulation
    //  - h: Pas de simulation le plus grand pour lequel on calcule l'erreur globale
    //  - n: Nombre de pas selectionnés, on va calculer l'erreur globale pour les pas h/2^i, avec i allant de 0 à N
    HH = h./2**(0:n)
    Err = []
    for h = HH
        Err = [Err,Erreur_Euler(X0,T,h)/N] // Se rapproche d'une erreur relative
    end
    [ord,ord2] = reglin(log(HH),log(Err))
endfunction


function [ord,ord2,HH,Err] = Ord_Conv_RK2(X0,T,h,n)
    // Fonction donnant l'ordre de convergence de la méthode de Runge-Kutta 2, ainsi que la liste des pas de temps choisis et les erreurs correspondantes (globales)
    // Paramètres d'entrée:
    //  - X0: Condition initiale (vecteur de format 3 x 1)
    //  - T: Durée de la simulation
    //  - h: Pas de simulation le plus grand pour lequel on calcule l'erreur globale
    //  - n: Nombre de pas selectionnés, on va calculer l'erreur globale pour les pas h/2^i, avec i allant de 0 à N
    HH = h./2**(0:n)
    Err = []
    for h = HH
        Err = [Err,Erreur_RK2(X0,T,h)/N] // Se rapproche d'une erreur relative
    end
    [ord,ord2] = reglin(log(HH),log(Err))
endfunction


[ord_Euler,ord2_Euler,HH_Euler,Err_Euler] = Ord_Conv_Euler(X0,T,h,10)
[ord_RK2,ord2_RK2,HH_RK2,Err_RK2] = Ord_Conv_RK2(X0,T,h,10)

scf();
plot("ll",HH_Euler,Err_Euler,"rs",'thickness',2)
plot("ll",HH_RK2,Err_RK2,"gs",'thickness',2)
plot("ll",HH_Euler,exp(ord2_Euler)*HH_Euler.^ord_Euler,'color',color(150,0,0),'thickness',2,'LineStyle','--')
plot("ll",HH_RK2,exp(ord2_RK2)*HH_RK2.^ord_RK2,'color',color(0,150,0),'thickness',2,'LineStyle','--')
xlabel("h")
ylabel("Erreur relative")
title("Ordres de convergence observés pour les méthodes d Euler et de Runge-Kutta 2")
legend("Erreurs (Euler)","Erreurs (RK2)","Ordre observé (Euler): "+string(ord_Euler),"Ordre observé (RK2): "+string(ord_RK2),'pos',"in_lower_right")
gca().grid=[1 1 1]









