dimanche 3 février 2019

Un solveur pour la FX-602P

Dans mon article précédent, je vous avais proposé de développer un solveur numérique pour la Casio FX-602P. Ce programme lui permet en effet de combler une de ses lacunes par rapport à la HP-15C qui possède la fonction [solve] d'origine. Il n'y a probablement pas de grand secret mais, malgré les explications données en annexe D du manuel de la 15C, je ne sais pas exactement quel est l'algorithme implémenté dans la HP. Alors j'ai choisi la méthode numérique de la sécante. Elle est proche de la méthode de Newton, mais plus facile à programmer, car on n'aura pas à calculer de dérivée. On obtiendra néanmoins une bonne approximation d'un zéro de f(x) en quelques itérations seulement.

Recherche d'un zéro de f(x) = x^3 - 2x -5

J'ai repris le vieux programme en Basic pour PB-1000 qui m'avait rendu de bons et loyaux services dans les années 90. Pour la toute première version, je m'étais inspiré d'un algorithme nommé "Newton" donné dans un vieux bouquin d'analyse numérique dont j'ai oublié le titre, et que j'avais emprunté à la bibliothèque de la Fac des Sciences Montpellier II. Andreas Wichmann avait publié mon programme pour PB-1000 sur son site quelques années plus tard, avec mon accord. Le voici de retour en 2019 dans une version écrémée :

10 ANGLE1:GOTO40
20 Y=X^3-2*X-5
30 RETURN
40 PRINT"X initial="X;:INPUTX:I=1E-8
50 A=X:GOSUB20:Z=Y:X=X+I:GOSUB20
60 X=A-Z*I/(Y-Z):PRINTA;X:IFABS(X-A)>I THEN50
70 PRINT"F(X)="Z;"X="X

Usage :
  1. éditer le programme pour entrer la fonction sous la forme Y=F(X) à la ligne 20
  2. si besoin, modifier la précision qui est codée en dur à la fin de la ligne 40 (variable I)
  3. lancer le programme
  4. répondre à la question "X initial= ... ?" par exemple avec la valeur 1.
  5. le programme pour PB-1000 affiche les valeurs de A et X à chaque itération
  6. affichage du résultat à la fin (contrôle avec F(X) qui doit être proche de zéro). Attention, dans certains cas particuliers, il se pourrait qu'il n'y ait pas convergence. Dans ce cas, on utilisera la touche [BRK] pour interrompre le programme. 
Pour illustrer le déroulement de ce programme, j'ai créé l'image ci-dessous avec un montage Gimp de trois copies d'écran PB-1000 :

Convergence de l'algorithme de la sécante sur PB-1000 (montage)

Pour le portage sur la FX-602P qui nous intéresse plus particulièrement dans cet article, j'utiliserai les quatre mémoires listées ci-dessous avec les variables équivalentes dans le listing en Basic :

Basic --- FX-602P
    X --- 00
    A --- 01
    Z --- 02
    I --- F

Ce qui nous donne le programme en 38 pas ci-dessous. Il est moyennement lisible, je vous l'accorde ! Vous pouvez le saisir par exemple en P1, toujours en une seule ligne de programme sur la 602P :

Min00 1 EXP 8 +/- MinF 
LBL0 MR00 Min01 GSBP0 Min02 MR00 + MRF = 
GSBP0 - MR02 = 1/x * MRF * MR02 = +/- + MR01 = PAUSE 
Min00 - MR01 = ABS x>=F GOTO0 MR00

Ce programme fait appel au sous-programme P0 ci-dessous, qui doit contenir la fonction f(x) ; ici [ x^3 - 2x - 5 ] en 10 pas :

Min00 x^y 3 - 2 * MR00 - 5 =

Usage :
  1. entrer la fonction f(x) dans P0
  2. si besoin, modifier la précision, codée en dur au pas 4 (variable F)
  3. taper le X initial, par exemple 1
  4. lancer le programme P1
  5. pause et affichage de la valeur de X à chaque itération ; le résultat à la fin
Sur la HP-15C, pour comparer les performances de son solveur intégré, j'ai utilisé le programme suivant :

001-42,21,11   LBL A
002-      36   ENTER
003-      36   ENTER
004-       3   3
005-      14   y^x
006-      34   x<>y
007-       2   2
008-      20   *
009-      30   -
010-       5   5
011-      30   -
012-   43 32   RTN

Je vous rétrocède ci-dessous les vitesses d’exécution que j'ai mesurées sur les trois machines avec le calcul donné en exemple. Sur les versions Casio, l'affichage intermédiaire permet de suivre la convergence, mais prend évidemment du temps :

- Casio PB-1000 : environ 5 secondes pour afficher le résultat (calcul à partir de X=1)
- Casio FX-602P avec [pause] à chaque itération : 16 secondes (calcul à partir de X=1)
- HP-15C avec la fonction [solve] intégrée : temps très variables, qui dépassent parfois la minute, en fonction de l'intervalle de recherche initial (registres X et Y).

Pour conclure, j'ai réalisé une petite vidéo, à visionner de préférence en 1080p, qui permet de revivre les sensations vintage procurées par la FX-602P et la HP-15C en les faisant tourner en parallèle :  


La calcul est super long sur la HP-15C parce que l'intervalle de recherche est arbitrairement fixé à [1, 51] dans la vidéo. Pour plus d'équité dans la comparaison, j'aurais pu fixer une valeur de Y plus basse. Par exemple, en fixant l'intervalle initial à [2, 3] le temps de calcul peut être réduit à 38 secondes environ. Parmi les autres optimisations possibles, on peut désactiver le mode de calcul complexe (merci C.Ret) via la séquence de touche peu intuitive [g] [CF] [8]. Il est également intéressant de retirer les deux [ENTER] superflus en début de programme. L'inconvénient de cette méthode, c'est qui ne faudra pas oublier le double [ENTER] manuel préalable si on utilise le programme pour un calcul direct (par exemple : f(1) = -6). Enfin, on peut factoriser la formule via la méthode de Horner ; comme cela a été suggéré dans le premier commentaire à cet article (ci-dessous). Voici un nouveau programme qui prend en compte ces remarques :

001-42,21,11   LBL A
002-      20   *
003-       2   2
004-      30   -
005-      20   *
006-       5   5
007-      30   -
008-   43 32   RTN

Dernier truc : un [ENTER] de trop, quand on saisit l'intervalle de recherche initial, a également une influence sur la performance ! Par conséquent, il vaut mieux taper [2] [ENTER] [3] sans deuxième [ENTER], et directement [f] [SOLVE] [f] [A]. Grâce à ces optimisations, le temps de calcul pour afficher le résultat 2,094551481 est ramené à 12 secondes pour la HP-15C. La différence avec les calculs précédents est substantielle. Cependant, elle ne suffit pas à faire la différence face à une FX-602P bien plus réactive.