TP 3. Résolution numérique de systèmes linéaires¶
%pylab inline
Utilisation des méthodes de Gauss-Seidel et de Jacobi¶
Soit $n$ un entier naturel non nul et $r\in[0,1]$. On cherche à résoudre le système linéaire de taille $n$, $A_r\vec{x}=\vec{b}_\epsilon$ par les méthodes de Jacobi et Gauss-Seidel, où $$ A_r=\begin{pmatrix}1&r&r^2&0&\cdots&0\\r&1&r&\ddots&\ddots&\vdots\\r^2&r&\ddots&\ddots&\ddots&0\\0&\ddots&\ddots&\ddots&\ddots&r^2\\\vdots&\ddots&\ddots&\ddots&1&r\\0&\cdots&0&r^2&r&1\end{pmatrix}\mbox{ et }\vec{b}_r=A_r\begin{pmatrix}1\\1\\\vdots\\1\end{pmatrix}. $$ La matrice peut être créée en appelant la fonction $\texttt{getA(n,r)}$ qui est fournie. On pose dans un premier temps $n=5$.
(a) On sait que si $A_r$ est à diagonale strictement dominante par lignes, alors la méthode de Gauss-Seidel converge.
Vérifier que $A_r$ est bien à diagonale strictement dominante par lignes quand $r=0.3$.
Créer une fonction $\texttt{Jacobi(A,b,x0,eps)}$ qui résoud un système $Ax=b$ par la méthode de Jacobi en prenant $\texttt{x0}$ comme vecteur initial et qui itère jusqu'à ce que le résidu $r_k=\Vert A x_k-b\Vert\leq \texttt{eps}$. La fonction doit rendre le vecteur solution approché $x_k$, le nombre d'itérations k nécessaires pour atteindre convergence et la valeur du résidu.
Même question pour Gauss-Seidel.
Approcher la solution $A_r\vec{x}=\vec{b}_r$ avec $r=0.3$ avec ces méthodes. Laquelle est la plus rapide?
(b) Réprésentation graphique:
Soient $B^{(Jacobi)}_r$ et $B^{(GS)}_r$ les matrice d'itération des méthodes de Jacobi et Gauss-Seidel. Tracer $\rho(B_r)$ (le rayon spectral de la matrice d'itération) pour $r\in[0,1]$ dans chaque méthode.
Que dire de la convergence des deux méthodes en fonction de la valeur de $r$?
Quelle méthode choisir pour résoudre le système quand $\varepsilon=0.5$?
(c) On se concentre ici sur la méthode de Jacobi. Compléter la fonction $\texttt{Jacobi(A,b,x0,eps)}$ pour qu'elle itère un nombre maximal kmax de fois et pour qu'elle stocke tout l'historique des valeurs du résidu. On pose à présent $n=100$. Pour $r=0.3$ et $r=0.35$, tracer (en utilisant la commande $\texttt{semilogy}$) les graphes de la norme du résidu $\vec{r}_{k}=\vec{b}_r-A_r\vec{x}_{k}$ en fonction du nombre d'itérations $1\leq k\leq50$. Commenter en particulier la pente des courbes. Dans le cas o`u $r=0.35$, estimer à partir du graphe le nombre d'itérations nécessaires pour que la norme du résidu soit plus petite que $10^{-10}$.
def getA(n,r):
A=eye(n)
i,j = np.indices(A.shape)
A[i==j-1]=r
A[i==j+1]=r
A[i==j-2]=r**2
A[i==j+2]=r**2
return A
# Réponse question (a)
def Jacobi(A,b,x0,eps):
# Initialisations
k=0
x=x0 # x_{k}
x_precedent=x0 # x_{k-1}
residu=norm(dot(A,x)-b)
i,j = np.indices(A.shape)
while(residu>eps):
# On incrémente indice
k+=1
# Actualisation x_precedent
x_precedent=x
# Nouvelle solution x utilisant x_precedent
x=(b-dot(A-diag(diag(A)),x_precedent))/diag(A)
# Actualisation residu
residu=norm(dot(A,x)-b)
return [x,k,residu]
def GS(A,b,x0,eps):
# Initialisations
k=0
x=x0 # x_{k}
residu=norm(dot(A,x)-b)
i,j = np.indices(A.shape)
while(residu>eps):
# On incrémente indice
k+=1
# Actualisation x. IMPORTANT: On n'a plus besoind de déclarer x_precedent dans GS!
for i in range(n):
row_i_of_A_without_coef_i=A[i,arange(len(A[i,:]))!=i]
x_without_coef_i=x[arange(len(x))!=i]
x[i]=(b[i]-dot(row_i_of_A_without_coef_i,x_without_coef_i))/A[i,i]
# Actualisation residu
residu=norm(dot(A,x)-b)
return [x,k,residu]
n=5
r=0.3
eps=0.001
x0=zeros(n)
A=getA(n,r)
b=dot(A,ones(n))
[x_Jacob,k_Jacob,residu_Jacob]=Jacobi(A,b,x0,eps)
[x_GS,k_GS,residu_GS]=GS(A,b,x0,eps)
print(x_Jacob,k_Jacob,residu_Jacob)
print(x_GS,k_GS,residu_GS)
# Réponse question (b)
def rhoJacobi(A):
D=diag(diag(A))
B=dot(inv(D),D-A)
[valPropre,vectPropre]=eig(B)
return max(abs(valPropre))
def rhoGS(A):
i,j = np.indices(A.shape)
N=A*(i<j)
M=A*(i>=j)
B=dot(inv(M),N)
[valPropre,vectPropre]=eig(B)
return max(abs(valPropre))
nValues=100
rhoValues=linspace(0,1,num=nValues,endpoint=True)
rhoJacobiValues=[]
rhoGSvalues=[]
for r in rhoValues:
A=getA(n,r)
rhoJacobiValues.append(rhoJacobi(A))
rhoGSvalues.append(rhoGS(A))
fig = figure()
plot(rhoValues,rhoJacobiValues,label='Jacobi')
plot(rhoValues,rhoGSvalues,label='Gauss-Seidel')
plot(rhoValues,ones(nValues),label='Seuil convergence')
legend()
Réponses question (b):
On observe sur la figure ci-dessus que le rayon spectral de la matrice d'itération de la méthode de Gauss-Seidel est toujours plus petit que celui de la matrice d'itération de la méthode de Jacobi. On confirme bien ce qui a été constaté dans la section (a): la convergence de la méthode de Gauss-Seidel est plus rapide que celle de la méthode de Jacobi pour ce problème. Par ailleurs, on voit encore que la méthode de Jacobi ne converge que pour $r\in[0,0.4]$ (environ) car ce n'est que pour cet intervalle que le rayon spectral est inférieur à 1. Gauss-Seidel converge pour une plage plus grande de valeurs $r\in[0,0.7]$.
Lorsque $r=0.5$, Jacobi ne converge pas. Il faut donc obligatoirement utiliser GS.
# Réponse question (c): On modifie légèrement la fonction Jacobi pour que residu soit un vecteur qui
# stocke l'historique des valeurs et que l'on itère kmax fois.
def Jacobi_residu(A,b,x0,eps,kmax):
# Initialisations
k=0
x=x0 # x_{k}
x_precedent=x0 # x_{k-1}
residu=[ norm(dot(A,x)-b) ]
i,j = np.indices(A.shape)
while(residu[-1]>eps and k<kmax): # residu[-1] ---> Le dernier élément du vecteur residu
# On incrémente indice
k+=1
# Actualisation x_precedent
x_precedent=x
# Nouvelle solution x utilisant x_precedent
x=(b-dot(A-diag(diag(A)),x_precedent))/diag(A)
# Actualisation residu
residu.append(norm(dot(A,x)-b))
return [x,k,residu]
n=100
eps=1E-3
x0=zeros(n)
kmax=50
# r=0.3
r=0.3
A=getA(n,r)
b=dot(A,ones(n))
[x_Jacob_30,k_Jacob_30,residu_Jacob_30]=Jacobi_residu(A,b,x0,eps,kmax)
# r=0.35
r=0.35
A=getA(n,r)
b=dot(A,ones(n))
[x_Jacob_35,k_Jacob_35,residu_Jacob_35]=Jacobi_residu(A,b,x0,eps,kmax)
# Figure
fig = figure()
semilogy(residu_Jacob_30,label='r=0.3, Jacobi')
semilogy(residu_Jacob_35,label='r=0.35, Jacobi')
legend()
Réponse question (c):
On observe sur la figure que la convergence est beaucoup moins rapide pour $r=0.3$ car la pente du segment de droite bleu est moindre. Dans ce cas, la norme du résidu est à peu près divisée par 10 toutes les 40 itérations itérations (à l'itération 10, la norme du résidu vaut environ 10 et elle vaut environ 1 à l'itération 50). Il faudra donc environ 450 itérations pour obtenir un résidu inférieur à $10^{-10}$.