Commit 92249c39 authored by Conor McCoid's avatar Conor McCoid

Reorganized folder, added assistant work, should now be easier to manage

parent 79f8e72a
function A=FDLaplacian(n,d)
% FDLAPLACIAN compute a finite difference Laplacian
% A=FDLaplacian(n,d) computes a finite difference Laplacian on the
% unit interval/square/cube (d=1,2,3) using n interior points
h=1/(n+1);
D=1/h^2*spdiags(ones(n,1)*[1 -2*d 1],[-1 0 1],n,n);
if d==1
A=D;
elseif d==2
X=spdiags(ones(n,1),0,n,n);
I2=1/h^2*spdiags(ones(n^2,1)*[1 1],[-n n],n^2,n^2);
A=kron(X,D)+I2;
elseif d==3
X=spdiags(ones(n^2,1),0,n^2,n^2);
I2=1/h^2*spdiags(ones(n^3,1)*[1 1],[-n n],n^3,n^3);
I3=1/h^2*spdiags(ones(n^3,1)*[1 1],[-n^2 n^2],n^3,n^3);
A=kron(X,D)+I2+I3;
end
% avec blkdiag
%if d==1
% A=D;
%elseif d==2
% A=[];
% for i=1:n
% A=blkdiag(A,D);
% end
% I2=1/h^2*spdiags(ones(n^2,1)*[1 1],[-n n],n^2,n^2);
% A=A+I2;
%elseif d==3
% C=[];
% for i=1:n
% C=blkdiag(C,D);
% end
% A=[];
% for i=1:n
% A=blkdiag(A,C);
% end
% I2=1/h^2*spdiags(ones(n^3,1)*[1 1],[-n n],n^3,n^3);
% I3=1/h^2*spdiags(ones(n^3,1)*[1 1],[-n^2 n^2],n^3,n^3);
% A=A+I2+I3;
%end
\documentclass[12pt,a4paper]{article}
\usepackage{amsmath,amssymb,amsthm,epsfig,amscd,verbatim}
\usepackage[french]{babel}
\usepackage{fancyvrb}
%\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\addtolength{\textwidth}{3cm}
\addtolength{\textheight}{4cm}
\addtolength{\hoffset}{-2cm}
\addtolength{\voffset}{-2cm}
%% commandes utiles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\tr}{\mbox{tr}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\noindent UNIVERSIT\'E DE GEN\`EVE \hfill Section de Math\'ematiques \\
\noindent Facult\'e des Sciences \hfill \`A rendre le mercredi 4 mars 2020 \\[-3mm] \hrule
\bigskip
\begin{center}
\textbf{\Large{Méthodes Itératives}} \\[1mm]\bigskip
\textbf{Série 1}
\end{center}
\bigskip
\begin{enumerate}
\item (8 pts) On rappel que pour $x\in\C^m$ et $1\leq p<\infty$, on définit
les normes suivantes,
\[
\|\vec{x}\|_p=\left(\sum_{i=1}^{m}|x_i|^p\right)^{\frac{1}{p}},\quad
\|\vec{x}\|_\infty=\max_{1\leq i\leq m}\{|x_i|\}.
\]
Et pour une matrice $A\in\C^{n\times m}$, nous utiliserons les normes,
\[
\|A\|_{pq}=\sup_{\vec{x}\neq0}\frac{\|A\vec{x}\|_p}{\|\vec{x}\|_q},
\]
pour $1\leq p,q\leq\infty$. Dans le cas particulier $p=q$, nous abrègerons
l'écriture par $\|\cdot\|_{pp}=\|\cdot\|_p$. On utilisera aussi la norme de {\it
Frobenius},
\[
\|A\|_F^2=\sum_{i=1}^n\sum_{j=1}^m|a_{ij}|^2.
\]
Montrer les quatres relations suivantes, où $\rho(A):=\max\{|\lambda_i|\}$ et
$\lambda_i$ sont les valeurs propres de $A$:
\[
\begin{array}{ll}
\displaystyle{\|A\|_1=\max_{1\leq j\leq m}\sum_{i=1}^n|a_{ij}|}, &
\displaystyle{\|A\|_\infty=\max_{1\leq i\leq n}\sum_{j=1}^m|a_{ij}|},\\
\displaystyle{\|A\|_2=(\rho(A^*A))^{\frac{1}{2}}}, &
\displaystyle{\|A\|_F=(\tr(A^*A))^{\frac{1}{2}} = (\tr(AA^*))^{\frac{1}{2}}}. \\
\end{array}
\]
\item (2 pts) Montrer la propriété de submultiplicité
\[
\|AB\|_F\leq \|A\|_F \|B\|_F.
\]
\item (8 pts) Montrer que pour tout $\epsilon>0$, il existe une norme de $\C^{n \times n}$ telle que
\[
\rho(A)\leq \|A\| \leq \rho(A)+\epsilon.
\]
Montrer que pour les normes $\|\cdot\|_p$ avec $1\leq p\leq\infty$, on a
\[
\rho(A)=\lim_{k\to\infty}\|A^k\|_p^{\frac{1}{k}}.
\]
%\item ( pts) Ecrire les méthodes classiques {\it Jacobi}, {\it Gauss-Seidel},
% {\it SOR} et {\it Richardson} dans la forme de correction vu au cours.
\item (2 pts) Considérer un maillage uniforme $\Omega_h$ dans $\mathbb{R}^2$ de taille $h$ et tel que
$(x_i,y_j) \in \Omega_h$ avec $x_{i+1} = x_i+h$ et $y_{j+1} = y_j+h$.
A l'aide de la série de Taylor, montrer que
%pour une fonction $u : \mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}$ suffisamment differentiable on a
%\begin{equation*}
%\begin{split}
%\Delta u(x,y) = \partial_{xx} u(x,y) + \partial_{yy} u(x,y) &= \frac{u(x+\Delta x,y)-2u(x,y)+u(x-\Delta x,y)}{\Delta x^2} + %O(\Delta x^2) \\
%&+ \frac{u(x,y+\Delta y)-2u(x,y)+u(x,y-\Delta y)}{\Delta y^2} + O(\Delta y^2), \\
%\end{split}
%\end{equation*}
%et donc
le schéma $D^2(u)$ des différences finies du Laplacien discret défini sur $\Omega_h$ est donné par
\begin{equation*}
D^2(u) = \frac{u_{i+1,j}+u_{i,j+1}-4u_{i,j}+u_{i-1,j}+u_{i,j-1}}{h^2},
\end{equation*}
et que ce schéma est d'ordre deux. % quelle que soit la dimension du domaine $d=1,2,3$.
\item (4 pts) \'Ecrire une fonction Matlab qui prend comme arguments le nombre de mailles et
la dimension du domaine et qui retourne la matrice du Laplacien discret avec
l'en-tête suivante,
\VerbatimInput[firstline=0,lastline=5]{FDLaplacian.m}
\end{enumerate}
{\bf \'Evaluation:}
\begin{itemize}
\item[$\bullet$] Les exercices
\item[$\bullet$] Un examen oral durant la session d'examens sur le cours.
\end{itemize}
La note finale est de : $30 \%$ exercices et $70\%$ examen oral.\\
Les exercices sont notés $15\%$ pour la présentation de séries et
$15\%$ pour la présentation au tableau. Il faut passer au tableau au
moins 2 fois dans le semestre. % ??
\end{document}
function u=Heat2D(f,g,p)
% HEAT2D solves the stationary heat equation over a square
% u=Heat2D(f,g,p) solves the heat equation over a square of size one
% by one. The source function is given by f, the boundary condition is
% given by b, the number of grid points for each direction is a quarter
% of the length(g) and one chooses the iterative solver with p,
% i.e. p=0 uses Jacobi, p=1 uses Gauss-Seidel.
\documentclass[12pt,a4paper]{article}
\usepackage{amsmath,amssymb,amsthm,epsfig,amscd,verbatim}
\usepackage[french]{babel}
\usepackage{fancyvrb}
%\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\addtolength{\textwidth}{3cm}
\addtolength{\textheight}{4cm}
\addtolength{\hoffset}{-2cm}
\addtolength{\voffset}{-2cm}
\usepackage{hyperref}
%% commandes utiles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\tr}{\mbox{tr}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\noindent UNIVERSIT\'E DE GEN\`EVE \hfill Section de Math\'ematiques \\
\noindent Facult\'e des Sciences \hfill \`A rendre le mercredi 22 mars \\[-3mm] \hrule
\bigskip
\begin{center}
\textbf{\Large{Méthodes Itératives}} \\[1mm]\bigskip
\textbf{Série 2}
\end{center}
\bigskip
\begin{enumerate}
\item (4 pts) Montrer les deux égalités suivantes pour $p$ un entier
positif.
\[
\lim_{k\to\infty}\binom{k}{p}^{\frac{1}{k}} = 1,
\]
et pour $0<\rho<1$,
\[
\lim_{k\to\infty}\binom{k}{p}\rho^k = 0.
\]
\item (4 pts) Soit $\vec{B}\in\C^{n\times n}$ une matrice carrée dont le
rayon spectral $\rho(\vec{B})$ vérifie $\rho(\vec{B})<1$. Démontrer que
$\vec{I}-\vec{B}$ est inversible et que
\[
(\vec{I}-\vec{B})^{-1} = \sum_{j=0}^\infty\vec{B}^j.
\]
Cette série est un cas particulier des séries de Neumann (généralisation aux
matrices carrées des séries entières).
\item (6 pts) Démontrer que la matrice $A$ du $-\Delta$ discret en dimensions un et deux
est une M-matrice.\\
\underline{Indice pour dimension un}: vous pouvez considérer un splitting $A = M-N$, démontrer que $\rho(M^{-1}N)<1$
et utiliser les résultats démontrés en cours.
Notons que pour une matrice tridiagonale il existe une formule pour calculer les valeurs propres.\\
\underline{Indice pour dimension deux}:
\begin{itemize}\itemsep0em
\item Soit $A_x$ et $A_y$ les matrices de $-\frac{d^2}{dx^2}$ et $-\frac{d^2}{dy^2}$ discrètes.
Alors, la matrice $A$ peut \^etre obtenue comme $A = A_x \otimes I + I \otimes A_y$,
o\`u $\otimes$ est le produit de Kronecker
(voir \url{https://fr.wikipedia.org/wiki/Produit_de_Kronecker} pour une explication détaillée).
\item Soit $(\lambda_x,v_x)$ et $(\lambda_y,v_y)$ deux valeurs/vecteurs propres de $A_x$ et $A_y$.
Demontrer que les valeurs propres de $A$ sont données par $\lambda = \lambda_x + \lambda_y$.
Pour montrer cette propriété, considérer le vecteur $w:=v_x \otimes v_y$ et étudier le produit
$A w = \bigl( A_x \otimes I + I \otimes A_y \bigr) \bigl( v_x \otimes v_y \bigr)$ en utilisant
la formule $(A \otimes B)(C \otimes D) = (AC) \otimes (BD)$.
\item Considérer un splitting $A = M-N$, démontrer que $\rho(M^{-1}N)<1$
et utiliser les résultats démontrés en cours.
\end{itemize}
\item (3 pts) Ecrire les méthodes itératives de {\it Jacobi}, {\it Gauss-Seidel}
et {\it SOR} sous les deux formes standard, de correction, de résidu et de différences.
\item (6 pts) Implémenter une fonction Matlab qui résout l'équation de la
chaleur stationnaire
\begin{equation*}
\begin{split}
\Delta u &= f \, \text{ dans $\Omega = (0,1)\times(0,1)$}, \\
u &= g \; \text{ sur $\partial \Omega$}, \\
\end{split}
\end{equation*}
en dimension deux à l'aide des méthodes de {\it Jacobi} et {\it
Gauss-Seidel}. Les solveurs itératifs seront implémentés dans des fonctions
externes au programme résolvant l'équation de la chaleur. Faire des en-têtes
claires et précises dans la description suivant le format ci-dessous.
La fonction Heat2D qui résout l'équation prendra en argument la source $f$,
le choix du solveur et le nombre de points de grille et résoudra la chaleur
sur un domaine carré de côtés de longueur un. Son en-tête doit ressembler à
la suivante,
\VerbatimInput[firstline=0,lastline=7]{Heat2D.m}
\end{enumerate}
{\bf \'Evaluation:}
\begin{itemize}
\item[$\bullet$] Les exercices
\item[$\bullet$] Un examen oral durant la session d'examens sur le cours.
\end{itemize}
La note finale est de : $30 \%$ exercices et $70\%$ examen oral.\\
Les exercices sont notés $15\%$ pour la présentation de séries et
$15\%$ pour la présentation au tableau. Il faut passer au tableau au
moins 3 fois dans le semestre.
\end{document}
\documentclass[12pt,a4paper]{article}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage[french]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{bm}
\usepackage{paralist}
\usepackage{graphicx}
\usepackage{lmodern}
\usepackage{verbatim}
\usepackage{fancyhdr}
\pagestyle{empty}
\setdefaultenum{\bfseries Corr{i}gé exerc{i}ce 1\mdseries}{a)}{i}{\textbullet}
\begin{document}
\noindent {\large UNIVERSITE DE GENEVE} \hfill Section de math\'ematiques \\
\noindent Facult\'e des sciences \hfill \\
\hrule
\bigskip
\bigskip
\begin{center}{\bf Méthodes itératives --
Corrigé série 3}\end{center}
\bigskip
\parskip 12pt
\begin{enumerate}
\item \textit{Sur 5 points} On considère un maillage régulier et rectangulaire
à $n_x$ lignes et $n_y$ colonnes.
Pour une fonction réelle $u$ définie sur
le domaine rectangulaire, on note par $u_{i,j}$ la
valeur $u(ih,jh)$ pour $1\leq i\leq n_x$ et
$1\leq j\leq n_y$. Si $u$ est suffisamment
régulière alors par la formule d'approximation de Taylor:
\begin{align*}
u_{i+1,j}&=u_{i,j}+h(\partial_xu)_{i,j}
+\frac{h^2}{2}(\partial_x^2u)_{i,j}
+\frac{h^3}{6}(\partial_x^3u)_{i,j}
+O(h^4),
\\
u_{i-1,j}&=u_{i,j}-h(\partial_xu)_{i,j}
+\frac{h^2}{2}(\partial_x^2u)_{i,j}
-\frac{h^3}{6}(\partial_x^3u)_{i,j}
+O(h^4),
\\
u_{i,j+1}&=u_{i,j}+h(\partial_yu)_{i,j}
+\frac{h^2}{2}(\partial_y^2u)_{i,j}
+\frac{h^3}{6}(\partial_y^3u)_{i,j}
+O(h^4),
\\
u_{i,j-1}&=u_{i,j}-h(\partial_yu)_{i,j}
+\frac{h^2}{2}(\partial_y^2u)_{i,j}
-\frac{h^3}{6}(\partial_y^3u)_{i,j}
+O(h^4),
\end{align*}
Il suffit alors de sommer ces inégalités puis de
diviser par $h^2$ pour obtenir
la formule aux différences finis:
\begin{equation*}
(\triangle u)_{i,j}=(A_hu)_{i,j}
+O(h^2).
\end{equation*}
\begin{equation}
(A_hu)_{i,j}=
u_{i+1,j}+
u_{i-1,j}+
u_{i,j+1}+
u_{i,j-1}+
-4u_{i,j}
\end{equation}
et avec les conventions
$u_{i,0}=u_{i,n_y+1}=u_{0,j}=u_{n_x+1,j}=0$. $A_h$ est l'opérateur du
Laplacien discret, c'est un opérateur linéaire et a donc une
représentation matricielle lorsque on choisit de représenter
le vecteur $u$ sous forme de vecteur colonne, \textit{i.e}
par renumérotation du maillage. Si la numérotation s'effectue en
colonnes, cette matrice est la matrice $\mathbf{A}$ de l'énoncé.
\item\textit{Sur 5 points}
La matrice $\mathbf{A}$ a la propriété $A$ donc par l'indication
(SOR avec $\omega=1$ est Gauss-Seidel)
$\lambda^2=\lambda\mu^2$ où les $\lambda$ sont
les valeurs propres de la matrice d'itération de Gauss-Seidel
$\mathbf{B}_{GS}$ et $\mu$ celles de la matrice d'itération de Jacobi
$\mathbf{B}_J$ donc $\rho(\mathbf{B}_{GS})=\rho(\mathbf{B}_J)$.
\item\textit{Sur 10 points dont 8 pour le a), 1 pour le b), 1 pour le c)}
\begin{enumerate}
\item Il faut s'inspirer du cas continu. En effet, sur le rectangle
$(O,L)\times(0,H)$, les fonctions
$(x,y)\mapsto\sin(k\pi x/L)\sin(l\pi y/L)$
sont des vecteurs propres du Laplacien. Cela conduit à
considérer les vecteurs $\bm{u}^{k,l}\in\mathbb{R}^{nx\times n_y}$,
$1\leq k\leq n_x$, $1\leq l\leq n_y$, définis par
\begin{equation*}
\bm{u}^{k,l}_{i,j}=
\sin(\frac{ik\pi}{n_x+1})\sin(\frac{jl\pi}{n_y+1}).
\end{equation*}
Utilisons maintenant la formule aux différences finies du corrigé de
l'exercice 1. On obtient, après applications des formules
trigonométriques usuelles:
\begin{align*}
(A_h\bm{u}^{k,l})_{i,j}&=
\frac{1}{h^2}
\left(
\sin(\frac{(i+1)k\pi} {n_x+1})
+\sin(\frac{(i-1)k\pi}{n_x+1})
-2\sin(\frac{ik\pi}{n_x+1})
\right)
\sin(\frac{jl\pi}{n_y+1})\\
&\phantom{=}+\frac{1}{h^2}
\left(
\sin(\frac{(j+1)l\pi}{n_y+1})
+\sin(\frac{(j-1)l\pi }{n_y+1})
-2\sin(\frac{jl\pi}{n_y+1})
\right)
\sin(\frac{ik\pi}{n_x+1})
\\
&=\frac{2}{h^2}
(\cos(k\pi/(n_x+1)-1)
+\cos(l\pi/(n_y+1)-1)
)
u^{k,l}_{i,j}\\
&=-\frac{4}{h^2}
\big(\sin^2(\frac{k\pi}{2(n_x+1)})
+\sin^2(\frac{l\pi}{2(n_y+1)})
\big)
u^{k,l}_{i,j}.
\end{align*}
Nous avons $n_x\times n_y$ vecteurs propres formant
une partie libre. Nous avons donc comptabilisé toutes les valeurs propres qui sont
bien les $\lambda_{i,j}$.
\item La matrice d'itération de Richardson
est donnée par $\mathbf{I}-\alpha\mathbf{A}$.
Si on note par $\lambda_{\mathrm{max}}$ la plus grande valeur propre
et $\lambda_{\mathrm{min}}$ la plus petite alors d'après le cours:
\begin{equation*}
\alpha_{\mathrm{opt}}=\frac{2}{\lambda_{\mathrm{min}}+\lambda_{\mathrm{max}}}.
\end{equation*}
Utilisant le a), on obtient:
\begin{multline*}
\lambda_{\mathrm{min}}+\lambda_{\mathrm{max}}=\frac{4}{h^2}
\big(\sin^2(\frac{\pi}{2(n_x+1)})
+\sin^2(\frac{\pi}{2(n_y+1)})\\
+\sin^2(\frac{n_x\pi}{2(n_x+1)})
+\sin^2(\frac{n_yl\pi}{2(n_y+1)})
\big)
=\frac{8}{h^2}.
\end{multline*}
D'où $\alpha_{\mathrm{opt}}=\frac{h^2}{4}$.
\item
La matrice d'itération de Jacobi est
\begin{equation*}
\mathbf{B}_J=\mathbf{I}+\frac{h^2}{4}\mathbf{A}.
\end{equation*}
Ses valeurs propres sont donc les
$\frac{1}{2}
(\cos(k\pi/(n_x+1))
+\cos(l\pi/(n_y+1))
)$ pour $1\leq i\leq n_x$,
$1\leq i\leq n_x$. Donc,
$\rho(\mathbf{B}_J)=\frac{1}{2}
(\cos(\pi/(n_x+1))
+\cos(\pi/(n_y+1))
)$.
Comme la matrice du Laplacien discrétisé a la propriété $A$,
il suffit maintenant d'appliquer la formule:
\begin{equation*}
w_{\mathrm{opt}}=
\frac{2}{1+\sqrt{1-\rho(\mathbf{B}_J)^2}}
\end{equation*}
\end{enumerate}
\item\textit{Sur 20 points dont 10 (5 pour SOR et 5 pour Richardson)
pour le a) , 5 pour le b) et 5 pour le c)}
Fonction Richardson:
\verbatiminput{Richardson.m}
Fonction SOR:
\verbatiminput{SOR.m}
Script d'exécution:
\verbatiminput{serie3.m}
\noindent\\
Figures obtenues:
\begin{center}
\includegraphics[width=0.85\textwidth]{serie3fig1.eps}
\includegraphics[width=0.85\textwidth]{serie3fig2.eps}
\end{center}
\item \textit{Sur 5 points} Il suffit de construire une suite de matrice
$(\mathbf{A}_k)_{0\leq k\leq n-2}$ et une suite de matrices de
Househölder $(\mathbf{P}_k)_{1\leq k\leq n-2}$ tels que
\begin{enumerate}
\renewcommand{\labelitemi}{\bullet}
\item $\mathbf{A}_0=\mathbf{A}$.
\item $\mathbf{A}_k=\mathbf{P}_k\mathbf{A}_{k-1}\mathbf{P}_k^{\mathrm{T}}$.
\item \label{item:FormePk} $\mathbf{P}_k=\mathbf{I}-2\bm{w}_k\bm{w}_k^{\mathrm{T}}$
avec $\lVert\bm{w}_k\rVert_2=1$ et tels que les
$k$ premières composantes de
$\bm{w}_k$ soient nulles.
\item
$\mathbf{A}_k$ est de Hessenberg sur ces $i$ premières colonnes,
\textit{i.e.} $(\mathbf{A}_k)_{i,j}=0$ pour $i\geq j+1$, $i\leq k$.
\end{enumerate}
En effet, dans ce cas $\mathbf{A}_{n-2}$ est de Hessenberg et
$\mathbf{A}_{n-2}=\mathbf{Q}\mathbf{A}\mathbf{Q}^{\mathrm{T}}$
avec $\mathbf{Q}=\mathbf{P}_1\ldots\mathbf{P}_{n-2}$.
Pour construire cette suite, il suffit de procéder par récurrence.
On suppose $\mathbf{A}_{k-1}$ construit. Remarquons d'abord qu'étant donné