algo recherche racine f(x) sur HP9100

[OlidaBel]

C'est très exactement ce qu'il ne faut pas faire ! C'est le mal ! ...
.. mais moi j'aime déroulé mes codes dans le vice et la luxure

Boulot impressionnant.
oui, tu es plein de vices hélas ! ("et c'est là qu'il ose")
Aurais-tu été g33k en Perl ou en C dans une autre vie ?

[C.Ret]

Comprenez bien que le précèdent code est inutilement compliqué afin de mimer l'affichage d'une HP-9100.

Il peut être perfectionné et simplifié :

RROF9100 on HP-28S (numérique).gif (17.25 Kio) Vu 2229 fois

Les couleurs servent à bien indiquer l'emboitement des structures de contrôle du programme et surtout la correspondance exacte nécessaire dans la pile opérationnelle des arguments. Cette nouvelle version moins profonde "tricote" moins que la précédente.

Outre le fait qu'elle rationnalise un peu l'ordre des arguments qui sont maintenant identiques entre la saisie et les résultats, cette dernière version a aussi l'avantage d'exécuter la recherche numérique des racines quelque soit le mode d'évaluation des constantes et fonctions du HP-28S (respectivement FLAG 35 et FLAG 36). En effet, avec le premier code, si l'on laisse les constantes symboliques dans la définition de la fonction, on peut avoir un fonctionnement perturbé ou impossible.

Je ne m'en étais pas rendu compte initialement, ayant exécuté un 35 CF afin que, par exemple, [■][ PI ] retourne 3.14159265359 et non pas 'π'. Ce dernier code fonctionne indépendamment du mode d'évaluation symbolique et quelque soit la façon dont est implémentée la function

J'oublie de dire que l'instruction RND permet d'arrondir x en fonction du format d'affichage et ainsi diminuer les temps de recherche en fonction de la précision souhaitée. Avec un petit n FIX on peut se limiter à une précision en rapport avec le sujet étudier. par contre, au format d'affichage STD la précision maximale (en fonction de l'amplitude des racines) sera automatiquement utilisée.


Je viens de recevoir un nouvel équipement qui a en outre l'avantage de pouvoir utiliser un module enfichable contenant une bibliothèque de programmes dont l'un a attiré mon attention :

ML 08 - Méthode utilisée.gif (7.25 Kio) Vu 2236 fois

Il semblerait donc, qu'à part un petit détail concernant la façon dont les Δx sont calculés, on ai ici un code fort similaire à celui que l'on trouve sur la carte magnétique du HP-9100 !

Pour trouver les trois racines de la fonction ci-dessus, il me suffit sur mon nouvel équipement muni de son Module de base (aka Master Library)de saisir le code suivant:

Code : Tout sélectionner

Lbl A'   (( STO 09 lnx 1/x * 1 exp * PI lnx ) y^x RCL 09 - PI exp )   INV SBR 

Ce qui pour un novice comme moi peu habitué à son clavier un peu capricieux et dont la disposition des touches m'est inconnue (et pleine de pièges) nécessite un peu de concentration.
Je vérifie donc un à un à l'aide de la touche SST chacun des codes mémorisés:

Code : Tout sélectionner

--- --   --- --   --- --   --- --   --- --   --- --   --- --   --- --   --- --
000 76   001 16   002 53   003 53   004 42   005 09   006 23   007 35   008 65
009 01   010 22   011 23   012 65   013 89   014 23   015 54   016 45
017 43   018 09   019 75   020 89   021 22   022 23   023 54   024 92
--- --   --- --   --- --   --- --   --- --   --- --   --- --   --- --

Je vérifie par quelques valeurs connues (en fait vérifiée sur un HP Prime) que cette fonction est convenablement programmée:
2nd Fix 4
1.8584 2nd A' affiche comme attendu -3.0786

ML 08 - Zeros of function (Testing Function Response).gif (69.35 Kio) Vu 2220 fois

5.8522 2nd A' affiche effectivement 4.3069
2nd INV Fix

Puis, je lance le programme de recherche des zéros d'une fonction depuis la biliothèque de programmes :
2nd Pgm 08
Et je saisis les paramètres :
1.1 A (limite inférieur - départ de la recherche)
9.1 B (limite supérieure - fin de la recherche)
2 C (Δx minimal existant entre chaque racine de la fonction)
11 +/- 2nd INV log D qui affiche 1.-11 (précision minimale demandée - je pousse l'engin dans ces retranchements - de toutes façon la seconde page du manuel indique que la garantie du constructeur ne s'applique qu'au premier propriétaire que je ne suis pas)
Enfin, je lance chaque recherche de racine à l'aide de la touche E :
E affiche après un temps assez long 1.366657987
E qui affiche 3.141592654
E qui affiche 8.907148774
E qui affiche rapidement "9.9999999 99" (affichage clignotant)

Marche bien malgré son grand age ! Juste parfois quelque rares rebonds de touches sur tout en tapant vite est mal (mauvaise habitude prise sur des claviers plus jeunes).

[Schraf]

@C.Ret : Même si le premier algorithme ressemble à de la dichotomie, il a l'avantage (ou l'inconvénient) de ne mémoriser qu'une seule valeur pour x, à chaque changement de signe du produit on divise le pas par 10. Pour la dichotomie traditionnelle on est obligé de mémoriser les 2 bornes puis de calculer le milieu pour choisir dans quel intervalle continuer. Par contre je suis étonné qu'en entrant 1.1 et 9.1 la machine arrive à te renvoyer 3 valeurs... on aurait pu imaginer qu'il prenne le milieu (1.1+9.1)/2 = 5.1 puis qu'il travaille dans [1.1 ; 5.1] avec à nouveau changement de signe dans [1.1 ; 3.1] et que finalement il ne trouve que la première solution 1.366.

