Merge branch 'master' of ssh://ssh.hesge.ch:10572/orestis.malaspin/math_tech_info

Makefile

@@ -37,6 +37,10 @@ deploy: all
 	make -C tpIntegrales
 	cp tpIntegrales/*.pdf mti/tpIntegrales/
 	cp tpIntegrales/tp_integrales_conv.html mti/tpIntegrales/index.html
+	mkdir -p mti/tpEdo
+	make -C tpEdo
+	cp tpEdo/*.pdf mti/tpEdo/
+	cp tpEdo/tpEquadiffs.html mti/tpEdo/index.html
 
 clean:
 	rm -rf *.html *.pdf

covid/python/covid.py

@@ -2,7 +2,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 swiss = np.array([2.100350058, 3.150525088, 4.900816803, 6.534422404, 10.501750292, 13.302217036, 24.970828471, 31.271878646, 39.323220537, 43.640606768, 57.292882147, 76.079346558, 100.116686114, 131.271878646, 158.576429405, 256.709451575, 256.709451575, 268.378063011, 309.218203034, 353.325554259, 453.675612602])
-swiss = np.array([18, 27, 42, 56, 90, 114, 214, 268, 337, 374, 491, 652, 858, 1125, 1359, 2200, 2200, 2300, 2650, 3028, 3888])
+swiss = np.array([18, 27, 42, 56, 90, 114, 214, 268, 337, 374, 491, 652, 858, 1125, 1359, 2200, 2200, 2300, 2650, 3028, 3888]) / 0.02
 
 days = np.array(range(1,len(swiss)+1))
 
@@ -27,18 +27,21 @@ def timestep(S0, I0, R0, dt, beta, lamb, N):
 S0 = 8000000
 I0 = swiss[0]
 R0 = 0
-max_t = 50*days[len(swiss)-1]
-n_steps = 10000
+max_t = 3*days[len(swiss)-1]
+n_steps = 100000
 dt = max_t / n_steps
 
 
 N = compute_n(S0, I0, R0) 
 
 lamb = 1.0 / 14.0
-beta_1 = 0.34
-beta_2 = beta_1 / 100
+beta_1 = 0.35
+beta_2 = beta_1 / 10
 beta = beta_1
 
+R0 = beta / lamb
+print(R0)
+
 s_list = [S0]
 r_list = [R0]
 i_list = [I0]
@@ -49,20 +52,21 @@ for i in range(0, n_steps):
     i_list.append(I1)
     r_list.append(R1)
     t_list.append(t_list[i]+dt)
-    if (t_list[i+1] >= 27 and t_list[i+1] <= 127) :
-        beta = beta_2
-    else:
-        beta = beta_1
+    # if (R1 >= 0.9 * N):
+    #     print(t_list[i]+dt)
+    #     break
+    # if (t_list[i+1] < 27) :
+    #     beta = beta_1
+    # elif (I1 > 0.01*N) :
+    #     beta = beta_2
     # else:
-    #     beta = 1
+    #     beta = beta_1
 
 s = np.array(s_list) 
 r = np.array(r_list)
 ii = np.array(i_list)
 t = np.array(t_list)
-
-print(t)
-
+s
 plt.semilogy(t, s, 'b')
 plt.semilogy(t, r, 'r')
 plt.semilogy(t, ii, 'k')

tpEdo/Makefile

+STYLES := ../css/tufte-css/tufte.css \
+	../css/pandoc.css \
+	../css/pandoc-solarized.css \
+	../css/tufte-extra.css
+
+OPTIONS = --filter=pandoc-numbering
+OPTIONS += --filter=pandoc-crossref
+
+PDFOPTIONS = --pdf-engine pdflatex
+PDFOPTIONS += --highlight-style kate
+PDFOPTIONS += --number-sections
+PDFOPTIONS += --template=./default.latex
+
+
+HTMLOPTIONS += -t html5
+HTMLOPTIONS += -c ../css/styling.css
+HTMLOPTIONS += --self-contained
+HTMLOPTIONS += --mathjax=../MathJax.js
+
+MD=$(wildcard *.md)
+HTML=$(MD:%.md=%.html)
+PDF=$(MD:%.md=%.pdf)
+
+
+all: $(HTML) $(PDF)
+
+%.pdf: %.md Makefile
+	pandoc -s $(OPTIONS) $(PDFOPTIONS) -o $@ $<
+
+%.html: %.md Makefile
+	pandoc -s $(OPTIONS) $(HTMLOPTIONS) -o $@ $<
+
+clean:
+	rm -rf *.html *.pdf

tpEdo/tpEquadiffs.md

+---
+# author:
+# - Orestis Malaspinas
+title: Travail pratique sur les équations différentielles
+autoSectionLabels: true
+autoEqnLabels: false
+eqnPrefix: 
+  - "éq."
+  - "éqs."
+chapters: false
+numberSections: false
+chaptersDepth: 1
+sectionsDepth: 3
+lang: fr
+documentclass: article
+papersize: A4
+cref: false
+pandoc-numbering:
+  - category: exercice
+urlcolor: blue
+---
+\newcommand{\real}{\mathbb{R}}
+
+# Rappel théorique : approximation numérique d'équations différentielles ordinaires 
+
+En général il n'existe pas de solutions analytiques pour les équations
+différentielles d'un intérêt pratique en ingénierie, c'est pourquoi on
+est obligé d'utiliser des méthodes numériques pour approximer la
+solution. Nous présentons ici une famille de méthodes numériques : les
+méthodes de Runge-Kutta qui permettent d'approximer les solutions
+d'équations différentielles ordinaires
+$$
+\frac{\mathrm{d}y}{\mathrm{d}t} = F(y,t), \quad y(t_0)=y_0,
+$$
+où $F$ est une fonction qui dépend de $y$ et de $t$, et où $y_0\in\real$ est la condition initiale de l'équation différentielle.
+
+Par exemple, dans le cas de l'équation différentielle
+$$
+y'(t)=y(t)+8t^2+1,
+$$
+$F(y,t)=y(t)+8t^2+1$.
+
+Afin d'avoir une solution unique, il est nécessaire de donner une condition
+initiale à une équation différentielle, de la forme $y(t_0)=y_0$. Ici on pourrait par exemple avoir
+$$
+y(t=0)=0.
+$$
+
+## Les méthodes de Runge-Kutta
+
+Les méthodes de Runge--Kutta sont une famille de méthodes explicites
+pour résoudre des équations différentielles ordinaires. Elles ont la
+forme général suivante 
+$$
+y_{n+1}=y_0+\delta t\sum_{i=1}^s b_ik_i,
+$$ où
+$\delta t$ est le pas de temps, et 
+$$
+y_n\equiv y(t_n),
+$$
+avec $t_n=t_0+n\delta t$. La valeur des $k_i$ dépend de la précision de
+l'algorithme (et donc de la valeur du nombre "d'étages" $s$) et ils sont
+donnés de façon itérative par $$\begin{aligned}
+ k_1&=F(y_n,t_n),\\
+ k_2&=F(y_n+a_{21}\delta t k_1,t_n+c_2\delta t),\\
+ &\cdots\\
+ k_s&=F(y_n+a_{s1}\delta t k_1+a_{s2}\delta t k_2+\cdots +a_{s,s-1}\delta t k_{s-1},t_n+c_s \delta t).\end{aligned}$$
+et où les $a_{ij}$, $b_i$ et $c_i$ sont propres de nouveaux au nombre
+d'étages $s$ de la méthode. Pour $s=1$, on a simplement $c_1=0$, $b_1=1$
+et donc il vient $$y_{n+1}=y_n+\delta t F(y_n,t_n).$$ Cette méthode
+s'appelle également méthode d'Euler explicite (vous l'avez peut-être
+reconnue). Pour $s=2$, on obtient $c_1=0$, $c_2=1/2$, $a_{21}=1/2$,
+$b_1=0$, et $b_2=1$. On a donc
+$$
+y_{n+1}=y_n+\delta t F\left(y_n+\frac{\delta t}{2}F(y_n,t_n),t_n+\frac{\delta t}{2}\right).
+$$
+Cette méthode s'appelle également méthode de point du milieu (elle est
+proche de la méthode du point du milieu pour les intégrales).
+
+Finalement, une méthode bien connue pour les méthode de Runge--Kutta est celle d'ordre 4
+$$
+y_{n+1}=y_n+\frac{1}{6}(k_1+2(k_2+k_3)+k_4),
+$$
+où
+\begin{align}
+k_1&=\delta tF(y_n, t_n),\\
+k_2&=\delta tF\left(y_n+\frac{k_1}{2}, t_n+\frac{\delta t}{2}\right),\\
+k_3&=\delta tF(y_n+\frac{k_2}{2}, t_n+\frac{\delta t}{2}),\\
+k_4&=\delta tF(y_n+k_3, t_n+\delta t).
+\end{align}
+
+## L'équation de Lorenz
+
+L'équation de Lorenz est un système d'équations différentielles
+ordinaires agissant sur une variable vectorielle à trois composantes
+$\vec y(t)=(y_x(t),y_y(t),y_z(t))$ 
+$$\begin{aligned}
+ \frac{\mathrm{d}y_x}{\mathrm{d}t}&=\sigma(y_y-y_x),\\
+ \frac{\mathrm{d}y_y}{\mathrm{d}t}&=y_x(\rho-y_z)-y_y,\\
+ \frac{\mathrm{d}y_z}{\mathrm{d}t}&=y_yy_x-\beta y_z,
+ \end{aligned}$$ 
+où les constantes sont données par $\beta=8/3$, $\sigma=10$, et $\rho=28$.
+
+![Solution dans le plan $x-z$ de l'équation de Lorenz pour une solution
+initiale donnée (image tirée de
+Wikipedia).](../figs/1024px-Lorenz_system_r28_s10_b2-6666.png){#fig_lorenz
+width="50%"}
+
+La solution de l'équation de Lorenz est dite *chaotique*. Cela signifie
+que si nous considérons deux conditions initiales $\vec y^1$ et $\vec y^2$ qui
+sont très proches. Elles s'éloigneront l'une de l'autre
+exponentiellement vite au cours du temps.
+
+Par ailleurs, comme on peut le constater sur la figure
+[1](#fig_lorenz){reference-type="ref" reference="fig_lorenz"} la
+solution a une structure de "papillon". Cette structure est retrouvée
+indépendemment de la condition initiale (on aura toujours cette forme de
+solution) bien que la position à un temps donné ne soit pas la même pour
+deux conditions initiales différentes.
+
+# Travail à effectuer
+
+## Écriture des solveurs
+
+Écrire un solveur pour la méthode d'Euler explicite (méthode de
+Runge--Kutta à un étage), un pour la méthode de Runge--Kutta à deux
+étages et un pour la méthode de Runge--Kutta d'ordre 4.
+
+## Équation de Lorenz
+
+Afin de valider les solveurs, appliquer les trois solveurs à l'équation de Lorenz et reproduire une
+figure approchant celle de l'énoncé (faire un graphique éventuellement
+tri-dimmensionnel de la trajectoire obtenue pour une condition initiale
+quelconque).
+
+## Orbites périodiques
+
+L'attracteur de Lorenz contient des trajectoires périodiques sur
+lesquelles le système revient au point initial après un certain temps.
+Ces trajectoires sont dites "instables", comme toute autre trajectoire
+de ce système : il suffit qu'on dévie un tout petit peu de la
+trajectoire, et on finit par s'en éloigner complètement. Le point
+suivant se trouve (dans la limite de la précision numérique fournie) sur
+une trajectoire périodique de période $t_{max}=1.5586522$:
+$y_0=(-0.9101673912,-1.922121396,18.18952097)$. Prendre ce point comme
+valeur initiale, et faites évoluer le système avec un pas de temps donné
+(par exemple $\delta t = 0.001$). Combien de temps arrivez-vous à rester
+sur l'orbite périodique avec chacun des solveurs? Tracez les deux
+trajectoires sur une figure 3D pour les comparer.
+
+# Modélisation d'épidémies
+
+Une fois les solveurs validés vous vous attaquerez à la tâche d'utiliser vos solveurs
+pour résoudre un vrai problème d'actualité.
+
+## Le modèle SEIR
+
+Afin de modéliser la propagation d'une épidémie, nous allons utiliser le modèle SEIR, qui est un modèle compartimental où $S$ est la population susceptible d'être infectée, $E$ la population exposée (pas encore infectieuse) qui correspond à la population où la population est en incubation, $I$ la population infectieuse, et finalement $R$ la population rétablie.
+
+La dynamique de ce système est caractérisée par la population initiale de chaque catégorie:
+$$
+S(t=0)=S_0,\quad E(t=0)=E_0,\quad I(t=0)=I_0,\quad R(t=0)=R_0,
+$$
+et par quatre équation différentielles ordinaires
+\begin{align}
+S'(t)&=-\frac{\mathcal{R}_0}{T_{inf}}I(t)\frac{S(t)}{N},\\
+E'(t)&=\frac{\mathcal{R}_0}{T_{inf}}I(t)\frac{S(t)}{N}-\frac{1}{T_{inc}}E(t),\\
+I'(t)&=\frac{1}{T_{inc}}E(t)-\frac{1}{T_{inf}}I(t),\\
+R'(t)&=-\frac{1}{T_{inf}}I(t),
+\end{align}
+où $\mathcal{R}_0$ est taux de reproduction de base, $T_{inf}$ est le temps où un individu est infectieux, et $T_{inc}$ est le temps d'incubation de la maladie. La taux de reproduction de base peut être remplacé par
+le taux de reproduction  effectif, $\mathcal{R}_t=\mathcal{R}_0\frac{S(t)}{N}$ qui représente
+la susceptibilité effective au temps $t$.
+
+Ce modèle peut être amélioré pour inclure d'autres termes (par exemple des termes reliés aux trajets de gens entrant ou sortant dans la population, ainsi que ternir compte des morts naturelles, etc). Pour le moment nous nous contentons de cette version simplifiée qui pourra être rendue encore plus réaliste par la suite.
+
+## Travail à effectuer
+
+L'utilisation que vous ferez de votre modèle ne dépendra que de vous. Essayez de concevoir des scénarios différents de réaction des pouvoirs publics et simulez les!
+
+# Évaluation
+
+Me cassez pas les pieds avec vos évaluations!
+
+
+