Remonter un système triangulaire, sans al.solve

Ce qu'on montre. Que la dernière équation d'un système triangulaire ne contient qu'une inconnue, et qu'on remonte ensuite ligne par ligne — c'est l'étape finale de toute résolution par le pivot.

Le principe. La $i$-ième ligne de $UX=B$ s'écrit $$u_{i,i}{x}_i+\sum _{j>i}{u}_{i,j}{x}_j=b_i.$$ Si l'on connaît déjà $x_{i+1},\dots ,x_{n-1}$, alors $x_i$ s'en déduit : $$x_i=\frac{1}{u_{i,i}}\left(b_i-\sum _{j>i}{u}_{i,j}{x}_j\right),$$ division licite puisque $u_{i,i}\ne 0$. On calcule donc les inconnues de la dernière vers la première.

Le programme.

import numpy as np
import numpy.linalg as al

def remontee(U, B):
    """Resout U X = B, U triangulaire superieure a diagonale non nulle."""
    n = np.shape(U)[0]
    X = np.zeros(n)
    for i in range(n - 1, -1, -1):     # de la derniere ligne a la premiere
        s = B[i]
        for j in range(i + 1, n):      # ce qui est deja connu
            s = s - U[i, j] * X[j]
        X[i] = s / U[i, i]
    return X

U = np.array([[2., 1., -1.],
              [0., 3.,  2.],
              [0., 0.,  4.]])
B = np.array([5., 8., 12.])
print(remontee(U, B))
print(al.solve(U, B))

La correction. L'invariant est : avant le tour d'indice $i$, les cases $X[i+1],\dots ,X[n-1]$ contiennent les valeurs correctes. Il est vrai avant le premier tour ($i=n-1$, aucune case à garantir). S'il est vrai avant le tour $i$, la boucle interne calcule exactement $b_i-{\sum }_{j>i}{u}_{i,j}{x}_j$, et la division par $u_{i,i}$ donne la valeur correcte de $x_i$ : l'invariant vaut avant le tour $i-1$. Après le tour $i=0$, toutes les cases sont correctes.

Le coût. Le tour $i$ effectue $n-1-i$ multiplications et une division. Au total $$\sum _{i=0}^{n-1}(n-i)=\frac{n(n+1)}{2}\text{ opeˊrations, soit de l’ordre de }\frac{n^2}{2}.$$ C'est bien moins que la résolution d'un système quelconque, qui exige d'abord la triangularisation.

Le contrôle à la main. La troisième ligne donne $4x_3=12$, donc $x_3=3$. La deuxième donne $3x_2+2\times 3=8$, donc $x_2=\frac{2}{3}$. La première donne $2x_1+\frac{2}{3}-3=5$, donc $2x_1=\frac{22}{3}$ et $x_1=\frac{11}{3}$. Le programme affiche $$X\approx (3,6667,0,6667,3),$$ et al.solve affiche la même chose.

Vérification directe. $2\times \frac{11}{3}+\frac{2}{3}-3=\frac{22+2-9}{3}=5$ ; $3\times \frac{2}{3}+2\times 3=8$ ; $4\times 3=12$. Les trois équations sont satisfaites.

Point délicat. range(n - 1, -1, -1) parcourt bien $n-1,n-2,\dots ,0$ : la borne $-1$ est exclue. Écrire range(n - 1, 0, -1) oublierait la première ligne, et la case X[0] resterait nulle — erreur silencieuse, car le programme ne signalerait rien.