Il doit manquer un bout dans la description de leur algorithme car visiblement la machine repart dans une recherche entre 3.1 et 5.1 puis entre 5.1 et 9.1.

PS. Je viens de lire page 26 du manuel et j'ai compris ! C'est le Delta qui permet de découper l'intervalle et il fait une recherche sur chaque intervalle, cool...

[C.Ret]

Il doit manquer un bout dans la description de leur algorithme car visiblement la machine repart dans une recherche entre 3.1 et 5.1 puis entre 5.1 et 9.1.

Oui, je n'ai pas scanné toutes les explications. En 1977, chaque programme du module prend autant de pages que nécessaire.
Je vois que tu as découvert de quelle engin il s'agit.

Oui, le Δx donné par l'utilisateur n'est pas exactement identique à celui de l'algorithme du HP-9100. Il y a une sophistication, l'idée directrice est que ce soit la machine qui cherche tous les zéros dans l'intervalle [ a b ). D'où aussi l'astuce du "9.9999999 99" clignotant qui indique qu'il n'y a plus de racines dans l'intervalle étudié.

La seul similitude est que s'il est pris trop grand, ce Δx va masquer des racines trop proches.

Evidemment, l'avantage de l'algo du HP-9100 est qu'il est concis et n'utilise qu'un minimum de ressources de la machine et l'on peut suivre les opération en voyant défiler en "temps réel" les résultats de la recherche !

[Schraf]

@C.Ret, il manque le pas 007 35 dans ton programme mais c'est un détail. J'ai une vraie Ti-58C avec le module mais bon, pour filmer c'est aussi bien avec l'appli sur téléphone :

En vidéo : Lancement du prog 08 + paramètres + résultats

[C.Ret]

Oups! Post corrigé ! Merci Eric.

Sur la vidéo on voit surtout que le simulateur va bien plus vite que la vraie. Rien qu'à l'insertion de la bande au-dessu des touches. je ne suis pas capable d'aller aussi vite !
Effectivement, pour une vidéo c'est mieux.

En vrai, ma TI-58c trouve les trois racines après plusieurs dizaines de secondes comptées à partir de chaque pression sur la touche utilisateur E.

Edit:
2'20" - 2'27" - 2'35" - 2" à 1.E-11
En fait il faut à chaque fois à peu près le même temps à une dizaine de secondes près; environ 2'30" pour une précision minimale demandée de 1E-11.
1'12" - 1'27" - 1'17" - 2" à 1.E-05
0'45" - 0'45" - 0'53" - 2" à 1.E-03

[Schraf]

Voici une adaptation en APL de la version proposée pour les Ti-58/59 avec le module Master Library.

Code : Tout sélectionner

  S←FIND;A;B;D;H;X;R
[1] 'Parametres : a b delta ?'
[2] A B D←⎕
[3] X←A+(B-A)×(¯1+⍳1+N)÷N←⌈(B-A)÷D
[4] X←0 (D←D÷2) D∘.+(0>,2×⌿FNC X)/,¯1↓X
[5] →(1E¯6 < D)/4
[6] S←S[⍋S←,1↑X]

C'est une dichotomie à ma façon, dans le sens où on continue à couper en 2 tous les intervalles tant qu'ils ont une largeur > 10^-6. Donc au lieu de faire une dichotomie sur chaque intervalle, en APL je traite tous les intervalles en même temps !

Une traduction assez fidèle en Python serait :

Code : Tout sélectionner

from math import exp,log,pi

def f(x):
  return (exp(1) * log(pi) / log(x)) ** x - exp(pi)

def find(a,b,d):
  n = int((b - a) / d)
  x = [a + (b - a) * k / n for k in range(1 + n)]
  while d > 1e-6:
    suiv = []
    for i,v in enumerate(x):
      if i < len(x) - 1 and f(v) * f(x[i+1]) < 0:
        suiv += [v, v + d/2, v + d]
    x = list(suiv)
    d /= 2    
  return sorted([x[0],x[3],x[6]])
  
>> find(1.1,9.1,2)
[1.3666568756103517, 3.1415916442871095, 8.907147979736328] 

# Valeurs successives de x :

