\documentclass{article} 
\usepackage{pstricks,color} 
\usepackage[latin1]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage[francais]{babel}
\usepackage{graphicx}
\usepackage{ifpdf}
\ifpdf
 \usepackage[colorlinks,pdftex]{hyperref}
\fi
\usepackage{times}
\usepackage{xspace}
\newcommand{\bs}{\symbol{92}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newtheorem{theorem}{Theorem}
\newtheorem{defn}{Definition}
\title{Mathématiques, climat et \'energie.}
\author{Bernard Parisse\\Institut Fourier\\UMR 5582 du CNRS\\
Université de Grenoble I}
\begin{document} 

\maketitle 
\begin{abstract}
Ce texte reprend certains thèmes des sciences du climat
où interviennent des mathématiques au niveau Licence. 
On abordera la définition des climats,
les causes naturelles de variation du climat (influence des paramètres orbitaux
de la Terre), et les causes anthropiques (effet de serre).
\end{abstract}

\tableofcontents


\section{Introduction}
La température du globe dépend d'une part
de l'intensité et de la répartition du rayonnement solaire,
d'autre part de la chaleur réémise vers l'espace.

Lorsque le Soleil est
à la verticale on reçoit en effet beaucoup plus d'énergie que s'il
est bas sur l'horizon. 
La position du Soleil et la durée de l'ensolleillement
journalier dépend d'une part de l'endroit de la Terre où l'on se situe
(la latitude: près des poles, des zones tempérées, des tropiques ou
de l'équateur), et d'autre part de l'inclinaison de l'axe de rotation
de la Terre et de la position de la Terre sur son orbite,
responsables des saisons, cf. figure \ref{fig:saisons}.
\begin{figure} \label{fig:saisons}
%\includegraphics[width=\textwidth]{1saison}
\includegraphics[width=10cm]{2obliquite}
\caption{Saisons et obliquit\'e (source
http://www.mnhn.fr/mnhn/lop/DIOGENE/morphologie/orbite.html)}
\end{figure}
Le calcul de l'inclinaison des
rayons du Soleil en fonction de la date et de la latitude 
fait intervenir principalement de la trigonométrie et un peu de géométrie 
dans l'espace.

La puissance du rayonnement solaire est inversement proportionnelle 
à la distance au carré entre la Terre et le Soleil (en admettant que
la puissance du Soleil soit constante à la source), la forme de l'orbite
de la Terre influe donc sur le climat. Si on néglige les autres planètes
et on assimile que la Terre est un point matériel,
l'orbite de la Terre est une ellipse (c'est-à-dire un cercle
aplati dans une direction) dont
le Soleil est un foyer, cf. figure \ref{fig:ellipse}
Au cours de son orbite la Terre se rapproche ou s'\'eloigne du Soleil
ce qui fait l\'eg\`erement varier le contraste et la dur\'ee des
saisons (actuellement les saisons sont moins contrast\'ees dans
l\'h\'emisph\`ere Nord que Sud et l'hiver bor\'eal est plus court
que l'\'et\'e bor\'eal ce qui explique que f\'evrier n'ait que 28
jours, plus pr\'ecis\'ement l'hiver bor\'eal dure 88.99 jours, le printemps
92.75 jours, l'\'et\'e 93.65 jours, l'automne 89.85 jours).
Le calcul précis de la puissance de l'ensolleillement nécessite donc
la connaissance
d'un peu de géométrie de l'ellipse, ainsi que de la résolution numérique
d'équations.

En réalité l'ellipse décrite par la Terre autour du Soleil
se déforme au cours du temps, son excentricité varie
ainsi que les dates de passage au plus près et au plus loin du Soleil,
cf. figure \ref{fig:precession} et voir pourquoi en annexe
\begin{figure}[h] \label{fig:precession}
\includegraphics[width=9cm]{51excentricite}
\caption{Pr\'ecession des \'equinoxes}
\end{figure}

De plus l'angle d'inclinaison de l'axe de rotation de la Terre
varie lentement ce qui répartit l'énergie solaire différemment,
plus l'inclinaison est grande plus les poles en reçoivent, cf. figure
\ref{fig:var_obliquite}.
\begin{figure}[h] \label{fig:var_obliquite}
\includegraphics[width=9cm]{32var_obliquite}
\caption{Variation de l'obliquit\'e (source Cyril Langlois)}
\end{figure}


Ces variations sont à la base de la théorie astronomique
des climats de Milankovitch. L'étude numérique des enregistrements 
du climat du passé sur plusieurs centaines de milliers d'années
(en particulier les forages en Antarctique, par exemple Vostok, Dome C, 
projet Epica)
est en bon accord avec cette théorie, on peut ainsi
mettre en évidence (par transformée de Fourier)
les périodicités principales de ces enregistrements et les comparer
avec celles calculées par les astronomes pour l'excentricité, la précession
des équinoxes et l'inclinaison de l'axe de rotation de la Terre,
cf. figure \ref{fig:temp_co2_methane}.
\begin{figure}[h] \label{fig:temp_co2_methane}
\includegraphics[width=10cm]{4temp_co2_methane}
\caption{Temp\'erature et gaz \`a effet de serre (Vostok)}
\end{figure}

L'arrivée au sol et la réémission de chaleur par la Terre est elle
influencée par les nuages et la présence dans l'air de molécules
ayant des fréquences de vibration propres proches de celles du
rayonnement émis par le sol (qui peuvent rendre l'air opaque à ce
type de rayonnement) provoquant l'effet de serre. 
Les molécules principalement responsables de l'effet de serre
sont le dioxyde de carbone, le méthane, les oxydes d'azote, l'ozone
ainsi que les gaz de type halofluorocarbone. Ainsi un doublement
de concentration de dioxyde de carbone équivaut à une augmentation
d'environ 4 watt par mètre carré au niveau du sol (à comparer aux 1360/4 watt
par mètre carré reçu au sommet de l'atmosphère). Leur concentration a varié
dans le passé de manière naturelle en phase avec les enregistrement 
de températures relevés dans les glaces, cf. ci-dessus, on peut le mettre
en évidence de manière statistique, ici par une régression linéaire de
la temp\'erature en fonction du CO2, figure \ref{fig:regression}
\begin{figure} \label{fig:regression}
\includegraphics[width=9cm]{Temperature_vs_CO2}
\caption{Temp\'erature en Antarctique et taux de CO2}
\end{figure}

Depuis le
début de l'ère industrielle, cette concentration a beaucoup augmenté,
largement au-delà des limites naturelles atteintes au cours des 750 mille
dernières années (par exemple 380 parties par million pour le CO2 a
Mauna Loa, largement
en-dehors de l'intervalle 180-290, en hausse de 2 ppm par an, figure
\ref{fig:co2})
\begin{figure} \label{fig:co2}
\includegraphics[width=9cm]{co2actuel}\\
\includegraphics[width=9cm]{Augmentation_co2_an}
\caption{\'Evolution r\'ecente du taux de CO2 (absolue et moyenne
  mobile de la hausse)}
\end{figure}
faisant craindre un réchauffement de la planète au
21ème siècle et au cours des siècles suivants, r\'echauffement
tr\`s probablement d\'ej\`a commenc\'e, figure
\ref{fig:ipcc}.
\begin{figure} \label{fig:ipcc}
\includegraphics[width=9cm]{6temp_actuel}
\caption{Temp\'erature actuelle (source IPCC)}
\end{figure}

La cause principale
de l'augmentation de ces concentrations est l'utilisation d'énergie
fossiles pour le CO2 (pétrole, gaz naturel et charbon) et d'engrais azotés
(agriculture intensive). Une simple r\`egle de trois permet de
faire la correspondance 
\begin{figure} \label{fig:tep2CO2}
\includegraphics[width=9cm]{80tep2CO2}
\includegraphics[width=9cm]{94scenario_b2_energie}
\includegraphics[width=9cm]{Bpopulation}
\caption{Carbone \'emis par tonne \'equivalent p\'etrole (source J.M. Jancovici) et exemple de sc\'enario \'energ\'etique et de population
pour le 21i\`eme si\`ecle (IPCC et ONU)}
\end{figure}
Ainsi, en un an (décembre 2004-novembre 2005), 
la France a consommé 12.9tonnes
équivalent pétrole (tep) de charbon, 95.1 tep de pétrole et 40.5 tep
de gaz, dont la combustion a émis dans l'atmosphère 122 millions 
de tonnes de carbone (sous forme de CO2),
soit 1.8 tonne de carbone par habitant, si le reste du monde consommait
au même rythme cela correspondrait à 6 parties par millions de CO2
émises par an.

La quantification
de l'ampleur du réchauffement (et des autres paramètres climatiques
en particulier la répartition et l'intensité des précipitations et
des vents) se fait en deux temps, d'une part
par des scénari d'emission de gaz à effet de serre, d'autre part
par modélisation numérique du climat en fonction des scénari d'émission.
La modélisation numérique donne des résultats relativement précis 
pour la température (3 à 4 degré pour un doublement du CO2), 
plus incertains pour les précipitations. 
%Les modèles sont validés en modélisant le climat 
%du 20ème siècle et en comparant
%avec les températures moyennes calculées pour le globe. 

Par contre
la modélisation des émissions de gaz à effet de serre est beaucoup
plus difficile, d'une part parce que les réserves de combustibles
fossiles utilisables sont controversées (par exemple sur la date
du pic de production pétrolier), d'autre part parce que le mix
énergétique du 21ème siècle est largement inconnu (fossile, nucléaire,
renouvelables?), mais aussi parce que la réponse de la biosphère
au réchauffement est encore mal quantifiée (par exemple on estime que
l'océan absorbe actuellement environ la moitié des émissions de CO2,
mais la solubilité du gaz carbonique dans l'eau diminue lorsqu'on
la réchauffe). Les mathématiques interviennent ici d'abord pour
définir le climat (par exemple pourquoi on utilise une série
de 30 années pour définir une température moyenne en un lieu donné),
dans la réalisation des modèles numériques (nous ne parlerons pas de
ce sujet qui dépasse le niveau licence scientifique), 
dans les modèles statistiques de réchauffement dus aux gaz
à effet de serre (corrélation), ainsi que dans les modèles d'émission 
de gaz à effet de serre, par
l'intermédiaire des modèles de production et 
consommation de combustibles fossiles.

\begin{figure} \label{fig:petrole_us}
\includegraphics[width=9cm]{A0petrole_us}
\includegraphics[width=9cm]{Curr_petrole}
\includegraphics[width=\textwidth]{Durr_fossiles}
\caption{Production de p\'etrole aux \'Etats-Unis (source
EIA et Oildrum) mettant en \'evidence le ph\'enom\`ene de pic
p\'etrolier et recherche du meilleur mod\`ele de pr\'ediction
1976-2004 connaissant l'avant 1976. En millions de barils par jour
(1 baril=159 litres).
Mod\'elisation des ressources ultimes en p\'etrole 
en milliards de barils et de production d'\'energie en combustibles
fossiles en milliards de tonne \'equivalent
p\'etrole par an (source Lah\'err\`ere)}
\end{figure}

\pagebreak 


\section{Le rayonnement solaire.}

\subsection{L'insolation au cours de l'ann\'ee.} \label{sec:soleil}
Pour connaitre la quantit\'e d'\'energie recue \`a un moment
donn\'e, il faut calculer l'angle entre la verticale du lieu
et la direction du Soleil. Plus g\'en\'eralement, on va calculer
les composantes du vecteur Terre-Soleil et les composantes
des vecteurs de la base locale (verticale locale, direction du Sud et
direction du parall\`ele). On choisit d'abord comme r\'ef\'erence
le plan $Txy$ de l'\'ecliptique (plan de l'orbite de la Terre
autour du Soleil), avec $Ty$ orthogonal \`a l'axe de rotation 
de la Terre (donc $Tx$ la projection de l'axe de rotation
de la Terre sur ce plan). Soit $\theta$ l'angle que fait la Terre
avec la direction du passage au p\'erih\'elie, et $\theta_0$
l'angle de la position de la Terre au solstice d'hiver avec la
direction du p\'erih\'elie, l'angle entre la direction Terre-Soleil
et $Tx$ est donc $\theta-\theta_0$\\
\includegraphics[width=\textwidth]{insol}\\
Dans $Txyz$ le vecteur unitaire $s$
de la direction Terre-Soleil a pour coordonn\'ees~:
\[ s=-(\cos(\theta-\theta_0),\sin(\theta-\theta_0),0)\]
On effectue ensuite une rotation autour
de $Ty$ d'angle $i$ l'inclinaison de l'axe de rotation de la Terre.
On obtient ainsi un rep\`ere $TXyZ$ ($TX$ et $TZ$ se déduisent de 
$Tx$ et $Tz$ par rotation d'angle $i$). Dans ce repère 
le vecteur untaire $s$ a pour coordonnées~:
\[ s=-(\cos(\theta-\theta_0)\cos(i),\sin(\theta-\theta_0),
\cos(\theta-\theta_0)\sin(i)) \]
Calculons maintenant dans ce rep\`ere $TXyZ$ les coordonn\'ees
des vecteurs de la base locale. On se place en un point de latitude
$l$ et de longitude $\phi$, on note $J$ la dur\'ee d'une p\'eriode
de r\'evolution de la Terre sur elle-m\^eme (23 heures 56 minutes,
c'est un peu moins d'un jour car il faut encore en moyenne 4 minutes
pour compenser le d\'eplacement de la Terre sur son orbite autour
du Soleil). La verticale locale a pour coordonn\'ees~:
\[ v=(\cos(l)\cos(\phi+2\pi t/J),\cos(l)\sin(\phi+2\pi t/J),\sin(l))
\]
L'\'energie solaire recue au lieu donn\'e (sur une surface horizontale ;
pour un panneau solaire, il faudrait calculer les coordonn\'ees d'un vecteur
perpendiculaire au panneau) est proportionnelle \`a
\[ \frac{s.v}{\rho^2} \]
o\`u $\rho(\theta)=a(1-e^2)/(1+e\cos(\theta))$ 
d\'esigne la distance Terre-Soleil.
Le calcul de $s.v$ donne, en notant $\varphi=\phi+2\pi t/J$~:
\[
-s.v  =  \cos(l) \cos(\varphi) \cos(\theta-\theta_0)\cos(i)
+ \cos(l)\sin(\varphi) \sin(\theta-\theta_0) 
+ \sin(l) \cos(\theta-\theta_0)\sin(i) \]
On rassemble les deux premiers termes qui d\'ependent rapidement
du temps par l'interm\'ediaire de $\varphi$
(le 3\`eme terme n'en d\'epend que par $\theta$ qui ne varie
que d'environ 1 degr\'e pendant une journ\'ee)
et on applique la formule de trigonom\'etrie:
\[ A \cos \alpha + B \sin \alpha = \sqrt{A^2+B^2}
\cos(\alpha-\alpha_0), \quad \left\{ \begin{array}{ccc}
\cos(\alpha_0) &=&
\frac{A}{\sqrt{A^2+B^2}} \\
\sin(\alpha_0)&=&\frac{B}{\sqrt{A^2+B^2}} \end{array} \right. \] 
Ici, apr\`es avoir factoris\'e $\cos(l)$, on a~:
\[ \sqrt{A^2+B^2}=\sqrt{
  \cos(\theta-\theta_0)^2\cos(i)^2+\sin(\theta-\theta_0)^2 }
= \sqrt{1 - \sin(i)^2 \cos(\theta-\theta_0)^2 } \]
On peut aussi calculer 
\[ \tan(\alpha_0)=\frac{B}{A} =  \frac{\tan(\theta-\theta_0)}{\cos(i)}\]
qui donne $\alpha_0$ modulo $\pi$ et compl\'eter en regardant 
le quadrant o\`u se trouve $(A,B)$, ici $\alpha_0$ et $\theta-\theta_0$ 
sont tous deux dans $[0,\pi]$ ou tous deux dans $[-\pi,0]$.
Finalement, on obtient le
\begin{theorem} \label{thm:soleil}
L'\'energie solaire recue au sol est proportionnelle \`a
\[ (1+e\cos(\theta))^2 s.v \]
o\`u $s.v$ est donn\'e par~:
\[ s.v = - \cos(l) \sqrt{1 - \sin(i)^2 \cos(\theta-\theta_0)^2 }
\cos(\varphi-\varphi_0)
- \sin(l) \cos(\theta-\theta_0)\sin(i)
\]
et~
\begin{itemize}
\item $e$ est l'excentricit\'e de l'orbite elliptique (environ
0.0167 actuellement)
\item $i$ est l'obliquit\'e (inclinaison de l'axe de rotation de la Terre,
environ 23 degr\'e 27 minutes actuellement)
\item $\theta$ est l'angle fait par la direction Terre-Soleil avec la 
direction du demi grand axe (Soleil-p\'erih\'elie), 
$\theta_0$ le m\^eme angle au solstice d'hiver
de l'h\'emisph\`ere Nord (environ -13 degr\'es). En premi\`ere 
approximation, on peut faire varier $\theta$ proportionnellement au temps,
voir la fin de la section \ref{sec:dureesaison} pour un calcul plus pr\'ecis.
\item $l$ est la latitude, $\varphi$ la longitude tenant compte
de la rotation de la Terre (somme de la longitude g\'eographique $\phi$ et 
du terme d\'ependant du temps $2\pi t/J$)
\item  $\varphi_0 \in [-\pi,\pi]$ est de m\^eme signe que $\theta-\theta_0$
et est d\'efini par~:
\begin{equation} \label{eq:culmination}
\tan \varphi_0 = \frac{\tan(\theta-\theta_0)}{\cos(i)} 
\end{equation}
\end{itemize}
\end{theorem}
{\bf Variations de $s.v$ au cours d'une journ\'ee dans l'approximation o\`u
$\theta$ ne varie pas}~:\\
On obtient une sinusoide entre les deux valeurs extr\^emes~:
\[ \pm \cos(l) \sqrt{1 - \sin(i)^2 \cos(\theta-\theta_0)^2 } -
 \sin(l) \cos(\theta-\theta_0)\sin(i) \]
Le maximum est atteint pour 
\[ \varphi=\varphi_0+ \pi \ \Rightarrow \ 
2 \pi t/J = \varphi_0 - \phi + \pi, \quad J=23h56m \]
le moment correspondant est appel\'e culmination (c'est le midi
solaire si le maximum est positif) et ne d\'epend pas de la latitude
(bien entendu la valeur du maximum en d\'epend).
Si le maximum est n\'egatif ou nul, la nuit dure 24h.
Si le minimum est positif ou nul, le jour dure 24h.
Par exemple au solstice d'hiver, $\theta=\theta_0$, selon la latitude
on obtient un maximum n\'egatif pour $l=\pi/2$ (pole Nord), positif 
pour $l=-\pi/2$ (pole Sud), le minimum et le maximum croissent 
entre ces 2 valeurs.
Si le maximum est positif et le minimum est n\'egatif, il y a
2 instants ou $s.v=0$ (lever et coucher du soleil).

L'\'energie solaire recue pendant une journ\'ee par une surface
horizontale est proportionnelle \`a l'int\'egrale entre le lever
et le coucher de $s.v/\rho^2$. S'il n'y a pas de lever/coucher, 
soit on ne recoit rien (nuit polaire),
soit on recoit l'int\'egrale entre 0 et 24h de $s.v$ (jour polaire).

L'intervalle entre 2 culminations n'est pas constant au cours de
l'ann\'ee, car $\varphi_0$ n'est pas une fonction lin\'eaire de
$\theta$ (qui lui m\^eme n'est pas lin\'eaire en fonction du temps
sauf en premi\`ere approximation avec une orbite terrestre
circulaire). On peut le calculer en d\'erivant (\ref{eq:culmination}).
Par exemple dans l'approximation d'une excentricit\'e nulle,
au solstice d'hiver ($\theta=\theta_0$), on obtient
\[ \varphi_0= 0, \quad (1+0^2) d\varphi_0 = (1+0^2)\frac{d\theta} {\cos(i)} \]
avec $d\theta$ qui correspond \`a 4 minutes, on trouve $d\varphi_0$
correspondant \`a 4.36 minutes.
L'\'ecart entre 2 culminations
est donc d'environ 24h 20secondes. 
Au moment du solstice, le Soleil se l\`eve et se couche donc environ
20 secondes plus tard entre un jour et son lendemain, dans 
l'hypothèse d'un mouvement circulaire de la Terre autour du Soleil. 
En réalité, l'orbite terrestre étant faiblement elliptique,
l'écart est un peu moins de 30 secondes en hiver et de 15 secondes en été, 
le mouvement de la Terre autour du Soleil étant plus rapide d'environ
3\% au solstice d'hiver et moins rapide d'environ 3\% au
solstice d'été. Comme 3\%
de l'écart moyen entre 2 culminations (4 minutes=240 secondes) 
correspond à 7 secondes cela explique la différence.

\subsection{Les saisons}
Les deux figures (\ref{fig:nord}, \ref{fig:sud})
qui suivent représentent la Terre en 4 points de son orbite, les 
2 solstices et les 2 \'equinoxes.
Les solstices sont définis par les 2 points de l'orbite où la projection
de l'axe de rotation terrestre est parallèle à l'axe Terre-Soleil.
Les équinoxes sont définis par les 2 points où il y a perpendicularité.
\begin{figure}[h] \label{fig:nord}
\includegraphics[width=10cm]{saison_nord}
\caption{Saisons (hémisphère Nord), source wikipedia}
\end{figure}
\begin{figure}[h] \label{fig:sud}
\includegraphics[width=10cm]{saison_sud}
\caption{Saisons (h\'emisph\`ere Sud), source wikipedia}
\end{figure}
Au solstice d'hiver, on voit que les parallèles situés aux hautes latitudes
Nord ne sortent jamais de l'obscurité. Aux latitudes intérmédiaires,
le morceau de parallèle situé au jour est nettement plus petit que
celui situé dans l'obscurité.
\`A l'équinoxe de printemps, chaque parallèle est à moitié au jour
et à moitié dans l'obscurité (derrière la grille). Au printemps,
la situtation est analogue.
Au solstice d'été, on est dans la situation inverse de l'hiver.

\pagebreak

\subsection{L'orbite de la Terre.}
En première approximation, l'orbite de la Terre est uniquement
influencée par la force de gravitation entre la Terre et le Soleil,
ce dernier pouvant être considéré comme fixe en raison de sa masse
(on peut éviter cette approximation en remplaçant le Soleil par
le centre de gravité du système Terre-Soleil).
La force de gravitation qui dérive d'un potentiel inversement
proportionnel à la distance Terre-Soleil est de la forme
\[ {\bf F'}=\frac{{\bf F}}{m_T}= K \frac{\bf r}{r^3}, \quad K=-\mu<0\]
o\`u ${\bf r}$ d\'esigne le vecteur Terre-Soleil, $m_T$ est la masse
de la Terre, $\mu=Gm_S$ est le produit de la constante
de gravitation universelle par la masse du Soleil.
Le moment cinétique de la rotation de la Terre autour du Soleil
est d\'efini par~:
\[ {\bf L} = {\bf r} \wedge \frac{d \bf r}{dt} \]
On v\'erifie que sa d\'eriv\'ee est nulle, donc ${\bf L}$ garde
une direction fixe {\bf k}, orthogonale \`a ${\bf r}$,
l'orbite de la Terre reste donc dans le plan défini \`a un instant
donn\'e par l'axe Terre-Soleil et le vecteur vitesse de la Terre.
De plus la conservation de {\bf L} entraine 
la loi des aires, l'aire balayée par le rayon Soleil-Terre
est proportionnelle au temps.

On utilise un repère en coordonnées polaires centré au Soleil,
$\rho$ désignant la distance Terre-Soleil et $\theta$ l'angle
fait par rapport à une direction fixe, on a alors
\[ {\bf L} = \rho^2 \frac{d\theta}{dt} {\bf k}\]
car si on calcule en coordonn\'ees polaires $d {\bf r}/dt$,
la composante sur le vecteur radial ${\bf e_r}$ est $d\rho/dt$, et
la composante sur le vecteur perpendiculaire ${\bf e_\theta}$
est $\rho d\theta/dt$.

\subsubsection{Calcul en utilisant le vecteur excentricit\'e.}
Montrons que le vecteur
\[ {\bf E}= \frac{1}{\mu} \frac{\bf dr}{dt} \wedge {\bf L} - 
\frac{\bf r}{\rho}\]
est aussi conserv\'e (o\`u on rappelle que $\mu$ provient
de la force de gravitation ${\bf F'}={\bf F}/m_T=-\mu {\bf r}/r^3$).
Le deuxi\`eme terme est proportionnel au vecteur
radial $-{\bf e_r}$, dont la d\'eriv\'ee est le vecteur orthogonal
$-d\theta/dt {\bf e_\theta}$.
Comme ${\bf L}$ est constant, la d\'eriv\'ee du premier terme est
\[ \frac{1}{\mu} {\bf F'} \wedge {\bf L} = \frac{\bf -e_r}{\rho^2} \wedge L {\bf k}
= \frac{L}{\rho^2} {\bf e_\theta} = - \frac{d\theta}{dt} \rho^2\]

Notons que ${\bf E}$ est dans le plan de l'orbite, prenons comme origine
des angles pour rep\'erer la Terre par rapport au Soleil la direction
de ${\bf E}$. En faisant le produit scalaire de ${\bf E}$ avec ${\bf r}$, 
on obtient en notant $e$ la norme de $E$
\begin{eqnarray*} e \rho \cos(\theta) &=&  {\bf E}.{\bf r} \\
& = &  (\frac{1}{\mu} \frac{\bf dr}{dt} \wedge {\bf L}- 
\frac{\bf r}{\rho}) .{\bf r} 
\\
& = & (\frac{1}{\mu} \frac{\bf dr}{dt} \wedge {\bf L}).{\bf r}
- \rho \\
& = & \frac{1}{\mu} ({\bf r} \wedge  \frac{\bf dr}{dt} ).{\bf L} - \rho \\
&= & \frac{1}{\mu} {\bf L}.{\bf L}- \rho \\
&= &\frac{L^2}{\mu} - \rho
\end{eqnarray*}
d'o\`u~:
\[ \rho = \frac{L^2}{\mu(1+e \cos(\theta))}\]

\subsubsection{Calcul par l'\'equation diff\'erentielle.}
On a les équations de conservation de l'\'energie et du moment cin\'etique~:
\[ \frac{K}{\rho}+\frac{m}{2} \left( \left(\frac{d\rho}{dt}\right)^2
+ \left(\rho \frac{d\theta}{dt} \right)^2 \right)
=C_1, \quad \rho^2 \frac{d\theta}{dt} = L, \quad K<0, m>0
\]
On change de variable dépendante pour $\rho$, 
en prenant $\theta$ au lieu de $t$,
comme $\frac{d\rho}{dt}=\rho' \frac{d\theta}{dt}$ 
(où $\rho'=\frac{d\rho}{d\theta}$), on a~:
\[
 \frac{K}{\rho}+\frac{m}{2}\left( (\rho'^2+\rho^2) 
\left(\frac{L}{\rho^2}\right)^2 \right) = C_1
\]
On effectue ensuite le changement de variable $\rho=1/u$, $\rho'=-u'/ u^2$,
d'où~:
\[
K u + \frac{m}{2} \left( 
\left(\frac{u'^2 }{u^4}+\frac{1}{u^2}\right) L^2 u^4 \right)
= C_1
\]
soit~:
\[ K u + \frac{mL^2}{2} ( u'^2+u^2) = C_1 \]
donc~:
\[ K' u + u^2 + u'^2 = C_3, \quad 
K'=\frac{2K}{mL^2} <0, C_3=\frac{2C_1}{mL^2}
\]
On pose maintenant $v=u+K'/2$, d'où~:
\[ v^2 + v'^2 = C_3 + \frac{K'^2}{4} = C_4 \]
On montre (en exprimant $v'$ en fonction de $v$ puis en séparant les
variables), que cette équation différentielle a pour solution générale
\[ v= \sqrt{C_4} \cos(\theta-\theta_0) \]
D'où~:
\[ \rho = \frac{1}{u} =
\frac{1}{v-\frac{K'}{2}} = 
\frac{1}{\sqrt{C_4} \cos(\theta-\theta_0)-\frac{K'}{2}}
\]
Comme $K'<0$ et comme la trajectoire de la Terre autour du Soleil
passe par tous les angles (donc $\rho$ est défini pour tout $\theta$,
le dénominateur ne peut pas s'annuler), on a~:
\[ \rho = \frac{\frac{2}{-K'}}{1+e \cos(\theta-\theta_0)}, \quad e \in [0,1[ \]
On définit ensuite $a$ par $\frac{2}{-K'}=a(1-e^2)$, et on obtient finalement
l'équation d'une ellipse dont l'origine (le Soleil) est un des
foyers~:
\[  \rho = \frac{a(1-e^2)}{1+e \cos(\theta-\theta_0)} \]
On suppose maintenant quitte à faire pivoter l'axe des $x$ que
$\theta_0=0$.

\subsubsection{Lois de K\'epler.}
L'orbite de la Terre est donc une ellipse dont le Soleil occupe un des foyers
(1\`ere loi de K\'epler).
On a aussi vu que $L=\rho^2 d\theta/dt$ est constant, ceci entraine la loi
des aires, infinit\'esimalement on a~:
\[ \frac{1}{2}\rho^2 d\theta = \frac{1}{2} L dt \]
ce qui se traduit par
l'aire balay\'ee par le rayon vecteur Soleil-Terre est proportionnelle
au temps (2\`eme loi de K\'epler). 
Au cours d'une p\'eriode $T$, l'aire parcourue est celle
de l'ellipse, donc 
\[ \pi a^2 \sqrt{1-e^2}=\frac{1}{2}LT \]
En prenant le carr\'e, et en appliquant
\[ \frac{L^2}{\mu}= a(1-e^2)\]
on en d\'eduit la troisi\`eme loi de K\'epler~:
\[ 4\pi^2 a^3 = \mu T^2 \Leftrightarrow \frac{a^3}{T^2} = \frac{\mu}{4\pi^2} \]
o\`u on rappelle que $\mu$ est le produit de la constante de gravitation
universelle par la masse du Soleil. (On peut \'evidemment faire
le m\^eme calcul pour la Lune autour de la Terre).

\subsection{Quelques propriétés de l'ellipse}
{\bf D\'efinition}\\
L'ellipse $E$ de foyers $F_1$ et $F_2$ de demi-grand axe
$a$ est l'ensemble des points $M$ du plan tels que 
\[ MF_1+MF_2=2a\]
Cf. figure \ref{fig:ellipse}.
\begin{figure}[h] \label{fig:ellipse}
\includegraphics[width=\textwidth]{30ellipse}
\caption{ellipse (source Florence Kalfoun)}
\end{figure}

On note $2c=F_1F_2$ la distance entre les deux foyers, qui doit \^etre
plus petite que $2a$ pour que l'ellipse soit non vide.
L'excentricit\'e de l'ellipse est d\'efinie par $e=c/a < 1$. Si $e=0$,
on obtient un cercle de centre $F_1=F_2$ et de rayon $a$. Si $e\neq 0$,
on va voir qu'il s'agit d'un cercle contracté
selon l'axe perpendiculaire \`a $F_1F_2$ dans un rapport de 
$\sqrt{1-e^2}$. On va \'egalement calculer l'\'equation en
coordonn\'ees polaires de $E$ pour montrer que l'\'equation obtenue
ci-dessus est bien celle d'une ellipse dont le Soleil occupe un foyer.

Soit $O$ le milieu de $F_1$ et $F_2$, on se place dans le rep\`ere orthonorm\'e
dont le premier axe $Ox$ contient $F_1$ et $F_2$ donc les
coordonn\'ees de $F_1$ sont $(c,0)$ et celles de $F_2$ sont $(-c,0)$. 
Soit $M(x,y)$ un
point de l'ellipse, on a d'une part~:
\[ MF_1^2 - MF_2^2 = (x-c)^2-(x+c)^2 = -4cx \]
et d'autre part~:
\[ MF_1^2 - MF_2^2 = (MF_1 + MF_2)(MF_1 - MF_2 ) = 2a (MF_1 - MF_2 )\]
donc~:
\[ MF_1 - MF_2 = \frac{-2cx}{a} \]
en additionnant avec $MF_1+MF_2=2a$ et en appliquant $c=ea$, on en d\'eduit~:
\begin{equation} \label{eq:MF1}
 MF_1 = a - \frac{cx}{a} = a-ex 
\end{equation}
En prenant le carr\'e, on a~:
\[ (x-ea)^2 + y^2 = (a-ex)^2\]
d'o\`u~:
\[ y^2 + x^2 (1-e^2) = a^2(1-e^2) \]
finalement~:
\[ x^2 + \frac{y^2}{1-e^2} = a^2 \]
qui est bien la contraction selon $Oy$ de rapport $\sqrt{1-e^2}$ du
cercle de centre $O$ et de rayon $a$ (appel\'e grand cercle de
l'ellipse).

En coordonn\'ees polaires, on note $\rho$ la distance de $F_1$ \`a
$M$, et $\theta$ l'angle entre l'axe $Ox$ et $F_1M$. L'abscisse de $M$
est donc~:
\[ x= ea + \rho \cos(\theta)\]
que l'on combine avec (\ref{eq:MF1}) pour obtenir~:
\[ \rho = a-ex =a(1-e^2) - e \rho \cos(\theta) \]
donc~:
\[ \rho = \frac{a(1-e^2)}{1+e\cos(\theta)} \]
ce qui nous permet d'affirmer que l'orbite de la Terre dans
l'approximation du point mat\'eriel soumis uniquement au Soleil
suppos\'e fixe est une ellipse dont le Soleil occupe un foyer.

\subsection{Influence de l'ellipse sur les saisons}
Il faut prendre garde \`a ne pas confondre les solstices et
\'equinoxes avec le moment o\`u la Terre coupe le grand axe
de son ellipse autour du Soleil. Il n'y a aucune raison que la
projection de l'axe de rotation de la Terre sur le plan de l'ellipse
soit parall\`ele ou perpendiculaire au grand axe de l'ellipse, et
actuellement ce n'est pas le cas, le solstice d'hiver a lieu le
21 d\'ecembre alors que le passage au plus proche du Soleil a lieu
vers le 3 janvier (donc pendant l'hiver de l'h\'emisph\`ere Nord)
et le passage au plus loin du Soleil a lieu d\'ebut juillet
(pendant l'\'et\'e). C'est pour cette raison que les saisons
sont moins marqu\'ees dans l'h\'emisph\`ere Nord que dans 
l'h\'emisph\`ere Sud. De plus la loi des aires oblige la Terre
a se d\'eplacer plus vite lorsqu'elle est proche du Soleil que
lorsqu'elle en est \'eloign\'ee ce qui diminue la dur\'ee de l'hiver
bor\'eal et augmente la dur\'ee de l'\'et\'e bor\'eal (c'est pour
cette raison que f\'evrier n'a que 28 jours alors que juillet et
aout ont 31 jours).

\subsection{L'équation du temps, la durée des saisons.} \label{sec:dureesaison}
\begin{figure}[h] \label{fig:ell}
\includegraphics[width=\textwidth]{ellipse}
\caption{Ellipse et \'equation du temps}
\end{figure}
La trajectoire elliptique $E$ de la Terre autour du Soleil est
représentée sur la figure \ref{fig:ell}
en bleu, l'excentricité de l'orbite a été énormément exagérée, il
s'agit d'une ellipse de foyers $S$ (le Soleil) et $S'$. Le point $A$ désigne
le p\'erihélie de l'orbite (passage de la Terre au plus proche du Soleil), 
qui a lieu
vers le 4 janvier. En noir, on a dessiné le grand cercle de l'ellipse
(l'ellipse s'obtient par contraction du grand cercle de rapport $\sqrt{1-e^2}$
où $e$ est l'excentricité de l'orbite).
L'aire décrite par le rayon Soleil-Terre ($ST$)
est proportionnelle au temps (loi des aires qui découle de la conservation
du moment cinétique), 
il en est donc de même de l'aire (en vert) du
décrite par le rayon $SM$. Si on ajoute à cette aire verte l'aire en 
rouge du triangle $OSM$, on obtient l'aire de l'arc de cercle $OAM$.
Donc
\[ \frac{1}{2} V \times OA^2 - \frac{1}{2} OS \times HM \]
est proportionnel au temps écoulé depuis le passage au p\'erihélie.
Comme $HM=OM \sin(V)$ et $OS= e \times OA$, on en déduit que 
\begin{equation} \label{eq:tempsV}
 V - e \sin(V) = C t = 2\pi \frac{t}{T}
\end{equation}
où la constante $C$ s'obtient en faisant varier $V$ de 0 à $2\pi$
ce qui correspond à la durée $T$ d'une révolution de la Terre 
autour du Soleil (1 an).

La relation entre $\theta$ (not\'e $t$ sur la figure)
et $V$ s'obtient par exemple en calculant 
l'abscisse de $M$ 
\begin{eqnarray*} 
x &= &a \cos(V) \\
 &= & ea + \rho \cos(\theta) \\
 &= & ea + a\frac{1-e^2}{1+e\cos(\theta)} \cos(\theta)
\end{eqnarray*}
Les angles $V$ et $\theta$ sont de même signe et
\begin{equation} \label{eq:Vt}
 \cos(V) = \frac{\cos(\theta)+e}{1+e\cos(\theta)}
\end{equation}
et réciproquement~:
\begin{equation} \label{eq:tV}
 \cos(\theta)=\frac{\cos(V)-e}{1-e \cos(V)} 
\end{equation}

{\bf Durée des saisons~}:\\
Il suffit de connaitre l'angle $\theta$ lors du solstice d'hiver et
de lui ajouter $k\pi/2$ pour $k=1,2,3$ pour connaitre l'angle $\theta$
au printemps, en été et à l'automne, on en déduit $V$ par (\ref{eq:Vt}) 
puis le temps écoulé depuis le p\'erihélie avec (\ref{eq:tempsV}).

{\bf Calcul de $\theta$ en fonction du temps écoulé depuis le passage
au p\'erihélie}~:\\
Il faut calculer $V$ par des méthodes numériques (point fixe ou 
méthode de Newton) en appliquant (\ref{eq:tempsV}), on en déduit
$\theta$ avec (\ref{eq:tV}).
En résumé, on a le~:
\begin{theorem}
Soit $\theta$ l'angle entre le demi grand axe de l'ellipse et 
la direction Soleil-Terre, $t\in [-T/2,T/2]$ le temps écoulé depuis le passage
au périhélie ($t=0$ lorsque $\theta=0$, $T=1$ an). Soit $V\in [-\pi,\pi]$ 
la solution de 
\[  V - e \sin(V) = 2\pi \frac{t}{T} \]
où $e$ est l'excentricité de l'ellipse.
Alors $\theta$ est donné par
\[  \cos(\theta)=\frac{\cos(V)-e}{1-e \cos(V)}  \]
\end{theorem}

\subsection{Les variations des paramètres orbitaux}
La Terre n'est pas une sphère idéale, elle a un renflement au niveau
de l'équateur, due à rotation de la Terre sur elle-même (la force
centrifuge y est plus importante). Ce renflement est dans un plan
qui fait un angle avec le plan de l'écliptique, le Soleil exerce
donc un couple sur ce renflement. Ce phénomène est à l'origine
de la précession des équinoxes, le passage au p\'erihélie de la Terre
se décale dans le temps.
De plus, la Terre n'est pas seulement soumise à l'influence du Soleil,
mais aussi des autres planètes, en particulier Jupiter. Cela
modifie sur de très longues périodes tous les paramètres
de l'orbite terrestre, en particulier l'excentricité, la précession
des équinoxes, mais aussi l'obliquité (inclinaison de l'axe de
rotation terrestre par rapport à la perpendiculaire au plan de
l'écliptique).
Le calcul de ces variations est bien au-delà des prétentions de ce texte,
le lecteur intéressé pourra se référer par exemple aux publications
de Laskar (chercher ce mot-clef ou des mots comme orbite,
perturbation, symplectique, hamiltonien, ...). On se bornera ici
\`a indiquer que le demi-grand axe ne varie pas, ce qui donne
une relation entre les variations de la constante des aires et 
de l'excentricit\'e 
\[ L \mbox{ est proportionnel \`a } \sqrt{1-e^2}\]

Les variations des paramètre orbitaux
modifient à long terme l'ensolleillement de la Terre (la valeur
de l'énergie reçue en un lieu sur une surface horizontale $s.v/\rho^2$
dépend de la latitude, de la position de la Terre sur son orbite
mais aussi de l'excentricité de l'orbite, de l'obliquité
et de la date du p\'erihélie par rapport aux saisons) et
sa répartition sur le globe par latitude, il est naturel de supposer
qu'elles influent sur le climat de la Terre.
Par exemple, l'\'energie moyenne recue par la Terre au cours d'une
p\'eriode $T$ de une ann\'ee est donn\'ee par
\[ \frac{1}{T}\int_0^{T} \frac{dt}{\rho^2} = 
\frac{1}{T}\int_0^{2\pi} \frac{d\theta}{L} = \frac{2 \pi}{T L}  \]
est proportionnelle \`a $\frac{1}{(T L)}$ donc \`a $(1-e^2)^{-1/2}$ (car
$T$ est aussi constant d'apr\`es la 3\`eme loi de K\'epler).
Au premier ordre, la variation de $e$ entraine donc une variation de 
l'ensolleillement global de
\[ \frac{1}{2} e^2\]
Pour la Terre, cela repr\'esente au plus 2.5 pour mille (la p\'eriode
la plus favorable aux glaciations \'etant celle o\`u l'orbite est
circulaire), soit, sans r\'etroactions, une variation globale
de 0.2 degr\'es Kelvin.

\section{Les climats du passé}
Les carottages de la glace accumulée à Vostok (Antarctique) et 
ailleurs permettent
de reconstituer la température en Antarctique (ou ailleurs) jusqu'à
700 mille ans. On peut effectuer une étude numérique des valeurs
obtenues pour faire apparaitre les principales périodicités
de ce signal numérique, et comparer avec les périodicités de
l'orbite terrestre. On utilisera aussi les données de ces carottes
pour mettre en évidence les corrélations entre gaz à effet de serre
et température dans le passé.

\subsection{La transformée de Fourier discrète.}
Soit $N$ un entier fixé. Une suite $x$ p\'eriodique de p\'eriode $N$ est
d\'etermin\'ee par le vecteur $x=[x_0,x_1,...x_{N-1}]$.
La transform\'ee de Fourier discr\'ete (DFT) notée $F_N$ fait correspondre 
à une suite $x$ p\'eriodique de p\'eriode $N$ une autre suite $y$
p\'eriodique de p\'eriode $N$, d\'efinie pour $k=0..N-1$ par :
\[ 
{(F_N(x))}_k=y_k=\sum_{j=0}^{N-1} x_j \omega_N^{-k\cdot j},
\]
où $\omega_N$ est une racine $N$-i\`eme primitive de l'unit\'e,
on prend $\omega=e^{\frac{2i\pi}{N}}$ si $x$ est à coefficients réels
ou complexes.
Cette transformation est linéaire, la transformée de la somme de 2
suites est la somme des transformées, et la transformée du produit
par une constante d'une suite 
est le produit par cette constante de la transformée
de la suite.

La transform\'ee de Fourier discr\'ete  $F_N$ est une transformation 
bijective dont la réciproque est donnée par~:
\[ F_N^{-1}=\frac{1}{N} \overline{F_N}, \quad
{(F_N^{-1}(y))}_k=\frac{1}{N}\sum_{j=0}^{N-1}y_j\omega_N^{k\cdot j}
\]
On le prouve en remplaçant $y$ par sa valeur~:
\begin{eqnarray*}
{(F_N^{-1}(y))}_k &=&
\frac{1}{N}\sum_{j=0}^{N-1} \sum_{l=0}^{N-1} x_l \omega_N^{-j\cdot l}
\omega_N^{k\cdot j} \\
&=&
\frac{1}{N}\sum_{j=0}^{N-1} \sum_{l=0}^{N-1} x_l \omega_N^{j (k-l)} \\
&=& 
\frac{1}{N} \sum_{l=0}^{N-1} x_l \sum_{j=0}^{N-1} (\omega_N^{(k-l)})^j \\
&=&
\frac{1}{N} \sum_{l=0}^{N-1} x_l 
\left\{ \begin{array}{lcl}
\frac{1-(\omega_N^{(k-l)})^N}{1-\omega_N^{(k-l)}} & \mbox{si} &
\omega_N^{(k-l)}\neq 1 \\
N & \mbox{si} & \omega_N^{(k-l)} = 1 
\end{array} \right.
\end{eqnarray*}
Or si $\omega_N^{k-l}=e^{2i\pi(k-l)/N}=1$ si et seulement si $k=l$
d'où le résultat.

{\bf Propriété}\\
La transformée de Fourier discrète d'une suite réelle vérifie
$y_{N-k}=\overline{y_k}$.\\
La preuve est immédiate en appliquant la définition.

Un des intérêts de la DFT est de mettre en évidence rapidement
d'éventuelles périodicités de $x$ divisant $N$. Plus précisément
soit $j$ est un entier divisant $N$. Considérons une suite réelle $x$ dont
la DFT $y$ est nulle sauf $y_l$ et $y_{N-l}$. Par linéarité, on
peut se ramener à 2 cas $y_l=y_{N-l}=1$ et $y_l=i, y_{N-l}=-i$. Dans
le premier cas, on obtient $x_k=\omega_N^{lk}+\omega_N^{-lk}=2\cos(2\pi kl/N)$,
dans le deuxième cas, on obtient $x_k=-2\sin(2\pi kl/N)$, qui sont périodiques
de période $N/l$.

Réciproquement, si $x$ a comme période $T=N/l$, alors en posant
$j=T m + r$ avec $m\in[0,l[$ et $r\in[0,T-1]$, on a
$x_j=x_r$ donc~:
\begin{eqnarray*}
 y_k &=&\sum_{j=0}^{N-1} x_j \omega_N^{-k\cdot j} \\
&=&\sum_{m=0}^{l-1} \sum_{r=0}^{T-1} x_r \omega_N^{-k (T m+r)} \\
&=&\sum_{r=0}^{T-1} x_r \sum_{m=0}^{l-1} \omega_N^{-k (T m+r)} \\
&=& \sum_{r=0}^{T-1} x_r \omega_N^{-kr} \sum_{m=0}^{l-1}
(\omega_N^{-kT})^m \\
&=& \sum_{r=0}^{T-1} x_r \omega_N^{-kr} 
\frac{1-(\omega_N^{-kT})^l}{1-\omega_N^{-kT}}
\end{eqnarray*}
si $\omega_N^{-kT} \neq 1$. Comme $(\omega_N^{-kT})^l=\omega_N^{-klT}=1$,
$y_k=0$ si $kT=kN/l$ n'est pas un multiple de $N$. 
Finalement si $k$ n'est pas un multiple de $l$, alors $y_k=0$.

Voyons maintenant le cas de ``pseudo-périodes'', supposons donc
que $x$ est périodique de période $N$ mais que de plus pour un $T>0$ 
quelconque (ne divisant pas forcément $N$), on ait
\[ x_{j+T}=x_j, \quad \forall j \in[0,N-T] \]
On peut refaire le raisonnement ci-dessus, modulo des erreurs.
plus précisément~:
\[ y_k - \sum_{j=N}^{\mbox{ceil}(N/T)T}x_j \omega_N^{-k\cdot j} 
= \sum_{m=0}^{\mbox{ceil}(N/T)} \sum_{r=0}^{T-1} x_r \omega_N^{-k (T m+r)} \]
On calcule donc $y_k$ à une erreur de $\mbox{ceil}(N/T)T-N$ termes majorés
par $|x_j|$ près.
Et le membre de droite vaudra~:
\[ \sum_{r=0}^{T-1} x_r \omega_N^{-kr} 
\frac{1-\omega_N^{-k\mbox{ceil}(N/T) T}}{1-\omega_N^{-kT}}
\]
Le module de la fraction est égal à
\[ |\frac{\sin(\pi k \mbox{ceil}(N/T) T/N)}{\sin(\pi k T/N)} |
= |\frac{\sin(\pi k (\mbox{ceil}(N/T) T/N-1))}{\sin(\pi k T/N)} |
\]
il est petit si $k$ n'est pas proche d'un multiple de $\mbox{ceil}(N/T)$.
Par exemple, prenons $N=2^16=65536$ et $T \approx N/10 =6554$. Dans
ce cas $\mbox{ceil}(N/T)T=10 \times 6554=65540$, il y a donc une erreur
de 4 termes sur le calcul de $y_k$. Si $k$ n'est pas proche d'un multiple
de 10, on doit trouver $y_k$ proche de 0 relativement à la valeur des $|x_j|$.

Les périodes et pseudo-périodes de $x$ correspondent donc aux
valeurs de $y_k$ grandes par la règle $k$ * période =$N$.

\subsection{La transformée de Fourier rapide}
Le calcul de la DFT est relativement lent, il nécessite de l'ordre
de $N^2$ opérations, car il revient à calculer la valeur du polynôme
de degré $N-1$:
\[ P(X)=\sum_{j=0}^{N-1} x_j X^j \]
aux $N$ points $1$, $\omega_N$, ..., $\omega_N^{N-1}$ 
(on a $y_k=P(\omega_N^k)$).
Mais si $N$ est une puissance de 2, on peut calculer de manière plus
astucieuse et réduire le nombre d'opérations à un ordre $N \ln(N)$.
En effet $N=2M$,
on d\'ecoupe $P$ en 2 parties de m\^eme longueur~:
\[ P(x)=x^M Q(X) + R(X) \]
on a alors
\[ P(\omega^{2k})= (Q+R) ((\omega^2)^{k}), \quad
P(\omega^{2k+1})= (-Q +R)_\omega ((\omega^2)^k) \quad
S_\omega(x)=\sum s_k w^k x^k
\]
On est donc ramen\'e \`a deux additions de 2 polyn\^omes de degr\'e $M$,
une multiplication coefficient par puissances ($4M$ op\'erations),
et au calcul des deux DFT de $Q+R$ et $R-Q$. Si $N=2^n$ on
v\'erifie que cela n\'ecessite $O(n2^n)$ op\'erations, donc $N\ln(N)$
opérations.
On appelle alors FFT cette méthode de calcul (DFT=FFT si $N2^n$).

\subsection{Application à la température}
Les enregistrements ne sont pas faits à des dates régulièrement espacées.
Pour se ramener à une suite, il faut interpoler la valeur de la
température à une date donnée en fonction de la température aux
dates les plus proches de l'échantillon. Le calcul de la FFT
de l'écart de température avec la température actuelle sur les 400 000
dernières années donne le graphe figure \ref{fig:vostok_fft2}
\begin{figure}[h] \label{fig:vostok_fft2}
\includegraphics[width=\textwidth]{vostok_fft2}
\caption{Norme de la FFT de la temp\'erature de Vostok}
\end{figure}
On observe que la FFT $y_k$ est de module relativement grand pour $k=4$,
$k=10$ et $k$ proche de 20. La valeur de $k=4$ correspond à une périodicité
de 100 000 ans, celle pour $k=10$ à 40 000 ans. Il est plus difficile
de tirer des conclusions du maximum relatif pour $k$ proche de 20
car 20 est un multiple des périodes précédentes 4 et de 10.

On peut comparer ces 2 périodes tirées de l'enregistrement numérique
des périodes de variation des paramètres orbitaux de la Terre~:
\begin{itemize}
\item l'excentricité de la Terre varie avec des périodes principales
de l'ordre de 100 000 et 400 000 années
\item l'obliquité de la Terre varie avec une période de 41 000 ans.
\item la précession des équinoxes a des périodicités principales
de 19 000 et 23 000 ans.
\end{itemize}
L'excentricité et l'obliquité ont donc des périodicités que l'on
retrouve dans les enregistrements de la température du passé
ce qui confirme l'influence supposéee de ces paramètres
sur le climat de la Terre.

La figure \ref{fig:vostok_fft3} montre la reconstruction par FFT
inverse en tronquant la FFT aux 30 premiers entiers (et sym\'etriquement)
\begin{figure}[h] \label{fig:vostok_fft3}
\includegraphics[width=\textwidth]{Temperature_vs_fft_30}
\caption{Reconstruction de la temp\'erature en tronquant la FFT aux 30 premiers entiers}
\end{figure}

\pagebreak

\section{L'effet de serre}
Nous nous sommes intéressés jusqu'à présent à la quantité d'énergie
qui arrive du Soleil, sans regarder ce qui peut se passer entre le sommet
de l'atmosphère et le sol. Or la composition chimique de l'atmosphère
a une influence sur le piégeage des radiations solaires, donc
sur le climat.
Les molécules des gaz à effet de serre ont des fréquences propres de vibrations
qui sont proche des fréquences de radiation émises par la Terre
et les rend plus ou moins opaques (alors qu'ils sont quasiment
transparents dans le domaine de fréquence des radiations provenant du Soleil).
Il s'agit sur Terre principalement du CO2, CH4, des chlorofluorocarbone,
de l'ozone (O3) et d'oxydes d'azote. 
Leur concentration a varié dans le passé, de manière naturelle
corrélée avec la température 
(cf. la section \ref{sec:corr}). Elle varie aujourd'hui du fait
de l'homme, et les scientifiques élaborent des modèles pour prévoir
la réaction du climat à ces variations.

\subsection{La définition du climat.}
Le domaine des mathématiques concerné sont les probabilités
et les statistiques.
Il s'agit de savoir quel sens attribuer à une phrase comme ``la hausse
de température de la Terre au 20ème siècle a été de 0.6 degré''.
Il s'agit aussi de savoir pourquoi on considère que 30 ans
est un nombre raisonnable d'années
pour définir le climat en un lieu donné. On peut faire
l'hypothèse que le temps un jour donné n'est pas indépendant
de celui du jour précédent (mémoire du temps), mais qu'il est
indépendant de celui du même jour l'année précédente~:
le temps suit une loi qui dépend de la saison. Cette loi subit une
lente dérive au cours du temps (variations du climat).
Un des résultats importants qui justifie la période de 30 années
est le théorème central-limite. Nous allons rappeler quelques
définitions avant de l'énoncer.

\subsubsection{Lois}
Lorsqu'on n'est pas capable de prédire de manière déterministe
plusieurs événements du même type, on les modélise souvent en disant 
qu'ils dépendent
du hasard mais pas de manière anarchique, ils suivent une loi. Ainsi,
le résultat du lancer d'un dé à 6 faces non pipé ne peut pas être
prédit, mais en revanche il suit une loi (uniforme), chaque face
a autant de chances d'apparaitre qu'une autre. La valeur d'un lancer
est une variable aléatoire à valeurs dans $\{1,2,3,4,5,6\}$ 
suivant la loi uniforme (la probabilité
d'observer chacune des faces de 1 à 6 est de 1/6). Il s'agit ici
d'une loi discrète, on ne peut observer que 6 valeurs. L'analogue
existe pour une variable aléatoire dont la valeur est réelle
(par exemple la température), on utilise alors une densité
de probabilité, c'est une fonction $f$ de $\R$ à valeurs
dans $\R^+$ telle que la probabilité d'observer la variable
aléatoire dans un intervalle $[a,b]$ soit $\int_a^b f(x) dx$.
L'intégrale sur $\R$ de $f$ est donc égale à 1.
Par exemple pour la loi uniforme sur $[1,6]$, $f$ est constante
et vaut $1/5$ sur $[1,6]$ et 0 ailleurs. 
Un autre exemple important est la loi normale développée plus bas.

\subsubsection{Ind\'ependance}
Lorsqu'on effectue plusieurs observations successives, une observation
peut d\'ependre de l'observation pr\'ec\'edente ou non. Par exemple,
la temp\'erature observ\'ee un jour donn\'e est en
g\'en\'eral proche de celle qui sera observ\'ee le lendemain. Mais ell
est ind\'ependante de celle observ\'ee l'ann\'ee pr\'ec\'edente \`a
la m\^eme date. De m\^eme qu'un jet de d\'e donne un r\'esultat
ind\'ependant du jet de d\'e pr\'ec\'edent. Math\'ematiquement
parlant, l'ind\'ependance se d\'efinit par le fait que la probabilit\'e
d'observer les \'ev\'enements $A$ et $B$ est le produit des
probabilit\'es d'observer $A$ et d'observer $B$.

\subsubsection{Moyenne et écart-type}
Il y a en fait deux notions de moyenne, la première
consiste à faire la moyenne d'un échantillon ($n$ observations
indépendantes) de la variable aléatoire, la deuxième, appelée espérance
de la variable aléatoire, est la moyenne
de la valeur de la variable aléatoire pondérée par la probabilité
d'observer chacune des valeurs (intuitivement on pense
s'en rapprocher lorsqu'on fait la moyenne d'un échantillon assez grand).
\begin{eqnarray*}
 \overline{X} &= &\frac{1}{n}\sum_{k=1}^n X_k (\mbox{\'echantillon}),
 \\
\mu=E(X)&=&= \sum x_i P(X=x_i) (\mbox{loi discrète}) \\
\mu=E(X) &=& \int x f(x) dx (\mbox{loi continue}) 
\end{eqnarray*}
Le lien entre les deux notions de moyenne est~:
\[ E(\overline{X}) = E(X) \]
L'espérance de la moyenne est égale à la moyenne, lorsqu'on ne
connait pas la loi d'une variable aléatoire, on peut estimer
son espérance de manière non biaisée par la moyenne de l'échantillon.
Intuitivement, si les observations sont ind\'ependantes et plus il
y en a, plus la moyenne devrait s'approcher de l'esp\'erance.
(on verra avec le th\'eor\`eme central-limite qu'on en sait plus
sur la distribution de la moyenne lorsque le nombre d'observations
ind\'ependantes pour calculer la moyenne est grand). 
Pour quantifier cela, il faut pouvoir \'evaluer l'\'ecart \`a la
moyenne.
La moyenne ne donne en effet pas d'informations sur les
fluctuations de part et d'autre de la moyenne. Pour cela, on 
utilise l'\'ecart au carr\'e \`a la moyenne, dont on fait la moyenne
puis on en prend la racine carr\'ee~: c'est l'écart-type.
On peut calculer l'écart-type pour un échantillon
(par rapport à la moyenne de cet échantillon ou par rapport à
l'espérance de la loi si elle est connue) 
et l'écart-type d'une loi (N.B. on appelle
variance le carr\'e de l'\'ecart type).
\begin{eqnarray*}
 \sigma^2(X) \ (\mbox{\'echantillon})&= 
&\frac{1}{n}\sum_{k=1}^n (X_k-\overline{X})^2
= \frac{1}{n}\sum_{k=1}^n X_k^2 - \overline{X}^2 \\
\sigma^2 \ (\mbox{loi discrète}) &=& E( (X-E(X))^2 )
= \sum (x_i-E(X))^2 P(X=x_i) \ dx \\
&=& E(X^2)-E(X)^2 =
= \sum x_i^2  P(X=x_i)\  - E(X)^2 \\
\sigma^2 \ (\mbox{loi continue}) &=& E( (X-E(X))^2 ) 
= \int (x-E(X))^2 f(x) \ dx \\
&=& E(X^2)-E(X)^2= \int x^2  f(x)\  dx - E(X)^2 
\end{eqnarray*}
On observe que si $m$ est une constante~: 
\[ \sigma^2(X-m)=\sigma^2(X) \]
C'est faux si $m$ n'est pas constant, on a
tout de m\^eme la relation
\[ \sigma^2(X+Y)=\sigma^2(X)+\sigma^2(Y) \]
lorsque $X$ et $Y$ sont ind\'ependants.

Contrairement au cas de la moyenne, il y a un ajustement à faire
lorsqu'on estime l'écart type d'une loi à partir d'un échantillon
si on ne connait pas la moyenne de la loi. En effet, on estime alors
la moyenne de la loi \`a partir de la moyenne
de l'\'echantillon, ce qui conduit \`a sous-estimer la somme des
carr\'es des \'ecarts \`a la moyenne.
Lorsque les $n$ valeurs de $X_i$ de l\'echantillon sont
ind\'ependants, on a~:
\[ E(\sigma^2 \mbox{ echantillon })= \frac{n}{n-1} E(\sigma^2(X)) \]
En effet~:
\begin{eqnarray*}
E(\sigma^2(X)) &=& 
E( \frac{1}{n}\sum_{k=1}^n X_k^2 - \overline{X}^2 ) \\
&=&
\frac{1}{n}\sum_{k=1}^n E(X_k^2) 
- E( (\frac{1}{n}\sum_{k=1}^n X_k)^2 )\\
&=&
\frac{1}{n} \sum_{k=1}^n E(X^2)
- \frac{1}{n^2} E( \sum_{k=1}^n X_k^2 + 
\sum_{k=1}^n \sum_{j=1, j \neq k}^n X_k X_j ) \\
&=&
E(X^2) - \frac{1}{n^2} (n E(X^2) + n(n-1) E(X)^2) \\
&=&
\frac{n-1}{n} (E(X^2)-E(X)^2) \\
&=&\frac{n-1}{n} \sigma^2
\end{eqnarray*}
On utilise donc
\[ \frac{1}{n-1}\sum_{k=1}^n (X_i-\overline{X})^2 \]
au lieu de 
\[ \frac{1}{n}\sum_{k=1}^n (X_i-\overline{X})^2 \]
pour estimer l'écart-type d'une loi \`a partir d'un \'echantillon.


La définition de ces deux paramètres 
va nous permettre d'énoncer deux r\'esultats sur la dispersion~:
\begin{itemize}
\item
\[ \sigma \neq 0 \Rightarrow P(|X-\overline{X}| \geq \mu \sigma) 
\leq \frac{1}{\mu^2}\]
En effet, on minore la variance de $X$ en restreignant aux $x_i$
v\'erifiant la condition ci-dessus~:
\[ \sigma^2 \geq \sum_{x_i/ |x_i-\overline{X}| \geq \mu \sigma} (x_i-\overline{X})^2
P(X=x_i) \geq \mu^2 \sigma^2 P(|X-\overline{X}| \geq \mu \sigma)
 \]
\item Si on a $n$ variables ind\'ependantes, alors la variance
de la somme de ces $n$ variables est la somme des variances de
chaque variable. En particulier si elles ont le m\^eme \'ecart-type $\sigma$,
l'\'ecart-type de la moyenne est $\sigma/\sqrt{n}$.
\end{itemize}

La d\'efinition de la loi normale (cf. infra) permet de donner un
sens quantitatif au fait qu'une moyenne d'observations ind\'ependantes
de m\^eme loi tend vers l'esp\'erance de la loi (c'est 
le théorème central-limite qui sera admis).

\subsubsection{La loi normale}
Il s'agit d'une loi continue, définie par la fonction
\[ N_{\mu,\sigma^2}(x)= \frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{(x-\mu)^2}{2\sigma^2}} \]
o\`u la constante devant l'exponentielle est déterminée par la condition 
$\int_{-\infty}^{+\infty} N_{\mu,\sigma^2}(x) dx=1$. 

On montre que la loi normale de paramètres $\mu$ et $\sigma^2$ a pour espérance
$\mu$ et variance $\sigma^2$. Pour l'esp\'erance, on fait le changement de
variable $x=\mu +y$, on obtient
\[ E(X) = \int   \frac{1}{\sigma \sqrt{2 \pi}} x
e^{-\frac{(x-\mu)^2}{2\sigma^2}} 
= \frac{1}{\sigma \sqrt{2 \pi}} \int (y+\mu)e^{-\frac{y^2}{2\sigma^2}}
= \mu +  \frac{1}{\sigma \sqrt{2 \pi}} \int y
  e^{-\frac{y^2}{2\sigma^2}}
= \mu
\]
la derni\`ere \'egalit\'e provenant de l'imparit\'e de la fonction
\`a int\'egrer.
Pour la variance, il faut faire une int\'egration par parties
pour enlever le terme en $x^2$.

Plage de normalité au niveau de confiance $\alpha \in ]0,1[$: il s'agit
d'un intervalle centr\'e autour de la moyenne 
$I=]\mu - t \sigma, \mu +t \sigma[$, 
tel que la probabilit\'e d'\^etre dans $I$ soit de $\alpha$.
Le calcul pratique se fait avec la fonction sp\'eciale erf (error function)
d\'efinie par 
\[ \mbox{erf}(x)= \frac{2}{\sqrt{\pi}} \int_0^x e^{-u^2} \ du\]
On montre par un changement de variable que $t$ v\'erifie l'\'equation
\[ \mbox{erf}(\frac{t}{\sqrt{2}})=\alpha\]
Ainsi, comme erf$(\sqrt(2))=0.954...$, on utilise souvent la plage
de normalit\'e au niveau de confiance de 95\% qui est
approximativement $[\mu-2\sigma,\mu+2\sigma]$ (la probabilit\'e
d'\^etre dans cet intervalle est l\'eg\`erement sup\'erieure
\`a 95\%)

\subsubsection{Théorème central-limite}
Ce théorème dit que la variable
aléatoire obtenue en faisant la moyenne de $N$ variables
aléatoires indépendantes et identiquement distribuées (de même loi inconnue)
suit une loi gaussienne de moyenne la moyenne de la loi inconnue
et d'écart type, l'écart type de la loi inconnue divisé par $\sqrt{N}$.
On estime alors la moyenne de la gaussienne par la moyenne $M$ des
observations, et l'écart type $S$ de la gaussienne par l'écart type
des observations divisé par $\sqrt{N}$. Pour une variable aléatoire
suivant une loi gaussienne
de moyenne $m$ et d'écart type $\sigma$, la probabilité de tirer
une valeur en-dehors de l'intervalle $[m-2\sigma,m+2\sigma]$ est
de 5\%, on peut estimer que la moyenne ``réelle'' pour le climat
est dans l'intervalle $[M-2S/\sqrt{N},M+2S/\sqrt{N}]$. Si $S/\sqrt{N}$
est suffisamment petit, on a une définition raisonnable du climat.
Plus $N$ est grand, plus $S/\sqrt{N}$ sera petit, mais d'une
part on ne dispose pas de séries de mesure du temps sur des périodes
très grandes, d'autre part la loi définissant le climat d'un lieu
évolue au cours du temps, on ne peut donc pas prendre $N$ trop grand.

On peut faire le calcul de $S/\sqrt{N}$ pour certaines stations météo
(les chiffres sont disponibles gratuitement sur Internet pour de
nombreuses stations météo américaines), on trouve pour $N=30$ des
valeurs de l'ordre de 0.1 degré. Lorsqu'on moyenne sur le globe,
les écarts-type diminuent au moins d'un ordre de grandeur, ce qui permet
de parler de manière raisonnable d'augmentation de température
(qui est largement supérieure aux incertitudes dues aux $S/\sqrt{N}$).

\subsection{Corrélation et régression}

\subsubsection{Corrélation entre deux variables aléatoires.}
Si deux variables $X$ et $Y$ sont indépendantes, alors 
\[ E(X Y ) = E(X) E(Y)\]
Dans le cas général, cette relation n'est pas vraie, on définit la
covariance entre $X$ et $Y$ par
\[ C(X,Y) = E(X Y ) - E(X) E(Y)\]
Cette valeur n'est pas ``normalisée'' (si $X$ et $Y$ ont des unités
physiques par exemple ce nombre a une dimension), 
on définit donc la corrélation par
\[ c(X,Y) = \frac{C(X,Y)}{\sigma(X) \sigma(Y)} \]
ce nombre est sans unité. On montre qu'il appartient à l'intervalle
$[-1,1]$. Si $Y=X$ 
alors $c=1$, si $Y=-X$ alors $c=-1$.
Si les variables sont indépendantes, alors $c=0$. Mais
la réciproque est fausse.  On utilise toutefois ce calcul
pour estimer si deux variables sont corrélées ou indépendantes.
Plus précisément, si deux variables paraissent corrélées, on peut
essayer d'exprimer l'une en fonction de l'autre, c'est ce qu'on
va faire dans la section suivante.

\subsubsection{Régression linéaire avec un mod\`ele lin\'eaire}
\label{sec:corr}
Soit $(X_i)_{1 \leq i \leq n}$ et $Y_i$ deux séries statistiques.
On cherche à savoir s'il est raisonnable de prévoir $Y_i$ en fonction
de $X_i$, avec par exemple une relation linéaire $Y_i = a + b X_i$.
Pour cela on minimise $f(a,b)=\sum (a + b X_i -Y_i)^2$ par rapport à $a$
et $b$, ce qui donne les équations :
\[ \left\{ \begin{array}{rcl}
\frac{\partial f}{\partial a} =0&=&2 \sum_{i=1}^n (a+bX_i-Y_i)\\
\frac{\partial f}{\partial b} =0&=&2 \sum_{i=1}^n X_i(a+bX_i-Y_i)
\end{array} \right. \]
équivalentes à :
\begin{equation} \label{eq:reg}
 \left\{ \begin{array}{llcl}
a \sum_{i=1}^n + & b \sum_{i=1}^n X_i &=& \sum_{i=1}^n Y_i\\
a \sum_{i=1}^n X_i + & b \sum_{i=1}^n X_i^2 &=& \sum_{i=1}^n X_i Y_i
\end{array} \right. 
\end{equation}
Si on note $U$ la matrice de première colonne remplie de 1, et de
deuxième colonne remplie par $X_i$ (en ligne $i$),
\[ U=\left( \begin{array}{cc}
1 & X_1 \\
1 & X_2 \\
... & ... \\
1 & X_n \end{array} \right) \]
et si on note $U^*$ la transposée de $U$,
\[ U^*=\left( \begin{array}{cccc}
1 & 1 & ... & 1 \\
X_1 & X_2 & ... & X_n
\end{array} \right) \]
on observe que 
\[ U^* U = \left( \begin{array}{cc}
\sum_{i=1}^n 1 & \sum_{i=1}^n X_i \\
\sum_{i=1}^n X_i & \sum_{i=1}^n X_i^2
\end{array} \right) \]
le système (\ref{eq:reg}) devient
\[ U^* U \left( \begin{array}{c} a \\ b \end{array} \right)
= \left( \begin{array}{c} \sum_{i=1}^n Y_i \\ \sum_{i=1}^n X_iY_i 
\end{array} \right) =
U^* Y, \quad \mbox{où }
Y=\left( \begin{array}{c} Y_1 \\ Y_2 \\ ... \\ Y_n \end{array}
\right)
\]
Donc la solution est donnée par
\[ A=\left( \begin{array}{c} a \\ b \end{array} \right)
= \mbox{inv}(U^* U) (U^* Y)\]
En pratique avec Xcas, on place dans un tableur deux colonnes,
une avec les $X_i$, une avec les $Y_i$, on rajoute une colonne
de 1 avant la colonne des $X_i$, 
on sélectionne les colonnes de 1 et de $X_i$, on retourne
dans l'historique (hist), on écrit \verb|U:=| puis on copie
(bouton du milieu de la souris ou touche copie) puis Entrée. On revient
dans mtrw, on sélectionne la colonne des $Y_i$, on repasse dans
hist, on écrit \verb|Y:=| et on recopie puis Entrée.
On effectue alors \verb|A:=inv(tran(U)*U)*(tran(U)*Y)|.

On peut ensuite juger de la régression en comparant $f(a,b)/n$
avec la variance de $Y_i$ ou en comparant $\sqrt{f(a,b)/n}$ avec
l'\'ecart-type de $Y_i$. On voit que $f(a,b)$
est la somme des carrés des composantes du vecteur $U A - Y$,
donc le calcul pratique de $f(a.b)/n$ se fait par la formule
\verb|l2norm(U*A-Y)^2/size(Y)|, \`a comparer avec \verb|stddev(Y)^2|. 
Plus ce rapport est proche de 0, plus on explique les variations de $Y_i$ par
la régression linéaire.

On peut montrer que les coefficients de la droite sont donn\'es par~:
\[ b=\frac{E(XY)-E(X)E(Y)}{\sigma(X)^2}, \quad
a = \overline{Y} - b \overline{X} \]
et
\[ f(a,b) = n (1-r^2) \sigma(Y)^2, 
\quad r=\frac{E(XY)-E(X)E(Y)}{\sigma(X)\sigma(Y)} \]


On montre que cette méthode se généralise (dans sa formulation
matricielle), si on veut prévoir par
exemple $Y_i$ à partir de deux séries statistiques $x_i$ et $X_i$
avec une relation linéaire de la forme $Y_i=a+bx_i+cX_i$, alors
$(a,b,c)$ est obtenu en calculant $\mbox{inv}(U^* U) U^* Y$,
où cette fois la matrice $U$ possède 3 colonnes, la colonne de 1,
la colonne des $x_i$ et la colonne des $X_i$. Les calculs
pratiques se font de manière identique.

\subsubsection{Autres régressions lin\'eaires}
On peut aussi supposer que $Y$ dépend de manière non linéaire de $X$,
par exemple logarithmique, exponentielle, polynomiale, etc. mais
toujours de mani\`ere lin\'eaire des param\`etres \`a d\'eterminer et 
calculer ces param\`etres en utilisant la même méthode
(minimiser l'erreur quadratique). Ceci peut servir par exemple
à modéliser la consommation d'une denrée non renouvelable (par
exemple les combustibles fossiles) par une fonction logistique
ou par une gaussienne.

\subsubsection{Corrélations entre gaz à effet de serre et température}
Si on fait une corrélation température-CO2 sur les 400 000 dernières
années, on trouve un rapport $f(a,b)/n$ sur 
variance de 0.35, soit une qualité de corrélation de 0.65 
(65\%\ de la variance
de la température est expliquée par le taux de CO2) et un
coefficient d'augmentation de température de 9 degrés pour 100 ppm
de CO2. On peut aussi
faire une corrélation température en fonction du CO2 et du CH4,
on trouve alors 6 degrés pour 100 ppmv de CO2 et 1.5 degré pour 100
ppbv de CH4 avec une qualité de corrélation de 0.68.

Bien entendu, une corrélation statistique ne suffit pas à conclure
qu'une des variables dépend (au moins partiellement) de l'autre,
mais ici la physique nous montre que l'augmentation du taux de CO2
augmente les radiations solaires capturées au niveau du sol, la
corrélation permet de quantifier cet effet (calcul que l'on pourra comparer
à celui obtenu par d'autres méthodes).


\section{Modélisation numérique (climat, \'energie)}
Dans cette section, on pr\'esente diff\'erents types de mod\`eles 
math\'ematiques qui id\'ealisent la situation r\'eelle. Les principaux
types s'appliquant ici sont les mod\`eles discrets, les mod\`eles
continus se ramenant \`a des \'equations ou syst\`emes diff\'erentiels
et les mod\`eles de type \'equations aux d\'eriv\'ees partielles (qui
ne seront pas pr\'esent\'es ici).

\subsection{Exemple~: un modèle très simplifié du climat} \label{sec:simplifie}
On suppose ici que la Terre est assimilable à un océan (pas de glace
modifiant l'albédo) recevant un rayonnement solaire constant et émettant
un rayonnement dont la puissance est composée d'une partie
proportionnelle à $T^4$ (comme un corps noir) dont on retranche
le piégage du au gaz carbonique (proportionnel au log de la concentration
de CO2). La différence entre les deux rayonnements incident et renvoyés
sert à réchauffer l'eau. On obtient l'équation (o\`u le coefficient
6 et le ln pour le taux de CO2 proviennent de calculs de physique)
\[ \frac{dT}{dt} = k ( 6 \ln (CO2/280) - \sigma (T^4-T_0^4) )\]
Le calcul de $k$ donne $(86400 \times 365)/(3e6 \times 4.18e3) =0.0025$ 
degrés par an
(on chauffe l'eau sur 3000 mètres en moyenne avec une puissance de $1W/m^2$).
On peut prendre $T_0=288$ (15 degré celsius), et pour $\sigma$ la constante
de Stefan-Boltzman ($5.67e-8$, ceci aura tendance à sous-estimer
le phénomène car l'eau de mer n'est quand même pas un corps noir).
On peut prendre par exemple une variation du taux de 
CO2 constante et égale à 2 ppm par an, soit si $t$ est exprimée en années
\[ CO2=380+2(t-2005) \] 

On ne sait pas résoudre exactement ce type d'équation différentielle,
on utilise alors des méthodes de résolution approchée. La plus simple
consisterait dans notre cas à supposer la température du globe constante
pendant une année et à calculer la différence entre l'année en cours
et l'année suivante
\begin{equation} \label{eq:euler}
 T(t+1)=T(t) + k ( 6 \ln (CO2/280) - \sigma (T^4-T_0^4) )
\end{equation}
Il s'agit de la méthode d'Euler (avec un pas de 1 an), graphiquement
cela revient \`a tracer le graphe de la solution d'une équation différentielle
en suivant les tangentes.

L'avantage de cette discrétisation est qu'on peut modéliser de manière
plus réaliste la variation du taux de CO2.

\subsection{Exemple~: scénarios d'émission de CO2}

\subsubsection{CO2 et consommation d'énergie}
L'homme émet dans l'atmosphère des gaz à effet de serre principalement
par combustion des pétrole, gaz, charbon, également par la déforestation
et les émissions de méthane de l'agriculture. Un petit calcul
élémentaire permet d'établir la correspondance entre 1 ppm de CO2
(une partie de CO2 par million en volume) et 2 Gt (gigatonnes) de carbone.

Le mètre a été pendant longtemps défini comme 1/40 millionième
de la circonférence de la Terre, on a donc $2\pi R=4e7 m$ où
$R$ désigne le rayon de la Terre. On en déduit la surface de la Terre
\[ 4 \pi R^2= (4e7)^2/\pi=5e14 m^2 \]
La pression de l'air au sol est équivalente à une colonne d'eau
de 10m (ou encore 76cm de mercure dont la densité est 13). La masse
de l'atmosphère terrestre est donc égale à celle de 5e15m$^3$ d'eau
soit $5e18$kg. L'air est pour trois quarts 
composé d'azote ($N_2$ masse molaire $2\*14=28$), et pour le reste
d'oxygène ($O_2$ masse molaire $2\*16=32$), alors que le carbone
$C$ a une masse molaire de 12.
Comme il y a 1 atome de carbone par molécule de $CO_2$, on en déduit
que la masse de carbone correspondant à une
concentration de 1ppm de $CO_2$ est de~:
\[ \frac{12}{3/4*28+1/4*32} * 5e18 * 1e-6= 2.1 e12 kg \]
soit 2.1Gt de carbone.

La consommation mondiale de combustibles fossiles s'établit pour 2003
à 3.5 Gtep (giga-tonne équivalent pétrole) de pétrole, 2.3 Gtep de gaz
naturel et 2.4 Gtep de charbon. Avec un taux moyen en tonne
de carbone par tep de 0.856 pour le pétrole, 0.651 pour le gaz
et 1.123 pour le charbon, on a émis
\[ 3.5 \times 0.856 + 2.3 \times 0.651 + 1.123 \times 2.4
= 7.1 Gt C\]
ce qui correspond à 3.5 ppm de carbone. Auquel il faut ajouter
le carbone émis par déforestation, retrancher le CO2 absorbé par
l'océan et les écosystèmes, on a mesuré une augmentation de 2ppm.

\subsubsection{Sc\'enarios de consommation d'\'energie et d'\'emissions de CO2}
Les incertitudes sur la sensibilité du climat à une augmentation
fixée des taux de gaz à effet de serre sont assez faibles. Par contre
les prévisions de taux de gaz à effet de serre sont très incertaines
puisque par exemple la concentration de CO2 à la fin du siècle dépend
de la quantité de combustibles fossiles que nous aurons brulée.
Celle-ci est très incertaine, car~:
\begin{itemize}
\item nous ne connaissons pas précisément la quantité de pétrole,
gaz, charbon à notre disposition. Si comme certains géologues le
pensent, nous avons consommé la moitié des pétroles conventionnels,
la consommation de pétrole va diminuer dans les années à venir,
au profit d'autres énergies fossiles plus ou moins polluantes
en CO2 par tep, d'énergies renouvelables ou d'économies d'énergie.
\item nous ne connaissons pas l'avenir des filières nucléaires 
\item nous ne pouvons pas prédire les conditions économiques
et démographiques qui prévaudront dans 20 ou 50 ans
\end{itemize}
Plusieurs modèles existent donc, par exemple les scénarios du GIEC
ou d'experts (par exemple ceux dérivés des projections de l'ASPO). 
Quelques principes mathématiques sous-tendent ces modèles~:
\begin{itemize}
\item Des modèles de type suite géométrique lorsque la croissance
n'est pas contrainte ou pour certains modèles simples de dépletion
d'une quantité (on suppose que la
quantité donnée va évoluer comme une suite géométrique
de raison >1 dans le premier cas, de raison < 1 dans le
deuxième cas).
\item Des modèles malthusiens lorsqu'on s'intéresse à une ressource
dont la quantité va limiter l'accroissement (de population, 
de consommation, ...), par exemple logistique~: on suppose
que la variation relative de la quantité est linéaire, avec
un coefficient directeur négatif pour limiter la croissance. Ce type de
modèle est utilisé par exemple pour estimer la quantité ultime récupérable
de pétrole dans le monde.
\item Des modèles linéaires lorsque plusieurs quantités sont
en interaction, on suppose que la variation de consommation
d'une quantité donnée dépend linéairement de toutes les quantités.
On résout alors un système différentiel linéaire.
\item etc.
\end{itemize}

\subsection{Mod\`eles discrets.}
Les mod\`eles les plus utilis\'es sont les suites arithm\'etiques,
les suites g\'eom\'etriques, et plus g\'en\'erallement les
suites r\'ecurrentes. Historiquement, c'est Malthus qui soutint la th\`ese
que la production croit comme une suite
arithm\'etique, alors que la population croit comme
une suite g\'eom\'etrique, pr\'edisant de gros probl\`emes
lorsque la suite g\'eom\'etrique d\'epassait la
suite arithm\'etique. C'est Verhulst qui proposa alors le
mod\`ele logistique pour am\'eliorer la mod\'elisation (cf. infra
dans le cas continu).

Plus g\'en\'erallement, les suites r\'ecurrentes $u_{n+1}=f(u_n)$
permettent de mod\'eliser des ph\'enom\`enes o\`u la relation
ne d\'epend pas explicitement du temps, il est alors int\'eressant
d'\'etudier les positions d'\'equilibre aussi appel\'es points fixes,
et leur stabilit\'e (si $u_0$ est proche d'un point fixe $l$,
la suite $u_n$ converge-t-elle vers $l$~?). On peut aussi \'etudier
la sensibilit\'e de l'\'equilibre \`a une variation des param\`etres
d\'efinissant la fonction $f$ (d\'eplacement de l'\'equilibre,
passage d'un \'equilibre stable \`a un \'equilibre instable).
Le r\'esultat principal sur la stabilit\'e est que si $|f'(l)|<1$ 
alors l'\'equilibre est stable (on peut le comprendre intuitivement
en faisant une repr\'esentation graphique de la suite en dimension 1).

Si $f$ d\'epend explicitement du temps (i.e. de $n$), par exemple
pour (\ref{eq:euler}), il n'y a pas
de r\'esultats simples. On peut r\'esoudre explicitement
quelques cas particuliers, par exemple $u_{n+1}=Au_n+P(n)$, o\`u $A$ est
une matrice constante et $P$ un polyn\^ome.

\subsection{Mod\`eles continus~: \'equations et systèmes différetiels}
On a vu dans les sections précédentes que de nombreux modèles se
ramenent à la résolution d'équations différentielles ou de
systèmes différentiels~:
\[ \frac{dX}{dt} = F(X(t),t) \]
o\`u $X$ est la quantit\'e recherch\'ee et $t$ le temps.
La plupart du temps ces équations
et systèmes ne peuvent être résolues que de manière numérique, mais
certains cas particuliers importants peuvent être résolus exactement.
Nous verrons le cas du modèle logistique comme exemple d'équation
r\'esoluble (qui est utilis\'e comme mod\`ele de production d'\'energie non
renouvelable) et les systèmes différentiels affines (qui aident \`a
comprendre une am\'elioration du mod\`ele tr\`es simple, \ref{sec:simplifie},
tenant compte des interactions entre temp\'erature, CO2 et alb\'edo). 
Pour les cas non résolubles, on essaie de déterminer le comportement
qualitatif de l'équation à partir de cas connus, en général en
linéarisant.
On mettra en évidence sur ces systèmes
la notion d'équilibre stable et instable, et la notion de 
rétroaction.

\subsubsection{Le modèle logistique}
Le modèle exponentiel, historiquement propos\'e par Malthus
et encore tr\`es (trop?) utilisé par les économistes, correspond
 à l'équation différentielle~:
\[ \frac{dy}{y dt} = b\]
o\`u $y(t)$ d\'esigne la quantit\'e totale extraite au temps $t$.
La solution de cette \'equation est une exponentielle qui tend
vers $+\infty$ lorsque $t$ tend vers l'infini, ce qui n'est pas
possible pour une ressource non renouvelable.
Pour avoir un modèle plus réaliste qui tienne compte d'une ressource finie,
Verhulst (1840) passe d'un taux d'accroissement constant à un taux linéaire
et qui diminue lorsque $y$ augmente
\[ \frac{dy}{y dt} = -ay+b =-a(y-c), \quad a>0, c=\frac{b}{a}\]
La solution de cette \'equation est appel\'ee fonction logistique.
On s\'epare les variables~:
\[ dt=-\frac{dy}{ay(y-c)} \]
on d\'ecompose en \'el\'ements simples~:
\[ a dt = \frac{dy}{c} \left( \frac{1}{y} - \frac{1}{y-c} \right) \]
soit~:
\[ b dt =  \left( \frac{dy}{y} - \frac{dy}{y-c} \right) \]
d'o\`u en int\'egrant~:
\[ bt+C= \ln(|\frac{y}{y-c}|)\]
puis, en enlevant les valeurs absolues et en posant $K=\pm e^{C}$~:
\[ \frac{y}{c-y}=Ke^{bt} \Rightarrow y(1+Ke^{bt})=Kce^{bt}\]
Donc $K>0$, on peut poser $K=e^{-bt_0}$ et finalement~:
\[ y = c \frac{e^{b(t-t_0)}} {1+e^{b(t-t_0)}}=c\frac{1}{1+e^{-b(t-t_0)}}\]
On peut \'etudier le comportement de la fonction $y$ en fonction du temps,
en $t=-\infty$, la limite de $y$ est nulle, la fonction croit pour $t$
n\'egatif comme une exponentielle
\[ y \approx_{t \rightarrow -\infty} c e^{b(t-t_0)} \]
(avec le m\^eme taux de croissance que lorsque $a=0$) et est convexe.
Pour $t=t_0$, $y=c/2$, la fonction a un point d'inflexion et devient concave
pour rejoindre sa limite $c$ en $t=+\infty$.
Si on calcule $y'$ qui mod\'elise la production instantann\'ee, on
obtient~:
\[ y' = bc \frac{e^{-b(t-t_0)}} {(1+e^{-b(t-t_0)})^2} = 
\frac{bc}{e^{b(t-t_0)}+2+e^{-b(t-t_0)}} \]
On observe que $y'$ a un graphe en forme de cloche avec deux limites
nulles aux extr\'emit\'es et un maximum en $t=t_0$~:\\
\includegraphics[width=\textwidth]{logistique}


En pratique, pour mod\'eliser $y$ on effectue une r\'egression lin\'eaire
du quotient de la production annuelle par la production cumul\'ee
(par exemple pour le p\'etrole depuis 1900) en ordonn\'ee (ce
qui repr\'esente $dy/ydt$) en fonction du temps.
Un grand int\'er\^et de ce mod\`ele est qu'il permet de pr\'edire la 
quantit\'e totale que l'on pourra extraire.
Sa principale faiblesse est la sensibilit\'e de cette estimation en fonction
de la qualit\'e des donn\'ees initiales et des irr\'egularit\'es
(par exemple chocs p\'etroliers). Pour la production
de p\'etrole, on peut un peu l'am\'eliorer en
mod\'elisant plusieurs cycles de p\'etrole donc en additionnant
plusieurs courbes logisitiques (p\'etrole on-shore facile,
p\'etrole offshore, deepwater, ultra-deepwater, polaire, p\'etroles
lourds).

\subsubsection{Systèmes différentiels linéaires}
On s'intéresse plus générallement à des systèmes affines:
\begin{equation} \label{eq:difflin} 
\frac{dX}{dt}=AX+C
\end{equation}
où $X(t)=(x_1,...,x_n)$ est un vecteur de dimension $n$, $A$
une matrice $n,n$ et $C=(c_1,...,c_n)$ est un vecteur constant
Plus g\'en\'erallement, on peut résoudre exactement ce type de syst\`eme
lorsque $C$ ne dépend que du temps.

Pour les syst\`emes qui nous int\'eressent, on suppose en
g\'eneral que les coefficients sur la diagonale sont n\'egatifs,
de sorte que si les interactions entre param\`etres \'etaient nulles,
la solution du syst\`eme convergerait.
Pour une mod\'elisation r\'ealiste, $A$ n'est pas diagonale, 
les coefficients non diagonaux traduisent
des r\'etroactions entre les param\`etres, ces r\'etroactions peuvent
\^etre positives si ces coefficients sont positifs
(l'interaction amplifie l'action de la variation originale du param\`etre)
ou n\'egatives dans le cas contraire.

Si la matrice $A$ est inversible, il existe une unique position
d'équilibre, en effet si $dX/dt=0$, on a $X=X_e=-A^{-1} C$ et réciproquement
cette valeur de $X$ est solution de (\ref{eq:difflin}).
Pour une condition initiale $X_0$ quelconque, on  peut alors
se demander si $X(t)$ va se rapprocher de la position d'équilibre $X_e$
(\'equilibre stable) ou non (instabilit\'e si on s'\'eloigne de
l'\'equilibre). Si le syst\`eme est stable, on peut aussi se demander 
comment la position d'équilibre varie lorsque $C$ varie 
(par exemple pour les param\`etres du climat si le forçage radiatif 
est modifi\'e par variation des param\`etres orbitaux de la Terre).

Il ne faut pas confondre la question de la stabilit\'e du syst\`eme
et l'existence de r\'etroactions positives ou n\'egatives. Un syst\`eme
avec des r\'etroactions positives peut tr\`es bien \^etre stable.

On se limitera ici au cas de matrices $A$ réelles et
diagonalisables (sur $\C$). Soient $\lambda_1,...,\lambda_n$ les
valeurs propres de $A$ (réelles ou par paires de complexes conjugués)
et $P$ la matrice de passage vers une base propre correspondante
\[ P^{-1} A P=D = \mbox{diag}(\lambda_1,...,\lambda_n)\]
On pose $X=PY$ ($Y=P^{-1}X$), on a alors
\[ \frac{dY}{dt}=P^{-1}\frac{dX}{dt} = P^{-1}(AX+C)  = P^{-1}APY + P^{-1}C
= D Y + P^{-1}C \]
Soit $B=P^{-1}C$, 
on est ainsi ramené à $n$ équations différentielles linéaires d'ordre 1
à coefficients constants~:
\[ \frac{dy_k}{dt} = \lambda_k y_k + b_k \]
Lorsque $b_k=0$, l'équation admet pour solution générale
\[ y_k(t)= c e^{\lambda_kt}\]
En faisant varier la constante $c=c(t)$, et en remplaçant dans l'équation
vérifiée par $y_k$, on obtient~:
\[ c' e^{\lambda_kt} = b_k \ \Rightarrow \ c' = e^{-\lambda_kt}  b_k \]
donc $c=\int e^{-\lambda_kt}  b_k$, et si $b_k$ est constant,
\[ c(t)=-\frac{b_k}{\lambda_k} e^{-\lambda_kt} + K \]
où la valeur de $K$ se calcule en posant $t=0$,
finalement la solution est donnée par~:
\[ y_k(t)=-b_k/\lambda_k + (y_k(0)+\frac{b_k}{\lambda_k}) e^{-\lambda_kt} \]
c'est la somme d'un terme constant ($k$-ième coordonnée de la position
d'équilibre dans la base propre) et d'un multiple d'une exponentielle
dont le comportement lorsque $t$ tend vers $+\infty$ permet de
déterminer si le système est stable.

Si $\Re(\lambda_k)>0$ pour au moins une des valeurs de $k$, alors 
la composante correspondante $y_k(t)$ tend vers l'infini, et le
système est instable. Si $\Re(\lambda_k)=0$, la composante reste bornée,
il n'y a pas convergence vers l'équilibre mais le système fait des tours
autour de la position d'équilibre. Si $\Re(\lambda_k)<0$ pour toutes
les valeurs de $k$, il y a convergence (exponentiellement rapide)
vers l'état d'équilibre, le système est stable,
la vitesse de convergence est caractérisée
par l'inverse de la valeur propre dont la partie réelle est la plus
proche de 0.

On prend souvent comme exemple de cette théorie le modèle proie-prédateur,
on peut ici dans le contexte énergétique proposer un modèle
baril (prix)-demande-production (sans contrainte géologique). 
Plus le prix du baril est élevé, plus la production augmente et
la demande baisse. 
Plus la production est élevée, plus le baril baisse.
Plus la demande est élevée, plus le baril monte. 
On est donc en dimension 3.
\[ \frac{d}{dt} \left(\begin{array}{c}
B \\
P \\
D
\end{array}\right)  =
 \left(\begin{array}{ccc}
0  & -c & d\\
a  & 0  & 0\\
-b & 0 & 0
\end{array}\right)
\left(\begin{array}{c}
B \\
P \\
D
\end{array}\right) + C
\]
On remarque que 0 est valeur propre, la trace de la matrice est nulle,
il y a 2 valeurs propres imaginaires pures $\pm i \sqrt{bd+ac}$
conjuguées et le système va effectuer des cycles (dont la période
est $2\pi/\sqrt{bd+ac}$).

\subsubsection{Les modèles non linéaires}
Un système différentiel général d'ordre 1 peut s'écrire
\[ \frac{dX}{dt} = F(X(t),t) \]
Si $F$ ne dépend pas explicitement du temps (systèréme autonome), 
on peut commencer par
chercher les positions d'équilibre possibles, i.e. les solutions
de $F(X)=0$, puis si $X_e$ est une solution,
linéariser, c'est-à-dire résoudre
\[ \frac{dX}{dt} = F'(X_e) (X-X_e) \]
pour voir si en partant de conditions initiales proches de l'\'equilibre,
le syst\`eme s'en approche (stable) ou s'en \'eloigne (instable).

Un autre cas est intéressant, celui d'une perturbation dépendant du temps
\[ F(X(t),t)=G(X)+p(t) \]
par exemple pour étudier la perturbation introduite par l'ajout de CO2 
à un système linéaire modélisant le système climatique non perturbé.
Dans le cas du CO2, on peut par exemple supposer un taux d'\'emission
constant ($p(t)=p$) et voir comment l'\'equilibre du syst\`eme
climatique est affect\'e. Plus g\'en\'erallement, pour trouver
l'\'equilibre on souhaitera trouver $X_e(t)$ tel que
\[ G(X_e(t))+p(t)=0\] 
c'est le théorème des fonctions implicites qui permet
de déterminer $X_e(t)$ en fonction de $t$, sous réserve que $\partial_X G$
soit inversible.
Si le système différentiel
linéarisé correspondant à $G(X_e(t))$ est stable et que la
perturbation $p(t)$ varie lentement par rapport aux constante
de temps caractérisant la stabilité ou est d'amplitude faible, 
alors la position d'équilibre va suivre.
Pour le système climatique soumis aux variations externes des param\`etres
orbitaux, on peut consid\'erer que c'est le cas. Par contre 
ce n'est pas vrai pour la perturbation anthropique du CO2, 
on peut toutefois \'etudier le syst\`eme
lin\'earis\'e, en supposant que $G'(X_e)$ reste \`a peu pr\`es constant.

\subsubsection{Exemples de mod\`eles de ce type}
On consid\`ere 3 param\`etres pour mod\'eliser le climat, 
la temp\'erature moyenne $T$, la masse de glace 
$G$ (qu'on suppose reli\'e \`a l'alb\'edo par la surface de glace)
et le taux de CO2 dans l'air $C$. 
La temp\'erature moyenne influe sur la puissance
rayonn\'ee (proportionnelle \`a $T^4$), sur le CO2 dissous
dans l'oc\'ean (si $T$ augmente, $C$ augmente) et sur la glace (si
$T$ augmente, $G$ diminue). Le CO2 influe sur la puissance rayonn\'ee
en fonction de $-\ln(C)$ et sur le taux de CO2 dissous dans l'eau de mer. 
La glace influe sur la puissance rayonn\'ee,
en fonction de la surface qu'on peut mod\'eliser par $G^(2/3)$. Enfin
la puissance rayonn\'ee est reli\'ee \`a la temp\'erature par
la capacit\'e calorifique de l'oc\'ean.

Il y a \'evidemment encore beaucoup de simplifications dans un tel
mod\`ele \`a 3 param\`etres, la temp\'erature de l'oc\'ean
n'est pas homog\`ene et la temp\'erature des eaux profondes ne varie
que tr\`es lentement, il y a d'autres gaz \`a effet de serre, ...

On a alors 
\[ \frac{d}{dt} \left(\begin{array}{c}
T \\
G \\
C
\end{array}\right) =
F = \left(\begin{array}{c}
c_c (\sigma(T_0^4-T^4)+\alpha ln(C) - \beta G^{2/3} ) \\
f(T) \\
g(C,T)
\end{array}\right)
\]
o\`u $c_c$ est proportionnel \a la capacit\'e calorifique de l'oc\'ean,
$f$ est une fonction d\'ecroissante, $\partial_T g$ est positif.

On peut aussi utiliser ce type de mod\'ele pour obtenir le pourcentage
de CO2 \'emis par l'homme restant dans l'atmosph\`ere au bout de $t$ ann\'ees
\[ CO_2(t)=18+26 e^{-t/3.4} + 24 e^{-t/21} + 18 e^{-t/70} + 14 e^{-t/420}\]

\subsection{Les modèles de type EDP}
Les fluides de l'atmosphère obéissent aux équations physique
de la dynamique des fluides. Contrairement au problème à 2 corps,
on ne peut pas résoudre ces équations de manière exacte. On utilise
alors des méthodes numériques, en discrétisant l'état de l'atmosphère.
On obtient alors une solution approchée des équations, c'est le
principe des prévisions numériques météorologiques. Ce principe
s'applique aussi au climat. On peut s'étonner que le même principe
qui ne permet pas des prévision fiables à plus de quelques jours
puisse permettre de faire des prévisions à l'échéance du siècle. La
contradiction n'est qu'apparente, parce que les modèles ne prévoient
pas le temps précis un jour donné, mais le temps moyen, en fonction
des grandes contraintes énergétiques (ensoleillement, effet de serre)
de même qu'on peut prévoir qu'il fera en moyenne plus chaud en été
qu'en hiver parce que l'ensoleillement est plus important en été
qu'en hiver. L'intérêt des modèles est de pouvoir quantifier
le ``il fera plus chaud si l'effet de serre est renforcé''.

Les problèmes mathématiques liés à la modélisation du climat
sont liés à la stabilité des méthodes de résolution numérique.
On peut en avoir une petite idée en étudiant les méthodes
numériques d'intégration ou la méthode
d'Euler pour résoudre des équations différentielles de type
$y'=f(x,y)$ (mais ici on a 3 dimensions d'espace et 1 de temps
et beaucoup plus qu'une variable, il y a la pression, la température,
l'humidité, ...). On peut aussi s'inspirer d'exemple de suites itératives
pour essayer de comprendre certains phénomènes.

\section{Energies renouvelables.}

\subsection{Le solaire thermique}
Le calcul du rendement d'un panneau solaire thermique découle
du calcul de l'ensolleillement fait au théorème \ref{thm:soleil}, mais
il faut tenir compte du fait que la normale $n$ au panneau solaire
n'est plus la verticale locale. On exprime les coordonnées de $n$
dans le repère local orthonormé
formé de la verticale locale, du Sud local, et de la direction
obtenue par produit vectoriel des deux vecteurs précédents.
La verticale locale $v$ est connue (cf. section \ref{sec:soleil}),
le vecteur normé porté par la direction du Sud local est égal
à l'opposé de la dérivée de $v$ par rapport à la latitude, après 
normalisation. 

Lorsque l'angle avec la verticale est fixé, on a bien sur
intérêt à aligner le panneau avec le Sud local.
On peut calculer l'angle optimal pour recevoir le maximum
d'énergie solaire tout au long de l'année, mais il est plus
intéressant de calculer l'angle optimal pour maximiser le
rayonnement reçu pendant la période de chauffe pour minimiser
l'utilisation d'un appoint (par exemple
sur la période automne-hiver, on trouve un angle voisin de 60 degrés
à nos latitudes sans tenir compte de la météo). On peut aussi
chercher à optimiser l'angle au solstice d'hiver si on souhaite
éviter l'utilisation d'un appoint (on utilisera par contre plus
de surface de panneaux, ceux-ci devant alors être toujours 
excédentaires).

\subsection{Le solaire à concentration.}
Il s'agit de générer de la chaleur (four solaire) ou 
du courant électrique comme avec une centrale
thermique classique, la chaleur étant fournie par les rayons solaires
au lieu de la combustion de charbon, gaz ou fuel.
Les rayons solaires ne sont bien sur pas suffisamment chauds pour
obtenir une chaleur ou un rendement suffisant sans concentration. 
La concentration
se fait à partir d'un ou de miroirs qui réfléchissent les rayons en une
zone la plus petite possible. 

Etude idéalisée~: on suppose que le
miroir peut prendre la forme d'une courbe r\'eguli\`ere
(ou que les miroirs sont suffisament
petits pour être assimilable à des tangentes à une courbe).
Quelle forme doit prendre la courbe pour focaliser les rayons
solaires incidents en un seul point~?
Une réponse est connue depuis longtemps, il s'agit d'une parabole (ou
plus précisément d'un paraboloide de révolution) qui focalise
les rayons perpendiculaires à la directrice de la parabole
vers le foyer de la parabole. Montrons 
qu'effectivement cette contrainte impose une parabole, en effectuant
une coupe par un plan passant par le foyer.

Rappel~: une parabole est d\'efinie g\'eom\'etriquement par une
droite $D$ appel\'ee directrice et un point $O$ appel\'e foyer,
c'est le lieu des points \'equidistants \`a $D$ et $O$.
\begin{figure} \label{fig:parabole}
\includegraphics[width=\textwidth]{parabole}
\caption{Parabole}
\end{figure}
Si on calcule l'\'equation de ce lieu, en prenant $O$ comme origine,
$M(x,y)$ et l'\'equation de $D$ sous la forme $y=-m$, on a
\[ x^2+y^2=OM^2=d(M,D)^2=(y+m)^2=y^2+2ym+m^2\]
donc
\[ y=\frac{x^2-m^2}{2m}\]
c'est bien l'\'equation usuelle d'une parabole.
La normale \`a la parabole est obtenue en prenant le gradient de l'\'equation
qui la d\'efinit
\[ OM-d(D,M)=0\]
c'est donc
\[ \frac{\overrightarrow{OM}}{OM} - \overrightarrow{j} \]
qui est le vecteur directeur de la bissectrice de la verticale et de $OM$
(plus pr\'ecis\'ement la bissectrice int\'erieure du vecteur 
$\overrightarrow{OM}$ et de $-\overrightarrow{j}$
ou de $\overrightarrow{MO}$ et de $\overrightarrow{j}$). 
La tangente est donc la bissectrice ext\'erieure
de cet angle, ou encore la bissectrice int\'erieure de $(MO,MK)$. On en
d\'eduit qu'un rayon vertical venant du haut se r\'efl\'echit sur
la parabole sym\'etriquement \`a la normale donc en passant par le foyer 
$O$. 

Consid\'erons maintenant une courbe ayant cette propri\'et\'e: les rayons
provenant du haut verticalement se r\'efl\'echissent dessus
en passant par $O$. Soit $M(x,y)$ un point de cette courbe. Prolongeons
le rayon vertical passant par $M$ vers le bas de la longueur $OM$,
on obtient un point $K$ d'ordonn\'ee 
\[ y_k=  y - \sqrt{x^2+y^2}\]
Montrons que lorsque $x$ varie, $y_k$ a une d\'eriv\'ee nulle. On a~:
\[ y_k'= y' - \frac{x+yy'}{\sqrt{x^2+y^2}} \]
D'autre part, la tangente en $M$ \`a la courbe est orthogonale \`a
$OK$ puisque c'est la bissectrice de l'angle $(MO,MK)$ (hypoth\`ese
que le rayon r\'eflechi passe par $O$) et que $MO=MK$ (par construction
de $K$), et le vecteur directeur de la tangente est $(1,y')$, on a donc
\[ (1,y').(x,y_k)=0 \Rightarrow x+y'y-y'\sqrt{x^2+y^2}=0\]
donc $y_k'$ est bien nul et $K$ est donc sur une droite horizontale fixe $D$,
d'o\`u l'on d\'eduit que $M$ est sur la parabole de foyer $O$ et
directrice $D$, la courbe est donc forc\'ement une parabole.

Rendement d'un tel système~: on peut le définir comme le rapport
de l'énergie focalisée sur la surface de miroir. 
Pour un miroir cylindrique avec une coupe parabolique, la surface
est proportionnelle \`a la longueur de la courbe en coupe. On peut
faire le calcul sur la parabole $y=x^2/2$ pour d\'eterminer par exemple
quel portion de parabole on peut raisonnablement utiliser. Le cout
de l'\'energie est proportionnel \`a
\[ \frac{1}{l} \int_0^l \sqrt{1+x^2} dx 
= \frac{1}{2} \left(
\sqrt{l^{2}+1}-\frac{1}{l}\ln(\sqrt{l^{2}+1}-l)
\right)
\]
Il ne faut pas prendre $l$ trop grand (sinon le cout augmente trop vite),
mais pas trop petit non plus (sinon on ne concentrera pas assez d'\'energie
pour avoir par exemple un bon rendement de Carnot).
On peut aussi jouer sur le param\`etre de la parabole (pour l'aplatir pr\`es de
son point le plus bas, ce qui permet d'augmenter $l$ \`a rendement
constant), mais il faut alors augmenter l'altitude du foyer.

Ce système est bien adapté s'il reste petit, par exemple
pour un four solaire pour faire la cuisson
d'un repas, mais n'est pas adapté à de grandes installations, la
taille du miroir parabolique à réaliser étant trop grande et trop
difficile à manoeuvrer (toute l'installation doit bouger pour que
la directrice de la parabole reste perpendiculaire aux rayons solaires).
Pour une centrale solaire à concentration, 
le foyer doit être à une position fixe, et 
les miroirs doivent pouvoir être
inclinés pour faire suivre les rayons du Soleil vers le foyer
au cours de la journée. La direction miroir-foyer est donc fixe, et la 
direction miroir-Soleil imposée par la position du Soleil,
le plan du miroir est déterminé par le point du miroir, la bissectrice
extérieure aux 2 directions (dans le plan les contenant) 
et la perpendiculaire à ce plan. Ceci impose bien d'avoir des miroirs
amovibles. L'énergie captée par ce système est
proportionnel au cosinus de l'angle entre la normale au miroir 
et la direction miroir-Soleil. Comme la normale au miroir est la
bissectrice intérieure de l'angle miroir-Soleil, miroir-foyer,
l'énergie captée (par miroir) est proportionnelle au cosinus de la
moitié de l'angle miroir-Soleil, miroir-foyer. On est donc ramené à un
calcul analogue à celui du rendement de panneaux solaires thermiques,
la normale au panneau étant remplacée par la direction miroir-foyer,
mais on divise ensuite l'angle par 2.

\subsection{Les aspects économiques.}
Il y a actuellement deux grands freins 
au d\'eveloppement des \'energies renouvelables~: l'intermittence
(conjointement avec la possibilit\'e de stocker l'\'energie)
et la facon dont le co\^ut est calcul\'e. Le premier frein
n\'ecessite des avanc\'ees technologiques pour pouvoir monter
en puissance, le deuxi\`eme ne me semble pas intrins\`eque mais
uniquement fond\'e sur la facon de financer un projet \'energie
renouvelables. En effet, un projet d'installation de panneaux solaires
ou d'\'eoliennes n\'ecessite un investissement important au d\'epart
(qui peut repr\'esenter jusqu'\`a plus de 90\% du co\^ut total
par exemple pour des panneaux solaires thermiques) contrairement
\`a un projet de centrale thermique au charbon ou au gaz.
Comment doit-on r\'epartir ce cout de mani\`ere juste sur la
dur\'ee de fonctionnement de l'installation~? Les \'economistes
utilisent un taux d'actualisation, ils comparent avec ce qu'une somme
d'argent \'equivalente aurait pu rapporter en euros constants. 
Plus pr\'ecis\'ement le cout annuel estim\'e correspond \`a
la moyenne des couts de fonctionnement plus le cout annuel qu'il faudrait
payer pour rembourser en euros constants l'investissement initial.
Le choix du taux d'actualisation est alors un param\`etre tr\`es important
qui peut faire varier du simple au double (voire plus) le co\^ut estim\'e
des \'energies renouvelables. Les taux actuels de 5\% \`a 10\%, 
calcul\'es sur la base d'une croissance soutenue
de l'\'economie mondiale rendue possible par la hausse de la 
consommation mondiale d'\'energie fossiles et l'am\'elioration de 
l'efficacit\'e \'energ\'etique, paraissent peu compatibles avec une 
p\'eriode de ``plateau fossiles'' voire de d\'epl\'etion.

%\subsection{L'éolien.}



\section{Conclusions}
Les variations de concentration de gaz à effet de serre au cours du passé
agissent sur le bilan énergétique de la Terre de manière logarithmique
de l'ordre de quelques
$W/m^2$ (Watt par mètre carré), une valeur qui est
du même ordre de grandeur que celles dues aux variations orbitales 
de la Terre ou que celles dues aux variations des volumes de glace
(modification de l'albédo terrestre). 
La corrélation entre gaz à effet de serre
et température calculée sur les 400 000 dernières années en Antarctique
donne une sensibilité de 9 degrés pour 100
ppm de CO2 au niveau de l'Antarctique. En divisant ce chiffre par 2 pour
tenir compte de la plus grande variabilité de la température en
Antarctique puis par 3 pour tenir compte des 3 causes principales
de variation de même importance (gaz à effet de serre, variations
orbitales, volume des glaces), on obtient une sensibilité empirique
de 1.5 degré pour 100ppm de CO2. Ce qui veut dire que si nous conservions
la concentration actuelle en CO2 (100 ppm de plus qu'avant le
début de l'ère industrielle en 2005), ce modèle prévoit une hausse de 1.5
degrés de la température, ce qui est parfaitement plausible avec l'élévation
de 0.6 degrés du 20ème siècle (l'océan jouant le role de régulateur
thermique de par son inertie
sur des échelles de temps entre le siècle et le millénaire, par
exemple le forçage radiatif actuel est estimé à plus de 2$W/m^2$,
or pour réchauffer l'océan de 1K avec une telle puissance il faut
attendre plus de 200 ans, pour faire ce calcul, on considère
qu'il y a une épaisseur moyenne de 3000 mètres d'eau, donc chaque
heure, 2$Wh$ servent à réchauffer 3000 mètre cubes d'eau, or il faut
un peu plus d'un $kWh$ pour réchauffer un mètre cube d'eau d'un degré,
il faut donc environ 1,5 millions d'heures pour réchauffer l'océan 
à cette puissance).

Bien évidemment, une corrélation seule ne suffit pas à justifier
l'existence d'une relation de cause à effet entre gaz à effet
de serre et température, ce sont les lois de la physique  qui
justifient cette relation, la corrélation nous donne toutefois
des informations quantitatives sur l'effet que l'on peut
rapprocher des calculs issus de la physique. Ceux-ci donne
un forçage radiatif de 4W par doublement du CO2 soit 0.3\% de
l'énergie reçue au sommet de l'atmosphère perpendiculairement
au Soleil, ou encore 1.2\% de l'énergie reçue au sommet
de l'atmosphère en moyenne (pour tenir compte de la surface
de la Terre $4\pi R^2$ qui est 4 fois plus grande que sa section 
perpendiculaire au Soleil $\pi R^2$). En tenant compte de la partie
du rayonnement arrivant effectivement dans les basses couches de
l'atmosphère, c'est un forçage de 2 à 3\%. Si la Terre
était un corps noir à une température de 288K, ce forçage se
traduirait par une hausse de plus de 0.5\% de la température pour compenser
(car la puissance émise par un corps noir est proportionnelle à $T^4$,
il faut donc diviser par 4 le forçage), donc de 2K environ. Ce chiffre
est raisonnablement en accord avec les modèles du climat qui prévoient 
une hausse de 3 à 4K pour un doublement du CO2, en tenant compte
des diverses rétroactions (modification de l'albédo par exemple).

Si nous voulons éviter une variation trop brusque de la température
et plus encore des précipitations, il faut donc agir dès maintenant,
ce qui ne sera pas facile car 1ppm de CO2 correspond à 2Gt de C
soit à peine plus de 300kg de C par Terrien (ou un peu plus
de 400 litres de pétrole)~: rappelons qu'un Français consomme 
(en gros pour moitié
directement par les transports et le chauffage, et l'autre moitié
indirectement par la consommation de biens) en moyenne 1.5 tep
de pétrole, 0.67 tep de gaz et 0.2 tep de charbon par an, il émet donc 2 tonnes
de C par an (ce qui généralisé à l'échelle de la planète correspondrait
à une émission de 6 ppmv de CO2 par an).

\section{R\'ef\'erences sur le web}
\begin{itemize}
\item Astronomie
\begin{itemize}
\item calcul de l'insolation en fonction de la latitude
au cours des millénaires écoulés (Laskar et al)\\
\verb|http://www.imcce.fr/fr/presentation/equipes/ASD/insola/earth/earth.html|
\item La théorie de Milankovitch: rechercher Milankovitch sur google
par exemple sur le site de l'ENS Lyon\\
\verb|http://www.ens-lyon.fr/Planet-Terre/Infosciences/Histoire/Paleoclimats/|
\end{itemize}

\item Climat
\begin{itemize}
\item IPCC \verb|http://www.ipcc.ch/|
\item Realclimate \verb|www.realclimate.org|
\item Reconstructions climatiques (nombreux liens)\\ 
\verb|http://www.ncdc.noaa.gov/paleo/recons.html|,
\item Données numériques de Vostok \\
\verb|ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/vostok/|, 
le fichier ch4nat.txt contient les données pour le méthane, co2nat.txt
pour le CO2, deutnat.txt pour la température. 
(remonter d'un cran dans l'arborescence du site pour d'autres lieux).
\item données du projet EPICA les plus récentes\\
\verb|ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/epica_domec/|
\item Données climatiques au 20 ième siècle \\
\verb|http://data.giss.nasa.gov/|\\
\verb|http://hydrolab.arsusda.gov/nicks/nicks.htm|
\end{itemize}

\item \'Energie, émission de gaz à effet de serre: 
\begin{itemize}
\item Concentration de CO2 à Mauna Loa (Hawai)\\
\verb|http://cdiac.esd.ornl.gov/ftp/trends/co2/maunaloa.co2|
\item Statistiques du gouvernement français\\
\verb|http://www.industrie.gouv.fr/energie/sommaire.htm|
\item The Oil Drum \verb|http://www.theoildrum.com/|
\item ASPO-France \verb|http://aspofrance.org/|, 
\item Institut Francais du pétrole (IFP) \verb|http://www.ifp.fr|
\item site de J.M. Jancovici (énergie et effet de serre),
 \verb|http://www.manicore.com/|
\end{itemize}

\item
En complément, vous pouvez aussi faire des recherches sur wikipedia, google
et d'autres moteurs de recherche.
\end{itemize}

\appendix

\section{Pr\'ecession.}
(D'apr\`es Fr\'ed\'eric Faure).
La rotation de la Terre a pour cons\'equence qu'elle n'est pas sph\'erique.
En effet, la force centrifuge vaut $m\omega^2 r$, o\`u
$r$ est le vecteur issu de la projection sur l'axe de rotation 
et $\omega=2\pi/1$ jour est la vitesse angulaire. 
Cette force d\'erive du potentiel $-m \omega^2 r^2/2$ qu'il faut
ajouter au potentiel d\^u \`a la gravit\'e, en premi\`ere
approximation $mgh$ o\`u $h+R$ est la distance au centre de la Terre,
($R$ \'etant la distance moyenne de la surface de la Terre \`a son centre).
Si l'on ne tient pas compte de la force centrifuge, $h=0$ pour avoir
un potentiel identiquement nul sur toute la surface, on a une sph\`ere.
Si l'on en tient compte, ce ne sera pas le cas sauf aux poles (o\`u $r=0$).
Ainsi \`a l'\'equateur, pour avoir un potentiel nul, il faut prendre
$h= \omega^2 R^2/(2g)$, soit $h=11$km environ, il y a donc un
bourrelet un peu plus \'epais au niveau de l'\'equateur.
En r\'ealit\'e, on observe un aplatissement de l'ordre du double,
l'approximation du potentiel de gravitation n'est en effet pas
tr\`es bonne, parce qu'elle suppose une Terre sph\'erique, ce
qui n'est pas le cas comme on vient de le voir.

L'existence d'un bourrelet \'equatorial et l'inclinaison de l'axe
de rotation de la Terre par rapport au plan de l'orbite terrestre entraine
l'existence d'un couple exerc\'e par le Soleil et la Lune sur la Terre.
Ce couple s'exerce dans des plans contenant l'axe de rotation de la Terre
(qui est axe de sym\'etrie)
et a donc un moment perpendiculaire \`a l'axe de rotation. Les
\'equations du mouvement (d\'eriv\'ee du moment cin\'etique
donc de la vitesse angulaire = moment du couple) entrainent alors que la
d\'eriv\'ee du vecteur rotation est perpendiculaire
\`a l'axe de rotation
(et cette d\'eriv\'ee est petite), le vecteur rotation conserve
donc sa norme au cours du temps mais sa direction bouge lentement
perpendiculairement \`a lui-m\^eme. Cela entraine en un an une petite
rotation du vecteur rotation de la Terre sur un cone autour
de la perpendiculaire au plan de l'orbite. On peut faire
le calcul en mod\'elisant la Terre par un ellipsoide de r\'evolution
(cf. par exemple
\verb|www-fourier.ujf-grenoble.fr/~faure/enseignement/|
puis \verb|meca_analytique_L3/index.html#Travaux_dirig%E9s|),
et on montre que l'action du Soleil et de la Lune s'additionne
pour faire faire un tour complet \`a l'axe de rotation de la Terre
autour de la perpendiculaire au plan de l'orbite en un peu
moins de 26 000 ans. C'est la pr\'ecession astronomique.

Le probl\`eme se complique si on consid\'ere les autres plan\`etes
du syst\`eme solaire, en particulier Jupiter. En effet l'orbite
elliptique de la Terre se d\'eforme au cours des mill\'enaires, 
en particulier la position de l'axe
des foyers, donc la direction de l'axe aph\'elie-p\'erih\'elie,
or c'est l'angle entre cet axe et la projection de l'axe de rotation
de la Terre qui importe pour les variations des saisons. C'est
pour cette raison que la p\'eriodicite de la pr\'ecession
``climatique'' est un peu plus courte que la pr\'ecession
astronomique.

\end{document}
