...
 
Commits (17)
\documentclass{article}
\usepackage{/home/mccoid/LaTeX/preamble}
\usepackage{C:/Users/conmc/Documents/LaTeX/preamble}
\begin{document}
......
function M=InterfaceMatrixAF(Na,Ta,Nb,Tb)
% INTERFACEMATRIX projection matrix for nonmatching triangular grids
% M=InterfaceMatrix(Na,Ta,Nb,Tb); takes two triangular meshes Ta
% and Tb with associated nodal coordinates in Na and Nb and
% computes the interface projection matrix M
na=size(Ta,1); nb=size(Tb,1); % number of triangles in each list
cand=zeros(na,nb); cand(1)=1; % init. candidate matrix
bl=1; % list of triangles in Tb to treat
bd=zeros(nb+1,1); bd([1,end])=1; % flag triangles already treated
M=sparse(size(Nb,2),size(Na,2));
while ~isempty(bl)
bc=bl(1); bl=bl(2:end); % pick out 1st tri. in bl
al=find(cand(:,bc)==1)'; % list of candidates intersecting bc
ad=zeros(na+1,1);
ad(end)=1; ad(al)=1; % flag triangles of Ta already treated
while ~isempty(al)
ac=al(1); al=al(2:end);
[P,nc,Mc]=Intersect(Nb(:,Tb(bc,1:3)),Na(:,Ta(ac,1:3)));
if ~isempty(P)
M(Tb(bc,1:3),Ta(ac,1:3))=M(Tb(bc,1:3),Ta(ac,1:3))+Mc;
ta=Ta(ac,4:6); % nhbrs of ac
al=[al ta(ad(ta)==0)]; % add those nhbrs not treated to al
ad(ta)=1; % flag new entries of al
cand(ta,bc)=1; % indicate candidate (superfluous?)
tb=Tb(bc,4:6); % nhbrs of bc
bl=[bl tb(bd(tb)==0)]; % add those nhbrs not treated to bl
bd(tb)=1; % flag new entries of bl
cand(ac,tb)=1; % ac is candidate for new bl entries
end
end
end
end
function [P,n,M]=Intersect(X,Y)
% INTERSECT intersection of two triangles and mortar contribution
% [P,n,M]=Intersect(X,Y); computes for the two given triangles X and
% Y (point coordinates are stored column-wise, in counter clock
% order) the points P where they intersect, in n the indices of
% which neighbors of X are also intersecting with Y, and the local
% mortar matrix M of contributions of the P1 elements X on the P1
% element Y. The numerical challenges are handled by including
% points on the boundary and removing duplicates at the end.
[P,n]=EdgeIntersections(X,Y);
P1=PointsOfXInY(X,Y);
if size(P1,2)>1 % if two or more interior points
n=[1 1 1]; % the triangle is candidate for all
end % neighbors
P=[P P1];
P=[P PointsOfXInY(Y,X)];
P=SortAndRemoveDoubles(P); % sort counter clock wise
M=zeros(3,3);
if size(P,2)>0
for j=2:size(P,2)-1 % compute interface matrix
M=M+MortarInt(P(:,[1 j j+1]),X,Y);
end
patch(P(1,:),P(2,:),'m') % draw intersection for illustration
% H=line([P(1,:) P(1,1)],[P(2,:),P(2,1)]);
% set(H,'LineWidth',3,'Color','m');
pause(0)
end
end
function P=PointsOfXInY(X,Y)
% POINTSOFXINY finds corners of one triangle within another one
% P=PointsOfXInY(X,Y); computes for the two given triangles X
% and Y (point coordinates are stored column-wise, in counter clock
% order) the corners P of X which lie in the interior of Y.
k=0;P=[];
v0=Y(:,2)-Y(:,1); v1=Y(:,3)-Y(:,1); % find interior points of X in Y
d00=v0'*v0; d01=v0'*v1; d11=v1'*v1; % using baricentric coordinates
id=1/(d00*d11-d01*d01);
for i=1:3
v2=X(:,i)-Y(:,1); d02=v0'*v2; d12=v1'*v2;
u=(d11*d02-d01*d12)*id; v=(d00*d12-d01*d02)*id;
if u>=0 && v>=0 && u+v<=1 % also include nodes on the boundary
k=k+1; P(:,k)=X(:,i);
end
end
end
function [P,n]=EdgeIntersections(X,Y)
% EDGEINTERSECTIONS computes edge intersections of two triangles
% [P,n]=EdgeIntersections(X,Y) computes for the two given triangles X
% and Y (point coordinates are stored column-wise, in counter clock
% order) the points P where their edges intersect. In addition,
% in n the indices of which neighbors of X are also intersecting
% with Y are given.
P=[]; k=0; n=[0 0 0];
for i=1:3 % find all intersections of edges
for j=1:3
b=Y(:,j)-X(:,i);
A=[X(:,mod(i,3)+1)-X(:,i) -Y(:,mod(j,3)+1)+Y(:,j)];
if rank(A)==2 % edges not parallel
r=A\b;
if r(1)>=0 && r(1)<=1 && r(2)>=0 && r(2)<=1 % intersection found
k=k+1; P(:,k)=X(:,i)+r(1)*(X(:,mod(i,3)+1)-X(:,i)); n(i)=1;
end
end
end
end
end
function P=SortAndRemoveDoubles(P)
% SORTANDREMOVEDOUBLES sort points and remove duplicates
% P=SortAndRemoveDoubles(P); orders polygon corners in P counter
% clock wise and removes duplicates
ep=10*eps; % tolerance for identical nodes
m=size(P,2);
if m>0
c=sum(P')'/m; % order polygon corners counter
for i=1:m % clockwise
d=P(:,i)-c; ao(i)=angle(d(1)+sqrt(-1)*d(2));
end
[tmp,id]=sort(ao);
P=P(:,id);
i=1;j=2; % remove duplicates
while j<=m
if norm(P(:,i)-P(:,j))>ep
i=i+1;P(:,i)=P(:,j);j=j+1;
else
j=j+1;
end
end
P=P(:,1:i);
end
end
function M=MortarInt(T,X,Y)
% MORTARINT computes mortar contributions
% M=MortarInt(T,X,Y); computes for triangles X and Y with nodal coordinates
% stored columnwise the integrals of all P1 element shape function combinations
% between triangles X and Y on triangle T.
% The result is stored in the 3 by 3 matrix M
Jd=-T(1,1)*T(2,3)-T(1,2)*T(2,1)+T(1,2)*T(2,3)+T(1,1)*T(2,2)+...
T(1,3)*T(2,1)-T(1,3)*T(2,2);
a=T(1,1);
d=T(2,1);
e=-T(2,1)+T(2,2);
b=-T(1,1)+T(1,2);
c=-T(1,1)+T(1,3);
f=-T(2,1)+T(2,3);
T1=1/2*(X(1,1)*(X(2,2)-X(2,3))-X(1,2)*(X(2,1)-X(2,3))+X(1,3)*(X(2,1)-X(2,2)));
T2=1/2*(Y(1,1)*(Y(2,2)-Y(2,3))-Y(1,2)*(Y(2,1)-Y(2,3))+Y(1,3)*(Y(2,1)-Y(2,2)));
A11=1/(2*abs(T1))*((X(2,2)-X(2,3))*b-(X(1,2)-X(1,3))*e);
B11=1/(2*abs(T1))*((X(2,2)-X(2,3))*c-(X(1,2)-X(1,3))*f);
C11=1/(2*abs(T1))*((X(2,2)-X(2,3))*a-(X(1,2)-X(1,3))*d+(X(1,2)*X(2,3)-X(2,2)*X(1,3)));
A12=1/(2*abs(T1))*((X(2,3)-X(2,1))*b+(X(1,1)-X(1,3))*e);
B12=1/(2*abs(T1))*((X(2,3)-X(2,1))*c+(X(1,1)-X(1,3))*f);
C12=1/(2*abs(T1))*((X(2,3)-X(2,1))*a+(X(1,1)-X(1,3))*d-(X(1,1)*X(2,3)-X(2,1)*X(1,3)));
A13=1/(2*abs(T1))*((X(2,1)-X(2,2))*b-(X(1,1)-X(1,2))*e);
B13=1/(2*abs(T1))*((X(2,1)-X(2,2))*c-(X(1,1)-X(1,2))*f);
C13=1/(2*abs(T1))*((X(2,1)-X(2,2))*a-(X(1,1)-X(1,2))*d+(X(1,1)*X(2,2)-X(2,1)*X(1,2)));
A21=1/(2*abs(T2))*((Y(2,2)-Y(2,3))*b-(Y(1,2)-Y(1,3))*e);
B21=1/(2*abs(T2))*((Y(2,2)-Y(2,3))*c-(Y(1,2)-Y(1,3))*f);
C21=1/(2*abs(T2))*((Y(2,2)-Y(2,3))*a-(Y(1,2)-Y(1,3))*d+(Y(1,2)*Y(2,3)-Y(2,2)*Y(1,3)));
A22=1/(2*abs(T2))*((Y(2,3)-Y(2,1))*b+(Y(1,1)-Y(1,3))*e);
B22=1/(2*abs(T2))*((Y(2,3)-Y(2,1))*c+(Y(1,1)-Y(1,3))*f);
C22=1/(2*abs(T2))*((Y(2,3)-Y(2,1))*a+(Y(1,1)-Y(1,3))*d-(Y(1,1)*Y(2,3)-Y(2,1)*Y(1,3)));
A23=1/(2*abs(T2))*((Y(2,1)-Y(2,2))*b-(Y(1,1)-Y(1,2))*e);
B23=1/(2*abs(T2))*((Y(2,1)-Y(2,2))*c-(Y(1,1)-Y(1,2))*f);
C23=1/(2*abs(T2))*((Y(2,1)-Y(2,2))*a-(Y(1,1)-Y(1,2))*d+(Y(1,1)*Y(2,2)-Y(2,1)*Y(1,2)));
M(1,1)=1/24*Jd*(2*A11*A21+B11*A21+A11*B21+2*B11*B21+4*C11*A21+4*A11*C21+4*C11*B21+4*B11*C21+12*C11*C21);
M(1,2)=1/24*Jd*(2*A11*A22+B11*A22+A11*B22+2*B11*B22+4*C11*A22+4*A11*C22+4*C11*B22+4*B11*C22+12*C11*C22);
M(1,3)=1/24*Jd*(2*A11*A23+B11*A23+A11*B23+2*B11*B23+4*C11*A23+4*A11*C23+4*C11*B23+4*B11*C23+12*C11*C23);
M(2,1)=1/24*Jd*(2*A12*A21+B12*A21+A12*B21+2*B12*B21+4*C12*A21+4*A12*C21+4*C12*B21+4*B12*C21+12*C12*C21);
M(2,2)=1/24*Jd*(2*A12*A22+B12*A22+A12*B22+2*B12*B22+4*C12*A22+4*A12*C22+4*C12*B22+4*B12*C22+12*C12*C22);
M(2,3)=1/24*Jd*(2*A12*A23+B12*A23+A12*B23+2*B12*B23+4*C12*A23+4*A12*C23+4*C12*B23+4*B12*C23+12*C12*C23);
M(3,1)=1/24*Jd*(2*A13*A21+B13*A21+A13*B21+2*B13*B21+4*C13*A21+4*A13*C21+4*C13*B21+4*B13*C21+12*C13*C21);
M(3,2)=1/24*Jd*(2*A13*A22+B13*A22+A13*B22+2*B13*B22+4*C13*A22+4*A13*C22+4*C13*B22+4*B13*C22+12*C13*C22);
M(3,3)=1/24*Jd*(2*A13*A23+B13*A23+A13*B23+2*B13*B23+4*C13*A23+4*A13*C23+4*C13*B23+4*B13*C23+12*C13*C23);
end
\ No newline at end of file
......@@ -4,39 +4,70 @@ function M=InterfaceMatrixRF(Na,Ta,Nb,Tb)
% and Tb with associated nodal coordinates in Na and Nb and
% computes the interface projection matrix M
bl =1; % bl: list of triangles of Tb to treat
bil=1; % bil: list of triangles Ta to start with
bd=zeros(size(Tb,1)+1,1); % bd: flag for triangles in Tb treated
bd(end)=1; % guard, to treat boundaries
bd(1)=1; % mark first triangle in b list.
na=size(Ta,1); nb=size(Tb,1); % number of triangles in each list
cand=zeros(na,nb); cand(1)=1; % init. candidate matrix
bl=1; % list of triangles in Tb to treat
bd=zeros(nb+1,1); bd([1,end])=1; % flag triangles already treated
M=sparse(size(Nb,2),size(Na,2));
while ~isempty(bl)
bc=bl(1); bl=bl(2:end); % bc: current triangle of Tb
al=bil(1); bil=bil(2:end); % triangle of Ta to start with
ad=zeros(size(Ta,1)+1,1); % same as for bd
ad(end)=1;
ad(al)=1;
n=[0 0 0]; % triangles intersecting with neighbors
bc=bl(1); bl=bl(2:end); % pick out 1st tri. in bl
al=find(cand(:,bc)==1)'; % list of candidates intersecting bc
ad=zeros(na+1,1);
ad(end)=1; ad(al)=1; % flag triangles of Ta already treated
while ~isempty(al)
ac=al(1); al=al(2:end); % take next candidate
ac=al(1); al=al(2:end);
[P,nc,Mc]=Intersect(Nb(:,Tb(bc,1:3)),Na(:,Ta(ac,1:3)));
if ~isempty(P) % intersection found
if ~isempty(P)
M(Tb(bc,1:3),Ta(ac,1:3))=M(Tb(bc,1:3),Ta(ac,1:3))+Mc;
t=Ta(ac,3+find(ad(Ta(ac,4:6))==0));
al=[al t]; % add neighbors
ad(t)=1;
n(nc>0)=ac; % ac is starting candidate for neighbor
ta=Ta(ac,4:6); % nhbrs of ac
al=[al ta(ad(ta)==0)]; % add those nhbrs not treated to al
ad(ta)=1; % flag new entries of al
cand(ta,bc)=1; % indicate candidate (superfluous?)
tb=Tb(bc,4:6); % nhbrs of bc
bl=[bl tb(bd(tb)==0)]; % add those nhbrs not treated to bl
bd(tb)=1; % flag new entries of bl
cand(ac,tb)=1; % ac is candidate for new bl entries
end
end
tmp=find(bd(Tb(bc,4:6))==0); % find non-treated neighbors
idx=find(n(tmp)>0); % take those which intersect
t=Tb(bc,3+tmp(idx));
bl=[bl t]; % and add them
bil=[bil n(tmp(idx))]; % with starting candidates Ta
bd(t)=1;
end
end
% bl =1; % bl: list of triangles of Tb to treat
% bil=1; % bil: list of triangles Ta to start with
% bd=zeros(size(Tb,1)+1,1); % bd: flag for triangles in Tb treated
% bd(end)=1; % guard, to treat boundaries
% bd(1)=1; % mark first triangle in b list.
% M=sparse(size(Nb,2),size(Na,2));
% while ~isempty(bl)
% bc=bl(1); bl=bl(2:end); % bc: current triangle of Tb
% al=bil(1); bil=bil(2:end); % triangle of Ta to start with
% ad=zeros(size(Ta,1)+1,1); % same as for bd
% ad(end)=1;
% ad(al)=1;
% n=[0 0 0]; % triangles intersecting with neighbors
% while ~isempty(al)
% ac=al(1); al=al(2:end); % take next candidate
% [P,nc,Mc]=Intersect(Nb(:,Tb(bc,1:3)),Na(:,Ta(ac,1:3)));
% if ~isempty(P) % intersection found
% M(Tb(bc,1:3),Ta(ac,1:3))=M(Tb(bc,1:3),Ta(ac,1:3))+Mc;
% t=Ta(ac,3+find(ad(Ta(ac,4:6))==0));
% al=[al t]; % add neighbors
% ad(t)=1;
% n(nc>0)=ac; % ac is starting candidate for neighbor
% end
% end
% tmp=find(bd(Tb(bc,4:6))==0); % find non-treated neighbors
% idx=find(n(tmp)>0); % take those which intersect
% t=Tb(bc,3+tmp(idx));
% bl=[bl t]; % and add them
% bil=[bil n(tmp(idx))]; % with starting candidates Ta
% bd(t)=1;
% end
% end
% try same without advancing front - is the problem the intersection or the
% advancing front?
function [P,n,M]=Intersect(X,Y)
% INTERSECT intersection of two triangles and mortar contribution
% [P,n,M]=Intersect(X,Y); computes for the two given triangles X and
......
% U = [-1,2,-0.5;2,2,-1];
U = [0.5,-1,5;0.5,5,-1];
V = [0,3,0;0,0,3];
triU = [U,U(:,1)];
triV = [V,V(:,1)];
plot(triU(1,:),triU(2,:),'-',triV(1,:),triV(2,:),'-')
[P,Q,R,n] = RefFree(U,V);
P = SortAndRemoveDoubles([P,Q,R]);
patch(P(1,:),P(2,:),'m')
\ No newline at end of file
function P=SortAndRemoveDoubles(P)
% SORTANDREMOVEDOUBLES sort points and remove duplicates
% P=SortAndRemoveDoubles(P); orders polygon corners in P counter
% clock wise and removes duplicates
ep=10*eps; % tolerance for identical nodes
m=size(P,2);
if m>0
c=sum(P')'/m; % order polygon corners counter
for i=1:m % clockwise
d=P(:,i)-c; ao(i)=angle(d(1)+sqrt(-1)*d(2));
end
[tmp,id]=sort(ao);
P=P(:,id);
i=1;j=2; % remove duplicates
while j<=m
if norm(P(:,i)-P(:,j))>ep
i=i+1;P(:,i)=P(:,j);j=j+1;
else
j=j+1;
end
end
P=P(:,1:i);
end
end
\ No newline at end of file
......@@ -18,7 +18,8 @@ set(fig,'DoubleBuffer','on');
PlotMesh(N1,T1,'b');
PlotMesh(N2,T2,'r');
M = InterfaceMatrixNew(N1,T1,N2,T2);
% M = InterfaceMatrixNew(N1,T1,N2,T2);
M = InterfaceMatrixRF(N1,T1,N2,T2);
%%
tic
......
[N,T,Np]=CircleMesh(10,eps);
[N,T,Np]=CircleMesh(20,eps);
%%
figure(1)
clf
PlotMesh(N,T,'b')
PlotMesh(Np,T,'r')
M = InterfaceMatrix(N,T,Np,T);
% M = InterfaceMatrix(N,T,Np,T);
M = InterfaceMatrixAF(N,T,Np,T);
% M = QuadComplex(N,T,Np,T);
pause(1)
% savefig('CirclePANG.fig')
......@@ -22,7 +23,5 @@ pause(1)
clf
PlotMesh(N,T,'b')
PlotMesh(Np,T,'r')
hold on
M = InterfaceMatrixRF(N,T,Np,T);
hold off
% savefig('CircleNew.fig')
\ No newline at end of file
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 8 mars 2017 \\[-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 les normes $\|\cdot\|_p$ avec $1\leq p\leq\infty$,
on a
\[
\rho(A)=\lim_{k\to\infty}\|A^k\|_p^{\frac{1}{k}}.
\]
De plus, pour tout $\epsilon>0$, il existe une norme de $\C^n$ telle que
\[
\rho(A)\leq \|A\| \leq \rho(A)+\epsilon.
\]
%\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 3 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{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}
\usepackage{float}
%% 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 5 avril \\[-3mm] \hrule
\bigskip
\begin{center}
\textbf{\Large{Méthodes Itératives}} \\[1mm]\bigskip
\textbf{Série 3}
\end{center}
\bigskip
\begin{enumerate}
\item (3 pts) On parle de matrice à diagonale strictement dominante $A$ si ses
coefficients diagonaux satisfont
\[
|a_{i,i}|>\sum_{j\neq i}|a_{i,j}|,\quad\forall i \in \{1,\ldots,n\}.
\]
Démontrer, pour une telle matrice, que la méthode de Jacobi converge.
\item (6 pts) Soit $A$ la matrice du moin Laplacien ($-\Delta$) discret en dimension 2:
\begin{equation*}
A = -
\begin{bmatrix}
T & I\\
I & T & \ddots\\
&\ddots & \ddots & I\\
& & I & T \\
\end{bmatrix}
\quad \text{avec} \quad
T =
\begin{bmatrix}
-4 & 1\\
1 & -4 & 1 \\
& & \ddots & \\
& & 1 & -4\\
\end{bmatrix}.
\end{equation*}
Considérez le splitting de bloc-Jacobi $A = D + L + L^\top = M - N$ avec
\begin{equation*}
M = D = -
\begin{bmatrix}
T \\
& T \\
& & \ddots\\
& & & T \\
\end{bmatrix}
\quad \text{et} \quad
N = - \bigl( L + L^\top \bigr) =
\begin{bmatrix}
0 & I \\
I & 0 & \ddots\\
&\ddots & \ddots & I\\
& & I & 0 \\
\end{bmatrix}.
\end{equation*}
En utilisant les résultats obtenus en cours, démontrer que la méthode de bloc-Jacobi converge.\\
\textit{Indication:} Vous devez démontrer que $-T$, $A$ et $2D - A$ sont définies positives.
Rappelez aussi la Série 2.
\item (3 pts) Pour la matrice du moin Laplacien discret en dimension 2, calculer
le paramètre optimal pour la méthode SOR. Faire de même pour la méthode de
Richardson.
\item (2 pts) Si une matrice $A$ a la propriété $A$, montrer que la
méthode de Gauss-Seidel converge deux fois plus vite que Jacobi.\\
\textit{Indication:} Utiliser le théorème
$(\lambda+\omega-1)^2=\lambda\omega^2\mu^2$, où $\lambda$ est une valeur
propre de la matrice d'itération de SOR, $\mu$ une valeur propre de la
matrice d'itération de Jacobi et $\omega$ le paramètre de relaxation dans
SOR.
\item (3 pts) Rappelons que la matrice d'itération de la méthode SOR est
\begin{equation*}
G_{\rm SOR} = \bigl( D + \omega L \bigr)^{-1} \Bigl( - \omega U + (1-\omega) D \Bigr).
\end{equation*}
Si $\lambda=0$ est une valeur propre de $G_{SOR}$, quelle est la valeur de $\omega$?
\item (6 pts) Implémenter une fonction Matlab qui résout un système linéaire
$A \vec{u} = \vec{f}$ à l'aide de la méthode {\it SOR}.
Testez votre code en résolvant l'équation suivante
\begin{equation*}
\begin{split}
-\Delta u &= 1 \, \text{ dans $\Omega = (0,1)\times(0,1)$}, \\
u &= 0 \; \text{ sur $\partial \Omega$}. \\
\end{split}
\end{equation*}
De plus, à l'aide des méthodes de {\it Jacobi}
et {\it Gauss-Seidel} (voir Série 2) essayez d'obtenir la figure suivante.
\begin{figure}[H]
\centering
\includegraphics[scale=0.6]{convergence.pdf}
\caption{La ligne bleue correspond à la méthode de Jacobi,
la ligne rouge correspond à la méthode de Gauss-Seidel
et la ligne violette correspond à la méthode SOR.
Les lignes noires pointillées représentent les estimations de la convergence des trois méthodes.}
\end{figure}
\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{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}
\usepackage{float}
%% 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 3 mai \\[-3mm] \hrule
\bigskip
\begin{center}
\textbf{\Large{Méthodes Itératives}} \\[1mm]\bigskip
\textbf{Série 4}
\end{center}
\bigskip
\textbf{Partie I: L'inégalité de Kantorovitch}\\
Dans cette partie nous vous demandons de démontrer l'inégalité de Kantorovitch au moyen de
l'inégalité de Wielandt.
Considérez une matrice hermitienne et définie positive $A \in \C^{n \times n}$ avec valeurs propres
$0<\lambda_1\leq\lambda_2\leq\cdots\leq\lambda_n$.
L'inégalité de Wielandt est
\begin{equation}\label{eq:Wielandt}
| \vec{u}^H A \vec{v} |^2 \leq \left( \frac{ \lambda_1 - \lambda_n }{ \lambda_1 + \lambda_n } \right)^2
(\vec{u}^H A \vec{u}) (\vec{v}^H A \vec{v}) \, \text{ pour tous $\vec{u},\vec{v} \in\C^n$ orthogonaux},
\end{equation}
et l'inégalité de Kantorovitch est
\begin{equation}\label{eq:Kantorovitch}
(\vec{u}^H A \vec{u}) (\vec{u}^H A^{-1} \vec{u}) \leq \frac{ ( \lambda_1 + \lambda_n )^2 }{ 4 \lambda_1 \lambda_n }
\| \vec{u} \|_2^4 \, \text{ pour tous $\vec{u} \in\C^n$}.
\end{equation}
L'inégalité de Kantorovitch peut être obtenue au moyen des étapes suivantes:
\begin{enumerate}\itemsep0em
\item[(a)] (3 pts) Démontrer que \eqref{eq:Wielandt} et \eqref{eq:Kantorovitch}
sont tous les deux satisfaits si $\vec{u} \in \C^n$ est un vecteur propre de $A$.
\textit{Indication :} utiliser l'inégalité arithmético-géométrique $\sqrt{ab} \leq \frac{1}{2} (a+b)$
pour tout $a,b>0$.
\item[(b)] (3 pts) Démontrer que si $\vec{u} \in \C^n$ n'est pas un vecteur propre de $A$,
alors $A^{-1} \vec{u} - (\vec{u}^H A^{-1} \vec{u}) \vec{u} \neq 0$
et $\vec{u} - (\vec{u}^H A^{-1} \vec{u}) A \vec{u} \neq 0$.
\item[(c)] (3 pts) Démontrer que si $\vec{u} \in \C^n$ est un vecteur unitaire, c.-à-d. $\|\vec{u}\|_2=1$,
alors $(\vec{u}^H A \vec{u}) (\vec{u}^H A^{-1} \vec{u}) \geq 1$, avec une inégalité stricte si $\vec{u}$
n'est pas un vecteur propre de $A$.\\
\textit{Indication:} commencer avec $1 = (\vec{u}^H \vec{u})^2 = (\vec{u}^H A^{1/2} A^{-1/2} \vec{u})^2$
et utiliser l'inégalité de Cauchy-Schwarz.
\item[(d)] (3 pts) Considérer un vecteur unitaire $\vec{u} \in \C^n$ qui n'est pas un vecteur propre de $A$.
Définir $\vec{v} := A^{-1} \vec{u} - (\vec{u}^H A^{-1} \vec{u}) \vec{u}$ et montrer que :
\begin{itemize}\itemsep0em
\item[(1)] $\vec{v}^H \vec{u} = 0$.
\item[(2)] $A\vec{v} \neq 0$.
\item[(3)] $\vec{u}^H A \vec{v} = 1 - (\vec{u}^H A \vec{u}) (\vec{u}^H A^{-1} \vec{u}) < 0$.
\item[(4)] $\vec{v}^H A \vec{v} = - (\vec{u}^H A^{-1} \vec{u}) (\vec{v}^H A \vec{u})$.
\end{itemize}
\textit{Indication :} utiliser les résultats démontrés dans (b) et (c).
\item[(e)] (4 pts) Utiliser l'inégalité de Wielandt \eqref{eq:Wielandt} et les résultats démontrés
ci-dessus pour obtenir l'inégalité de Kantorovitch \eqref{eq:Kantorovitch}.
\end{enumerate}
$\;$\\
\textbf{Partie II: Résumé des méthodes itératives stationnaires}
\begin{enumerate}
\item (2 pts) Démontrer les équivalences suivantes :
\begin{itemize}\itemsep0em
\item Jacobi:
\begin{equation*}
\begin{split}
&(\vec{u}_{k+1})_i = \frac{1}{a_{ii}}\left[\vec{f}_i - \sum_{j=1}^{i-1} a_{ij} (\vec{u}_{k})_j - \sum_{j=i+1}^{n} a_{ij} (\vec{u}_{k})_j \right] \\
&\Leftrightarrow \;
\vec{u}_{k+1} = - D^{-1} (L+U) \vec{u}_{k} + D^{-1} \vec{f}.
\end{split}
\end{equation*}
\item Gauss-Seidel:
\begin{equation*}
\begin{split}
&(\vec{u}_{k+1})_i = \frac{1}{a_{ii}}\left[\vec{f}_i - \sum_{j=1}^{i-1} a_{ij} (\vec{u}_{k+1})_j - \sum_{j=i+1}^{n} a_{ij} (\vec{u}_{k})_j \right] \\
&\Leftrightarrow \;
\vec{u}_{k+1} = - (D+L)^{-1} U \vec{u}_{k} + (D+L)^{-1} \vec{f}.
\end{split}
\end{equation*}
\item SOR:
\begin{equation*}
\begin{split}
&(\vec{u}_{k+1})_i = (1-\omega)(\vec{u}_{k})_i + \frac{\omega}{a_{ii}}\left[\vec{f}_i - \sum_{j=1}^{i-1} a_{ij} (\vec{u}_{k+1})_j - \sum_{j=i+1}^{n} a_{ij} (\vec{u}_{k})_j \right] \\
&\Leftrightarrow \;
\vec{u}_{k+1} = - (D+\omega L)^{-1} \left[ \bigl( (1-\omega)D - \omega U \bigr) \vec{u}_{k} + \omega \vec{f} \right].
\end{split}
\end{equation*}
\end{itemize}
\item (3 pts) Considérer le Laplacien discret $A$ de $-\Delta$ en dimension 2.
Démontrer que la méthode de Jacobi et la méthode de Richardson (avec le paramètre optimal $\alpha_{opt}$
obtenu dans la Série 3) ont le même comportement de convergence.
\item (3 pts) Considérer une matrice $A \in \R^{2 \times 2}$ et le splitting $A = D + L + U$ avec $D$ inversible.
Démontrer que si la méthode de Jacobi converge, alors la méthode de Gauss-Seidel converge aussi et plus rapidement.
\item (4 pts) Considérer le Laplacien discret $A\in \R^{n\times n}$ de $-\Delta$ en dimension 2
et soit $h=\frac{1}{n+1}$ la taille de la grille.
Dénoter avec $\rho(G_{\rm J};h)$, $\rho(G_{\rm GS};h)$, et $\rho(G_{\rm SOR};h)$ les rayons spectraux
des méthodes de Jacobi, Gauss-Seidel, et SOR en fonction de $h$.
Démontrer que :
\begin{itemize}\itemsep0em
\item $\rho(G_{\rm J};h) = 1 - \frac{\pi^2}{2} h^2 + O(h^4)$.
\item $\rho(G_{\rm GS};h) = 1 - \pi^2 h^2 + O(h^4)$.
\item $\rho(G_{\rm SOR};h) = 1 - 2 \pi h + O(h^2)$.
\end{itemize}
\end{enumerate}
$\;$\\
\\
\textbf{Partie III: Méthodes de descente maximale (ou du gradient) et du gradient conjugué}
\begin{enumerate}\itemsep0em
\item (4 pts) Considérer une fonction $f : \R^n \rightarrow \R$ et on dénote avec $\nabla f(\vec{v})$ le
gradient de $f$ en $\vec{v} \in \R^n$ obtenu au moyen du produit scalaire habituel pour $\R^n$.
Démontrer que $\nabla f(\vec{v})$ est orthogonal à la ligne de niveau de $f$ en $\vec{v}$ et que
$-\nabla f(\vec{v})$ est la direction de descente maximale de $f$ en $\vec{v}$.
\item (4 pts) Implémenter une fonction Matlab qui résout un système linéaire $A \vec{u} = \vec{f}$ à l'aide de la méthode
du gradient.
\item (4 pts) Implémenter une fonction Matlab qui résout un système linéaire $A \vec{u} = \vec{f}$ à l'aide de la méthode
du gradient conjugué.
\item (1 pts) Tester vos codes en résolvant l'équation de Laplace (comme dans la Série 3).
\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{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}
\usepackage{float}
%% 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 3 mai \\[-3mm] \hrule
\bigskip
\begin{center}
\textbf{\Large{Méthodes Itératives}} \\[1mm]\bigskip
\textbf{Série 4}
\end{center}
\bigskip
\textbf{Part I: The Kantorovitch inequality}\\
In this fist part you must prove the Kantorovitch inequality using the Wielandt inequality.
Consider a Hermitian positive definite matrix $A \in \C^{n \times n}$ having eigenvalues
$0<\lambda_1\leq\lambda_2\leq\cdots\leq\lambda_n$. The Wielandt inequality is
\begin{equation}\label{eq:Wielandt}
| \vec{u}^H A \vec{v} |^2 \leq \left( \frac{ \lambda_1 - \lambda_n }{ \lambda_1 + \lambda_n } \right)^2
(\vec{u}^H A \vec{u}) (\vec{v}^H A \vec{v}) \, \text{ for all orthogonal $\vec{u},\vec{v} \in\C^n$},
\end{equation}
and the Kantorovitch inequality is
\begin{equation}\label{eq:Kantorovitch}
(\vec{u}^H A \vec{u}) (\vec{u}^H A^{-1} \vec{u}) \leq \frac{ ( \lambda_1 + \lambda_n )^2 }{ 4 \lambda_1 \lambda_n }
\| \vec{u} \|_2^4 \, \text{ for all $\vec{u} \in\C^n$}.
\end{equation}
The Kantorovitch inequality can be obtained by the Wielandt inequality following the next steps.
\begin{enumerate}
\item[(a)] (3 pts) Show that \eqref{eq:Wielandt} and \eqref{eq:Kantorovitch} are both satisfied if $\vec{u} \in \C^n$
is an eigenvector of $A$.
Hint: use the arithmetic-geometric mean inequality $\sqrt{ab} \leq \frac{1}{2} (a+b)$
for any $a,b>0$.
\item[(b)] (3 pts) Show that if $\vec{u} \in \C^n$ is not an eigenvector of $A$,
then $A^{-1} x - (\vec{u}^H A^{-1} \vec{u}) \vec{u} \neq 0$
and $x - (\vec{u}^H A^{-1} \vec{u}) A \vec{u} \neq 0$.
\item[(c)] (3 pts) Show that if $\vec{u} \in \C^n$ is a unit vector, that is $\|\vec{u}\|_2=1$,
then $(\vec{u}^H A \vec{u}) (\vec{u}^H A^{-1} \vec{u}) \geq 1$, with strict inequality if $\vec{u}$ is
not an eigenvector of $A$.
Hint: begin with $1 = (\vec{u}^H \vec{u})^2 = (\vec{u}^H A^{1/2} A^{-1/2} \vec{u})^2$
and use the Cauchy-Schwarz inequality.
\item[(d)] (3 pts) Consider a unit vector $\vec{u} \in \C^n$ which is not an eigenvector of $A$.
Define $\vec{v} := A^{-1} \vec{u} - (\vec{u}^H A^{-1} \vec{u}) \vec{u}$.
Show that
\begin{itemize}\itemsep0em
\item[(1)] $\vec{v}^H \vec{u} = 0$.
\item[(2)] $A\vec{v} \neq 0$.
\item[(3)] $\vec{u}^H A \vec{v} = 1 - (\vec{u}^H A \vec{u}) (\vec{u}^H A^{-1} \vec{u}) < 0$.
\item[(4)] $\vec{v}^H A \vec{v} = - (\vec{u}^H A^{-1} \vec{u}) (\vec{v}^H A \vec{u})$.
\end{itemize}
Hint: use the results proved in (b) and (c).
\item[(e)] (4 pts) Use the Wielandt inequality \eqref{eq:Wielandt} and the results proved above
to obtain the Kantorovitch inequality \eqref{eq:Kantorovitch}.
\end{enumerate}
$\;$\\
\\
\textbf{Part II: Summary of stationary iterative methods}
\begin{enumerate}
\item (2 pts) Show the following equivalences:
\begin{itemize}
\item Jacobi:
\begin{equation*}
\begin{split}
&(\vec{u}_{k+1})_i = \frac{1}{a_{ii}}\left[\vec{f}_i - \sum_{j=1}^{i-1} a_{ij} (\vec{u}_{k})_j - \sum_{j=i+1}^{n} a_{ij} (\vec{u}_{k})_j \right] \\
&\Leftrightarrow \;
\vec{u}_{k+1} = - D^{-1} (L+U) \vec{u}_{k} + D^{-1} \vec{f}.
\end{split}
\end{equation*}
\item Gauss-Seidel:
\begin{equation*}
\begin{split}
&(\vec{u}_{k+1})_i = \frac{1}{a_{ii}}\left[\vec{f}_i - \sum_{j=1}^{i-1} a_{ij} (\vec{u}_{k+1})_j - \sum_{j=i+1}^{n} a_{ij} (\vec{u}_{k})_j \right] \\
&\Leftrightarrow \;
\vec{u}_{k+1} = - (D+L)^{-1} U \vec{u}_{k} + (D+L)^{-1} \vec{f}.
\end{split}
\end{equation*}
\item SOR:
\begin{equation*}
\begin{split}
&(\vec{u}_{k+1})_i = (1-\omega)(\vec{u}_{k})_i + \frac{\omega}{a_{ii}}\left[\vec{f}_i - \sum_{j=1}^{i-1} a_{ij} (\vec{u}_{k+1})_j - \sum_{j=i+1}^{n} a_{ij} (\vec{u}_{k})_j \right] \\
&\Leftrightarrow \;
\vec{u}_{k+1} = - (D+\omega L)^{-1} \left[ \bigl( (1-\omega)D - \omega U \bigr) \vec{u}_{k} + \omega \vec{f} \right].
\end{split}
\end{equation*}
\end{itemize}
\item (3 pts) Consider $A$ the discrete matrix of $-\Delta$ in dimension 2.
Show that the method of Jacobi and the method of Richardson (with the optimal parameter $\alpha_{opt}$
computed in the Série 3) have the same convergence behavior.
\item (3 pts) Consider a matrix $A \in \R^{2 \times 2}$ and the splitting $A = D + L + U$ with $D$ invertible.
Show that if the method of Jacobi converges, then the method of Gauss-Seidel converges as well with a faster
rate of convergence.
\item (4 pts) Consider the $A\in \R^{n\times n}$ the discrete matrix of $-\Delta$ in dimension 2
and let $h=\frac{1}{n+1}$ be the mesh size.
Denote by $\rho(G_{\rm J};h)$, $\rho(G_{\rm GS};h)$, and $\rho(G_{\rm SOR};h)$ the spectral radii
of the Jacobi, Gauss-Seidel, and SOR methods as a function of $h$.
Prove that:
\begin{itemize}\itemsep0em
\item $\rho(G_{\rm J};h) = 1 - \frac{\pi^2}{2} h^2 + O(h^4)$.
\item $\rho(G_{\rm GS};h) = 1 - \pi^2 h^2 + O(h^4)$.
\item $\rho(G_{\rm SOR};h) = 1 - 2 \pi h + O(h^2)$.
\end{itemize}
\end{enumerate}
$\;$\\
\\
\textbf{Part III: Steepest descent and conjugate gradient methods}
\begin{enumerate}\itemsep0em
\item (4 pts) Consider a function $f : \R^n \rightarrow \R$ and denote by $\nabla f(\vec{v})$ the
gradient of $f$ in $\vec{v} \in \R^n$ obtained by the usual scalar product for $\R^n$.
Prove that $\nabla f(\vec{v})$ is orthogonal to the level set of $f$ in $\vec{v}$ and that
$-\nabla f(\vec{v})$ defines the direction of steepest descent of $f$ in $\vec{v}$.
\item (4 pts) Implement the steepest descent method to solve a linear system $A \vec{u} = \vec{f}$.
\item (4 pts) Implement the conjugate gradient method to solve a linear system $A \vec{u} = \vec{f}$.
\item (1 pts) Test your codes to solve the Laplace equation.
\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{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}
\usepackage{float}
\usepackage{algorithm}
\usepackage{algorithmic}
\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}
%% commandes utiles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\tr}{\mbox{tr}}
\newcommand{\Span}{\mbox{span}}
\DeclareMathOperator*{\argmin}{arg\,min}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\noindent UNIVERSIT\'E DE GEN\`EVE \hfill Section de Math\'ematiques \\
\noindent Facult\'e des Sciences \hfill \`A rendre le mercredi 17 mai \\[-3mm] \hrule
\bigskip
\begin{center}
\textbf{\Large{Méthodes Itératives}} \\[1mm]\bigskip
\textbf{Série 5}
\end{center}
\bigskip
\textbf{Exercice 1: méthode du gradient conjugué.} (4 pts)\\
Considérer l'algorithme suivant qui donne la méthode du gradient conjugué pour résoudre
le système linéaire $A \vec{u} = \vec{f}$.
\begin{algorithm}\caption{(Algorithme du gradient conjugué)}
\begin{algorithmic}[1]
\REQUIRE Approximation initiale $\vec{u}_0$, nombre maximal d'itérations $k_{max}$, tolérance $tol$;
\ENSURE Solution approximative $\vec{u}_k$;
\STATE $\vec{r}_0 = \vec{f} - A \vec{u}_0$;
\STATE $\vec{p}_0 =\vec{r}_0$ et $k=0$;
\WHILE{$k < k_{max}$ \AND $\|\vec{r}_k\|_2 > tol$}
\STATE $\alpha_k = \frac{ \vec{r}_k^\top \vec{p}_k }{ \vec{p}_k^\top A \vec{p}_k }$;
\STATE $\vec{u}_{k+1} = \vec{u}_{k} + \alpha_k \vec{p}_{k}$;
\STATE $\vec{r}_{k+1} = \vec{f} - A \vec{u}_{k+1}$;
\STATE $\beta_k = \frac{ \vec{r}_{k+1}^\top A \vec{p}_k}{\vec{p}_{k}^\top A \vec{p}_k }$;
\STATE $\vec{p}_{k+1} = \vec{r}_{k+1} + \beta_k \vec{p}_{k}$;
\STATE $k = k + 1$;
\ENDWHILE
\RETURN $\vec{u}_k$;
\end{algorithmic}
\end{algorithm}
Noter que dans la boucle ``while'' il y a 4 produits matrice-vecteur.
Récrire l'algorithme tel qu'un seul produit matrice-vecteur est utilisé
dans la boucle. Essayer également de réduire le nombre de produits vecteur-vecteur.
\underline{Indication :} montrer que $\alpha_k = \frac{ \|\vec{r}_k\|_2^2}{ \vec{p}_k^\top A \vec{p}_k }$
et $\beta_k = \frac{ \|\vec{r}_{k+1}\|_2^2}{ \|\vec{r}_{k}\|_2^2 }$.
\bigskip
\textbf{Exercice 2: Optimalité du pas $\alpha_k$.} (3 pts)
Soit $A \in \mathbb{R}^{n\times n}$ symétrique et définie positive. Montrer que :
\begin{itemize}
\item le pas $\alpha_k = \frac{ \vec{r}_k^\top \vec{p}_k }{ \vec{p}_k^\top A \vec{p}_k }$ est optimal
dans le sens où il minimise la fonction
$\alpha \mapsto \phi(\alpha) := \frac{1}{2} (\vec{u}_k + \alpha \vec{p}_k )^\top A (\vec{u}_k + \alpha \vec{p}_k )
- \vec{f}^\top (\vec{u}_k + \alpha \vec{p}_k )$.
\item $\alpha_k$ minimise aussi $\alpha \mapsto \varphi(\alpha) := \| \vec{u}_k + \alpha \vec{p}_k - \vec{u} \|_A^2$.
\end{itemize}
\bigskip
\textbf{Exercice 3: Méthodes de Lanczos et du gradient conjugué.}\\
Dans cet exercice vous montrerez que la méthode du gradient conjugué,
pour résoudre $A \vec{u} = \vec{f}$ ($A\in \mathbb{R}^{n \times n}$ symétrique et définie positive),
peut être obtenue par la méthode de Lanzcos.
Considérer les étapes suivantes:
\begin{enumerate}\itemsep0em
\item (3 pts) Rappelons la décomposition $AQ_k = Q_{k+1} H_k$, où
$Q_k = \begin{bmatrix} \vec{q}_1 & \cdots & \vec{q}_k \end{bmatrix} \in \mathbb{R}^{n \times k}$
et supposons que $\vec{q}_1 = \frac{1}{\|\vec{f}\|_2} \vec{f}$.
La matrice $H_k$ a la forme
$H_k = \begin{bmatrix}
T_k \\
\vec{v} \\
\end{bmatrix}$
\begin{equation*}
T_k = \begin{bmatrix}
\alpha_1 & \beta_1 \\
\beta_1 & \alpha_2 & \beta_2 \\
& & \ddots \\
& & \beta_{k-1} & \alpha_k \\
\end{bmatrix},
\end{equation*}
et $\vec{v} = \begin{bmatrix}
0 & \cdots & 0 & \beta_k
\end{bmatrix}$.
\begin{itemize}\itemsep0em
\item[(a)] Montrer que $Q_k^\top A Q_k = T_k$.
\item[(b)] Montrer que
\begin{equation*}
\vec{u}_k := \argmin_{\vec{u}_k \in \Span\{ \vec{q}_1 \, \cdots \, \vec{q}_k \}} \frac{1}{2} \vec{u}_k^\top A \vec{u}_k-\vec{u}_k^\top \vec{f} = Q_k T_k^{-1} \bigl( \| \vec{f}\|_2 \vec{e}_1 \bigr) .
\end{equation*}
\underline{Indication :} Considérez un vecteur $\vec{y} \in \mathbb{R}^k$ tel que $\vec{u}_k = Q_k \vec{y}$
et montrez que $\vec{y} = T_k^{-1} Q_k^\top \vec{f}$.
\end{itemize}
\item (3 pts) Considérer la décomposition LU de $T_k$, c'est-à-dire $T_k = L_k U_k$
\begin{equation*}
L_k = \begin{bmatrix}
1 \\
\gamma_1 & 1 \\
& \ddots & \ddots \\
& & \gamma_{k-1} & 1 \\
\end{bmatrix}
\quad \text{et} \quad
U_k = \begin{bmatrix}
\eta_1 & \beta_1 \\
& \eta_2 & \beta_2 \\
& & \ddots & \\
& & & \eta_k \\
\end{bmatrix}.
\end{equation*}
\begin{itemize}\itemsep0em
\item[(a)] Montrer que $\eta_1 = \alpha_1$, et $\gamma_{j-1} = \frac{\beta_{j-1}}{\eta_{j-1}}$ et $\eta_j = \alpha_j - \gamma_{j-1} \beta_{j-1}$
pour $j=2,\dots,k$.
\item[(b)] Définir $P_k := Q_k U_k^{-1}$ et montrer que les colonnes $\vec{p}_j$ de $P_k$ satisfont
$\vec{p}_1 = \frac{1}{\eta_1} \vec{q}_1$ et
$\vec{p}_j = \frac{1}{\eta_j}\bigl( \vec{q}_j - \beta_{j-1} \vec{p}_{j-1} \bigr)$ pour $j=2,\dots,k$.
\item[(c)] Montrer que $\vec{u}_k = \vec{u}_{k-1} + \zeta_k \vec{p}_k$ pour un $\zeta_k \in \mathbb{R}$.
\underline{Indication :} commencer avec $\vec{u}_k = Q_k T_k^{-1} \bigl( \| \vec{f}\|_2 \vec{e}_1 \bigr)$
et utiliser la décomposition LU de $T_k$.
\end{itemize}
\item (4 pts) Montrer que:
\begin{itemize}\itemsep0em
\item[(a)] les vecteurs $\vec{p}_k$ sont $A$-orthogonaux.
\item[(b)] $\vec{r}_k = \sigma_k \vec{q}_{k+1}$ pour un $\sigma_k \in \mathbb{R}$,
$\vec{r}_k = \vec{f} - A \vec{u_k}$ est le résidu.
\item[(c)] $\vec{r}_k \perp \vec{q}_j$ pour $j=1,\dots,k$ et que $\vec{r}_k \perp \vec{r}_{k-1}$.
\end{itemize}
\item (3 pts) Considérons un vecteur $\tilde{\vec{p}}_{k-1}$ parallèle à $\vec{p}_{k}$. Alors on peut récrire
$\vec{u}_k = \vec{u}_{k-1} + \zeta_k \vec{p}_k$ (2.c) comme $\vec{u}_k = \vec{u}_{k-1} + \tilde{\alpha}_{k-1} \tilde{\vec{p}}_{k-1}$ pour un certain $\tilde{\alpha}_{k-1} \in \mathbb{R}$.
\begin{itemize}\itemsep0em
\item[(a)] Montrer que $\tilde{\vec{p}}_{k} = \vec{r}_{k} + \tilde{\beta}_{k-1} \tilde{\vec{p}}_{k-1}$
pour un certain $\tilde{\beta}_{k-1} \in \mathbb{R}$. \underline{Indication :} Rappeler 2.b.
\item[(b)] Calculer $\tilde{\alpha}_k$ en utilisant l'orthogonalité $\vec{r}_k \perp \vec{r}_{k-1}$.
\item[(c)] Calculer $\tilde{\beta}_{k-1}$ en utilisant la $A$-orthogonalité des vecteurs $\tilde{p}_k$.
\underline{Indication :} vous devez obtenir que $\tilde{\alpha}_k = \frac{ \| \vec{r}_k \|_2^2 }{ \tilde{\vec{p}}_k^\top A \tilde{\vec{p}}_k}$ et $\tilde{\beta}_k = \frac{ \| \vec{r}_k \|_2^2 }{ \| \vec{r}_{k-1} \|_2^2 }$.
\end{itemize}
\item (1 pts) Résumez les résultats obtenus et comparez-les avec l'algorithme de l'Ex. 1.
\end{enumerate}
\bigskip
\textbf{Exercice 4: Espace de Krylov.} (pts 6)\\
Considérer les vecteurs $\vec{r}_j$ et $\vec{p}_j$ pour $j=1,\dots,k$ générés
par les itérations de la méthode du gradient conjugué jusqu'à l'étape $k$.
Montrer que
\begin{equation*}
\Span \{ \vec{r}_0 , \dots, \vec{r}_k \} = {\cal K}_k(A,\vec{r}_0)
\quad \text{et} \quad
\Span \{ \vec{p}_0 , \dots, \vec{p}_k \} = {\cal K}_k(A,\vec{r}_0).
\end{equation*}
\bigskip
{\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{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}
\usepackage{float}
\usepackage{algorithm}
\usepackage{algorithmic}
\renewcommand{\algorithmicrequire}{\textbf{Input:}}