1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
|
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"lang": "fr"
},
"source": [
"__Programmation vectorielle__\n",
"\n",
"Le but des exercices est\n",
"\n",
"* d'avoir un programme qui donne la bonne réponse\n",
"* qui soit le plus rapide possible (et pour cela on utilise massivement Numpy)\n",
"\n",
"En règle général si vous avez des `for` imbriqués c'est mauvais signe."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"np.set_printoptions(precision=10, linewidth=150, suppress=True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"lang": "fr"
},
"source": [
"## Méthode du pivot de Gauss partiel\n",
"\n",
"L'ennoncé est dans le cours. On verifiera sur le cas qui génère des\n",
"erreurs d'arrondis."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x mieux = [1.0000019 0.99999905]\n",
"[1. 3.]\n"
]
}
],
"source": [
"def solve_partial_gauss(A, b):\n",
" A = A.astype(np.float64) # si A a des entiers, on va avoir des erreurs de calculs\n",
" for i in range(len(A)-1):\n",
" pivot = i + np.argmax(np.abs(A[i:,i]))\n",
" A[[i,pivot]] = A[[pivot, i]]\n",
" b[[i,pivot]] = b[[pivot, i]]\n",
" coefs = - A[i+1:,i] / A[i,i] # Normalement il faut tester que A[i,i] != 0\n",
" A[i+1:, i:] += np.outer(coefs, A[i, i:]) # ou np.einsum('i,j -> ij', coefs, A[i, i:])\n",
" b[i+1:] += coefs * b[i]\n",
" # A est maintenant triangulaire surpérieur\n",
" res = np.zeros(len(b), dtype=b.dtype)\n",
" res[-1] = b[-1] / A[-1,-1]\n",
" for i in range(len(A)-1)[::-1]:\n",
" res[i] = (b[i] - A[i,i+1:] @ res[i+1:]) / A[i,i]\n",
" return res\n",
"e = 1E-6\n",
"A = np.array([[e, 1], [1, 2]], dtype='float32')\n",
"b = np.array([1., 3.], dtype='float32')\n",
"x = solve_partial_gauss(A,b)\n",
"print(\"x mieux =\", x)\n",
"print(A @ x )"
]
},
{
"cell_type": "markdown",
"metadata": {
"lang": "fr"
},
"source": [
"## Factorisation de Choleski\n",
"\n",
"Il s'agit de décomposer A en $A = B\\, B^T$ avec B une matrice triangulaire inférieure. Cela n'est possible\n",
"que si la matrice A est symétrique et définie positive (c'est d'ailleurs un facon de vérifier qu'une\n",
"matrice est définie positive).\n",
"\n",
"Écrire l'algorithme de Choleski qui prend A et retourne B (pour deviner l'algorithme, essayez de trouver les \n",
"coefficients de B à partir des coefficients de A sur une matrice A 4x4)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"lang": "fr"
},
"source": [
"Rappel : pas de boucles `for` imbriquées mais des opérations vectorielles et matricielles (sur des sous-matrices).\n",
"\n",
"Créer une matrice symétrique définie positive et vérifier que sonprogramme marche bien."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"lang": "fr"
},
"source": [
"## Amméliorer Jacobi\n",
"\n",
"Lorsqu'on écrit une itération de la méthode de Jacobi avec l'ensemble des coefficients, on constate que\n",
"si on calcule la nouvelle valeur de **x** élément par élement alors lorsqu'on veut mettre à jour `x[1]`, \n",
"on connait déjà `x[0]`. Idem lorsqu'on met à jour `x[2]` on connait déjà `x[0]` et `x[1]`, etc.\n",
"\n",
"L'idée de la version amméliorée de Jacobi est d'utiliser la nouvelle valeur de `x[0]` et non pas l'ancienne\n",
"comme c'est le cas dans l'algorithme du cours. Ainsi en utilisant des valeurs mise à jour on peut espérer\n",
"converger plus vite.\n",
"\n",
"Dans ce cas, chaque itération demande un calcul ligne par ligne et donc une boucle `for` dans la boucle `for` sur\n",
"les itérations."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"lang": "fr"
},
"source": [
"#### Test d'arrêt\n",
"\n",
"On ajoutera un argument `error` à la fonction qui indique la précision désirée du résultat. Malheureusement pour connaitre la précision de notre calcul il faut connaitre la solution. On pourrait aussi calculer \n",
"${\\bf A \\, x}^t - {\\bf b}$ mais ce n'est pas non plus l'écart entre ${\\bf x}^t$ et la solution (en plus\n",
"c'est un calcul en n² donc lourd).\n",
"\n",
"Aussi on va regarder quand la\n",
"convergence ralentit et donc on s'arrête lorsque\n",
"$||{\\bf x}^{t+1} - {\\bf x}^t|| < \\texttt{error}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"lang": "fr"
},
"source": [
"#### Bonus\n",
"\n",
"Réécrire la méthode sans la boucle `for` mais en prenant bien le nouvelles valeurs de `x`. Pour cela on va découper la matrice **A** en deux matrices triangulaires."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# les méthodes de Numpy triu et tril seront utiles\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
|