>> find(1.1,9.1,2)
[1.1, 3.1, 5.1, 7.1, 9.1]				# Découpage initial avec delta = 2
[1.1, 2.1, 3.1, 3.1, 4.1, 5.1, 7.1, 8.1, 9.1]		# Découpage avec delta = 1 pour 1.1, 3.1 et 7.1
[1.1, 1.6, 2.1, 3.1, 3.6, 4.1, 8.1, 8.6, 9.1]		# Découpage avec delta = 0.5 pour 1.1, 3.1 et 8.6, etc.
[1.1, 1.35, 1.6, 3.1, 3.35, 3.6, 8.6, 8.85, 9.1]
[1.35, 1.475, 1.6, 3.1, 3.225, 3.35, 8.85, 8.975, 9.1]
[1.35, 1.4125, 1.475, 3.1, 3.1625, 3.225, 8.85, 8.9125, 8.975]
[1.35, 1.38125, 1.4125, 3.1, 3.13125, 3.1625, 8.85, 8.88125, 8.9125]
...

Quelques explications supplémentaires :

Code : Tout sélectionner

      FIND
Parametres : a b delta ?
⎕:
      1.1 9.1 0.1
      
1.366657257 3.141592407 8.907148743

[3] : X = 1.1 1.2 1.3 ... 8.8  8.9 9 9.1	Découpage de l'intervalle suivant delta
[4] : 2×⌿FNC X : On cherche toutes les images de X par la fonction puis on fait les produits 2 à 2 et on regarde ceux qui sont négatifs
Par exemple :
2×/ 3 4 5 6
12 20 30		(Càd 3*4, 4*5 et 5*6)
Ce qui permet de filtrer les X que l'on peut garder : Par exemple (0 1 1 0) / 1 2 3 4 masque 1 et 4 mais garde 2 et 3
[4] :       X
1.3  3.1  8.9 
1.35 3.15 8.95
1.4  3.2  9 

On ajoute aux X restants 3 valeurs : 0, D/2 et D (c'est le début de la dichotomie)
Exemple :
      1 2 3 ∘.+ 40 50 60
41 51 61
42 52 62
43 53 63

[5] : On continue tant que D > 10^-6
[6] : S←S[⍋S←,1↑X] permet de trier les valeurs par ordre croissant...
Exemple :
      S←8 3 5 1
      ⍋S
4 2 3 1	8 doit être à la 4e position, 3 à la 2e position, 5 à la 3e position et 1 à la 1ere position
      S[⍋S]
1 3 5 8

Exemple de déroulement du contenu de X :

Code : Tout sélectionner

      FIND
Parametres : a b delta ?
⎕:
      1.1 9.1 2
      
1.1 3.1 7.1
2.1 4.1 8.1
3.1 5.1 9.1

1.1 3.1 8.1
1.6 3.6 8.6
2.1 4.1 9.1

1.1  3.1  8.6 
1.35 3.35 8.85
1.6  3.6  9.1 
...
8.907128906 3.141564941 1.36663208 
8.907144165 3.1415802   1.366647339
8.907159424 3.141595459 1.366662598

8.907144165 3.1415802   1.366647339
8.907151794 3.14158783  1.366654968
8.907159424 3.141595459 1.366662598

1.366647339 3.1415802 8.907144165

[Schraf]

Une version Applesoft BASIC (On est sur une grosse machine mais le code proposé doit être facilement transposable sur les pockets SHARP, CASIO...). Vous pourrez tester le programme en le copiant-collant ici. Vous tapez 1.1 Entrée, 9.1 Entrée et 2 Entrée.

Démonstration en vidéo sur cet autre émulateur plus réaliste

Code : Tout sélectionner

5 DIM PR(100):DIM SU(100):P=3.1415926
10 DEF FN F(X) = (exp(1) * log(P) / log(X)) ^ X - exp(P)
15 HOME
20 INPUT "A,B,Delta ?";A,B,D
30 N = INT((B-A)/D): U = -1
40 FOR K=0 TO N
45 X = A + K * (B - A) / N
50 IF FN F(X) * FN F(X + D) <= 0 THEN GOSUB 200
55 NEXT K
60 V = -1
70 FOR K=0 TO U-1
80 IF FN F(PR(K)) * FN F(PR(K + 1)) <= 0 THEN GOSUB 300
100 NEXT K
110 IF D < 1E-6 GOTO 600
115 D = D / 2: PRINT "D=";D
120 FOR K = 0 TO V: PR(K) = SU(K): PRINT PR(K):NEXT K
140 GOTO 60
200 FOR J=0 TO 2
210 U = U + 1: PR(U) = X + J * D / 2
230 NEXT J
240 RETURN
300 FOR J=0 TO 2
310 V = V + 1: SU(V) = PR(K) + J * D / 2
330 NEXT J
340 RETURN
600 PRINT "--":FOR J=0 TO V / 3: PRINT SU(3 * J): NEXT J

Version Apple II

AppleII.jpg (21.98 Kio) Vu 2174 fois

Autre test :

Code : Tout sélectionner

10 DEF FN F(X) = X * X - X - 1

A,B,Delta ?-5
??5
??.1
-0.6180343627929689
1.6180335998535151