Je dois programmer sur Python la méthode de Jacobi, de Gauss-Seidel et SOR pour résoudre le système Ax=b.
On utilise une décomposition A=M-N avec :
• M = D et N = A-M pour Jacobi
• M = D + E et N = A-M pour Gauss-Seidel, avec L matrice triangulaire inférieure stricte de A
• $ M=\frac{1}{\omega} D +L $, N=A-M pour SOR
On doit écrire une fonction jacobi(A,b,Imax,eps) qui prend en entrée une matrice A, un vecteur b, un nombre d'itérations maximal Imax et une tolérance eps (epsilon). On utilisera un test d'arrêt du type $ [ k \leq Imax \text{ ou } ||r^k||_2 \leq \epsilon] $ où $ r^k = Ax^k - b $ est le "résidu" à l'étape k.
Pour la méthode de Gauss-Seidel, pareil et pour la méthode SOR, on rajoute un paramètre omega.
Pour les tests, on prendra Imax=1000 et eps=10^(-3).
Ensuite je dois faire un tableau montrant les temps de calcul de la méthode pour obtenir $ ||r^n|| \leq 10^{-3} $.
Je mets mon code pour SOR :
Code : Tout sélectionner
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
def SOR(A,b,Imax,eps,x0,omega):
D=np.diag(np.diag(A))
M=np.dot((1/omega),D)+np.tril(A,-1)
N=M-A
invM = la.inv(M)
r=np.dot(A,x0)-b
x=x0
i=0
err=1+eps
res=[]
while((i<Imax) and ((la.norm(r))>=eps)):
x=np.dot(np.dot((invM),N),x)+np.dot((invM),b)
r=np.dot(A,x)-b
err = la.norm(r,2)
res.append(err)
i=i+1
return (x,i,res)
# on peut bien sûr modifier les différentes matrices A, b et x0
A=np.array([[2,0,0],[4,5,0],[7,8,9]])
x0=np.array([[1],[1],[1]])
b=np.array([[20],[8],[7]])
print(SOR(A,b,1000,10**(-3),x0,0.5))