Calcul d'une racine cubique

Il n'y a pas de fonction système qui calcule la racine cubique. Il faut la définir avec exp() et log().
Appelons r3(x) = exp(log(x)/3) (avec r3(0) = 0) cette fonction racine cubique "système".

from math import *

def r3(x):
    """ Racine cubique système """
    if x == 0:
        return 0
    return exp(log(x) / 3)
    
print(r3(125), r3(970299), r3(10**300)) #renvoie

5.000000000000001 98.99999999999999 9.999999999999825e+99

  On doit pouvoir faire mieux...

  On veut écrire une fonction rac3 ( a ) calculant, pour un nombre flottant a , la racine cubique de a en respectant le cahier des charges suivant:
    (1) On ne doit faire que des calculs arithmétiques élémentaires (+, -, ×, /, puissance, partie entière)
    (2) Le calcul de rac3 ( a ) doit être faisable pour un flottant a 0 avec | a | < + = 2 1024 - 2 971 1.79 10 308
    (3) Le calcul de rac3 ( a ) doit être de complexité O(1) opérations arithmétiques élémentaires
    (4) Le résultat doit être correct: on doit avoir pour tout flottant a , rac3 ( a ) - a 3 ε a 3 , (avec ε = 2 - 53 dans l'idéal)

II) Stratégie

1) Idée

  Les multiplications par des puissances de 2 sont des calculs exacts sur les flottants. Pour calculer rac3 ( a ) , l'idée est d'écrire a = 8 n b , avec b [ 1 , 8 [ et n Z . Alors rac3 ( a ) = 2 n rac3 ( b ) .  On se ramène ainsi au calcul de rac3() sur l'intervalle I = [ 1 , 8 ] .
  Pour b I , on construit une suite ( u n ) qui converge très rapidement vers b 3 , en exploitant la méthode de Newton associée à la fonction g ( x ) = x 3 - b , pour laquelle b 3 est solution de g ( x ) = 0 .  

2) Calcul sur l'intervalle [1, 8[

Soient f ( x ) = x - g ( x ) g ' ( x ) = x - x 3 - b 3 x 2 = 2 3 x + b 3 x 2 = 1 3 ( 2 x + b x 2 ) et ( u n ) définie par u 0 = 3 2 et n N , u n + 1 = f ( u n ) avec f ( x ) = 1 3 ( 2 x + b x 2 ) .
Notons α = b 3 . On a f ( x ) - α = 1 3 ( 2 x + α 3 x 2 ) - α = 1 3 x 2 ( 2 x 3 - 3 α x 2 + α 3 ) = 1 3 ( x - α ) 2 ( 2 x + α x 2 ) = 1 3 ( x - α ) 2 ( 2 x + α x 2 )
Comme b [ 1 , 8 ] , α [ 1 , 2 ] . Alors: x I , 0 f ( x ) - α 4 3 ( x - α ) 2 (1)
Mais aussi 0 f ( 3 2 ) - α 1 3 ( 3 2 - α ) 2 ( 4 3 + 8 9 ) 1 3 1 4 20 9 < 1 5 ( car α [ 1 , 2 ] ) 1 10 . Donc 0 u 1 - α 1 5
Puis     avec (1)   0 u 2 - α 4 3 ( u 1 - α ) 2 4 3 1 25 = 4 75 .  Et par récurrence à partir de (1), 0 u 2 + n - α 1 M ( M ( u 2 - α ) ) 2 n avec M = 4 3 .
On en déduit que 0 u 2 + n - α 3 4 ( 16 225 ) 2 n pour n 0 . Avec n = 4 , 0 u 2 + 4 - α 3.2 * 10 - 19
Le réel u 6 approche   α = b 3 à 3.2 * 10 - 19 près pour tout réel b [ 1 , 8 ] . On a alors | u 6 - α | 3.2 * 10 - 19 α   (car α 1 ). Par ailleurs, 3.2 * 10 - 19 < 2 - 53 1.1 * 10 - 16 .
En posant rac3 ( b ) = u 6 , on définit une fonction qui (théoriquement) approche la fonction b 3 avec une précision relative de 3.2 * 10 - 19   sur [1,8] .

3) Bilan: calcul de rac3(a)

Soit a un flottant strictement positif.
    (1) On calcule l'entier n Z tel que 8 n a < 8 n + 1 puis le flottant b tel que   a = 8 n b ( b = a 8 n [ 1 , 8 [ ),
    (2) On calcule le terme u 6 de la suite ( u n ) définie par u 0 = 1.5 et n N , u n + 1 = 1 / 3 * ( 2 * u n + b / u n 2 )
    (3) Alors rac3 ( a ) = 2 n * u 6

