Version imprimable
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
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.Citation:
Envoyé par anapurna
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.
Où est le problème :D ? 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.Citation:
Envoyé par anapurna
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.
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.
À condition de disposer de la bonne origine ... On ne peut rien faire pour l'instant, faute de données.Citation:
Envoyé par PatDo10
1) Tu ferais mieux de produire les valeurs numériques (TIk , Sk).Citation:
Envoyé par PatDo10
2) Les ordonnées (Sk) 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 ...
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.
S'il en était ainsi, cela ôterait toute signification au nuage de points ... Voir le graphe du message #20 .Citation:
Envoyé par 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.
Donne les coordonnées des 8 points !
:rose: S'il te plaît ... :rose2:
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.
:fleur: :fleur2: Enfin des données tangibles, que l'on peut tester ! :D
Re-bonjour, :D
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 = xancien - x1ancien , y = yancien - y1ancien
pour les autres listes (4 à 6).
Pièce jointe 418094
Ce qui nous occupe, c'est la recherche du minimum de la somme des carrés des écarts:
S = (1/2)*Somk=27(yk - F(a, b, xk))2 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:
DaF(a, b, x) = (1 - Exp(-x/b) ,
DbF(a, b, x) = (-a*x/b^2)*Exp(-x/b) .
# Remarque: le premier terme, nul par définition, n'est pas compté:
y1 - F(a, b, x1) = 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') = DaS(a, b)*da + DbS(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)║ = ║DaS(a, b).ua + DbS(a, b).ub║ = (DaS(a, b)2 + DbS(a, b)2)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:
DaS(a, b) = -Somk=27(yk - F(a, b, xk))*DaF(a, b, xk) ;
DbS(a, b) = -Somk=27(yk - F(a, b, xk))*DbF(a, b, xk) .
La suite au prochain message, pour éviter la perte éventuelle d'images.
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:
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:Code:
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
Pièce jointe 418134
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):
Pièce jointe 418137_Pièce jointe 418138_Pièce jointe 418142
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):
Pièce jointe 418146_Pièce jointe 418149_Pièce jointe 418153_Code:
2
3
4
5
6
Pièce jointe 418156_Pièce jointe 418157
Pièce jointe 418160
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 (N2) 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 (h0) la plus grande puissance entière de 10 simultanément inférieure à (a) et (b):
puis on lance ensuite le processus de comparaison avec un pas (h) égal à 10-3*h0 (ordre 3, pour l'approximation liée au cadrage de la racine).Code:
2
3
4
5
6
7
8
9
10
11
12
13
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*h0 < 10-20*Min(a, b)
très inférieur au seuil de précision courante (10-19 pour Pascal).
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 (Aini, Bini) = (1500, 300):Code:
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
Pièce jointe 419031
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 :mrgreen: ... 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 :furieux:
Pièce jointe 419033
Cela se comprend dès que l'on observe la carte des variations de la fonction étudiée:
Pièce jointe 419035
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(1), 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:
# Pour résumer tout ce qui précède,Code:
2
3
4
5
6
- 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 :aie:
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 (1) des carrés des écarts:
S(a, b) = Somk=27(yk - F(a, b, xk))2 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):
Pièce jointe 420566_Pièce jointe 420570_Pièce jointe 420571
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:
Pièce jointe 420572_Pièce jointe 420576_Pièce jointe 420580
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:
Pièce jointe 420584_Pièce jointe 420590_Pièce jointe 420594
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] _ a0 = 1273.01888 _ b0 = 299.37390
# pour L6: [1505.9 ; 1506.1]×[242.1 ; 242.3] _ a0 = 1506.03815 _ b0 = 242.21727
# pour L5: [1741.9 ; 1742.1]×[199.9 ; 200.1] _ a0 = 1741.93775 _ b0 = 199.94739
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 + ha, b + hb) = S(a, b) + DaS(a, b)*ha + DbS(a, b)*hb + (1/2)*(D2aS(a, b)*ha2 + 2*DabS(a, b)*ha*hb + D2bS(a, b)*hb2) + (ha2 + hb2)*R(ha, hb) ,
dans lequel le reste R(ha, hb) 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:
DaS(a, b) = 0 ; DbS(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)*Smin ;
il vient à la limite:
S(a', b') - S(a, b) = Smin/6 = (1/2)*(D2aS(a, b)*ha2 + 2*DabS(a, b)*ha*hb + D2bS(a, b)*hb2) + (ha2 + hb2)*R(ha, hb)
et pour des écarts suffisamment petits pour que le reste devienne négligeable:
Smin/6 ~ (1/2)*(D2aS(a, b)*ha2 + 2*DabS(a, b)*ha*hb + D2bS(a, b)*hb2) .
On obtient finalement par le changement de coordonnées: X = ha = a' - a ; Y = hb = b' - b :
Smin/3 ~ D2aS(a, b)*X2 + 2*DabS(a, b)*X*Y + D2bS(a, b)*Y2 ,
é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 :aie: .
Pièce jointe 420620_Pièce jointe 420624_Pièce jointe 420628
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:
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).Code:
2
3
4
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]
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:
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.Code:
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Soit (L) une longueur caractérisant le domaine étudié [Amin ; Amax]×[Bmin ; Bmax] ;
on a pour la somme (S) des carrés des écarts S ~ L2 ~ Amax * Bmax (ou Amax2 + Bmax2) ;
pour la norme du gradient: N ~ L2/L = L soit (Amax2 + Bmax2)1/2.
La plus sévère condition d'arrêt (N2 <~ 1E-36 * (Amax2 + Bmax2) 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.
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é:Code:
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
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 = N2min/(Amax2 + Bmax2)1/2 vaut 10-36)
Pièce jointe 421351
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:
Voici ce qui apparaît dans le cas de 12 positions initiales (images 3 et 4), puis 40 (dernière image):Code:
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25PROCEDURE 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;
Pièce jointe 421426_Pièce jointe 421430
Pièce jointe 421492_Pièce jointe 421436
Pièce jointe 421440
