% article mod�e 2D avec porosit�variable et discontinu sans darcy et forchheimer \documentclass[12pt,a4 paper]{article} \usepackage[latin1]{inputenc} \usepackage[T1]{fontenc} \usepackage{graphicx} \usepackage{amsfonts,amssymb,amsmath} \usepackage{amstext} \usepackage{natbib} \newcommand{\n}{\boldsymbol n} \newcommand{\btau}{\boldsymbol \tau} %\usepackage{doublespace} \title{GMRES in FVLib} \author{ S. Clain $^a$} \begin{document} \maketitle \section{The theory} \subsection{Projection method and Krylov Spaces} Let consider a linear problem $Ax=b$ in $\mathbb R^n$ with $A$ an invertible matrix. For $k\leq n$, we consider two sub vectorial spaces $\mathcal K_k$, $\mathcal L_k$ of $\mathbb R^n$ of dimension $k$. Let $x_0\in\mathbb R^n$, the projection method consists in providing $x_k\in \mathbb R^n$ such that we both have $$ x_k\in x_0+\mathcal K_k, \qquad Ax_k-b \perp \mathcal L_k. $$ In other words, we seek $u_k\in\mathcal K_k$ such that \begin{equation} ==0,\qquad \forall v\in \mathcal L_k \label{Petrov_Galerkin} \end{equation} with $r_0=Ax_0-b$.\\ If $A$ is invertible, we have existence and uniqueness of the solution. We introduce the Krylov space $K_k(A,r_0,)$ given by $$ K_k(A,r_0,)=\textrm{Span}(r_0,Ar_0,...,A^{k-1}r_0) $$ and set $\mathcal K_k=K_k(A,r_0)$, $\mathcal L_k=AK_k(A,r_0,)=\textrm{Span}(Ar_0,...,A^{k-1}r_0,A^kr_0)$.\\ {\sc proposition} $x_k$ is the solution of (\ref{Petrov_Galerkin}) if and only if $x_k$ minimizes the functional \begin{equation} u\in\mathcal K_k \to E(x_0+u)= \label{minimize_functional} \end{equation} with corresponds to the minimization of the residu on $\mathcal K_k$.\\ Moreover we have $E(x_0+u_{k+1})\leq E(x_0+u_{k})$, thus the residual always decreases. $\blacksquare$\\ \noindent If $u_k\in \mathcal K_k$ is the minimizer, then we have using definition of space $\mathcal L_k$ \begin{equation} =0 \quad \forall v\in \mathcal K_k. \label{minimizer} \end{equation} We have the following result.\\ {\sc proposition} If the Krylov space is invariant by the composition with $A$ \begin{equation} AK_k(A,r_0) = K_k(A,r_0) \label{Krylov_invariant} \end{equation} then $x_k=x_0+u_k$ is the solution, namely $Ax_k-b=0$. $\blacksquare$\\ One can prove that there exists at least a $k\leq n$ such that (\ref{Krylov_invariant}) is achieved but $k$ coul be very large. In the iterative method, we only require that the residual part $E(x_0+u_{k})$ is small enough. \subsection{Orthogonal basis and the Arnoldi Algorithm} When property (\ref{Krylov_invariant}) is not achieved, vectors $r_0,Ar_0,...,A^{k-1}r_0$ define a basis but the major inconvenient is that the vector are merely parallel (remember that $A^k r_0$ converges to the first eigenvector). In order to provide a more suitable basis we build an orthogonal basis using the Arnoldi argorithm.\\ We start setting $q_1=r_0/\vert r_0\vert$. Assume now that $q_1,q_2,..,q_k$ is a base of $K_k(A,r_0)$. If property (\ref{Krylov_invariant}) is not achieved then we can prove that $Aq_k\notin K_k(A,r_0)$ and we write $$ A q_k=\sum_{i=1}^{k+1} h_{jk} q_j $$ where $q_{k+1}$ is a normalized vector, we have to determine. Using the orthogonal condition, we compute $h_{jk}=$, $j=1,...,k$. Then we set $$ p_{k+1}=Aq_k-\sum_{i=1}^k h_{jk} q_j, \qquad h_{k+1,k}=\vert p_{k+1}\vert, \qquad q_{k+1}=\frac{p_{k+1}}{h_{k+1,k}}. $$ From a matricial point of view, we have $AQ_k=Q_{k+1}H(k+1,k)$ where $Q_k$ is a $n\times k$ Stiefel matrix and $H(k+1,k)$ is a $k+1\times k$ Hessenberg matrix $$ Q=[q_1,...,q_k], \quad H(k+1,k)=\left [ \begin{array}{ccccc} h_{11} & \hdots & \hdots &\hdots &h_{1k}\\ h_{21} & h_{22} & \hdots &\hdots & \vdots \\ 0 & h_{32} & \ddots &\hdots & \vdots \\ \vdots & \hdots & \ddots &\ddots & \vdots \\ 0 & \hdots & 0 &h_{k,k.1}&h_{k,k} \\ 0 & \hdots & \hdots & 0 &h_{k+1,k} \end{array} \right ]. $$ Note that $h_{k+1,k}=0$ is equivalent to property (\ref{Krylov_invariant}). Hence we can proceed as long as $h_{k+1,k}>0$. \subsection{An equivalent minimizing condition} Vector $u_k$ is the unique element of $\mathcal K_k$ wich minimizes $$ E(x_0+u)==\Vert Au-r_0\Vert^2. $$ Using the orthogonal basis $Q$ of $\mathcal K_k$, we can write any $u \in\mathcal K_k$ as $$ u=\sum_{i=1}^k y^k_i q_i, \qquad y^k\in\mathbb R^k $$ Hence $Au=AQ_ky^k$. Note that $u$ has a fix dimension while dimension of vector $y^k$ changes with $k$. Moreover let us remark that $r_0=\vert r_0\vert Q_{k+1}e^{k+1}=(1,0,...,0)$ the first canonical vector of $\mathbb R^{k+1}$. We then have \begin{eqnarray*} E(x_0+u)&=&\Vert Au-r_0\Vert^2\\ &=&\Vert AQ_ky^k-\vert r_0\vert Q_{k+1}e^{k+1} \Vert^2\\ &=&\Vert Q_{k+1}H(k+1,k)y^k-\vert r_0\vert Q_{k+1}e^{k+1}\Vert^2\\ &=&\Vert Q_{k+1}\{H(k+1,k)y^k-\vert r_0\vert e^{k+1}\}\Vert^2 \end{eqnarray*} Since $Q_{k+1}$ is unitary, we deduce $$ E(x_0+u)=\widetilde E_k(y^k)=\Vert H(k+1,k)y^k-\vert r_0\vert e^{k+1}\Vert^2. $$ So we reduce the problem of minimization in $\mathbb R^k$ of the overdetermined system $H(k+1,k)y^k=\vert r_0\vert e^{k+1}$. \subsection{Orthogonalization of the Hessenberg matrix} To solve the minimizination problem, we use a $QR$ transformation of the $H(k+1,k)$ matrix using Givens rotations.\\ \noindent {\sc Step 1} One have $H(2,1)=[h_{11},h_{21}]^T$. We apply a Givens matrix (unitary matrix) of the form $$ G_1=\left [ \begin{array}{cc} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{array} \right ]. $$ If we choose $$ \cos(\theta)=\frac{h_{11}}{\sqrt{h_{11}^2+h_{21}^2}}, \quad \sin(\theta)=\frac{h_{21}}{\sqrt{h_{11}^2+h_{21}^2}} $$ We find $$ U(2,1)=G_1H(2,1)=\left [ \begin{array}{c} \sqrt{h_{11}^2+h_{21}^2}\\ 0 \end{array} \right ], \qquad c^1=\vert r_0\vert G_1 e^1_{2}=\vert r_0\vert \left [ \begin{array}{c} \displaystyle\frac{h_{11}}{\sqrt{h_{11}^2+h_{21}^2}}\\ \displaystyle\frac{-h_{21}}{\sqrt{h_{11}^2+h_{21}^2}} \end{array} \right ]= \left [ \begin{array}{c} c^1_1 \\c^1_2 \end{array} \right ]. $$ \noindent {\sc Step k} We first note that $$ H(k+1,k)=\left [ \begin{array}{cc} H(k,k-1) & h_{1k} \\ & \vdots \\ 0 & h_{k+1,k} \end{array} \right ]. $$ Now assume that we have determined the Givens matrizes $G_1$,...,$G_{k-1}$ such that we have a $k\times k-1$ upper triangular matrix $U(k,k-1)$ and a vector $c^{k-1} \in\mathbb R^k$ and $$ \Vert H(k,k-1)y^{k-1}-\vert r_0\vert e^{k-1}\Vert=\Vert U(k-1,k)y^{k-1}-\vert r_0\vert c^{k-1}\Vert $$ Problem of minimizing $\Vert H(k+1,k)y^k-\vert r_0\vert e^{k}\Vert$ is equivalent to minimize $\Vert \widetilde U(k+1,k) y^k -\vert r_0\vert \widetilde c^{k}\Vert$ with $$ \widetilde U(k+1,k)=\left [ \begin{array}{cc} U(k,k-1) & h_{1k} \\ & \vdots \\ 0 & h_{k+1,k} \end{array} \right ], \quad \widetilde c^{k}=\left [ \begin{array}{c} c^{k-1}_1 \\ \vdots \\c^{k-1}_k\\0 \end{array} \right ]. $$ We consider the $k+1\times k+1$ Givens matrix $$ G_k=\left [ \begin{array}{ccc} Id & 0 & 0 \\ 0 & cs_k & sn_k \\ 0 &-sn_k & cs_k \end{array} \right ], $$ if one chooses $$ cs_k=\frac{h_{kk}}{\sqrt{h_{kk}^2+h_{k+1,k}^2}}, \quad sn_k=\frac{h_{k+1,k}}{\sqrt{h_{kk}^2+h_{k+1,k}^2}} $$ then we have $$ U(k+1,k)=G_k\left [ \begin{array}{cc} & \widetilde h_{1k} \\ U(k,k-1) & \vdots \\ & \widetilde h_{k-1,k} \\ 0 & \widetilde h_{k,k} \\ 0 & h_{k+1,k} \end{array} \right ], \quad c_{k}=G_k\left [ \begin{array}{c} c^{k-1}_1 \\ \vdots \\c^{k-1}_{k-1} \\ c^{k-1}_{k}\\ 0 \end{array} \right ], $$ where we have compute the successive products $\widetilde {h}^k=G_{k-1}G_{k-2}...G_2G_1 h^k$. Note that we We get the solution $y^k$ using a backward substitution of system $$ U(k+1,k) y^k= c^k $$ starting at line $k$ since the last line (k+1) is degenerated. We evaluate the minimum with $$ x_k=x_0+\sum_{i=1}^y y^k_i q_i. $$ \end{document}