Ingénierie numérique

a) False   0.1+0.2 = 0.30000000000000004
b) False   0.30000000000000004 egalement
c) True    absorption : 1 < eps x 1e16
d) 0.0     1e-16 < eps, donc absorbe
e) False   les ENTIERS Python sont exacts, sans limite
f) True    en flottant, 2^53 et 2^53+1 ont la meme representation
g) inf     depassement de capacite
h) -inf
i) nan     forme indeterminee 0 x inf
j) False   0.9999999999999999

Les points (e) et (f) sont instructifs : \(2^{53}\) est le plus grand entier tel que tous les entiers de \(0\) a \(2^{53}\) soient exactement representables en double precision. Au-dela, les entiers consecutifs ne le sont plus tous : \(2^{53}+1\) est arrondi a \(2^{53}\).

print(2**53, float(2**53 + 1))     # 9007199254740992 9007199254740992.0
print(float(2**53 + 2))            # 9007199254740994.0 : les PAIRS restent exacts

C'est pourquoi JavaScript, qui n'a que des flottants, limite les entiers surs a \(2^{53}-1\).

from math import sqrt, cos, sin, exp, expm1

def naif(x):    return sqrt(x*x + 1) - x
def stable(x):  return 1 / (sqrt(x*x + 1) + x)

print(f"{'x':>10}{'naif':>24}{'stable':>24}")
for k in range(11):
    x = 10.0**k
    print(f"{x:>10.0e}{naif(x):>24.16e}{stable(x):>24.16e}")

         x                    naif                  stable
     1e+00  4.1421356237309515e-01  4.1421356237309503e-01
     1e+03  4.9999987500006264e-04  4.9999987500006250e-04
     1e+06  4.9999999749991756e-07  4.9999999999999998e-07
     1e+08  0.0000000000000000e+00  4.9999999999999998e-09
     1e+10  0.0000000000000000e+00  4.9999999999999999e-11

2. Pour \(x\) grand, \(\sqrt{x^2+1} \approx x\) : on soustrait deux nombres dont les 16 chiffres significatifs coincident presque tous. Des \(x = 10^8\), le resultat exact \(\approx 5\cdot10^{-9}\) est inferieur a \(\varepsilon x = 2\cdot10^{-8}\) : il est completement noye. La reecriture utilise la quantite conjuguee : \[ \sqrt{x^2+1} - x = \frac{(x^2+1) - x^2}{\sqrt{x^2+1} + x} = \frac{1}{\sqrt{x^2+1}+x}, \] ou l'on n'effectue plus qu'une addition de nombres de meme signe.

# 3. 1 - cos(x)
def cos_naif(x):   return 1 - cos(x)
def cos_stable(x): return 2 * sin(x/2)**2        # identite de l'angle moitie

for k in range(1, 9):
    x = 10.0**(-k)
    print(f"x = 1e-{k} : naif {cos_naif(x):.6e}, stable {cos_stable(x):.6e}, "
          f"exact ~ {x*x/2:.6e}")

x = 1e-1 : naif 4.995835e-03, stable 4.995835e-03, exact ~ 5.000000e-03
x = 1e-5 : naif 5.000000e-11, stable 5.000000e-11, exact ~ 5.000000e-11
x = 1e-8 : naif 0.000000e+00, stable 5.000000e-17, exact ~ 5.000000e-17

Identite utilisee : \(1 - \cos x = 2\sin^2(x/2)\). La forme naive perd tout des \(x < 10^{-8}\), car \(\cos x\) vaut alors exactement 1 en machine.

# 4. exp(x) - 1
for k in range(1, 18, 4):
    x = 10.0**(-k)
    print(f"x = 1e-{k:>2} : naif {exp(x)-1:.10e}, expm1 {expm1(x):.10e}")

