Semaine 3 : Systèmes linéaires

using LaTeXStrings
using LinearAlgebra
using Plots
using Polynomials
import Random

Interpolation et approximation

[Exercice 1] Approximation par moindres carrés

L'objectif de cet exercice est d'approximer des données x, y données par un polynôme de degré fixé, potentiellement bien inférieur au nombre de points de données.

1. Sans utiliser de bibliothèque, écrire une fonction polyfit(x, y, d) qui, étant donnés des points $(x_1, y_1), \dotsc, (x_N, y_N)$ et un entier $0 \leq d \leq N-1$, calcule l'unique polynôme $p \in \mathbb{R}[x]$ de degré au plus $d$ minimisant l'erreur totale

   $$ E := \frac{1}{2} \sum_{n=0}^{N} \bigl\lvert p(x_n) - y_n \bigr\rvert^2. $$

   La fonction doit renvoyer un vecteur contenant les $d+1$ coefficients de $p$, en commençant par le terme constant.

   Indice (cliquer pour afficher)

   Dans la fonction, on peut procéder comme suit. D'abord, créer la matrice et le vecteur suivants : $$ \mathbf{A} = \begin{pmatrix} 1 & x_0 & \dots & x_0^d \\ \vdots & \vdots & & \vdots \\ 1 & x_{N} & \dots & x_N^d \end{pmatrix}, \qquad \mathbf{b} = \begin{pmatrix} y_0 \\ \vdots \\ y_N \end{pmatrix}. $$ On remarque que l'erreur $E$ se réécrit comme suit : $$ E(\boldsymbol{\alpha}) = \frac{1}{2} \bigl\lVert \mathsf{A} \boldsymbol{\alpha} - \boldsymbol{b} \bigr\rVert^2, $$ où $\boldsymbol{\alpha}$ est un vecteur contenant les coefficients du polynôme, par ordre de degré croissant, et $\lVert \cdot \rVert$ est la norme euclidienne. En d'autres termes, la fonction $E$ est une fonction quadratique du vecteur des coefficients du polynôme. Écrire $\nabla E(\boldsymbol{\alpha}) = 0$ conduit aux dites équations normales : $$ \mathsf{A}^{\top} \mathsf{A} \boldsymbol{\alpha} = \mathsf{A}^{\top} \mathbf{b}. $$ C'est un système linéaire avec une matrice carrée inversible au membre de gauche, qui peut être résolu en utilisant l'opérateur backslash \ ; en fait on peut écrire simplement A\b au lieu de (A'A)\(A'b), car l'opérateur \ effectue une approximation par moindres carrés par défaut.

function polyfit(x, y, d)
    ### BEGIN SOLUTION
    A = (^).(x, transpose(0:d))
    return A\y
    ### END SOLUTION
end;

n = 10 ; x = 1:n ; Random.seed!(2025); y = randn(n)
@assert polyfit([0.], [0.], 0) ≈ [0.]
@assert polyfit(1:5, 1:5, 1) ≈ [0., 1.]
@assert polyfit(x, y, 0) ≈ [sum(y)/n]
@assert polyfit(x, y, 0) ≈ [sum(y)/n]
@assert polyfit(x, y, 2) ≈ fit(x, y, 2).coeffs

2. Écrire également une fonction polyval(α, X) pour évaluer le polynôme $$ p(x) = \alpha_0 + \alpha_1 x + \dotsc + \alpha_d x^d $$ en tous les points de X. La fonction doit renvoyer le résultat dans un vecteur.

function polyval(α, X)
    ### BEGIN SOLUTION
    d = length(α) - 1
    (^).(X, transpose(0:d)) * α
    ### END SOLUTION
end;

n = 10 ; Random.seed!(2025); α = randn(n)
@assert polyval([0.], [1., 2., 3.]) == [0., 0., 0.]
@assert polyval(α, [0., 1.]) ≈ [α[1], sum(α)]

Utiliser les données ci-dessous, donnant l'altitude d'une bille en chute libre en fonction du temps sur une planète lointaine, pour tester votre code. L'expérience a été réalisée sur une autre planète. Pouvez-vous identifier laquelle d'après le coefficient du terme quadratique ? Voir Accélération gravitationnelle.

# Time since dropping the marble
x = [0., 1., 2., 3., 4., 5.]

# Altitude of the marble
y = [100., 98.26, 93.56, 81.79, 71.25, 53.22]

# Fit using polyfit
α = polyfit(x, y, 2)

# Evalutate at X
X = LinRange(0, 5, 200)
Y = polyval(α, X)

plot(X, Y, label="My approximation")
scatter!(x, y, label="Data")

# Modify the planet, in lowercase English, then remove the exception
planet = "earth"
### BEGIN SOLUTION
planet = "mars"
### END SOLUTION

"mars"

@assert planet != "earth"
### BEGIN HIDDEN TESTS
@assert planet == "mars" || planet == "mercury"
### END HIDDEN TESTS

Approximation par moindres carrés avec `Polynomials.fit` (cliquer pour afficher)

Nous avons vu dans le cours précédent que Polynomials.fit peut être utilisée pour l'interpolation polynomiale. Cette fonction peut également être utilisée pour ajuster des données en minimisant la somme des résidus au carré, ce qui peut être fait comme suit :

# This returns structure of type `Polynomial`, with associated degree 2
p = fit(x, y, 2)
@show p

# The coefficients can be obtained as follows
@show p.coeffs

# The structure `p` behaves like a function
@show p(0)

X = LinRange(0, 5, 200)
plot(X, p.(X), label="Polynomials.jl approximation")

p = Polynomial(100.11035714285717 + 0.02374999999997361*x - 1.8716071428571388*x^2)
p.coeffs = [100.11035714285717, 0.02374999999997361, -1.8716071428571388]
p(0) = 100.11035714285717

Méthodes directes

[Exercice 2] Décomposition LU sans élimination de Gauss

L'objectif de cet exercice est de proposer un algorithme permettant de réaliser la décomposition LU d'une matrice réelle $\mathsf{A}\in\mathbb{R}^{n×n}$, non pas par élimination gaussienne mais par identification des entrées de $\mathsf A$ avec celles de $\mathsf L \mathsf U$. Il s'agit de trouver un matrice triangulaire inférieure $\mathsf L$ formée de 1 sur la diagonale et une matrice triangulaire supérieure $\mathsf U$ telles que : $$ \tag{LU} \mathsf{A}=\mathsf{L}\mathsf{U} $$

1. Écrire une fonction my_lu(A) qui prend comme argument une matrice A et qui renvoie les matrices L et U. Pour calculer ces matrices, s'appuyer sur une identification successive des éléments des deux membres de (LU), ligne par ligne de haut en bas, et de gauche à droite au sein de chaque ligne.

   Hint (click to display)

   Lorsqu'on suit l'ordre conseillé, la comparaison de l'élément $(i, j)$ fournit une équation pour $\ell_{ij}$ si $j < i$, et une équation pour $u_{ij}$ si $j \geq i$. Notons qu'il est possible de parcourir les éléments dans d'autres ordres.

function my_lu(A)
    ### BEGIN SOLUTION
    n = size(A, 1)
    L, U = zeros(n, n), zeros(n, n)
    for i in 1:n, j in 1:i
        U[j, i] = A[j, i] - L[j, 1:j-1]'U[1:j-1, i]
        L[i, j] = (A[i, j] - L[i, 1:j-1]'U[1:j-1, j]) / U[j, j]
    end
    ### END SOLUTION
    return L, U
end;

ref_lu(A) = LinearAlgebra.lu(A, NoPivot())
@assert my_lu(diagm([1; 2; 3])) == (diagm([1; 1; 1]), diagm([1; 2; 3]))
@assert my_lu([2 -1 0; -1 2 -1; 0 -1 2])[1] ≈ [1 0 0; -1/2 1 0; 0 -2/3 1]
@assert my_lu([2 -1 0; -1 2 -1; 0 -1 2])[2] ≈ [2 -1 0; 0 3/2 -1; 0 0 4/3]
@assert (C = [1 2 3 4; 4 3 2 1; 1 2 1 2; 1 5 4 1]; my_lu(C)[1] ≈ ref_lu(C).L)
@assert (C = [1 2 3 4; 4 3 2 1; 1 2 1 2; 1 5 4 1]; my_lu(C)[2] ≈ ref_lu(C).U)
@assert (C = randn(100, 100); my_lu(C)[1] ≈ ref_lu(C).L)
@assert (C = randn(100, 100); my_lu(C)[2] ≈ ref_lu(C).U)

2. On suppose maintenant que la matrice réelle définie positive A est à largeur de bande b supposée beaucoup plus petite que n. Réécrire la fonction de décomposition LU en exploitant la largeur de bande.

   Hint (click to display)

   Pour rappel, la largeur de bande est le plus petit nombre naturel $b$ tel que $a_{ij} = 0$ pour tout $i,j \in \{1, \dots, n\}$ tels que $|i-j| > b$.

function my_banded_lu(A, b)
    ### BEGIN SOLUTION
    n = size(A, 1)
    L, U = zeros(n, n), zeros(n, n)
    for i in 1:n, j in max(1, i-b):i
        U[j, i] = A[j, i] - L[j, max(j-b, 1):j-1]'U[max(j-b, 1):j-1, i]
        L[i, j] = (A[i, j] - L[i, max(j-b, 1):j-1]'U[max(j-b, 1):j-1, j]) / U[j, j]
    end
    ### END SOLUTION
    return L, U
end;

@assert begin C = randn(100, 100); my_banded_lu(C, 99)[1] ≈ ref_lu(C).L end
@assert begin C = randn(100, 100); my_banded_lu(C, 99)[2] ≈ ref_lu(C).U end

### BEGIN HIDDEN TESTS
_generate_banded(n, b) = [abs(i - j) ≤ b ? randn() : 0. for i in 1:n, j in 1:n]
@assert (C = _generate_banded(10, 2); my_banded_lu(C, 2)[1] ≈ ref_lu(C).L)
@assert (C = _generate_banded(10, 2); my_banded_lu(C, 2)[2] ≈ ref_lu(C).U)
@assert (C = _generate_banded(10, 3); my_banded_lu(C, 3)[1] ≈ ref_lu(C).L)
@assert (C = _generate_banded(10, 3); my_banded_lu(C, 3)[2] ≈ ref_lu(C).U)
@assert (C = _generate_banded(1000, 1); my_banded_lu(C, 1)[1] ≈ ref_lu(C).L)
@assert (C = _generate_banded(1000, 1); my_banded_lu(C, 1)[2] ≈ ref_lu(C).U)
### END HIDDEN TESTS

3. Construire une fonction generate_banded(n, b) permettant de générer une matrice carrée aléatoire de taille n à largeur de bande donnée b.

function generate_banded(n, b)
    ### BEGIN SOLUTION
    A = zeros(n, n)
    for i in 1:n
        for j in max(1, i-b):min(n, i+b)
            A[i, j] = randn()
        end
    end
    return A
    ### END SOLUTION
end;

@assert generate_banded(10, 2)[1, 5] == 0
@assert generate_banded(10, 2)[2, 5] == 0
@assert generate_banded(10, 2)[3, 5] != 0
@assert generate_banded(10, 2)[4, 5] != 0
@assert generate_banded(10, 2)[5, 5] != 0
@assert generate_banded(10, 2)[6, 5] != 0
@assert generate_banded(10, 2)[7, 5] != 0
@assert generate_banded(10, 2)[8, 5] == 0
@assert generate_banded(10, 2)[9, 5] == 0

4. En utilisant generate_banded, tester votre implémentation de my_banded_lu, pour n = 100 et des valeurs de b égales à 2, 3 et 10.

   Hint (click to display)

   Vous pouvez utiliser la fonction lu de la bibliothèque LinearAlgebra, avec l'argument NoPivot(), correspondant à la fonnction ref_lu définie ci-dessus, comme fonction de référence. Vous pouvez également utiliser la macro @assert pour vos tests.

### BEGIN SOLUTION
# See hidden tests above
### END SOLUTION

Méthodes itératives de base

[Exercice 3] Itération de Richardson

Considérer le système linéaire suivant: $$ \mathsf A \mathbf x := \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 11 \\ 9 \end{pmatrix} =: \mathbf b. $$

1. Illustrer à l'aide de la fonction Plots.contourf les lignes de niveau de la fonction

$$ f(\mathbf x) = \frac{1}{2} \mathbf x^T \mathsf A \mathbf x - \mathbf b^T \mathbf x. $$

A = [3. 1.; 1. 3.]
b = [11.; 9.]
sol = A\b

### BEGIN SOLUTION
ℓ, nplot = 4, 100
xplot = sol[1] .+ LinRange(-ℓ, ℓ, nplot)
yplot = sol[2] .+ LinRange(-ℓ, ℓ, nplot)

plot(aspect_ratio=:equal)
plot!(xlims=xplot[[1,end]], ylims=yplot[[1,end]], legend=:none)
contourf!(xplot, yplot, (x₁, x₂) -> [x₁; x₂]'A*[x₁; x₂]/2 - b'*[x₁; x₂], c=:viridis)
### END SOLUTION

2. Implémenter l'itération de Richardson avec $\omega = 0.1$ pour résoudre le système. Votre fonction devra renvoyer une liste contenant toutes les itérations. Initialiser l'algorithme à $\mathbf x = 0$ et, comme critère d'arrêt, utiliser $$ \lVert \mathsf A \mathbf x - \mathbf b \lVert \leq \varepsilon \lVert \mathbf b \lVert, \qquad \varepsilon = 10^{-50} $$

   Hint (click to display)

   Pour ajouter un élément à la fin d'une liste, utiliser la fonction push!:

   push!(xs, x)  # Adds x at the end of xs
   

function richardson(ω)
    ε = 1e-50
    x = zeros(BigFloat, 2)
    xs = [x]
    # BEGIN SOLUTION
    while norm(A*x - b)/norm(b) ≥ ε
        x -= ω*(A*x - b)
        push!(xs, x)
    end
    # END SOLUTION
    return xs
end

ω = .1
xs = richardson(ω)
scatter!(eachrow(hcat(xs...))...)
plot!(eachrow(hcat(xs...))...)

3. Faire un plot de la norme du résidu $\lVert r_k \rVert := \lVert \mathsf A \mathbf x^{(k)} - \mathbf b \rVert$ en fonction de $k$, en utilisant une échelle linéaire pour l'axe des abcisses et une échelle logarithmique pour l'axe des ordonnées, gràce à l'argument yscale=:log passé à la fonction Plots.plot.

### BEGIN SOLUTION
errors = [norm(A*x - b) for x in xs]
plot(0:length(xs) - 1, errors, label="Error")
plot!(yscale=:log, xlabel=L"k")
### END SOLUTION

4. En utilisant Polynomials.fit, calculer une approximation du type $$ \log(r_k) \approx α + βk \qquad \Leftrightarrow \qquad r_k \approx \exp(α) \times \exp(β)^k. $$ Comparer la valeur de $\exp(β)$ au rayon spectral $\rho(\mathsf I - \omega \mathsf A)$.

# Define β and ρ below
### BEGIN SOLUTION
p = fit(1:length(errors), log.(errors), 1)
β = p[1]  # p[1] == p.coeffs[2]
λs = eigvals(A)
ρ = maximum(abs.(1 .- ω*λs))
### END SOLUTION

exp(β), ρ

(0.7997825113938985902912009059015771785038565910611528260786746665011047405987195, 0.8)

5. Calculer le paramètre $\omega$ optimal et refaire le plot de la décroissance de la norme du résidu dans ce cas.

### BEGIN SOLUTION
ω = 2 / (minimum(λs) + maximum(λs))
xs = richardson(ω)
errors = [norm(A*x - b) for x in xs]
plot(yscale=:log, xlabel=L"k")
plot!(0:length(xs) - 1, errors, label="Error")
### END SOLUTION

La méthode des gradients conjugués

On considère un système linéaire de la forme $$ \tag{*} \mathsf A \mathbf x = \mathbf b, $$ où $\mathsf A \in \mathbb R^{n\times n}$ est une matrice symétrique définie positive. On note $\mathbf x_*$ la solution exacte du système. L'objectif de cet exercice est de dériver théoriquement puis implémenter la méthode des gradients conjugués. Afin de se concentrer sur la seconde partie, les solutions de la première partie sont données.

[Exercice 4a] Dérivation théorique de la méthode

On suppose qu'un vecteur $\mathbf x_0$ est donné et on définit, pour $k \in \mathbb N$, l'espace de Krylov $$ \mathcal K_k := \text{Span} \Bigl\{ \mathbf r_0, \mathsf A \mathbf r_0, \mathsf A^2 \mathbf r_0, \dotsc, \mathsf A^{k-1} \mathbf r_0 \Bigr\}, \qquad \mathbf r_0 := \mathsf A \mathbf x_0 - \mathbf b. $$ Soit $\mathbf x_k$ le minimiseur de $f$ sur l'espace affine $\mathbf x_0 + \mathcal K_k$, où $f$ est la fonction $$ f(\mathbf x) := \frac{1}{2} \mathbf x^T \mathsf A \mathbf x - \mathbf b^T \mathbf x. $$ Rappelons que le vecteur $\mathbf x_*$ est minimiseur global de cette fonction strictement convexe. Le but de cet exercice est de motiver la méthode des gradients conjugés, qui peut être vue comme une formule de récurence simple pour construire les $\mathbf x_k$ de manière itérative. Pour ce faire, nous allons suivre les étapes suivantes :

1. Montrer qu'il existe $k_* \in \{0, \dotsc, n\}$ tel que $\mathbf x_k = \mathbf x_*$ pour tout $k \geq k_*$.

   Solution (click to display)

   • La suite des sous-espaces $\{\mathcal K_k\}$ est croissante et ces sous-espaces sont de dimension au plus $n$, donc il existe $k_* \in \{0, \dotsc, n\}$ tel que $\mathcal K_{k_* + 1} = \mathcal K_{k_*}$. En particulier $\mathsf A \mathcal K_{k_*} \subset \mathcal K_{k_*}$.

   • Comme $\mathsf A$ est inversible, les deux sous-espaces $\mathsf A \mathcal K_{k_*}$ et $\mathcal K_{k_*}$ sont de même dimension. Par conséquent, l'inclusion $\mathsf A \mathcal K_{k_*} \subset \mathcal K_{k_*}$ implique qu'en fait $\mathsf A \mathcal K_{k_*} = \mathcal K_{k_*}$.

   • En particulier $\mathbf r_0 \in \mathsf A \mathcal K_{k_*}$, donc il existe $\mathbf z_* \in \mathcal K_{k_*}$ tel que $\mathsf A \mathbf z_* = -\mathbf r_0$.

   • On conclut car alors $\mathsf A (\mathbf x_0 + \mathbf z_*) = \mathbf b$, et donc $\mathbf x_* = \mathbf x_0 + \mathbf z_*$ appartient à l'espace affine $\mathbf x_0 + \mathcal K_{k_*}$.

   
   

2. Montrer que $f(\mathbf x) = \frac{1}{2} \lVert \mathbf x - \mathbf x_* \lVert_{\mathsf A}^2 + f(\mathbf x_*)$, où $\langle \mathbf x, \mathbf y \rangle_{\mathsf A} := \mathbf x^T \mathsf A \mathbf y$ et $\lVert \cdot \rVert_{\mathsf A}$ est la norme associée.

   Remarque (click to display)

   Ce point implique que $\mathbf x_k - \mathbf x_0$ est la projection $\mathsf A$-orthogonale de $\mathbf x_* - \mathbf x_0$ sur $\mathcal K_k$.

   
   

3. Soit $\mathbf r_k := \mathsf A \mathbf x_k - \mathbf b.$ Montrer que $\mathbf r_k = \mathbf r_0 + \mathsf A (\mathbf x_k - \mathbf x_0)$, et en déduire que $\mathbf r_k \in \mathcal K_{k+1}$.

4. Montrer que $\langle \mathbf r_{k}, \mathbf v \rangle = 0$ pour tout $\mathbf v \in \mathcal K_{k}.$

   Solution (click to display)

   Noter que $$ f(\mathbf x_k + \varepsilon \mathbf v) = \frac{1}{2} \lVert \mathbf x_k - \mathbf x_* \lVert_{\mathsf A}^2 + \varepsilon \langle \mathbf x_k - \mathbf x_*, \mathbf v \rangle_{\mathsf A} + \frac{\varepsilon^2}{2} \lVert \mathbf v \lVert_{\mathsf A}^2 + f(\mathbf x_*). $$ Par définition de $\mathbf x_k$, cette expression est minimale pour $\varepsilon = 0$. En imposant que la dérivée par rapport à $\varepsilon$ du membre de droite soit égale à 0 pour $\varepsilon = 0$, on obtient $\langle \mathbf x_k - \mathbf x_*, \mathbf v \rangle_{\mathsf A} = 0,$ ce qui permet de conclure car $\langle \mathbf x_k - \mathbf x_*, \mathbf v \rangle_{\mathsf A} = \langle \mathbf r_{k}, \mathbf v \rangle$.

   
   

5. Déduire des deux points précédents que si $\mathbf r_k \neq \mathbf 0$, alors $\mathbf r_k \in \mathcal K_{k+1} \setminus \mathcal K_k$.

6. Déduire aussi que $\langle \mathbf r_{k}, \mathbf v \rangle_{\mathsf A} = 0$ pour tout $\mathbf v \in \mathcal K_{k-1}$, pour $k \geq 1$.

   Solution (click to display)

   C'est clair car $\langle \mathbf r_{k}, \mathbf v \rangle_{\mathsf A} = \langle \mathbf r_{k}, \mathsf A \mathbf v \rangle$, et $\mathsf A \mathbf v \in \mathcal K_{k}$ pour tout $\mathbf v \in \mathcal K_{k-1}$, par définition des espaces de Krylov $\{\mathcal K_k\}$.

   
   

7. Conclure que, pour tout $k$ tel que $\mathbf r_k \neq \mathbf 0$, les vecteurs $(\mathbf d_0, \dotsc, \mathbf d_k)$ définis par la relation de récurrence $$ \tag{CG1} \mathbf d_{\ell} = \mathbf r_{\ell} - \frac{\mathbf r_{\ell}^T \mathsf A \mathbf d_{\ell-1}}{\mathbf d_{\ell-1}^T \mathsf A \mathbf d_{\ell-1}} \mathbf d_{\ell-1}, \qquad \mathbf d_0 = \mathbf r_0 $$ forment une base $\mathsf A$-orthogonale de $\mathcal K_{k+1}$.

   Solution (click to display)

   La propriété est trivialement vérifiée pour $k = 0$. Prenons maintenant $k > 0$ et supposons par récurrence que $(\mathbf d_0, \dotsc, \mathbf d_{k-1})$ forment une base $\mathsf A$-orthogonale de $\mathcal K_k$. Par construction $\left\langle \mathbf d_k, \mathbf d_{k-1} \right\rangle_{\mathsf A} = 0$. On a aussi que $\left\langle \mathbf d_k, \mathbf d_{\ell} \right\rangle_{\mathsf A} = 0$ pour tout $\ell \in \{0, \dotsc, k-2\}$, car on a montré précédemment que $\langle \mathbf r_k, \mathbf v \rangle_{\mathsf A} = 0$ pour tout $\mathbf v \in \mathcal K_{k-1}$.

   
   

8. Montrer que $\mathbf x_{k+1} - \mathbf x_{k}$ est $\mathsf A$-orthogonal à $\mathcal K_k$, d'où $\mathbf x_{k+1} = \mathbf x_{k} - \omega_k \mathbf d_k$ pour un $\omega_k \in \mathbb R$ quand $\mathbf r_k \neq 0$.

   Solution (click to display)

   Noter que $\langle \mathbf x_{k+1} - \mathbf x_k, \mathbf v \rangle_{\mathsf A} = \langle \mathbf r_{k+1} - \mathbf r_k, \mathbf v \rangle$, et utiliser un des résultats prouvés ci-dessus.

   
   

9. Pour terminer, montrer que $$ \tag{CG2} \mathbf x_{k+1} = \mathbf x_k - \frac{\mathbf r_k^T \mathbf d_k}{\mathbf d_k^T \mathsf A \mathbf d_k} \mathbf d_k. $$

   Solution (click to display)

   Il suffit de trouver le paramètre $\omega_k$ du point précédent par minimisation de $f(\mathbf x_k - \omega_k \mathbf d_k)$.

[Exercice 4b] Implémentation de la méthode

Les relations (CG1) et (CG2) permettent de construire les vecteurs $(\mathbf x_k)$ de manière récursive, à un coût par itération indépendent de $k$; c'est la méthode des gradients conjugués.

• Sur base de ces équations, implémenter une fonction Julia pour résoudre le système (*). Votre fonction devra renvoyer un vecteur de vecteur (de type Vector{Vector{Float64}}) contenant toutes les itérations $(\mathbf x_0, \mathbf x_1, \dotsc, \mathbf x_K)$. Utiliser comme vecteur initial $\mathbf x_0$ un vecteur de 1, qui pourra être construit en utilisant la fonction ones. Comme critère d'arrêt, utiliser la condition $$ \lVert \mathsf A \mathbf x_k - \mathbf b \rVert \leq ε \lVert \mathbf b \rVert. $$

function conjugate_gradients(A, b; ε=1e-10)
    ### BEGIN SOLUTION
    n = length(b)
    x = ones(n)
    d = r = A*x - b
    xs = [x]
    while √(r'r) ≥ ε * √(b'b)
        ω = d'r / (d'A*d)
        x -= ω*d
        r = A*x - b
        β = r'A*d/(d'A*d)
        d = r - β*d
        push!(xs, x)
    end
    xs
    ### END SOLUTION
end;

@assert typeof(conjugate_gradients([1. 0.; 0. 1.], [2.; 3.])) == Vector{Vector{Float64}}
@assert conjugate_gradients([1. 0.; 0. 1.], [2.; 3.])[end] ≈ [2.; 3.]
@assert conjugate_gradients([2. 1.; 1. 2.], [3.; 3.])[end] ≈ [1.; 1.]
@assert conjugate_gradients([4. 1. 0.; 1. 4. 1.; 0. 1. 4.], [4.; 1.; 0.])[end] ≈ [1.; 0.; 0.]

• Tester votre implémentation sur la matrice $\mathsf A$ et le vecteur $\mathbf b$ du code ci-dessous, qui correspondent à une discrétisation par différences finies de l'équation de Poisson sur le domaine $\Omega = (0, 2) \times (0, 1)$, avec condition de Dirichlet homogène à la frontière: $$ \begin{aligned} - \Delta u(x, y) &= b(x, y), \qquad (x, y) \in \Omega, \\ u(x, y) &= 0, \quad \qquad \quad (x, y) \in \partial \Omega. \end{aligned} $$ Le membre de droite qu'on prend est $$ b(x, y) = \sin(4\pi x) + \sin(2\pi y). $$

using LaTeXStrings
using LinearAlgebra
using Plots
import SparseArrays

# Domain size
Lx, Ly = 2, 1

# Number of discretization points along x and y, including the boundary points
nx, ny = 101, 101

function discretize(nx, ny)
    hx, hy = Lx/(nx - 1), Ly/(ny - 1)
    Dxx = (1/hx^2) * SymTridiagonal(2ones(nx-2), -ones(nx-3))
    Dyy = (1/hy^2) * SymTridiagonal(2ones(ny-2), -ones(ny-3))
    A = kron(Dxx, I(ny-2)) + kron(I(nx-2), Dyy)
    xgrid = Lx/(nx-1) * (1:nx-2)
    ygrid = Ly/(ny-1) * (1:ny-2)
    x_2d = reshape([x for y in ygrid, x in xgrid], (nx-2)*(ny-2))
    y_2d = reshape([y for y in ygrid, x in xgrid], (nx-2)*(ny-2))
    b = sin.(4π*x_2d) + sin.(2π*y_2d)
    return SparseArrays.SparseMatrixCSC(A), b
end

function plot_solution(f)
    f = reshape(f, ny-2, nx-2)
    z = [zeros(nx)'; zeros(ny-2) f zeros(ny-2); zeros(nx)']  # Add boundary
    xgrid = Lx/(nx-1) * (0:nx-1)
    ygrid = Ly/(ny-1) * (0:ny-1)
    heatmap(xgrid, ygrid, z, c=:viridis, levels=50)
end

A, b = discretize(nx, ny)
xs = conjugate_gradients(A, b)
x₊ = A\b  # Exact solution
plot_solution(x₊)
plot_solution(xs[end])

• Faire un graphe illustrant l'évolution du résidu $\lVert \mathbf r_k \rVert$, de l'erreur $\lVert \mathbf x_k - \mathbf x_* \rVert$, et de l'erreur $\lVert \mathbf x_k - \mathbf x_* \rVert_{\mathsf A}$ en fonction de $k$. Choisissez des échelles appropriées pour le graphe (linéaire-linéaire, log-log, linéaire-log ou log-linéaire).

### BEGIN SOLUTION
rs = [A*x - b for x in xs]
es = [x - x₊ for x in xs]
plot(xlabel=L"k", yaxis=:log)
plot!([√(r'r) for r in rs], label=L"\Vert Ax - b \Vert")
plot!([√(e'r) for (e, r) in zip(es, rs)], label=L"\Vert x - x_{*} \Vert_A")
plot!([√(e'e) for e in es], label=L"\Vert x - x_{*} \Vert")
### END SOLUTION