(Python et mathématiques) Comment changer la taille de ma matrice ?

Bonjour,
Je suis débutante en Python et je dois faire un code qui programme les différentes méthodes de relaxation (Jacobi, Gauss-Seidel et SOR) pour calculer Ax=b.

Pour la méthode de Jacobi, on décompose A=M-N avec M=diagonale(A) et N=A-M. Partant d'un x0 donné, on écrit Mx^(k+1) = Nx^(k) + b.
On cherche à écrire une fonction qui prend en entrée une matrice A, un vecteur colonne b, un nombre d'itérations maximal Imax, une tolérance eps (pour epsilon). On utilisera un test d'arrêt [k <= Imax ou ||r^k|| < eps, où r^k = Ax^k - b est le "résidu" à l'étape k.

Voici mon code :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

def jacobi(A,b,Imax,eps,x0): D=np.diag(np.diag(A)) N=D-A 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((la.inv(D)),N),x)+np.dot((la.inv(D)),b) r=np.dot(A,x)-b err = la.norm(r,2) res.append(err) i=i+1 return (x,i,res)

---

Puis j'ai fait un test avec une matrice A de taille 3x3 pour voir si mon code est correct.

Code:

---

1 2 3 4

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(jacobi(A,b,1000,10**(-3),x0))

---

ça me semble pas trop mal, mais j'ai quelques questions.

Pensez-vous qu'il y ait des erreurs dans mon code ?
L'énoncé me demande de tester le code avec une matrice 20x20. Je me demande comment tester avec une telle matrice sans entrer les coefficients un par un (ce qui me semble long et fastidieux). Avez-vous une idée?

Merci d'avance, je suis débutante en Python, j'espère que mon sujet est bien placé...

Salut,

Tout d'abord, es-tu sûre du double np.diag ?

Code:

---

D=np.diag(np.diag(A))

---

Ici:

Code:

---

1 2	while ((i<Imax) and ((la.norm(r))>=eps)):

---

d'où provient la ?
Il est définit après mais pas avant, tu auras une NameError.

Autre chose qui paraîtra peut-être accessoire mais qui facilite la vie, aére ton code et nomme les choses de leur vrai nom.

En Python on appele un chat un chat. À la rigueur felix mais pas C ou fgh.
Cela évite de rendre le code aussi indigeste qu'un bon cours de math.

Pour tester avec des matrices de tailles différentes, ajoute à ton code une fonction qui les créera pour toi:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

>>> def get_matrix(max_, size): ... mat = [] ... in_mat = [] ... for row in range(size): ... for colomn in range(size): ... in_mat.append(random.randrange(max_)) ... mat.append(in_mat) ... in_mat = [] ... return mat ... >>> matrix = get_matrix(10, 20) >>> for i in matrix: ... print(i) ... [7, 5, 1, 5, 6, 5, 1, 6, 7, 6, 3, 9, 7, 6, 5, 3, 4, 3, 4, 9] [2, 6, 1, 5, 8, 6, 1, 1, 2, 3, 0, 5, 0, 0, 1, 4, 2, 2, 4, 0] .....

---

Quand tu écris dans le while :

Code:

---

x=np.dot(np.dot((la.inv(D)),N),x)+np.dot((la.inv(D)),b)

---

Imagine à ce moment que D soit une matrice très grosse. Donc tu exécute cette ligne à chaque itération, et à chaque itération tu calcules l'inverse de D. 2 fois en plus. Ca peut devenir couteux comme opération ca...
Pourquoi ne pas simplement le calculer une bonne fois pour toute avant le while en écrivant qqch du genre :

Code:

---

inv_D=la.inv(D)

---

On peut se le permettre ici vu que D ne change pas au fur et à mesure que tu itères.

Ensuite je trouve ca assez bizarre de retourner la liste des erreurs comises au fil des itérations. On peut retourner l'erreur à la dernière itération avec le x calculé, ca suffit amplement.

Enfin je metterais des valeurs par défaut au paramètre de ta fonction. L'utilisateur n'a pas forcément envie de s'embêter avec ces paramètres techniques lorsqu'il veut juste la solution d'un système. Donc proposer des valeurs par défaut pour Imax, eps et x0, permettrait à l'utilisateur de ne pas avoir à trop ce soucier de ceux-ci.


Pour faire des matrices de tailles arbitraires, tu peux utiliser np.random.
Avec ca tu peux générer une matrice A aléatoire (avec une taille que tu choisi) et un vecteur x au hasard aussi.
Là tu calcule b=Ax. Ensuite tu files A et b à ta fonction et tu vérifies que ta fonction te renvoie bien le x de départ.

Merci beaucoup ! :)