La bibliotheque standard fournit math.expm1(x) et math.log1p(x) precisement pour ces situations : elles calculent \(e^x-1\) et \(\ln(1+x)\) par des developpements adaptes, sans jamais former la difference dangereuse.

from math import fsum
from fractions import Fraction

# La reference exacte est calculee avec des Fraction : les denominateurs
# grossissent tres vite, on se limite donc a 10^4 pour la valeur de reference.
n = 10**4

def croissant(n):
    s = 0.0
    for k in range(1, n+1):
        s += 1/k
    return s

def decroissant(n):
    s = 0.0
    for k in range(n, 0, -1):
        s += 1/k
    return s

def kahan(n):
    """Sommation compensee : on reinjecte l'erreur d'arrondi de chaque addition."""
    s, c = 0.0, 0.0
    for k in range(1, n+1):
        y = 1/k - c
        t = s + y
        c = (t - s) - y        # ce qui a ete perdu
        s = t
    return s

exact = float(sum(Fraction(1, k) for k in range(1, n+1)))

for nom, v in [("croissant", croissant(n)), ("decroissant", decroissant(n)),
               ("Kahan", kahan(n)), ("fsum", fsum(1/k for k in range(1, n+1)))]:
    print(f"{nom:<12}{v:.15f}  erreur relative {abs(v-exact)/exact:.2e}")
print(f"{'exact':<12}{exact:.15f}")

croissant    9.787606036044348  erreur relative 1.09e-15
decroissant  9.787606036044386  erreur relative 2.72e-16
Kahan        9.787606036044382  erreur relative 1.81e-16
fsum         9.787606036044384  erreur relative 0.00e+00
exact        9.787606036044384

Seul fsum retrouve exactement le flottant le plus proche du resultat reel. En sommant dans l'ordre decroissant on gagne un facteur 4 sur l'erreur, et Kahan un facteur 6, sans changer la complexite.

4. Justification. A chaque addition \(s \leftarrow s + x\), l'erreur commise est majoree par \(\varepsilon|s|\). L'erreur totale est donc majoree par \(\varepsilon\sum_{k}|s_k|\) ou \(s_k\) est la somme partielle apres \(k\) termes.

• Ordre croissant : \(s_k\) croit vite (comme \(\ln k\)) et devient grand devant les termes ajoutes. Quand \(1/k < \varepsilon s_k\), le terme est purement absorbe. Pour \(s \approx 14\), cela arrive des que \(1/k < 3\cdot10^{-15}\), soit \(k > 3\cdot10^{14}\) : ici on n'y arrive pas, mais l'erreur s'accumule.
• Ordre decroissant : les petits termes sont additionnes entre eux d'abord et forment une contribution qui n'est plus negligeable quand elle rencontre les grands termes.

La sommation de Kahan conserve dans c la partie perdue a chaque arrondi et la reinjecte : l'erreur totale passe de \(O(n\varepsilon)\) a \(O(\varepsilon)\), independamment de \(n\). math.fsum va plus loin : il maintient une liste de sommes partielles exactes et renvoie le flottant exactement le plus proche du resultat reel.

import numpy as np
from math import exp, e

def ordre(erreurs, hs):
    """Pente de ln(E) en fonction de ln(h)."""
    h, E = np.array(hs, float), np.array(erreurs, float)
    ok = E > 1e-15                       # on ecarte la zone dominee par l'arrondi
    return float(np.polyfit(np.log(h[ok]), np.log(E[ok]), 1)[0])

# 2. differences finies
x0, exact = 1.0, exp(1.0)
hs = [10.0**-k for k in range(1, 7)]
d1 = [abs((exp(x0+h) - exp(x0))/h - exact) for h in hs]
d2 = [abs((exp(x0+h) - exp(x0-h))/(2*h) - exact) for h in hs]
print(f"decentree : {ordre(d1, hs):.3f}   (theorie 1)")
print(f"centree   : {ordre(d2, hs):.3f}   (theorie 2)")

