\NeedsTeXFormat{LaTeX2e}
\documentclass[10pt,a4paper]{amsart}
\usepackage{amsmath,amsfonts,amssymb,latexsym}
\usepackage{a4,mathptmx,helvet,courier}
\usepackage[francais]{babel}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}

\newcounter{exonum}[section]
\newenvironment{exo}{\begin{enumerate}\stepcounter{exonum}%
\renewcommand{\theenumi}{\thesection.\arabic{exonum}}%
\renewcommand{\labelenumi}{\bf\theenumi.}\item}{\end{enumerate}}

\newcommand{\N}{\mathbb{N}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\usepackage{graphicx}

\begin{document} 

\pagestyle{empty}
\thispagestyle{empty}

{\footnotesize\noindent
Universit\'e Joseph Fourier, Grenoble I   \hfill $\bullet$ \hfill
L3 Intro. mod\'elisation num\'erique                       \hfill $\bullet$ \hfill
Ann\'ee 2015/2016           \\[-2mm]\hrule}


\begin{center}
%  {\Large\bf Math\'ematique assist\'ees par ordinateur} \\
  Examen du 18 mai 2016, de 9h \`a 12h. \\
  \it Documents, calculatrices et ordinateurs ultraportables
  (netbooks) d\'econnect\'es du r\'eseau autoris\'es.\\
  Ce sujet comporte 2 pages. 
Bar\^eme donn\'e \`a titre indicatif et non contractuel.
\end{center}

\section{Valeurs singuli\`eres d'une matrice (11 points)}
Soit $A$ une matrice non nulle carr\'ee d'ordre $n$ et $A^*$ sa transpos\'ee
(ou sa transconjugu\'ee si $A$ est \`a coefficients complexes), on note $B=A^* A$
\begin{enumerate}
\item Montrer que $B$ est diagonalisable de valeurs propres r\'eelles
  positives.
\item Soit $v$ un vecteur non nul tel que $\| (B-\mu) v\| \leq \varepsilon \| v
  \|$, montrer que $B$ admet une valeur propre dans l'intervalle
$[\mu-\epsilon,\mu+\epsilon]$ (on pourra \'ecrire $v$
dans une base orthonormale bien choisie relativement \`a $B$).
\item Soit $\lambda>0$ la plus grande valeur propre de $B$, 
donner un algorithme permettant de trouver une valeur approch\'ee
$\tilde{\lambda}$ de $\lambda$ avec une pr\'ecision $\varepsilon$
en utilisant les suites $v_{k}=w_{k}/\|w_k\|, w_{k+1}=Bv_k $ pour un vecteur
$w_0$ \`a coefficients al\'eatoires. 
\item Programmer l'algorithme et 
l'appliquer \`a la matrice de Hilbert de taille 5 et 10. On
  pourra utiliser la commande \verb|hilbert(5), hilbert(10)|.
\item Quel est le cout d'une it\'eration de l'algorithme en fonction
de la taille $n$ de la matrice~? Comparer
avec le cout du calcul de $B$. Peut-on re\'ecrire l'algorithme pour
ne pas avoir \`a calculer $B$~?
\item Proposer une algorithme de calcul de la valeur
  approch\'ee de la plus petite valeur propre de $B$ en utilisant
  $B^{-1}$. Tester l'algorithme avec la matrice de Hilbert de taille
5. Pour une matrice de Hilbert de taille 10, il faut choisir une
tol\'erance relative, donc un test d'arr\^et du type $\|
w_{k+1}-\lambda v_k\| \leq \varepsilon |\lambda| \|v_k\|$,
expliquer pourquoi.
\item Quel est le cout de cet algorithme en fonction de $n$~? 
Est-il possible d'optimiser en utilisant la d\'ecomposition $LU$ de $B$
ou de $A$~?
\item On se propose d'appliquer l'algorithme de la question (3)
\`a $\tilde{\lambda}I_n-B$ o\`u $\tilde{\lambda}$ est une valeur approch\'ee
de la plus grande valeur propre de $B$. Montrer que la valeur propre
approch\'ee obtenue donne une estimation de la plus petite valeur
propre de $B$. Comparer le cout de l'algorithme avec celui de la
question pr\'ec\'edente et comparer la qualit\'e des approximations 
obtenues pour la plus petite valeur propre de $B$.
\item Rappeler l'expression du conditionnement en norme euclidienne
d'une matrice $A$ en fonction de ses valeurs singuli\`eres et
l'expression des valeurs singuli\`eres de $A$ en fonction des valeurs
propres de $B$.\\
Donner un algorithme de calcul du conditionnement
de $A$ relativement \`a la norme euclidienne.
\item Tester
avec une matrice de Vandemonde (\verb|vandermonde|) 
prise en 11 points \'equidistants
de $[-1,1]$ puis en 11 points de Tchebyshev de $[-1,1]$ (i.e. les 11
racines du 11-i\`eme polyn\^ome de Tchebyshev). 
\item  On voudrait calculer la deuxi\`eme plus grande valeur propre de
  $B$, pour cela on consid\`ere la matrice $B'=B-\tilde{\lambda} w
  w^*$ o\`u $w$ est un vecteur propre
colonne norm\'e approch\'e de $B$ correspondant \`a $\tilde{\lambda}$.
D\'eterminer les valeurs propres de $B'$ en faisant l'approximation que $w$ 
est un vecteur propre de $B$. En d\'eduire un algorithme permettant de
d\'eterminer la deuxi\`eme plus grande valeur propre de $B$.
Peut-on \'ecrire cet algorithme sans calculer explicitement $B$~?
\end{enumerate}


\section{\'Equation diff\'erentielle (9 points)}
%
Pour r\'esoudre de mani\`ere approch\'ee
une \'equation (ou un syst\`eme diff\'erentiel) $y'=f(t,y)$ avec $f$ de
classe $C^2$,
on va appliquer la m\'ethode explicit\'ee ci-dessous, appel\'ee
m\'ethode de Heun de pas constant $h$. 

Soit $y_0=y(t_0)$ la condition initiale, $t_n=t_0+nh$, on approche
$y(t_n)$ par la suite r\'ecurrente $y_n$ d\'efinie par~:
$$ \tilde{y}_{n+1}=y_n+hf(t_n,y_n), \quad
y_{n+1}=y_n+\frac{h}{2}(f(t_n,y_n)+f(t_n+h,\tilde{y}_{n+1}))$$
\begin{enumerate}
\item On suppose que $f(t,y)$ ne d\'epend pas de $y$,
quelle est la m\'ethode num\'erique d'int\'egration correspondante~?
En d\'eduire une majoration de l'erreur locale (ou erreur de consistance)
$|y_1-y(t_0+h)|$ en fonction de $h$ dans ce cas.
\item % Le montrer?
Dans toute la suite $f$ peut d\'ependre de $y$.
On pourra utiliser sans le justifier que $F: \ t \rightarrow f(t,y(t))$
est $C^2$ sur l'intervalle de d\'efinition maximal de $y$.\\
Donner une valeur approch\'ee de l'int\'egrale
$I=\int_{t_0}^{t_0+h} f(t,y(t))$ par la m\'ethode des rectangles
et une majoration de l'erreur $ \tilde{y_1}-y(t_0+h)$.\\
Donner une valeur approch\'ee de l'int\'egrale $I$ par la m\'ethode
de la question (1), 
puis une majoration de l'erreur locale $|y_1-y(t_0+h)|$.
\item On pose $y_{n+1}=y_n+h\Phi(t_n,y_n,h))$,
donner l'expression de $\Phi$ en fonction de $f$, justifier que $\Phi$
est lipschitzienne par rapport \`a $y$, et en d\'eduire une majoration
de l'erreur globale $|y_n-y(t_n)|$.
\item On suppose dans la suite que $f(t,y)=ty$ et que $h\leq 0.1, t
  \in [0,T]$.
Calculer une constante de Lispchitz pour $\Phi$.
\item Programmer la m\'ethode. D\'eterminer une valeur approch\'ee
de $y(1)$ obtenue par cette m\'ethode
pour l'\'equation $y'=ty, y(0)=1$ avec $h=0.1, 0.01, 0.001$.
Comparer avec la valeur exacte de $y(1)$, observe-t-on le comportement attendu
en puissance de $h$~?
\item D\'eterminer de m\^eme une valeur approch\'ee $y_N$ 
de $y(t_N)$ pour $t_N=4$
et $h=0.1,0.01, 0.001$ (donc $N=40, 400, 4000$), calculer l'erreur 
en utilisant la valeur exacte de $y(t_N)$, observe-t-on le
comportement attendu en fonction de $t_N$~?
\item Soit l'\'equation d'ordre 2 $y''+ty'+y=0$. L'\'ecrire comme
un syst\`eme diff\'erentiel d'ordre 1 en $Y=(y,y')$ et appliquer la m\'ethode
pr\'ec\'edente pour d\'eterminer une valeur approch\'ee
de $y(1)$ pour la solution de condition initiale $y(0)=1, y'(0)=0$
($Y(0)=(1,0)$).
\end{enumerate}

\vfill

\pagebreak

Notes de correction 1 (13.5 points=.5+1.5+1.5+.5+1.5+1.5+1+2+1+1+1.5)
\begin{enumerate}
\item car B est sym\'etrique positive
\item par l'absurde, on minore la norme euclidienne au carr\'ee de
  $(B-\mu)v$ par $\epsilon^2$ fois la norme au carr\'e de $v$
\item M\'ethode de la puissance (cf. TD/cours)
\begin{verbatim}
svd1(A,eps,N):={
  local w,v,j,At,l;
  At:=trn(A);
  w:=ranv(size(A),0..1);
  for j from 1 to N do
    v:=normalize(w);
    w:=A*v; w:=At*w;
    l:=dot(v,w);
    if (l2norm(w-l*v)<eps*abs(l)){ return v,l; }
  od;
  retourne "maxiter";
}:;
\end{verbatim}
\item En principe c'\'etait $A$ la matrice de Hilbert, et donc $B=A^*A$.
\item Le cout de l'it\'eration est domin\'ee par la multiplication matrice-vecteur en
  $O(n^2)$ (le reste en $O(n)$). N\'egligeable devant le cout
  multiplication matrice-matrice pour calculer $B$ en $O(n^3)$.
On peut \'eviter ce cout en faisant la multiplication du vecteur par
$A$ puis par $A^*$.
\item on prend l'inverse de la plus grande valeur propre de $B^{-1}$.
La plus grande valeur propre de $B^{-1}$ est grande pour \verb|hilbert(5)|
et m\^eme  tr\`es grande pour \verb|hilbert(10)| (onze chiffres),
obtenir une convergence avec un test absolu est impossible si on
travaille en flottants (pour $\epsilon<1e-5$ en tout cas, et 
serait beaucoup trop lent en flottants multi-pr\'ecision).
\item le cout est le m\^eme que pour la plus grande valeur propre de
  $B$ une fois $B^{-1}$ calcul\'e ce qui n\'ecessite $O(n^3)$
  op\'erations (calcul de $B$ et de $B^{-1}$). On peut un peu
  optimiser la constante dans le $O(n^3)$ en rempla\c{c}ant
\verb|w:=A*v; w:=At*w;| par \verb|w:=linsolve(p,l,u,v)|
o\`u $p,l,u$ est la d\'ecomposition $LU$ de $B$, voire en
\'economisant la multiplication $B=A^*A$, en
d\'ecomposant la r\'esolution de $A^*Aw=v$ en
plusieurs syst\`emes triangulaires. Mais cela reste en $O(n^3)$.
\item La plus grande valeur propre de $\tilde{\lambda} I_n-B$
est $\tilde{\lambda}$  moins la plus petite valeur propre de $B$.
On peut donc appliquer la m\'ethode de la puissance directement
sur $\tilde{\lambda} I_n-A^*A$ pour avoir une estimation de la plus
petite valeur propre de $B$, ce qui coutera $O(kn^2)$ si $k$ est
le nombre d'it\'erations qui d\'epend lin\'eairement de la
pr\'ecision souhait\'ee.\\
Mais pour les matrices de Hilbert, la plus petite valeur propre
est trop proche de 0 pour obtenir une estimation par cette m\'ethode,
on va trouver une estimation proche de $\tilde{\lambda}$ et 
on va donc soustraire deux nombres tr\`es proches l'un de l'autre,
en perdant tous les chiffres significatifs (sauf si on travaille en
multi-pr\'ecision mais le nombre d'it\'erations sera tr\`es grand).
\item conditionnement $L^2$ de $A$=racine carr\'ee du
rapport de la plus grande valeur singuli\`ere sur la plus petite de $A^*A$,
donc racine carr\'ee de l'estimation de la plus grande valeur propre
de $B$ sur la plus petite.
\item On applique \`a \verb|vandermonde(seq(-1,1,0.1)|
et \verb|vandermonde(proot(tchebyshev1(11)))| (on peut
aussi \'ecrire explicitement les racines).
\item L'orthogonal de $w$ est stable par $B$ et $B'$ et les valeurs propres
sont identiques. Pour $w$ on remplace la valeur propre $\lambda$
par 0. Donc la plus grande valeur propre de $B'$ est la seconde valeur
propre de $B$. On peut \'ecrire l'algorithme de la puissance
sans calculer $B$ en multipliant par $A$ puis $A^*$.
\end{enumerate}

Notes de correction 2: (10.5=1+2+1.5+1+3+1+1)
\begin{enumerate}
\item m\'ethode des trap\`ezes, donc erreur locale sur un pas en $M_2
  h^3/12$ ($M_2=$max$|f'{'}|$)
\item les rectangles \`a gauche donnent $\tilde{y}_{1}$ donc
  $|\tilde{y}_{1}-y(t_0+h)|$
est major\'e par $M_1h^2/2$ (ici $M_1$=max$|F'|$).
On applique les trap\`ezes \`a $F$, on en d\'eduit que
$$ | y(t_0+h)-y(t_0) - \frac{h}{2}(f(t_0,y_0)+f(t_0+h,y(t_0+h))) |
\leq M_2 h^3 \frac{1}{12}$$
avec ici $M_2$=max$|F''|$. Donc
\begin{eqnarray*}
 |y_1-y(t_0+h)| 
&\leq& \frac{h}{2} |f(t_0+h,\tilde{y}_1)-f(t_0+h,y(t_0+h))
| + M_2 h^3 \frac{1}{12} \\
&\leq &\frac{h}{2} \Lambda |\tilde{y}_{1}-y(t_0+h)| + M_2h^3\frac{1}{12} \\
&\leq &\frac{h^3}{12} (3 \Lambda M_1+M_2)
\end{eqnarray*}
\item $$ \Phi(t_n,y_n,h)=f(t_n,y_n)+f(t_n+h,y_n+hf(t_n,y_n))$$
est donc lipschitzienne en $y_n$~:
\begin{eqnarray*}
 |\Phi(t_n,z,h)-\Phi(t_n,y,h)| &\leq& |f(t_n,z)-f(t_n,y)|
+ |f(t_n+h,z+hf(t_n,z))-f(t_n+h,y+hf(t_n,y))| \\
&\leq &\Lambda (|z-y|+|z+hf(t_n,z)-(y+hf(t_n,y))| \\
&\leq &\Lambda (2+h\Lambda)|z-y|
\end{eqnarray*}
D'o\`u une erreur globale en $h^2/12 (3\Lambda M_1+M_2)
\frac1C(e^{Ct_n}-1)$
avec $C$ la constante de Lipschitz de $\Phi$.
\item On peut prendre $\Lambda=T$ donc $C=T(2+.1T)$
\item 
\begin{verbatim}
heun(f,t0,t1,N,y0):={
  local h,j,t,y1;
  h:=evalf(t1-t0)/N;
  pour j de 0 jusque N-1 faire
    t:=t0+j*h;
    y1:=y0+h*f(t,y0);
    y0:=y0+h/2*(f(t,y0)+f(t+h,y1));
  fpour; 
  retourne y0;
}:;
\end{verbatim}
\verb|heun(f,0,1,10,1.0)| renvoie 1.64788134551.
\verb|heun(f,0,1,100,1.0)| renvoie 1.6487142374
\verb|heun(f,0,1,100,1.0)| renvoie 1.64872120183
La solution exacte est $y'/y=t$ donc $y=Ke^{t^2/2}$ avec $K=1$
(condition initiale 1 en $t=0$). On compare donc avec
$e^{0.5}=-1.6487212707$. On a le comportement attendu d'une erreur
en $h^2$ (on gagne 2 d\'ecimales si on divise $h$ par 10, erreur de
\verb|6e-8| pour \verb|h:=1e-3|).
\item Si $t_n=4$ au lieu de 1, on a
\verb|heun(f,0,4,40,1.0)-exp(8)| qui vaut -242,
\verb|heun(f,0,4,400,1.0)-exp(8)| qui vaut -3,
\verb|heun(f,0,4,400,1.0)-exp(8)| qui vaut -0.03,
on a toujours le bon comportement en puissances de $h$,
mais une erreur environ 1 million de fois plus grande. 
La constante de Lipschitz qui passe de 2.1 \`a 9.6, qui modifie
le terme $\frac{1}{C}(e^{Ct_n}-1)$ de 3 \`a \verb|5e15| expliquant
un rapport d'erreur tr\`es grand (on aurait d'ailleurs pu s'attendre \`a beaucoup
plus que 1 million).
\item \verb|f(t,y):=[[0,1],[-1,-t]]*y|
et \verb|heun(f,0,1,10,[1,0])| renvoie \\
$[0.605257530382,-0.605257530382]$
%(valeur exacte exp(-t^2/2))
\end{enumerate}

\end{document}


