Discussion :

[anapurna]

salut

je crois que wiwaxia est parti un peu loin

sur l'un des deux axe rien ne changera
et sur l'autre trouver le centre entre les deux valeur bin c'est simplement l'axe qui as servit à faire le miroir

c'est c'est le milieu de la partie miroir c'est aussi l'inverse de la moitié positive

[wiwaxia]

... je crois que wiwaxia est parti un peu loin ...

Je n'ai fait que remédier au mieux à des données corrompues, inutilisables en l'état ... et les changements proposés conservent la forme de l'équation, ainsi que la valeur de la constante de temps (b) - probablement la plus importante.

L'erreur de compréhension des données est devenue manifeste au message #14 .
Le graphique donné plus loin (#17) parle de lui-même.

... c'est le milieu de la partie miroir c'est aussi l'inverse de la moitié positive

Où est le problème ? C'est le milieu des deux points images, ou si l'on veut l"image du milieu des deux premiers points ... Les ordonnées correspondantes sont opposées, et non pas inverses.

Je peux bien sûr me tromper ... Il serait intéressant que PatDo10 nous fournisse des éclaircissements sur ce point, et produise la liste des valeurs numériques.

[PatDo10]

Du coup la fonction à utiliser est bien S(t) = S0(1 - 2 exp(-t/T1)).

J'ai fait des recherches et comme il y a des valeurs aberrantes il faut utiliser la méthode des moindres carrés pondérés bisquare (Matlab utilise cette méthode : https://fr.mathworks.com/help/curvef...s-fitting.html) mais je n'arrive pas à comprendre ce que représente les hi dans le calcule.
Si quelqu'un à une idée.
Merci d'avance.

[wiwaxia]

Du coup la fonction à utiliser est bien S(t) = S0(1 - 2 exp(-t/T1)) ...

À condition de disposer de la bonne origine ... On ne peut rien faire pour l'instant, faute de données.

... J'ai fait des recherches et comme il y a des valeurs aberrantes il faut utiliser la méthode des moindres carrés pondérés bisquare ...

1) Tu ferais mieux de produire les valeurs numériques (TI_k , S_k).

2) Les ordonnées (S_k) ne sont-elles pas en réalité des valeurs absolues ? Les deux premières ordonnées n'étaient-elles pas négatives ? En quoi consistait le système étudié ? La notion de "temps d'inversion" implique un changement de signe ...

[PatDo10]

Les 4 premières valeurs sont négatives du coup.

Si j'ai bien compris c'est justement la méthode des moindres carrés pondérés bisquare qui permettra de déterminer les variables S0 et T1, en donnant au début de l'algorithme des données de départ pour S0 et T1 (ex S0 = 1000 et T1 = 1000). Et l'algorithme les modifiera si j'ai bien compris.

[wiwaxia]

Les 4 premières valeurs sont négatives du coup ...

S'il en était ainsi, cela ôterait toute signification au nuage de points ... Voir le graphe du message #20 .

[PatDo10]

Si les 4 premiers points sont négatifs justement la courbe seraient croissante (avec certaines valeurs aberrantes).

J'ai testé avec le solveur d'excel et celui si me donne bien les 4 premières valeurs négatives.

[wiwaxia]

Donne les coordonnées des 8 points !

S'il te plaît ...

[PatDo10]

Demandé si gentiment...

(100, 870), (180, 1357), (570, 207), (665, 154), (1137, 317), (1195, 346), (1688, 419), (2165, 441)
Du coup si on utilise la valeur absolue et sinon le y des 4 premiers points est négatifs.

[wiwaxia]

Enfin des données tangibles, que l'on peut tester !

[wiwaxia]

Re-bonjour,