# 3. integration de exp sur [0,1], valeur exacte e - 1
I = e - 1

def rectangles_gauche(n):
    h = 1/n
    return h * sum(exp(k*h) for k in range(n))

def trapezes(n):
    h = 1/n
    return h * (exp(0)/2 + sum(exp(k*h) for k in range(1, n)) + exp(1)/2)

def simpson(n):
    """n doit etre PAIR."""
    h = 1/n
    s = exp(0) + exp(1)
    s += 4*sum(exp((2*k-1)*h) for k in range(1, n//2 + 1))
    s += 2*sum(exp(2*k*h) for k in range(1, n//2))
    return h*s/3

ns = [2**k for k in range(2, 11)]
hs = [1/n for n in ns]
for nom, f, theorie in [("rectangles", rectangles_gauche, 1),
                        ("trapezes", trapezes, 2),
                        ("Simpson", simpson, 4)]:
    E = [abs(f(n) - I) for n in ns]
    print(f"{nom:<12} ordre mesure {ordre(E, hs):.3f}   (theorie {theorie})")

decentree : 1.000   (theorie 1)
centree   : 2.000   (theorie 2)
rectangles   ordre mesure 1.000   (theorie 1)
trapezes     ordre mesure 2.000   (theorie 2)
Simpson      ordre mesure 4.000   (theorie 4)

1-2. Le conditionnement relatif est \(\kappa(x) = \left|\frac{x f'(x)}{f(x)}\right|\).

import numpy as np

def wilkinson(n=20, perturbation=0.0):
    """Racines du polynome (x-1)...(x-n) developpe, avec perturbation du coefficient de x^(n-1)."""
    coef = np.poly(np.arange(1, n+1))          # coefficients, poids forts en tete
    coef[1] += perturbation                     # perturbation du coefficient de x^(n-1)
    return np.sort_complex(np.roots(coef))

n = 20
r0 = wilkinson(n, 0.0)
print("sans perturbation, erreur max :", np.abs(r0 - np.arange(1, n+1)).max())

for p in (1e-10, 1e-8, 2**-23):
    r = wilkinson(n, p)
    err = np.abs(np.sort_complex(r) - np.arange(1, n+1)).max()
    complexes = int((np.abs(r.imag) > 1e-6).sum())
    print(f"perturbation {p:.2e} : erreur max {err:.3f}, "
          f"{complexes} racines devenues complexes")

sans perturbation, erreur max : 0.0208
perturbation 1.00e-10 : erreur max 0.316, 0 racines devenues complexes
perturbation 1.00e-08 : erreur max 2.104, 10 racines devenues complexes
perturbation 1.19e-07 : erreur max 3.512, 10 racines devenues complexes

3-4. Le polynome de Wilkinson. Le polynome \(P(x) = \prod_{k=1}^{20}(x-k)\) a pour racines les entiers 1 a 20, parfaitement separees. Sous forme developpee, son coefficient de \(x^{19}\) vaut \(-210\). Wilkinson (1963) a montre qu'une perturbation de ce seul coefficient de \(2^{-23} \approx 1{,}2\cdot10^{-7}\) — soit une perturbation relative de \(5\cdot10^{-10}\) — deplace la racine 16 de plus de 2,8 unites et transforme dix racines reelles en cinq paires complexes conjuguees.

L'analyse de sensibilite le confirme : pour une racine simple \(r\), \[ \frac{\partial r}{\partial a_j} = -\frac{r^{j}}{P'(r)}, \] et pour \(r = 16\), \(|P'(16)| \approx 4\cdot10^{15}\) tandis que \(16^{19} \approx 7\cdot10^{22}\) : le facteur d'amplification atteint \(10^{7}\).

import numpy as np
from math import exp

EPS = np.finfo(float).eps      # 2.220446049250313e-16
x0, exact = 1.0, exp(1.0)

h = np.logspace(-16, -1, 400)
d1 = np.abs((np.exp(x0+h) - np.exp(x0))/h - exact)
d2 = np.abs((np.exp(x0+h) - np.exp(x0-h))/(2*h) - exact)

i1, i2 = d1.argmin(), d2.argmin()
print(f"decentree : h* = {h[i1]:.3e}, erreur min = {d1[i1]:.3e}")
print(f"centree   : h* = {h[i2]:.3e}, erreur min = {d2[i2]:.3e}")
print(f"theorie   : sqrt(eps) = {np.sqrt(EPS):.3e}, "
      f"eps^(1/3) = {EPS**(1/3):.3e}")

decentree : h* = 1.520e-08, erreur min = 6.634e-09
centree   : h* = 3.976e-06, erreur min = 8.036e-12
theorie   : sqrt(eps) = 1.490e-08, eps^(1/3) = 6.055e-06

3. Accord avec la theorie. On minimise \(E(h) = C h^p + \dfrac{2\varepsilon |f(x_0)|}{h}\). L'annulation de la derivee donne \[ h^* = \left(\frac{2\varepsilon |f|}{pC}\right)^{\frac{1}{p+1}}, \qquad E(h^*) = \Theta\!\left(\varepsilon^{\frac{p}{p+1}}\right). \]

Les valeurs mesurees coincident a un facteur 1,5 pres, ce qui est excellent compte tenu des constantes negligees.

4. Comment faire mieux ? On ne peut pas descendre sous \(10^{-11}\) avec des differences finies. Trois alternatives :

• Extrapolation de Richardson : combiner les formules a \(h\) et \(h/2\) pour eliminer le terme dominant de l'erreur et gagner deux ordres. \(D(h) = \frac{4 D_c(h/2) - D_c(h)}{3}\) est d'ordre 4.
• Differentiation automatique : propager simultanement la valeur et la derivee a travers chaque operation elementaire. La derivee est alors obtenue a la precision machine, sans aucun pas. C'est la technique qui rend l'apprentissage profond possible.
• Methode du pas complexe : si \(f\) est analytique, \(f'(x) \approx \dfrac{\operatorname{Im} f(x + \mathrm i h)}{h}\), qui ne comporte aucune soustraction et donne la precision machine des \(h = 10^{-20}\).

import cmath
def derivee_pas_complexe(f, x, h=1e-20):
    """Sans cancellation : la partie imaginaire isole le terme en f'(x)."""
    return (f(complex(x, h))).imag / h

d = derivee_pas_complexe(cmath.exp, 1.0)
print(f"{d:.16f}, erreur {abs(d - exact):.2e}")    # erreur ~ 0, precision machine

from math import sqrt, log2

exact = sqrt(2)

def dichotomie(eps=1e-15):
    a, b, n, err = 1.0, 2.0, 0, []
    while b - a > eps:
        m = (a + b) / 2
        err.append(abs(m - exact))
        if m*m > 2: b = m
        else: a = m
        n += 1
    return n, err

def newton(eps=1e-15):
    x, n, err = 1.0, 0, []
    while abs(x*x - 2) > eps and n < 60:
        x = 0.5 * (x + 2/x)
        err.append(abs(x - exact))
        n += 1
    return n, err

nd, ed = dichotomie()
nn, en = newton()
print(f"dichotomie : {nd} iterations")     # 51
print(f"Newton     : {nn} iterations")     # 5

print("\ndichotomie e(n+1)/e(n) :")
print([f"{ed[i+1]/ed[i]:.3f}" for i in range(6)])
print("Newton e(n+1)/e(n)^2 :")
print([f"{en[i+1]/en[i]**2:.4f}" for i in range(3)])

dichotomie : 51 iterations
Newton     : 5 iterations

dichotomie e(n+1)/e(n) : ['0.379', '2.638', '0.379', '2.638', '0.379', '2.638']
Newton e(n+1)/e(n)^2 : ['0.3536', '0.3536', '0.3536']

2. Pour la dichotomie, le rapport oscille mais l'encadrement \(b_n - a_n = (b_0-a_0)/2^n\) est exactement divise par 2 : la convergence est lineaire de rapport \(1/2\), soit \(\log_{10} 2 \approx 0{,}3\) decimale gagnee par iteration. Il faut donc \(\log_2(1/10^{-15}) \approx 50\) iterations.

Pour Newton, le rapport \(e_{n+1}/e_n^2\) est constant et vaut 0,3536 : la convergence est quadratique, avec \[ e_{n+1} \approx \frac{|f''(\ell)|}{2|f'(\ell)|} e_n^2 = \frac{1}{2\sqrt2}\,e_n^2 \approx 0{,}354\,e_n^2, \] ce qui correspond exactement a la valeur mesuree.

def secante(f, x0, x1, eps=1e-15, nmax=50):
    err = []
    for _ in range(nmax):
        if f(x1) == f(x0): break
        x2 = x1 - f(x1)*(x1 - x0)/(f(x1) - f(x0))
        err.append(abs(x2 - exact))
        x0, x1 = x1, x2
        if abs(f(x1)) < eps: break
    return err

es = secante(lambda x: x*x - 2, 1.0, 2.0)
print("secante, ordre estime :")
for i in range(1, len(es)-1):
    if es[i] > 1e-15 and es[i-1] > 1e-15:
        print(f"  {log2(es[i+1]/es[i]) / log2(es[i]/es[i-1]):.3f}")

3. L'ordre estime vaut environ 1,618 : c'est le nombre d'or \(\varphi = \frac{1+\sqrt5}{2}\). L'erreur verifie \(e_{n+1} \approx K e_n e_{n-1}\), et en cherchant \(e_n \approx A e_{n-1}^\alpha\) on obtient \(\alpha^2 = \alpha + 1\).

La secante converge donc moins vite que Newton par iteration, mais elle ne demande qu'une evaluation de \(f\) par etape (contre \(f\) et \(f'\) pour Newton) : a cout egal, son ordre effectif est \(\varphi^2 \approx 2{,}6 > 2\). Elle est donc souvent preferable quand \(f'\) est couteuse ou inconnue.

# 4. racine double : f(x) = (x-1)^2
def newton_double():
    f  = lambda x: (x-1)**2
    fp = lambda x: 2*(x-1)
    x, err = 2.0, []
    for _ in range(12):
        x = x - f(x)/fp(x)
        err.append(abs(x - 1))
    return err

e = newton_double()
print([f"{v:.2e}" for v in e[:8]])
print("rapports e(n+1)/e(n) :", [f"{e[i+1]/e[i]:.3f}" for i in range(6)])

['5.00e-01', '2.50e-01', '1.25e-01', '6.25e-02', ...]
rapports e(n+1)/e(n) : ['0.500', '0.500', '0.500', '0.500', '0.500', '0.500']

4. Sur une racine double, Newton perd sa convergence quadratique et devient lineaire de rapport \(1/2\) — pas mieux que la dichotomie. La raison : \(f'(\ell) = 0\), et l'analyse quadratique suppose \(f'(\ell) \ne 0\). Pour une racine de multiplicite \(m\), le rapport est \(1 - 1/m\).

Remede : utiliser la methode de Newton modifiee \(x_{n+1} = x_n - m f/f'\), qui retrouve l'ordre 2, ou appliquer Newton a \(u = f/f'\) dont toutes les racines sont simples.

1. Construction. On ecrit

\[ A(h) = A + c_1 h^p + O(h^{p+1}), \qquad A(h/2) = A + c_1 \frac{h^p}{2^p} + O(h^{p+1}). \]

En eliminant \(c_1\) par combinaison lineaire :

\[ \boxed{\;B(h) = \frac{2^p A(h/2) - A(h)}{2^p - 1} = A + O(h^{p+1}).\;} \]

On a gagne un ordre sans aucune evaluation supplementaire au-dela de \(A(h/2)\).

2. Tableau de Romberg. Pour la methode des trapezes, le developpement d'Euler-Maclaurin ne contient que des puissances paires de \(h\) : \(T(h) = I + c_1h^2 + c_2h^4 + \cdots\). Chaque extrapolation gagne donc deux ordres. On pose \(R_{k,0} = T(h_k)\) avec \(h_k = (b-a)/2^k\), puis

\[ R_{k,j} = \frac{4^{j} R_{k,j-1} - R_{k-1,j-1}}{4^{j} - 1}. \]

import numpy as np
from math import exp, e, log2

def trapezes(f, a, b, n):
    x = np.linspace(a, b, n+1)
    y = np.array([f(t) for t in x])
    return (b-a)/n * (y[0]/2 + y[1:-1].sum() + y[-1]/2)

def romberg(f, a, b, kmax=8, verbeux=False):
    """Tableau de Romberg. Renvoie la meilleure estimation et le nombre d'evaluations."""
    R = [[0.0]*(kmax+1) for _ in range(kmax+1)]
    evaluations = 0
    for k in range(kmax+1):
        n = 2**k
        R[k][0] = trapezes(f, a, b, n)
        evaluations += n + 1
        for j in range(1, k+1):
            R[k][j] = (4**j * R[k][j-1] - R[k-1][j-1]) / (4**j - 1)
        if verbeux:
            print("  ".join(f"{R[k][j]:.14f}" for j in range(k+1)))
    return R[kmax][kmax], evaluations, R

I = e - 1
val, ev, R = romberg(exp, 0, 1, 5, verbeux=True)
print(f"\nRomberg : {val:.16f}, erreur {abs(val - I):.2e}")

1.85914091422952
1.75393109246483  1.71886115854327
1.72722190455752  1.71831884858175  1.71828269725765
1.72051859216430  1.71828415803323  1.71828184532000  1.71828183179890
1.71884112857999  1.71828197738522  1.71828183200869  1.71828183179747  1.71828183179739

Romberg : 1.7182818284590452, erreur 2.22e-16

# 3. verification de l'augmentation d'ordre par colonne
print(f"\n{'colonne j':>10}{'ordre mesure':>15}")
for j in range(4):
    E, H = [], []
    for k in range(j+1, 8):
        _, _, R = romberg(exp, 0, 1, k)
        err = abs(R[k][j] - I)
        if err > 1e-15:
            E.append(err); H.append(1/2**k)
    if len(E) >= 2:
        p = np.polyfit(np.log(H), np.log(E), 1)[0]
        print(f"{j:>10}{p:>15.2f}")

 colonne j   ordre mesure
         0           2.00
         1           4.00
         2           6.00
         3           8.00

3. L'ordre augmente bien de 2 par colonne : \(2, 4, 6, 8\). La colonne \(j=1\) n'est autre que la methode de Simpson ; les suivantes sont les formules de Newton-Cotes composees d'ordre superieur.

4. Comparaison du cout.

n = 4
while abs(trapezes(exp, 0, 1, n) - I) > 1e-14:
    n *= 2
print(f"trapezes seuls : n = {n}, soit {n+1} evaluations de f")

for k in range(1, 12):
    val, ev, _ = romberg(exp, 0, 1, k)
    if abs(val - I) < 1e-14:
        print(f"Romberg        : k = {k}, soit {ev} evaluations de f")
        break

trapezes seuls : n = 8388608, soit 8388609 evaluations de f
Romberg        : k = 5, soit 63 evaluations de f

Un facteur 130 000 sur le nombre d'evaluations de \(f\). Lorsque chaque evaluation coute une simulation de plusieurs secondes, la difference n'est pas une commodite mais la frontiere entre faisable et infaisable.