III) Algorithmes

  L'entier n Z tel que 8 n a < 8 n + 1 est   n = ln ( a ) ln ( 8 ) . Mais l'usage du ln() n'est pas autorisé ici. Comme un flottant a est codé sous la forme a = s * m * 2 e , le système a un accès immédiat à la mantisse m et l'exposant e .
  Python propose la fonction frexp(a) du module math qui calcule pour un flottant a le couple ( x , p ) tel que a = x * 2 p avec x [ 0.5 , 1 [ et p Z .

Q1) En utilisant la fonction frexp(), écrire une fonction exposant_8(a) (avec a flottant > 0) calculant l'entier relatif n tel que 8 n a < 8 n + 1 . Par exemple, vérifier que:

print(exposant_8(8**10 - 0.01), exposant_8(8**11), exposant_8(8**13+0.1))→9 11 13
print(exposant_8(1/511.99), exposant_8(1/512), exposant_8(1/512.1))→-3 -3 -4

Q2) Ecrire une fonction rac3(a) (avec a ≥ 0) calculant rac3(a) en suivant la méthode exposée au II). Par exemple, vérifier que:

print(rac3(125), rac3(970299), rac3(10**300)) → 5.0 99.0 1e+100

IV) Tests de contrôle

Q4) Justifier rapidement que rac3() vérifie les points (1), (2) et (3) du cahier des charges.

  Dans la suite on vérifie au mieux le point (4) du cahier des charges, la correction. Un problème est qu'il n'y a pas de fonction système de référence pour calculer la racine cubique qui permettrait de contrôler la correction de la fonction rac(3). La solution est simple: on vérifiera que rac3 3 = Id

La fonction suivante controle_F(f, a, b, n) trace le graphe de la fonction f 3 Id - 1 sur [a, b] avec au total 2 n + 1 points de tracés équirépartis dans [a,b]. Si a et b sont des flottants, alors tous les points de la subdivision le sont aussi.

import numpy as np
import pylab as plt

def controle_F(f, a, b, n):
    """ graphe de f(x)^3/x - 1 sur [a,b] avec 2**n points"""
    plage_x = np.linspace(a, b, 2**n + 1)
    plage_y = [f(x)**3 / x - 1 for x in plage_x]
    plt.grid(True)
    titre = "Erreur relative de {0}^3 / Id".format(f.__name__)
    plt.title(titre)
    plt.plot(plage_x, plage_y)
    plt.show()

Q5) En faisant les tests suivants, constater que la fonction système r3() est un peu meilleure que la fonction rac3() aux alentours de 1, mais que ce n'est plus vrai sur un intervalle se rapprochant de 0 ou de l'infini.

controle_F(rac3, 1, 8, 13)
controle_F(r3, 1, 8, 13)
controle_F(rac3, 2**-20, 2**-10, 13)
controle_F(r3, 2**-20, 2**-10, 13)
controle_F(rac3, 2**10, 2**20, 13)
controle_F(r3, 2**10, 2**20, 13)

V) Amélioration de la fonction rac3()

  Lors des tests de la Q5), la plage des valeurs prises par h ( x ) = rac3 ( x ) 3 x - 1 est [-6ε, 6ε] avec ε = 2 - 53 1.1 * 10 - 16 .  C'est bien mais pas parfait: on doit pouvoir ramener cette plage à [-4ε, 4ε]. En effet, on peut démontrer que si f est une fonction à arrondi correct (c'est le cas de f ( x ) = x 3 ) et si A ( x ) x - 1 ε avec  ε= 2 - 53 , alors l'erreur err = A ( f ( x ) ) f ( x ) - 1 commise en évaluant f ( x ) est majorée par err ε ( 1 + x f ' ( x ) f ( x ) ) . Avec f ( x ) = x 3 , alors err 4 ε car f ' ( x ) = 3 x 2 .
  Ici, comme rac3 ( x ) devrait être évaluée  (relativement) à ε près, alors rac3 3 ( x ) doit être évaluée à 4ε près. Ce n'est pour l'instant pas toujours le cas.

Q6) En revenant au codage d'un nombre flottant suivant la norme IEEE754, expliquer pourquoi les fonctions pred(x) et succ(x) suivantes renvoient le prédécesseur et le successeur d'un flottant x, c'est à dire le plus grand flottant strictement plus petit que x et le plus petit flottant strictement plus grand que x.   (On pourra se ramener au cas où 1 ≤ x < 2 )