J'ai regardé ce que l'on pouvait tirer de tes résultats, en tenant compte:
a) de l'inversion quasi-certaine du signe des deux premières ordonnées (#14, #16, #17), et
b) de la mauvaise disposition des deux premiers points, à laquelle on est contraint de porter remède (#20).

Voici donc les listes de 7 paires de réels, résultant:
1°) de la suppression de l'un des deux premiers vecteurs (Listes 1 et 2), ou du remplacement de ces termes douteux par leur moyenne arithmétique (Liste 3);
2°) du changement de coordonnées résultant du choix du premier vecteur comme nouvelle origine, soit:
x = x_ancien - x_1ancien , y = y_ancien - y_1ancien
pour les autres listes (4 à 6).



Ce qui nous occupe, c'est la recherche du minimum de la somme des carrés des écarts:
S = (1/2)*Som_k=2⁷(y_k - F(a, b, x_k))² dans laquelle intervient la fonction:
F(a, b, x) = a*(1 - Exp(-x/b)) , dont tu as déjà livré les expressions des dérivées partielles:
D_aF(a, b, x) = (1 - Exp(-x/b) ,
D_bF(a, b, x) = (-a*x/b^2)*Exp(-x/b) .
# Remarque: le premier terme, nul par définition, n'est pas compté:
y₁ - F(a, b, x₁) = 0 - F(a, b, 0) = 0 .

Dans un repère orthonormé dont les deux axes représentent les valeurs des coordonnées (a, b), les variations élémentaires de la somme étudiée admettent pour expression le produit scalaire:
dS = (GradS)│(MM') = D_aS(a, b)*da + D_bS(a, b)*db ,
de sorte que tout extremum de (S) est caractérisé par l'annulation de son gradient:
Grad(S) = 0 (vecteur nul) ,
ou ce qui revient au même celle de sa norme:
║Grad(S)║ = ║D_aS(a, b).u_a + D_bS(a, b).u_b║ = (D_aS(a, b)² + D_bS(a, b)²)^1/2 = 0 .

La fonction (S(a, b) se calcule sur une liste donnée (L4, L5 ou L6); elle admet pour dérivées partielles:
D_aS(a, b) = -Som_k=2⁷(y_k - F(a, b, x_k))*D_aF(a, b, x_k) ;
D_bS(a, b) = -Som_k=2⁷(y_k - F(a, b, x_k))*D_bF(a, b, x_k) .

La suite au prochain message, pour éviter la perte éventuelle d'images.

[wiwaxia]

L'exploration graphique constitue un préliminaire prudent à toute recherche de solution dans un plan.
En conséquence, on entreprend de visualiser les variations de la grandeur
z = NgradS(a, b) = ║Grad(S)║
sur un domaine approprié D = [Amin ; Amax]×[Bmin ; Bmax]
en attribuant une couleur unique à chaque valeur du rapport (r = z/Zmax), dont le dénominateur correspond à la plus grande des normes calculables aux quatre coins de l'image; un écrêtage à l'unité est évidemment nécessaire, pour le cas où l'on trouverait de plus grandes valeurs à l'intérieur du domaine.

Lorsque (r) varie de 0 à 1, la couleur du fond évolue du noir au rouge, en passant par le bleu et le mauve; les 19 ou 20 stries à dominante verte discontinue constituent des courbes de niveau.
Voir ci-dessous le codage de la matrice de pixels:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 CONST Cm = 255; Cm4 = 4 * Cm; Km = 20 * Cm;
 
 PROCEDURE CalcMatIm(Va1, Va2, Vb1, Vb2, Zm: Reel;
                     La, Ha: Z_32; VAR Ma: Tab_Pix);
   VAR Kjv:  Word;Ha1, La1, Ix, Iy: Z_32;
       Kx, Ky, Kz, p, q, r, Va, Vb, z, Zmin: Reel; Px: Pixel;
   BEGIN
     La1:= La - 1; Kx:= (Va2 - Va1) / La1;
     Ha1:= Ha - 1; Ky:= (Vb2 - Vb1) / Ha1; Zmin:= Zm;
     FOR Ix:= 0 TO La1 DO
       BEGIN
         p:= Kx * Ix; Va:= Va1 + p;
         FOR Iy:= 0 TO Ha1 DO
           BEGIN
             p:= Ky * Iy;                Vb:= Vb1 + p;
             z:= N_GradS(Va, Vb, Liste);
             IF (z<Zm) THEN CalcMin(z, Va, Vb, Zm, Azero, Bzero);
             r:= z /Zmax;                IF (r>1) THEN r:= 1;
             Kjv:= Round(Km * r);        Px[2]:= Kjv MOD Cm;
             p:= 1 - r; q:= r * p;       Px[3]:= Round(Cm4 * q);
             q:= Sqr(p);                 Px[1]:= Round(Cm * (1 - q));
             Ma[Ix, Iy]:= Px
           END
       END
   END;

L'image ci-dessous livre une vue d'ensemble sur le domaine [0 ; 3000]×[0.001 ; 3000] dans le cas de la dernière liste:



La zone qui nous intéresse se situe près du bord inférieur, sur la médiane verticale; on la retrouve approximativement dans la même région dans le cas des 3 listes (L4, L5, L6), la dernière conduisant toujours à des résultats systématiquement intermédiaires par rapport aux deux autres.
Ici les mêmes domaines [1000 ; 2000]×[100 ; 400], et les images données par (L4, L6, L5):

__

On peut ainsi localiser le zéro de la fonction avec une précision croissante, en choisissant visuellement des domaines de plus en plus étroits.
Quelques instructions permettent de mémoriser les coordonnées du minimum, lors de la synthèse de l'image Bitmap; on trouve ainsi, dans le cas de la dernière liste (L6):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
                                                Azero           Bzero
A	0	3000	B	0.0001	500	1503.76	        240.60
A	1000	2000	B	200	400	1506.27	        242.61
A	1480	1520	B	230	250	1506.07	        242.23
A 	1505.6	1506.4	B	242.0	242.4	1506.04	        242.22
A	1506.02	1506.06	B	242.20	242.24	1506.038	242.217

___
_

[wiwaxia]

Le résultat issu du calcul concomitant à la précédente synthèse de l'image constitue une information utile; le procédé mis en oeuvre ne saurait cependant être retenu: calculer (N²) valeurs, pour n'en retenir qu'une seule et en déduire sa position avec une précision médiocre (1/N), c'est pour un programme une bien faible performance.

Une localisation du zéro de la norme du gradient est envisageable par la recherche systématique de la plus faible valeur de cette fonction en quatre points voisins d'un point donné (a, b); on compare à cette fin la valeur centrale NgradS(a, b) à celles observées aux points de coordonnées (a ± h, b) et (a, b ± h).

On repère d'abord (h₀) la plus grande puissance entière de 10 simultanément inférieure à (a) et (b):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
 
 PROCEDURE CalcHini(Ai_, Bi_: Reel; VAR Hi_: Reel);
   VAR i: Byte; h, p: Reel;
   BEGIN
     h:= 1;
     WHILE ((h<Ai_) OR (h<Bi_)) DO BEGIN
                                     p:= h * 10; h:= p
                                   END;
     WHILE ((h>Ai_) OR (h>Bi_)) DO BEGIN
                                     p:= h * 0.1; h:= p
                                   END;
     Hi_:= h
   END;

puis on lance ensuite le processus de comparaison avec un pas (h) égal à 10^-3*h₀ (ordre 3, pour l'approximation liée au cadrage de la racine).
On passe systématiquement au nouveau couple de coordonnées correspondant à la plus faible des nouvelles valeurs, à condition qu'elle soit aussi inférieure à la valeur centrale NgradS(a, b).
S'il n'en existe pas, on reste au même point, mais on prend un pas 10 fois plus faible: h' = h / 10 , et l'ordre augment d'une unité.
L'arrêt intervient lorsque l'on atteint l'ordre 20, donc pour un pas h = 10^-20*h₀ < 10^-20*Min(a, b)
très inférieur au seuil de précision courante (10^-19 pour Pascal).

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
 
 PROCEDURE Rech_Zero(InL:Byte; Ai, Bi: Reel;
                     VAR V_A, V_B: Reel; Liste: Tab7V);
   TYPE Lst5R = ARRAY[0..4] OF Reel;
   CONST m = 13; n = m - 3;
   VAR i, Ordre, y: Byte; k: Z_32; h, r, s, Va, Vb, Wa, Wb: Reel; Ls: Lst5R;
   BEGIN
     k:= 0;                   Va:= Ai;         Vb:= Bi;
     CalcHini(Aini, Bini, r); h:= 0.001 * r;   Ordre:= 3;
     y:= 28;                  Inc(y, 2 * InL); E(0008);
     REPEAT
       Ls[0]:=N_GradS(Va    , Vb    , Liste);
       Ls[1]:=N_GradS(Va + h, Vb    , Liste);
       Ls[2]:=N_GradS(Va    , Vb + h, Liste);
       Ls[3]:=N_GradS(Va - h, Vb    , Liste);
       Ls[4]:=N_GradS(Va    , Vb - h, Liste);
       s:= Ls[0]; i:= 0;
       IF (s>Ls[1]) THEN BEGIN
                           s:= Ls[1]; i:= 1; Wa:= Va + h; Wb:= Vb
                         END;
       IF (s>Ls[2]) THEN BEGIN
                           s:= Ls[2]; i:= 2; Wa:= Va;     Wb:= Vb + h
                         END;
       IF (s>Ls[3]) THEN BEGIN
                           s:= Ls[3]; i:= 3; Wa:= Va - h; Wb:= Vb
                         END;
       IF (s>Ls[4]) THEN BEGIN
                           s:= Ls[4]; i:= 4; Wa:= Va;     Wb:= Vb - h
                         END;
       IF (i>0) THEN BEGIN
                       Va:= Wa; Vb:= Wb
                     END
                ELSE BEGIN
                       r:= h * 0.1; h:= r; Inc(Ordre)
                     END;
       Inc(k); We(5, y, k, 9);
       Write(Ordre:7, '   ', h:7, '   ', Va:m:n, '   ', Vb:m:n)
     UNTIL ((Ordre>19) OR KeyPressed);
     E(0014);  We(5, y, k, 9); Write(Ordre:7, '   ', h:7, '   ');
     E(0010);  Write(Va:m:n, '   ', Vb:m:n);
     V_A:= Va; V_B:= Vb;
   END;

L'algorithme converge rapidement pour les 3 listes de points, vers les limites indépendantes des valeurs initiales; voici les résultats obtenus à partir de (A_ini, B_ini) = (1500, 300):



Les valeurs affichées sur les trois dernières lignes sont dépourvues de signification numérique, hormis leur ordre de grandeur (~ 10^-15); ce sont les restes résiduels de sommes algébriques théoriquement nulles, dont les termes sont <~200 (en valeur absolue). Le rapport (~ 10^-17) dépasse à peine le seuil de précision des calculs, de sorte que l'on peut considérer qu'au point limite la relation Grad(S) = 0 est bien vérifiée.

# Cependant tout n'est pas aussi simple dans le meilleur des mondes virtuels ... parce que le programme lancé à l'aveugle sur un choix arbitraire les valeurs initiales diverge très souvent, ce qui a entraîné une longue et pénible vérification des calculs



Cela se comprend dès que l'on observe la carte des variations de la fonction étudiée:



que l'on peut voir comme la vue plongeante d'une surface en relief, d'équation z = NgradS(a, b); le parcours codé est celui d'un randonneur muni d'une carte IGN, et qui s'orienterait en permanence vers les points les plus bas⁽¹⁾, donc les zones les plus sombres.
Dès que la position de départ s'éloigne trop du point limite recherché, la progression est attirée:
- soit vers le vallon rectiligne indéfiniment descendant (vers le haut, à droite),
- soit vers la dépression située en-dessous, sur l'axe des abscisses (b ~ 0, graphe quasi-rectiligne).
Par exemple, pour A = 1500, la bonne limite n'est atteinte que pour (162 < b < 573.49); le nombre d'étapes devient très grand lorsque l'on s'approche des bornes de l'intervalle:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
A = 1500
Bini    Niter
570.00	14457
573.30	76834
573.40	77046
573.45	77183

# Pour résumer tout ce qui précède,
- l'examen de l'image représentant les variations de la fonction à deux variables permet un choix approprié des valeurs initiales;
- les résultats expérimentaux restent utilisables, moyennant la transformation décrite (listes 4, 5 et surtout 6).

Sur ce dernier point, une confirmation décisive serait apportée par le calcul des écarts type à associer aux valeurs de (a) et (b); il faudrait poursuivre de ce côté-là.

(1) ce que je vous déconseille catégoriquement, dans le monde réel de la moyenne ou haute montagne

[wiwaxia]

Une recherche des moindres carrés sans évaluation d'écart-type, c'est comme un thé sans citron, ou un poulet à la crème sans morilles.

1°) Considérons les variations de la somme ⁽¹⁾ des carrés des écarts:
S(a, b) = Som_k=2⁷(y_k - F(a, b, x_k))² dans laquelle intervient toujours la fonction:
F(a, b, x) = a*(1 - Exp(-x/b)) .

(1) Il n'intervient plus ici le facteur (1/2), qui permettait auparavant de se débarrasser du facteur (2) dans les expressions des dérivées (#31); la rectification de ce message m'a échappé: j'aurais alors dû parler de demi-somme.

La représentation des variations de S(a, b) sur un large domaine [1000 ; 3000]×[1E-6 ; 2000] suggère que le minimum de cette fonction (zone noire) se situe sur l'axe des abscisses, dans le cas des trois nuages de points étudiés (soit, dans l'ordre, les listes 4, 6 et 5):

__

Une simple réduction d'un facteur 4 des dimensions des domaines ([1000 ; 1500]×[1E-6 ; 500]) permet de vérifier qu'il n'en est rien:

__

On retrouve à plus grande échelle - ici respectivement:
[1260 ; 1280]×[290 ; 310] ; [1495 ; 1515]×[230 ; 250] ; [1730 ; 1750]×[190 ; 210]
un faisceau d'ellipses coaxiales et de même forme:

__

Le recours à des domaines encore plus étroits conduit à une localisation précise du minimum, par quelques instructions analogues à celles déjà utilisées (#32):
# pour L4: [1272.9 ; 1273.1]×[299.3 ; 299.5] _ a₀ = 1273.01888 _ b₀ = 299.37390

# pour L6: [1505.9 ; 1506.1]×[242.1 ; 242.3] _ a₀ = 1506.03815 _ b₀ = 242.21727

# pour L5: [1741.9 ; 1742.1]×[199.9 ; 200.1] _ a₀ = 1741.93775 _ b₀ = 199.94739

[wiwaxia]

Une connexion défaillante m'a contraint d'interrompre.

2°) La forme des courbes observées est confirmée par l'analyse; la fonction étudiée admet en effet au voisinage de tout point (a, b) un développement au second ordre selon la formule de Taylor^(1,2):
S(a + h_a, b + h_b) = S(a, b) + D_aS(a, b)*h_a + D_bS(a, b)*h_b + (1/2)*(D_2aS(a, b)*h_a² + 2*D_abS(a, b)*h_a*h_b + D_2bS(a, b)*h_b²) + (h_a² + h_b²)*R(h_a, h_b) ,
dans lequel le reste R(h_a, h_b) tend vers zéro en même temps que ses arguments.

Plaçons-nous au voisinage du minimum de la fonction S(a, b), caractérisé par l'annulation de ses dérivées partielles premières:
D_aS(a, b) = 0 ; D_bS(a, b) = 0 ,
et limitons-nous aux écarts pour lesquels (S) varie d'une quantité au plus égale à la contribution moyenne des six points du nuage hors de son origine, soit:
S(a', b') <= S(a, b) + S(a, b)/6 = (7/6)*S_min ;
il vient à la limite:
S(a', b') - S(a, b) = S_min/6 = (1/2)*(D_2aS(a, b)*h_a² + 2*D_abS(a, b)*h_a*h_b + D_2bS(a, b)*h_b²) + (h_a² + h_b²)*R(h_a, h_b)
et pour des écarts suffisamment petits pour que le reste devienne négligeable:
S_min/6 ~ (1/2)*(D_2aS(a, b)*h_a² + 2*D_abS(a, b)*h_a*h_b + D_2bS(a, b)*h_b²) .

On obtient finalement par le changement de coordonnées: X = h_a = a' - a ; Y = h_b = b' - b :
S_min/3 ~ D_2aS(a, b)*X² + 2*D_abS(a, b)*X*Y + D_2bS(a, b)*Y² ,
équation d'une ellipse centrée en (a, b).

Les écarts-types associés aux valeurs de (a) et (b) peuvent être représentés par les demi-dimensions (Lx, Ly) du rectangle circonscrit à l'ellipse frontière, et dont les arêtes sont parallèles aux axes du repère. Les valeurs s'obtiennent facilement par voie graphique, sur une image d'échelle appropriée; on peut aussi envisager un calcul à partir de l'équation cartésienne, mais c'est moins évident puisqu'il faut alors connaître les valeurs des 3 dérivées secondes .

__

Il est sans doute réducteur de ramener aux dimensions d'un rectangle les influences simultanées des deux variables (a, b); mais cela donne au moins une idée de la dispersion des résultats.

On obtient finalement:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
A4 = 1273   Ea4 = 29   B4 = 299   Eb4 = 35
A6 = 1506   Ea6 = 30   B6 = 242   Eb6 = 29
A5 = 1742   Ea5 = 31   B5 = 200   Eb5 = 25

soit une dispersion relative (Eu/U) de l'ordre de 1.8 à 2.3 % pour (A), beaucoup plus forte (~ 12 %) dans le cas de (b).

La précision sur le second résultat est décevante, parce que l'influence du paramètre (b) ne se fait sentir que sur la partie gauche du graphe (à cause du terme exponentiel e^-x/b), justement là où les points sont très mal connus.

(*) Voir 1 et 2

# L'exploitation graphique a été faite sur des images de dimensions initiales 500 pixels, correspondant à des domaines de largeur égale à 100:
- pour L4: [1220 ; 1320]×[250 ; 350]
- pour L6: [1450 ; 1550]×[200 ; 300]
- pour L5: [1690 ; 1790]×[150 ; 250]

[wiwaxia]

Un processus relativement simple permet de localiser le minimum de la somme S(a, b), en envisageant une succession de petits déplacements, proportionnels au gradient local de la fonction étudiée et de sens opposé:
MM' = -Kf.Grad(S) .

Le nombre d'étapes nécessaires pour atteindre la limite souhaitée, dépend de la valeur de (Kf), ainsi que de la condition d'arrêt.

a) On observe pour la constante (Kf) et un ensemble donné de positions initiales une valeur optimale résultant d'un compromis entre de trop faibles valeurs (allongeant excessivement le temps de calcul), et d'autres trop élevées (conduisant à l'explosion numérique, donc à l'absence de limite); on obtient ainsi dans le cas de la dernière liste (L6), toujours reprise par la suite:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Kf      Nmin    Nmax
0.050    192    5374
0.060    157    4476
0.070    132    3835
0.080    122    3354 
0.082    110    3272
0.084    107    3193
0.086    104    3119
0.088    101    3048
0.089    100    3013
0.090    100    2980
0.091    101    2949
0.092    103    2926
0.100    201    2838

b) La condition d'arrêt est liée à la décroissance de la norme (N) du vecteur Grad(S) , qui tend vers 0 (vecteur nul) lorsque l'on se rapproche du minimum recherché; elle exige une évaluation de sa valeur initiale.
Soit (L) une longueur caractérisant le domaine étudié [A_min ; A_max]×[B_min ; B_max] ;
on a pour la somme (S) des carrés des écarts S ~ L² ~ A_max * B_max (ou A_max² + B_max²) ;
pour la norme du gradient: N ~ L²/L = L soit (A_max² + B_max²)^1/2.
La plus sévère condition d'arrêt (N² <~ 1E-36 * (A_max² + B_max²) correspondra par conséquent à la précision maximale des calculs effectués (~ 10^-18).

