Comment calculer la machine epsilon dans MATLAB?

D'abord et avant tout, il n'y a pas une telle chose comme un do mot-clé dans MATLAB, afin de l'éliminer de votre code. Aussi, n'utilisez pas de eps comme une variable réelle. C'est un pré-définis en fonction MATLAB qui calcule machine epsilon, c'est aussi ce que vous êtes en train de calculer. En créant une variable appelée eps, vous fait de l'ombre à la fonction elle-même, et ainsi de toutes les autres fonctions de MATLAB qui nécessitent son utilisation va se comporter de façon inattendue, et ce n'est pas ce que vous voulez.

Utiliser quelque chose d'autre à la place, comme macheps. Aussi, votre algorithme est légèrement incorrect. Vous avez besoin de vérifier pour 1.0 + (macheps/2) dans votre while boucle, pas 1.0 + macheps.

En d'autres termes, ce faire:

macheps = 1;

while 1.0 + (macheps/2) > 1.0
    macheps = macheps / 2;
end

Cela devrait vous donner 2.22 x 10^{ -16}, qui est d'accord avec MATLAB si vous tapez dans eps dans l'invite de commande. Pour le vérifier:

>> format long
>> macheps

macheps =

     2.220446049250313e-16

>> eps

ans =

     2.220446049250313e-16

---

Bonus

Dans le cas où vous ne le saviez pas, la machine epsilon est la limite supérieure de l'erreur relative due à l'arithmétique à virgule flottante. En d'autres termes, ce serait le maximum de la différence attendue entre un vrai nombre à virgule flottante, et qui est calculée sur un ordinateur en raison du nombre fini de bits utilisés pour stocker un nombre à virgule flottante.

Si vous vous souvenez, les nombres flottants, inévitablement, sont représentés comme des bits binaires sur votre ordinateur (ou à peu près tout numérique). En termes de La norme IEEE 754 virgule flottante standard, MATLAB suppose que toutes les valeurs numériques sont de type double, ce qui représente des nombres à virgule flottante 64 bits. Vous pouvez bien évidemment remplacer ce comportement par explicitement la conversion en un autre type. Avec la norme IEEE 754 virgule flottante standard, pour double type pour la précision des chiffres, il y a 52 bits qui représentent la fractionnaire partie du nombre.

Voici un beau diagramme de quoi je parle:

^{Source: Wikipédia}

Vous voyez qu'il est un peu réservé pour le signe du nombre, 11 bits sont réservés pour l'exposant de la base et enfin, 52 bits sont réservés pour la partie fractionnaire. Ce total s'ajoute à 64 bits. La partie fractionnaire est une collection ou la somme des nombres de la base 2, avec des exposants négatifs à partir de -1 en bas à -52. Le MSB du nombre à virgule flottante commence avec2^{ -1}, tout le chemin vers le bas pour 2^{ -52} que le LSB. Essentiellement, epsilon machine calcule la résolution maximale de la différence pour une augmentation de 1 bits en binaire entre deux nombres, étant donné qu'ils ont la même signe et la même exposant de la base. Techniquement parlant, la machine epsilon en fait égal à 2^{ -52} que c'est la résolution maximale d'un seul bit à virgule flottante, compte tenu de ces conditions dont j'ai parlé plus tôt.

Si vous regardez le code de plus près, la division par 2 est de décalage de bits votre nombre à la droite d'une position à chaque itération, en commençant à la valeur entière de 1, ou 2^{0}, et nous profitons de ce numéro et les ajouter à 1. Nous gardons de décalage de bits, et de voir ce que la valeur est égale à par l'ajout de cette peu décalée de la valeur de 1, et nous allons jusqu'au point où lorsque nous décalage de bits vers la droite, un changement n'est plus enregistré. Si vous avez de décalage de bit plus à droite, la valeur devient 0 en raison de dépassement de capacité, et donc 1.0 + 0.0 = 1.0, et ce n'est plus > 1.0, qui est ce que la while de la boucle de contrôle.

Une fois le while boucle se ferme, c'est ce seuil qui définit machine epsilon. Si vous êtes curieux, si vous coup de poing dans 2^{ -52} dans l'invite de commande, vous obtiendrez ce eps est égale à:

>> 2^-52

ans =

     2.220446049250313e-16

Ce sens que vous êtes décalage d'un bit vers la droite 52 fois, et le point avant que la boucle s'arrête serait à son LSB, qui est 2^{ -52}. Par souci d'être complet, si vous étiez à la place d'un compteur à l'intérieur de votre while boucle, et compter combien de fois le while boucle s'exécute, il va exécuter exactement 52 fois, représentant 52 bits se décale vers la droite:

macheps = 1;
count = 0;
while 1.0 + (macheps/2) > 1.0
    macheps = macheps / 2;
    count = count + 1;
end

>> count

count =

  52