Commit 813e9d2c authored by Conor McCoid's avatar Conor McCoid
Browse files

Asst: added TP06 for Anal. Num. 2021-2022

parent 4d53528d
This diff is collapsed.
\documentclass{article}
\usepackage{geometry}
\geometry{noheadfoot, margin=0.5in}
\usepackage{prerex}
\begin{document}
\thispagestyle{empty}
\begin{chart}
\end{chart}
\end{document}
\documentclass[11pt,a4paper]{article}
\usepackage{amsmath} % Pour les symboles complementaire comme les matrices !
\usepackage{amssymb}
\usepackage{verbatim}
\usepackage{epsfig}
%\usepackage[french]{babel}
\usepackage{color}
\usepackage[latin1]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[frenchb]{babel}
%\usepackage[T1]{fontenc} % Pour la bonne cesure du francais
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{example}{Example}
\newcommand{\qed}{\hfill$\qedsquare$\goodbreak\bigskip}
\advance\voffset by -35mm \advance\hoffset by -25mm
\setlength{\textwidth}{175mm} \setlength{\textheight}{260mm}
\pagestyle{empty}
\newdimen\figcenter
\newdimen\figlarge
\newdimen\figheight
\def\figinsert#1{
\figcenter=\hsize\relax
\advance\figcenter by -\figlarge
\divide\figcenter by 4
\vskip +1.truecm
\vskip\figheight\noindent\hskip\figcenter
\special{psfile=#1}
\vskip -.7truecm}
\begin{document}
\noindent {\large UNIVERSIT\'E DE GEN\`EVE} \hfill Section de math\'ematiques \\
\noindent Facult\'e des Sciences \hfill \`A rendre le jeudi 16 decembre 2021 \\
\null \hfill L'\'evaluation aura lieu le vendredi 17 decembre 2021 \\
\hrule
\bigskip
\bigskip
\begin{center} {\bf {\Large Analyse Num\'erique}}\end{center}
\begin{center} {\bf {\Large Travaux Pratiques~-~S\'erie 6}}\end{center}
\bigskip
\paragraph{Exercice 1 : Le pivot de Gauss.}
On aimerait utiliser l'algorithme d'\'elimination de Gauss
avec recherche partielle de pivot
pour r\'esoudre le syst\`eme lin\'eaire
\begin{equation} \label{eq:linsys}
Ax=b,
\end{equation}
o\`u $A$ est une matrice carr\'ee.
\begin{enumerate}
\item D'abord, il faut écrire des ``solvers'' pour résoudre des systèmes
triangulaires. Pour un système $Ux = b$, où $U$ est une matrice triangulaire supérieure,
utiliser l'en-tête suivante:
{\small
\begin{verbatim}
function y = usol(U,b)
% USOL solves the system Ux = b for upper triangular U
% y = usol(U,b) solves the system Ux = b, where U is a
% upper triangular matrix.
\end{verbatim}
}
Voici l'en-tête pour $Lx=b$, où $L$ est une matrice triangulaire inférieure:
{\small
\begin{verbatim}
function y = lsol(L,b)
% LSOL solves the system Lx = b for lower triangular L
% y = lsol(L,b) solves the system Lx = b, where L is an
% lower triangular matrix.
\end{verbatim}
}
Pour vérifier vos programmes, créer d'abord des matrices triangulaires
aléatoires $L$ et $U$. (Les commandes {\tt rand}, {\tt tril} et {\tt triu}
peuvent être utiles.) Ensuite, créer la solution exacte
{\tt x\_exact} et calculer {\tt b = L*x\_exact} et {\tt c = U*x\_exact}. Quelle est la
différence entre {\tt x\_exact} et les résultats de {\tt lsol}/{\tt usol}?
\item Ensuite, écrire un programme qui calcule la factorisation LU sans
recherche partielle de pivot:
{\small
\begin{verbatim}
function [L,U] = LUNoPivot(A)
% LUNOPIVOT computes LU factorization of A without pivoting
% [L,U] = LUNoPivot(A) computes the LU factorization of A using
% Gaussian elimination without pivoting. L is a unit lower triangular
% matrix and U is an upper triangular matrix.
\end{verbatim}
}
Appliquer à la matrice
\[
A=\begin{bmatrix}
4 & -3 & -1 \\
0.5 & 3 & -1 \\
0 & 1 & 0
\end{bmatrix},
\]
puis comparer avec les matrices $L$ et $U$ données par la commande Matlab {\tt lu(A)}. Pour la matrice $A$ proposée, la matrice de permutation $P$ est simplement l'identité (Pourquoi?). Ensuite, utiliser {\tt lsol} et {\tt usol} pour résoudre $Ax = b$ avec un vecteur $b$ aléatoire.
\item Rajouter la recherche partielle de pivot en modifiant le programme
ci-dessus (voir section 5.5). Utiliser l'en-tête suivante:
{\small
\begin{verbatim}
function [L,U,p] = LUPivot(A)
% LUPIVOT computes LU factorization of A with partial pivoting
% [L,U,p] = LUPivot(A) computes the LU factorization of A with partial
% pivoting. L is a unit lower triangular matrix, U is an upper triangular
% matrix, and p is a vector of permutations where p(i) is the row of A
% from which the i-th pivot is selected. In other words, p is computed such
% that solving A*x = b is equivalent to solving L*U*x = b(p).
\end{verbatim}
}
Pour ce faire, il faut faire très attention lors des échanges des équations!
Les facteurs $L$ et $U$, ainsi que le vecteur $p$, doivent être modifiés
avec chaque échange. Comment peut-on utiliser {\tt lsol}/{\tt usol} pour résoudre
$Ax = b$ dans ce cas général?
\item Vérifier le programme {\tt LUPivot} sur les exemples ci-dessous:
\begin{align*}
&\mbox{(a)}
\quad
A=
\left( \begin{array}{ccc}
1 & 5 & 2 \\
0 & 0 & 8 \\
2 & 4 & 1
\end{array}\right),
\quad
b=
\left( \begin{array}{c}
9 \\
8 \\
9
\end{array} \right)\\
&\mbox{(b)}
\quad
A=\left( \left( \frac{j}{n} \right)^{i-1} \right)_{i,j=1}^{n},
\quad
b=
\left( \begin{array}{c}
1 \\
1/2 \\
\vdots \\
1/n
\end{array} \right).
\end{align*}
Pour $n=8,10,12$, \`a quelle pr\'ecision
num\'erique du r\'esultat peut-on s'attendre?
\item Comparer les résultats obtenus en résolvant $Ax=b$ en utilisant d'abord {\tt LUPivot} puis la commande matlab " \textbackslash ". Qu'observez-vous ? Choisissez n'importe quels matrices $A$ et $b$.
\item V\'erifier que le pivot trouvé par {\tt LUPivot} est le même que celui utilisé
par l'algorithme {\tt LUNoPivot} pour la matrice suivante de taille $n\times
n$ (matrice impliquée dans le calcul de splines). Comment expliquer cela ?
\[
A=\frac{1}{h}\begin{bmatrix}
4 & 1 & & & & \\
1 & 4 & 1 & & & \\
& \ddots & \ddots & \ddots & & \\
& & 1 & 4 & 1 & \\
& & & 1 & 4 & \\
\end{bmatrix}
\]
o\`u $h=\frac{1}{n}$, pour certains $n$ grand.
\end{enumerate}
\paragraph{Exercice 2 : La règle de Cramer, une méthode peu efficace
numériquement.}
On souhaite maintenant résoudre le système linéaire (\ref{eq:linsys}) en utilisant la règle de Cramer
$$
x_i = \frac{\det(A_i)}{\det(A)}, \quad i=1\ldots n,
$$
$A_i$ est la matrice obtenue en remplaçant la $i$ème colonne de $A$ par le vecteur $b$.
La fonction déterminant $\det$ sera calculée récursivement en utilisant le développement de Laplace:
$$
\det(A)=\sum_{j=1}^n (-1)^{1+j} a_{1,j} \det{M_{1,j}}
$$
$M_{1,j}$ est la matrice obtenue à partir de la matrice $A$ en supprimant la ligne $i$ et la colonne $j$
(par convention, le déterminant de la matrice vide est égal à $1$).
\begin{enumerate}
\item
\emph{Coût des calculs.}
\begin{enumerate}
\item
Écrivez un programme pour calculer le déterminant par la formule de Laplace
%What is the number of arithmetic operations of the algorithm?
et testez-le en utilisant des matrices aléatoires (\texttt{rand} en \textsc{Matlab}) de dimension $n=3,5,7,9$.
\item
Écrivez un programme pour calculer $\det(A)$ pour une matrice $A$ raisonnablement grande, par exemple $n \geq 100$, sans utiliser la formule de Laplace.
\item
Estimer le coût des calculs pour résoudre un système linéaire en utilisant la règle de Cramer.
\end{enumerate}
\item
\emph{Stabilité numérique.}
Pour une matrice donnée $A$, vous pouvez suivre les instructions ci-dessous pour tester la stabilité numérique de la résolution de systèmes linéaires.
\begin{enumerate}
\item[(1)]
Générer une solution ``exacte'' $x=\mathtt{randn}(n,1)$, et le membre de droite correspondant $b:=Ax$.
\item[(2)]
On appelle $\widehat x$ la solution numérique de $Ax=b$
par la règle de Cramer. Calculer l'erreur forward et les normes des résidus:
\[
\mathtt{err} :=\frac{\|\widehat x - x\|_2}{\|x\|_2},\qquad
\mathtt{res}:=\frac{\|A\widehat x-b\|_2}{\|A\|_2\|x\|_2+\|b\|_2}.
\]
\end{enumerate}
\medskip
Pour simplicité, on considère le système linéaire $2\times 2$
$Ax=b$, pour lequel la règle de Cramer devient
\begin{align*}
x_1 = (b_1a_{22} - b_2a_{12})/d,\qquad x_2 = (b_2a_{11} - b_1a_{21})/d
\end{align*}
$d=a_{11}a_{22} - a_{21}a_{12}$.
Utiliser la fonction \texttt{matgen} (disponible sur Chamilo) pour générer des matrices $A$ de dimension $n=2$, et condition
$\kappa=10^2,10^6,10^{10},10^{18}$, respectivement.
Observer les erreurs $\texttt{err}$ et $\texttt{res}$ pour la règle de Cramer, et comparer le résultat avec l'algorithme de l'élimination de Gauss de l'Exercice~1, ainsi qu'avec la fonction de \uppercase{MATLAB} backslash
$\backslash$.
\end{enumerate}
\paragraph{Exercice 3 : Inversion-et-multiplication.} Dans cet exercice on examine l'idée de résolution d'un système linéaire $Ax=b$ par inversion directe: d'abord calculer la matrice inverse $M=A^{-1}$, puis appliquer la multiplication matrice-vecteur $x=Mb$.
\begin{enumerate}
\item
Écrire un programme pour résoudre efficacement
$x^{(1)},x^{(2)},\dots,x^{(k)}$ d'un ensemble de $k$ systèmes linéaires
\[
Ax^{(i)} = b^{(i)}\qquad\text{for}\qquad i=1,\dots,k,
\]
$b^{(i)}\in \mathbb R^n$ sont des membres de droite donnés.
\emph{Indice : Comme $A$ est la même pour tout les systèmes linéaires, on peut réutiliser sa factorisation LU.}
\item
Modifier votre programme pour calculer la matrice inverse d'une matrice $A$.
\emph{Indice: $A^{-1}$ est la solution de $AX = I$, où $I$ est une matrice identité.}
\item
%Solve the linear system by invert-and-multiply: $x=A^{-1}b$.
Étudier la stabilité de l'algorithme inversion-et-multiplication comme dans l'Exercice 2.2, en utilisant des matrices de dimension $n=100$ et condition
$\kappa=10^2,10^6,10^{10},10^{18}$.
\end{enumerate}
\paragraph{Exercice 4 : Équation de Poisson en dimension 1}
Prenons une fine barre métallique homogène de longueur $L$ et de section constante. On suppose que la barre est thermiquement isolée, de telle manière que la chaleur se propage uniquement suivant la barre. On applique une source thermique extérieure $H(x)$ tout au long de la barre. On pose $u(x,t)$ la température au point $x$ de la barre et au temps $t$. Alors en appliquant le principe de conservation de l'énergie, on obtient l'équation suivante, dite de la chaleur 1D:
\begin{equation}
\frac{\partial u(x,t)}{\partial t} =\alpha^2 \, \frac{\partial^2 u(x,t)}{\partial^2 x} +H(x), \qquad x\in [0,L], \quad t>0, \notag
\end{equation}
$\alpha^2$ est le coefficient de diffusion. On suppose aussi que les extrémités de la barre sont maintenues à une température nulle,
\begin{equation}
u(0,t) = u(L,t)=0. \notag
\end{equation}
Si la température initiale est donnée par $g(x)$, alors on a
\begin{equation}
u(x,0)= g(x) \qquad \mbox{with}\qquad g(0)=g(L)=0. \notag
\end{equation}
À présent, on s'intéresse aux solutions stationnaires de cette équation, c'est à dire les solutions vérifiant $\frac{\partial u(x,t)}{\partial t}=0$. L'équation de la chaleur stationnaire ainsi obtenue est appelée équation de Poisson. À présent, la température ne dépend plus du temps $t$. Si l'on appelle la température $\phi(x)$ et la source de chaleur extérieure $F(x)$, l'équation de Poisson est
\begin{eqnarray}\label{PE}
-\, \frac{\partial^2 \phi(x)}{\partial^2 x} &=& F(x), \qquad x\in [0,L], \\
\phi(0) &=& \phi(L)=0\notag.
\end{eqnarray}
Pour résoudre numériquement \eqref{PE}, on divise l'intervalle $[0,L]$ en $n$ sous-intervalles de longueur $h=L/n$ avec $0=x_{0}<x_{1}<...<x_{n}<x_{n}=L$. On utilise alors un schéma de différences finies pour approcher $\phi_{xx}:= \frac{\partial^2 \phi(x)}{\partial^2 x}$, et on obtient finalement
\[\frac{1}{h^2} \begin{pmatrix} 2 & -1\\-1 & 2 & -1\\ & -1 & \ddots & \ddots\\
& & \ddots & \ddots & -1\\
& & & -1 & 2 \end{pmatrix} \begin{pmatrix}
\phi_{1} \\
\phi_{2}\\
\vdots\\
\phi_{n-2}\\
\phi_{n-1}
\end{pmatrix} = \begin{pmatrix}
F_{1} \\
F_{2}\\
\vdots\\
F_{n-2}\\
F_{n-1}
\end{pmatrix},
\]
$F_{i}=F(x_{i})$.
\begin{enumerate}
\item Résoudre le système d'équations tridiagonal précédent en utilisant l'opération Matlab $"\textbackslash$" pour $F(x) = -2$, $L=1$ et $n=10,20,50,100$. Comparer avec la solution exacte.
\item Utiliser votre fonction {\tt LUNoPivot} pour observer la structure des matrices de la décomposition LU de la matrice de Poisson vue ci-dessus. Modifier votre fonction {\tt LUNoPivot.m} pour le rendre plus rapide pour ce système d'équations tridiagonal.
\end{enumerate}
\begin{figure}[ht]
\centering
\includegraphics[scale =0.6]{n10}
\label{fig1}
\end{figure}
\hspace{1cm}
\begin{figure}[ht]
\centering
\includegraphics[scale =0.6]{n100}
\label{fig2}
\caption{Solutions de l'équation de Poisson avec $F(x)=-2$ pour $n=10$ et $100$.}
\end{figure}
\end{document}
function [L,U] = LUNoPivot(A)
% LUNOPIVOT computes LU factorization of A without pivoting
% [L,U] = LUNoPivot(A) computes the LU factorization of A using
% Gaussian elimination without pivoting. L is a unit lower triangular
% matrix and U is an upper triangular matrix.
n = length(A(:,1));
L = diag(ones(n,1));
U = A;
for j = 1 : n-1
for i = j+1:n
l = U(i,j)/U(j,j);
L(i,j) = l;
U(i,j:n) = U(i,j:n) - U(j,j:n)*l;
end
end
function [L,U,p] = LUPivot_ibrahim(A)
% LUPIVOT computes LU factorization of A with partial pivoting
% [L,U,p] = LUPivot(A) computes the LU factorization of A with partial
% pivoting. L is a unit lower triangular matrix, U is an upper triangular
% matrix, and p is a vector of permutations where p(i) is the row of A
% from which the i-th pivot is selected. In other words, p is computed such
% that solving A*x = b is equivalent to solving L*U*x = b(p).
% ~this algorithm is the one found in the polycpie.
[m,n]=size(A);
if m~=n
disp('Error: Please enter a square matrix!');
return;
end
U=A; L=eye(m); p=1:m;
for k=1:m-1
[pivot,i]=max(abs(U(k:m,k))); % serching for pivot
i=i+k-1; % to get the correct index
if pivot~=0
U([k,i],k:m)=U([i,k],k:m); % swaping rows
L([k,i],1:k-1)=L([i,k],1:k-1);
p([k,i])=p([i,k]);
for j=k+1:m
L(j,k)=U(j,k)/U(k,k);
U(j,k:m)=U(j,k:m)-L(j,k)*U(k,k:m);
end
end
end
\ No newline at end of file
%% Part 1
U = triu(rand(5,5),0);
x_exact = rand(5,1);
b = U*x_exact;
y = usol(U,b);
x_exact-y
L = tril(rand(5,5),0);
x_exact1 = rand(5,1);
b1 = L*x_exact1;
y1 = lsol(L,b1);
x_exact1-y1
%% Part 2
A=[4 ,-3, 1;0.5, 3,-1;0,1,0];
[L,U] = LUNoPivot(A);
[L1,U1,P]=lu(A);
b2=rand(3,1);
x_exact2=A\b2;
y2=lsol(L,b2);
x2=usol(U,y2);
x2-x_exact2
%% Part 3
n=3;
A=rand(n,n);
[L,U,p] = LUPivot(A);
P=zeros(n,n);
for i=1:n
P(i,p(i))=1;
end
[L1,U1,P1]=lu(A);
b3=rand(n,1);
x_exact3=P*A\b3;
y3=lsol(L,b3);
x3=usol(U,y3);
x3-x_exact3
%% Part 4
% a)
A = [1 5 2; 0 0 8; 2 4 1];
b = [9;8;9];
[L,U,p] = LUPivot(A);
P=zeros(3,3);
for i=1:3
P(i,p(i))=1;
end
x_exact3=P*A\b;
y3=lsol(L,b);
x3=usol(U,y3);
x3-x_exact3
%b)
N=12;
for i=1:N
for j=1:N
A(i,j)=((j/N)^(i-1))^N;
end
b(i)=1/i;
end
[L,U,p] = LUPivot(A);
[L1,U1,P1]=lu(A);
P=zeros(N,N);
for i=1:N
P(i,p(i))=1;
end
x_exact4=P*A\b;
y4=lsol(L,b);
x4=usol(U,y4);
x4-x_exact4
%% Part 5
n=10;
A=rand(n,n);
b=randn(n,1);
tic
x_exact5=A\b;
toc
tic
[L,U,p] = LUPivot(A);
b3=b;
for i=1:n
b5(i)=b(p(i));
end
y5=lsol(L,b5);
x5=usol(U,y5);
toc
x5-x_exact5
%% Part 6
n=100;
h=1/n;
A=(1/h)*(gallery('tridiag',n,1,4,1));
[L,U,p] = LUPivot(A);
[L1,U1] = LUNoPivot(A);
%check if p=(1,2,......,n)
clear all;
%% PArt 1
n=9;
A=rand(n);
det(A);
d = mydet(A);
%% Part 2
n=100;
A=randn(n,n);
[L,U]=LUNoPivot(A);
prod(diag(L))*prod(diag(U))
det(A)
clear all;
n = 2;
format long
kappa = [1.0E2, 1.0E6, 1.0E10, 1.0E18];
for i = 1:length(kappa)
% Data
kappa(i)
A = matgen(n, kappa(i));
x = randn(n,1);
b = A*x;
% Cramer's rule
d=A(1,1)*A(2,2)-A(2,1)*A(1,2);
x_cramer = [(b(1)*A(2,2)-b(2)*A(1,2))/d ; (b(2)*A(1,1)-b(1)*A(2,1))/d];
% Gauss
[L,U,p]=LUPivot(A);
P=zeros(n,n);
for i=1:n
P(i,p(i))=1;
end
y=lsol(L,b);
x_Gauss=usol(U,y);
% Backslash
x_back = A\b;
R_cramer = [ norm(x_cramer-x)/norm(x), norm(A*x_cramer-b)/(norm(A)*norm(x)+norm(b))]
R_Gauss = [ norm(x_Gauss-x)/norm(x), norm(A*x_Gauss-b)/(norm(A)*norm(x)+norm(b))]
R_back = [ norm(x_back-x)/norm(x), norm(A*x_back-b)/(norm(A)*norm(x)+norm(b))]
end
clear all;
n = 100;
kappa = [1.0E2, 1.0E6, 1.0E10, 1.0E18];
for i = 1:length(kappa)
% Data
A = matgen(n, kappa(i));
x = randn(n,1);
b = A*x;
% Invert and multiply
M = inv(A);
x1 = M*b;
% LU fac
x2 = A\b;
R = [ norm(x1-x)/norm(x); norm(A*x1-b)/(norm(A)*norm(x1)+norm(b));