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
|
{
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 2,
"cells": [
{
"cell_type": "markdown",
"metadata": {
"lang": "en"
},
"source": [
"## Vector Programming\n",
"\n",
"The aim of the exercises is\n",
"\n",
"* to have a program that gives the correct answer\n",
"* which is as fast as possible (and for this we use massively Numpy)\n",
"\n",
"Generally if you have nested `for` it's a bad sign."
]
},
{
"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": "en"
},
"source": [
"## Partial Gaussian pivot method\n",
"\n",
"The announcement is in the course. We will check on the case which generates\n",
"rounding errors."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"lang": "en"
},
"source": [
"## Choleski factorization\n",
"\n",
"This is to decompose A into $A = B\\, B^T$ with B a lower triangular matrix. This is not possible\n",
"that if the matrix A is symmetric and positive definite (this is moreover a way of verifying that a\n",
"matrix is positive definite).\n",
"\n",
"Write Choleski's algorithm that takes A and returns B (to guess the algorithm, try to find the\n",
"coefficients of B from the coefficients of A on a 4x4 matrix A)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"lang": "en"
},
"source": [
"Reminder: no nested `for` loops but vector and matrix operations (on sub-matrices).\n",
"\n",
"Create a positive definite symmetric matrix and verify that its program works."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"lang": "en"
},
"source": [
"## Improve Jacobi\n",
"\n",
"When we write an iteration of the Jacobi method with all the coefficients, we find that\n",
"if we calculate the new value of **x** element by element then when we want to update `x[1]`,\n",
"we already know `x[0]`. Ditto when we update `x[2]` we already know `x[0]` and `x[1]`, etc.\n",
"\n",
"The idea of the improved version of Jacobi is to use the new value of `x[0]` and not the old one\n",
"as is the case in the algorithm of the course. So by using updated values we can expect\n",
"converge faster.\n",
"\n",
"In this case, each iteration requires a line-by-line calculation and therefore a `for` loop within the `for` loop on\n",
"the iterations."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"lang": "en"
},
"source": [
"#### Test stop\n",
"\n",
"We will add an `error` argument to the function which indicates the desired precision of the result. Unfortunately, to know the accuracy of our calculation, it is necessary to know the solution. We could also calculate\n",
"${\\bf A \\, x}^t - {\\bf b}$ but it's also not the gap between ${\\bf x}^t$ and the solution (besides\n",
"it is a calculation in n² therefore cumbersome).\n",
"\n",
"Also we will look when the\n",
"convergence slows down and therefore we stop when\n",
"$||{\\bf x}^{t+1} - {\\bf x}^t|| < \\texttt{error}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"lang": "en"
},
"source": [
"#### Bonuses\n",
"\n",
"Rewrite the method without the `for` loop but taking the new values of `x`. For this we will cut the matrix **A** into two triangular matrices."
]
},
{
"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": []
}
]
}
|