Les procédures qui suivent constituent le coeur du programme , et illustrent les considérations précédentes.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
 
 PROCEDURE CalcF(Va, Vb: Reel; Vxy: Ve2D; VAR DaF_, DbF_, Dy_: Reel);
   VAR ExpE, p, q, r: Reel;
   BEGIN
     IF (Vb=0) THEN Alarme(q);
     q:= Vxy.u / Vb; IF (q<100) THEN ExpE:= Exp(-q)
                                ELSE ExpE:= 0;
     r:= 1 - ExpE;   DaF_:= r;
     p:= Va * r;     Dy_:= Vxy.v - p;
     p:= Va * Vxy.u; q:= Sqr(Vb);
     r:= p / q;      DbF_:= -r * ExpE
   END;
 
 PROCEDURE CalcS_GradS(Wa, Wb: Reel; VAR L_: Tab7V;
                       VAR Sc_, DaS_, DbS_: Reel);
   VAR k: Byte; DaF, DaS, DbF, DbS, DeltaY, Dy2, p, q, Sy2: Reel;
   BEGIN
     CalcF(Wa, Wb, L_[2], DaF, DbF, DeltaY);
     DaS:= DeltaY * DaF; DbS:= DeltaY * DbF; Sy2:= Sqr(DeltaY);
     FOR k:= 3 TO 7 DO BEGIN
                         CalcF(Wa, Wb, L_[k], DaF, DbF, DeltaY);
                         p:= Sqr(DeltaY);  q:= Sy2 + p; Sy2:= q;
                         p:= DeltaY * DaF; q:= DaS + p; DaS:= q;
                         p:= DeltaY * DbF; q:= DbS + p; DbS:= q
                       END;
     Sc_:= Sy2; DaS_:= -2 * DaS; DbS_:= -2 * DbS
   END;
 
 PROCEDURE Evolution(Ia: Byte; Aini, Bini: Reel; Adeg: Z_16);
   CONST C1 = 59; C2 = C1 + 11; L1 = 40; o = 502;
         Epsilon = 1E-36; Kf = 0.090; m = '    ';
   VAR y: Byte; k: Word; DaS, DbS, Sc, Lim, N2, p, q, Va, Va1, Vb, Vb1, Vs: Reel;
   BEGIN
     AffKf(C1 - 6, L1 - 2, Kf, Adeg); Va:= Aini;     Vb:= Bini;
     p:= Sqr(Va);                     q:= Sqr(Vb);   Lim:= Epsilon * (p + q);
     CalcS_GradS(Va, Vb, Liste, Sc, DaS, DbS); k:= 0;
     Trace_Point(Amin, Bmin, Va, Vb, Lx, Ly, Matr_Image);
     REPEAT
       p:= Kf * DaS; Va1:= Va - p;
       q:= Kf * DbS; Vb1:= Vb - q;
       p:= Sqr(DaS); q:= Sqr(DbS);     N2:= p + q;
       Inc(k);       We(C1, L1, k, 5); Wr(C2, L1, N2, 10);
       Trace_Point(Amin, Bmin, Va1, Vb1, Lx, Ly, Matr_Image);
       CalcS_GradS(Va1, Vb1, Liste, Sc, DaS, DbS);
       Va:= Va1;     Vb:= Vb1
     UNTIL ((N2<Lim) OR (k>65000));
     y:= 2 * Ia;      Inc(y, 46); E(0010);
     We(2, y, Ia, 3); Write(m, k:6, m, N2:6, m);
     Write(Va:20:14, m, Vb:20:15); A_
   END;          
 
 PROCEDURE Trace_Courbe(La, Ha: Z_32; Va1, Va2, Vb1, Vb2: Reel;
                        VAR Lx_, Ly_: Reel);
   CONST Nmax = 12; D_Pi = 2 * Pi; DPiN = D_Pi / Nmax; Ad = 360 / Nmax;
   VAR k: Byte; Aini, Am, Bini, Bm, Ct, Da, Db, DeltaV, p, St, t: Reel;
   BEGIN
     DeltaV:= Va2 - Va1; Lx_:= (La - 1) / DeltaV; Da:= 0.45 * DeltaV;
     DeltaV:= Vb2 - Vb1; Ly_:= (Ha - 1) / DeltaV; Db:= 0.45 * DeltaV;
     Am:= 0.5 * (Va1 + Va2); Bm:= 0.5 * (Vb1 + Vb2);
     AffTC;
     FOR k:= 0 TO (Nmax - 1) DO
       BEGIN
         t:= DpiN * k;
         Ct:= Cos(t); p:= Da * Ct; Aini:= Am + p;
         St:= Sin(t); p:= Db * St; Bini:= Bm + p;
         Evolution(k, Aini, Bini, Round(Ad * k))
       END
   END;