def pred(x):
    """Prédécesseur du flottant x"""
    return x - 2**-53 * x
def succ(x):
    """Successeur du flottant x"""
    y = x + 2**-53 * x
    if y == x:
        return x + 2**-52 * x
    return y

  Lors du calcul de   a 3 , on va choisir, parmi les trois candidats p = pred(r), r = rac3(a) et s = succ(r) quel est le meilleur en mesurant l'écart | x 3 - a | pour chacun de ces trois candidats x et en choisissant celui pour lequel l'écart est minimum.  Le problème est d'éviter l'erreur d'élimination qui se produit en calculant x 3 - a . Voilà ce que l'on suggère:

  On se ramène déjà au cas où a ∈ [1,8[. On note p, r, s les nombres pred(r), r, succ(r) avec r = rac3(a). On veut calculer significativement   x 3 - a avec x { p , r , s } .
  La mantisse de x est codée sur 52 bits et comme x [ 1 , 2 [ , alors x = ( 1 , m 1 ... m 52 ) 2 * 2 0 .

  On écrit x = x 1 + x 2 avec x 1 = [ 2 16 * x ] 2 16 = ( 1 , m 1 ... m 16 0 ... 0 ) 2 et x 2 = x - x 1 = ( 0 , m 17 ... m 52 ) * 2 - 16 . Alors le calcul de x 1 3 se fait exactement, car ( 2 - 16 ) 3 = 2 - 48 et 48 < 52. De plus, remarquons que x 2 < 2 - 16 .
  Alors: x 3 - a = ( x 1 + x 2 ) 3 - a = x 1 3 - a + 3 x 1 2 x 2 + 3 x 1 x 2 3 + x 2 3 = A + B avec A = x 1 3 - a et B = 3 x 1 2 x 2 + 3 x 1 x 2 3 + x 2 3 = x 2 * ( 3 x 1 * x + x 2 * x 2 ) .
Le calcul de A est fait exactement . Le calcul de B se fait assez bien, malgré une erreur d'absorption car x 2 << x 1 .

D'après les test du IV), les valeurs de | x 3 - a | sont dans l'intervalle [0,6α=6×ε×a] avec α = 2 - 53 a le plus petit écart possible entre deux nombres proches de a .
Dans le calcul de x 3 - a = A + B il se produit une erreur d'élimination car A 3 x 1 2 x 2 B , mais la somme x 3 - a = A + B pourra valoir n β avec β 2 - 53 B et n { 0 , ... , N } avec N β 6 α = 6 2 - 53 a , c'est à dire N 6 a B 6 x 3 3 x 1 2 x 2 6 x 1 3 3 x 1 2 x 2 = 2 x 1 x 2 > 2 17 > 100000 . Il y a plus de 100 000 valeurs possibles pour x 3 - a ...
Le calcul de l'écart   x 3 - a est beaucoup plus fin: l'effet de l'erreur d'élimination initiale est considérablement  atténuée.

Q7) Ecrire une fonction ecart(x, a) calculant | x 3 - a | = | A + B | en suivant les explications précédentes .

Q8) Ecrire une fonction Rac3(a) calculant a 3 en suivant les consignes ci-dessous:

def Rac3(a):
    si a = 0 alors le résultat est 0
    n = exposant_8(a)
    b = a / 8**n
    r = 1.5
    for i in range(6):
        r = (2*r + b / (r * r))/3
    p, s = pred(r), succ(r)
    trouver x ∈ {p, r, s} tel que ecart(x, b) soit minimum
    le résultat est 2**n * x

Q9) En refaisant les tests de Q5 avecRac3(), vérifier que l'objectif semble atteint.

Q10) En utilisant le module decimal, vérifier que même sur [1, 8] la fonction Rac3() est meilleure que r3() en testant:

from decimal import *
getcontext().prec = 20  # Fixe le nombre de chiffres significatifs

def controle_D(f, a, b, n):
    """ graphe de f(x)^3/x - 1 sur [a,b] avec 2**n points"""
    plage_x = np.linspace(a, b, 2**n + 1)
    plage_y = [Decimal(f(x))**3 / Decimal(x) - 1 for x in plage_x]
    plt.grid(True)
    titre = "Erreur relative de {0}^3 / Id".format(f.__name__)
    plt.title(titre)
    plt.plot(plage_x, plage_y)
    plt.show()

controle_D(r3, 1, 8, 13)
controle_D(Rac3, 1, 8, 13)

© 2015 - Eric Obermeyer      Powered by      Corrigé