Voici ce que donne le lancement de l'itération à partir de 12 positions initiales régulièrement disposées sur une ellipse entourant le minimum cherché:
Aini = 1500 + 1350 * Cos(2*k*Pi/12) ; Bini = 250 + 225 * Sin(23*k*Pi/12) :
(Ntb : nombre d'étapes ; N2fin : valeur finale du carré de la norme considérée;
la constante Epsilon = N²_min/(A_max² + B_max²)^1/2 vaut 10^-36)



Les valeurs limites obtenues pour (a) et (b) apparaissent constantes sur les 18 chiffres affichés, donc bien indépendantes du point de départ.

# À ce stade, il n'y a pas pour les point successifs de trajectoire observable en raison de l'importance des premiers déplacements; on ne constate que leur accumulation dans la vallée quasi-rectiligne où se trouve la cuvette.
Les trajectoires n'apparaissent qu'à la condition de resserrer les positions successives, de telle sorte que chaque déplacement ne dépasse pas un pixel de l'image; cela peut s'obtenir en introduisant un nouveau facteur (Kf1) dépendant de la norme (N) du gradient:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 PROCEDURE Evolution01(Aini, Bini: Reel; Adeg: Z_16);
   CONST C1 = 59; C2 = C1 + 11; L1 = 39; o = 502;
         Epsilon = 1E-36; Kf = 0.090;
   VAR k: Word;
       DaS, DbS, Sc, Kf1, Lim, N1, N2, p, q, Va, Va1, Vb, Vb1, Vs: Reel;
   BEGIN
     AffKf(C1 - 6, L1 - 2, Kf, Adeg); Va:= Aini;     Vb:= Bini;
     p:= Sqr(Va);                     q:= Sqr(Vb);   Lim:= Epsilon * (p + q);
     CalcS_GradS(Va, Vb, Liste, Sc, DaS, DbS); k:= 0;
     Trace_PointF(Amin, Bmin, Va, Vb, Lx, Ly, Matr_Image);
     REPEAT
       p:= Sqr(DaS);  q:= Sqr(DbS);     N2:= p + q; N1:= Sqrt(N2);
       IF (N1>1) THEN Kf1:= Kf / N1
                 ELSE Kf1:= Kf;
       p:= Kf1 * DaS; Va1:= Va - p;
       q:= Kf1 * DbS; Vb1:= Vb - q;
       Inc(k);        We(C1, L1, k, 5); Wr(C2, L1, N2, 10);
       Trace_Point(Amin, Bmin, Va1, Vb1, Lx, Ly, Matr_Image);
       CalcS_GradS(Va1, Vb1, Liste, Sc, DaS, DbS);
       Va:= Va1;      Vb:= Vb1
     UNTIL ((N2<Lim) OR (k>65000));
     Trace_PointF(Amin, Bmin, Va1, Vb1, Lx, Ly, Matr_Image); E(0700);
     Wr(C1 - 1, L1 + 2, Va, o); Wr(C1 + 15, L1 + 2, Vb, o)
   END;

Voici ce qui apparaît dans le cas de 12 positions initiales (images 3 et 4), puis 40 (dernière image):

_

_