Code de regression polynomiale: [Résolu]

Messages postés
5
Date d'inscription
jeudi 2 juillet 2009
Dernière intervention
15 juillet 2009
- 7 juil. 2009 à 14:51 - Dernière réponse :
Messages postés
2
Date d'inscription
jeudi 22 avril 2010
Dernière intervention
25 mai 2010
- 25 mai 2010 à 09:09
Bonjour à tous !
je suis nouveau sur le fofo, on m'a dit qu'il y avait des pros qui trainaient dans le coin, alors je suis venu voir.

Voilà,
je suis ingénieur et j'essaye d'utiliser un programme de regression polynomiale pour trouver un polynome  approchant la position d'un ballon stratosphèrique en fonction du temps car la précision du capteur gps ne me convient pas.
Voici ce que j'aimerai approximer:
il se trouve que le programme me sort une approximation mais pas de tous les points, il s'arrête avant la fin.
Qu'en pensez vou ?
j'avoue que j'ai du mal à comprendre le programme.
(Le code est en bas du message)

<col style=\"width: 60pt;\" width=\"80\" /><col style=\"width: 62pt;\" width=\"83\" />----
<table x="" style="border-collapse: collapse; width: 120pt;" border="0" cellpadding="0" cellspacing="0" width="160"><col style="width: 60pt;" width="80" /><col style="width: 60pt;" width="80" />----, Temps
sec </td>Alt, ----
, , ----
0, 200, ----
60, 351, ----
120, 498, ----
240, 791, ----
300, 944, ----
360, 1094, ----
420, 1431, ----
600, 1764, ----
660, 1913, ----
720, 2070, ----
780, 2232, ----
840, 2395, ----
900, 2527, ----
960, 2659, ----
1020, 2820, ----
1080, 2990, ----
1140, 3156, ----
1260, 3469, ----
1320, 3621, ----
1380, 3768, ----
1440, 3918, ----
1500, 4072, ----
1560, 4220, ----
1620, 4363, ----
1680, 4513, ----
1740, 4672, ----
1800, 4856, ----
1860, 5060, ----
1920, 5257, ----
1980, 5443, ----
2040, 5619, ----
2100, 5791, ----
2160, 5967, ----
2220, 6159, ----
2280, 6359, ----
2340, 6550, ----
2400, 6734, ----
2460, 6915, ----
2520, 7100, ----
2580, 7305, ----
2700, 7781, ----
2760, 8028, ----
2820, 8248, ----
2880, 8427, ----
2940, 8585, ----
3000, 8775, ----
3060, 8998, ----
3120, 9219, ----
3180, 9421, ----
3240, 9617, ----
3300, 9844, ----
3360, 10101, ----
3420, 10395, ----
3480, 10673, ----
3540, 10906, ----
3600, 11112, ----
3660, 11329, ----
3720, 11550, ----
3780, 11772, ----
3840, 11940, ----
3900, 12062, ----
3960, 12247, ----
4020, 12526, ----
4080, 12826, ----
4140, 13102, ----
4200, 13358, ----
4260, 13613, ----
4320, 13875, ----
4380, 14127, ----
4440, 14345, ----
4500, 14555, ----
4560, 14825, ----
4620, 15135, ----
4680, 15437, ----
4740, 15707, ----
4860, 16180, ----
4920, 16270, ----
4980, 16284, ----
5040, 16275, ----
5100, 16277, ----
5220, 16254, ----
5340, 16324, ----
5400, 16304, ----
5460, 16295, ----
5520, 16279, ----
5580, 16297, ----
5640, 16338, ----
5700, 16316, ----
5820, 16281, ----
5880, 16292, ----
5940, 16319, ----
6000, 16327, ----
6120, 16271, ----
6180, 16258, ----
6360, 16285, ----
6420, 16270, ----
6480, 16299

</td><td style="width: 62pt;" width="83">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" height="17">
</td><td>
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr><tr style="height: 12.75pt;" height="17"><td style="height: 12.75pt;" x="" align="right" height="17">
</td><td class="xl22" x="">
</td></tr></tbody></table>
voici le code que j'ai trouvé sur ce site:

 Option Explicit

 Public x() As Double 'Données d'entrée (xi et yi)
 Public y() As Double
 Public ai_p() As Double 'Valeurs des coefficients du polynôme de tendance

 Sub main()

 Dim i As Integer, j As Integer
 Dim ligne As String
 Dim virg_deb As Long
 Dim virg_fin As Long

 Dim a As String

 i = 0
 'Boucle de récupération des données xi et yi à traiter
 Open "valeurs_init.txt" For Input As #1
 While Not EOF(1)
 i = i + 1
 ReDim Preserve x(1 To i)
 ReDim Preserve y(1 To i)

 Line Input #1, ligne

 'Remplacement des blancs par des séparateurs (chr(47) = "/")
 ligne = Replace(ligne, Chr(9), Chr(47))
 ligne = Replace(ligne, Chr(32), Chr(47))
 ligne = Chr(47) & ligne & Chr(47)

 'Récupération des données x(i) et y(i)
 'Recherche des séparateurs (/) et stockage des valeurs indiquées entre les séparateurs
 virg_deb = InStr(1, ligne, Chr(47))
 virg_fin = InStr((virg_deb + 1), ligne, Chr(47))
 x(i) = CDbl(Mid(ligne, (virg_deb + 1), (virg_fin - virg_deb - 1)))

 virg_deb = InStr((virg_fin), ligne, Chr(47))
 virg_fin = InStr((virg_deb + 1), ligne, Chr(47))
 y(i) = CDbl(Mid(ligne, (virg_deb + 1), (virg_fin - virg_deb - 1)))
 Wend

 Close #1

 'Appel de la procédure de régression polynômiale :
 'le premier argument est le degré du polynôme de tendance chosi, x() et y() sont les données à traiter
 'ai_p() est un argument de retour renvoyant les coef ai du polynôme
 Call regression_polynomiale(12, x(), y(), ai_p())

 Dim xf() As Double
 Dim yf() As Double
 ReDim xf(1 To UBound(x))
 ReDim yf(1 To UBound(y))

 'Boucle générant un fichier texte contenant les valeurs xf(i) et yf(i) approchées
 Open "valeurs_fin.txt" For Output As #1
 For i = 1 To UBound(x)
 xf(i) = x(i)
 yf(i) = ai_p(0)
 For j = 1 To (UBound(ai_p()))
 yf(i) = yf(i) + ai_p(j) * xf(i) ^ j
 Next j

 Print #1, xf(i) & Chr(9) & Format(yf(i), "####0.000")
 Next i

 End Sub

 Sub regression_polynomiale(deg_pol As Integer, xi() As Double, yi() As Double, coef_pol() As Double)
 
 Dim i As Integer, j As Integer
 Dim k As Integer, l As Integer
 Dim nb_pts As Integer
 Dim mat_A() As Double, mat_B() As Double
 Dim mat_At() As Double, mat_Bt() As Double
 Dim dim_mat As Integer

 dim_mat = deg_pol + 1

 'Nombre de points
 nb_pts = UBound(xi)
 If nb_pts <> UBound(yi) Then MsgBox "Problème, il n'y a pas le même nombre de valeurs pour les x et les y"

 ReDim mat_A(1 To dim_mat, 1 To dim_mat)
 ReDim mat_B(1 To dim_mat)
 ReDim mat_At(1 To dim_mat, 1 To dim_mat)
 ReDim mat_Bt(1 To dim_mat)
 ReDim coef_pol(0 To deg_pol)

 ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 '' Etape 1 :
 '' Mise en équation du problème de minimisation par la méthode des moindres carrés
 '' on cherche à minimiser : E(a0,..., ap) = somme(yi-somme(aj*xi^j)²) où p est le degré du pôlynome
 '' les dérivées partielles de E par rapport à chacun des coef du polynôme de tendance doivent être nulles,
 '' c'est une condition nécessaire !
 ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 ' les xi et yi sont transmis en paramètres
 ' la procédure renvoie alors en sortie 2 matrices mat_A et mat_B tq :
 ' mat_A * vect_des_aj = mat_B
 Call mise_en_éq(nb_pts, deg_pol, xi(), yi(), mat_A(), mat_B())

 ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 '' Etape 2 :
 '' On réalise alors la triangulation de la mat_A et les opérations équivalentes sur mat_B
 ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 ' mat_A et mat_B sont données en arguments d'entrée
 ' mat_At et mat_Bt sont renvoyées en arguments de sortie, après triangulation de mat_A
 Call triangulation(dim_mat, mat_A(), mat_B(), mat_At(), mat_Bt())

 ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 '' Etape 3 :
 '' Cette procédure résout alors le système et renvoie les coefficients aj du pôlynome de tendance
 ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 'les matrices mat_At et mat_Bt sont données en arguments d'entrée
 'les coef aj sont renvoyées en arguments de sortie sous le nom de variable coef_pol
 Call resolution(dim_mat, mat_At(), mat_Bt(), coef_pol())

 End Sub

 Sub mise_en_éq(n As Integer, p As Integer, vect1() As Double, vect2() As Double, mat1() As Double, mat2() As Double)

 Dim i As Integer, j As Integer, k As Integer

 For k = 0 To p ' on balaie les lignes
 For j = 0 To p ' on balaie les colonnes
 mat1(k + 1, j + 1) = 0
 For i = 1 To n ' on somme tous les xi^(k+j)
 mat1(k + 1, j + 1) = mat1(k + 1, j + 1) + vect1(i) ^ (k + j)
 Next i
 Next j
 Next k

 For k = 0 To p ' on balaie les lignes
 mat2(k + 1) = 0
 For i = 1 To n ' on somme tous les (yi*xi^(k))
 mat2(k + 1) = mat2(k + 1) + (vect2(i) * vect1(i) ^ (k))
 Next i
 Next k

 End Sub

 Sub triangulation(n As Integer, mat1() As Double, mat2() As Double, mat3() As Double, mat4() As Double)

 Dim i As Integer, j As Integer, k As Integer
 Dim mat3k() As Double, mat4k() As Double
 ReDim mat3k(1 To n, 1 To n, 1 To n) As Double
 ReDim mat4k(1 To n, 1 To n) As Double

 For i = 1 To n
 mat4k(i, 1) = mat2(i)
 For j = 1 To n
 mat3k(i, j, 1) = mat1(i, j)
 Next j
 Next i

 For k = 1 To (n - 1)
 For i = 1 To n
 For j = 1 To n
 If i <= k Then
 mat3k(i, j, (k + 1)) = mat3k(i, j, k) ' on conserve les mêmes valeurs
 mat4k(i, (k + 1)) = mat4k(i, k)
 Else
 mat4k(i, (k + 1)) = mat4k(i, k) - (mat3k(i, k, k) * mat4k(k, k) / mat3k(k, k, k))
 If j <= k Then
 mat3k(i, j, (k + 1)) = 0
 Else
 mat3k(i, j, (k + 1)) = mat3k(i, j, k) - (mat3k(i, k, k) * mat3k(k, j, k) / mat3k(k, k, k))
 End If
 End If
 Next j
 Next i
 Next k

 For i = 1 To n
 mat4(i) = mat4k(i, n)
 For j = 1 To n
 mat3(i, j) = mat3k(i, j, n)
 Next j
 Next i
 End Sub

 Sub resolution(n As Integer, mat3() As Double, mat4() As Double, coef_pol() As Double)

 Dim i As Integer, j As Integer, k As Integer
 Dim xsol() As Double
 Dim som_aijxj As Double
 ReDim xsol(1 To n)

 xsol(n) = mat4(n) / mat3(n, n)
 coef_pol(n - 1) = xsol(n)

 For i = (n - 1) To 1 Step (-1)
 som_aijxj = 0
 For j = i + 1 To n
 som_aijxj = som_aijxj + mat3(i, j) * xsol(j)
 Next j
 xsol(i) = 1 / mat3(i, i) * (mat4(i) - som_aijxj)
 coef_pol(i - 1) = xsol(i)
 Next i

 End Sub
Afficher la suite 

Votre réponse

10 réponses

Meilleure réponse
Messages postés
106
Date d'inscription
lundi 11 avril 2005
Dernière intervention
16 juillet 2010
- 10 juil. 2009 à 09:19
3
Merci
Bonjour,


Si je comprends bien ton problème, tu as un ballon sonde doté d'un récepteur GPS. A l'intant t, tu enregistres la date et l'altitude de ton ballon sonde.


Et c'est là que je ne suis pas sur de comprendre. Tu cherches à interpoler la position de ton ballon n'importe quand entre ton instant zéro et ta dernière mesure, où alors tu cherches à savoir où est ton ballon entre ta dernière mesure et la suivante que tu ne connais pas encore.


Pour la première solution, je ne m'embêterais même pas avec les moindres carrés pour chercher un polynôme. Je ferai une FFT. On s'en fout des fréquences et des longueurs d'onde. Tu auras simplement une fonction numérique correpsondant à tous tes points de mesures. Quelle que soit l'allure de ta courbe en fonction de la météo, etc, ton calcul fonctionnera, alors que tel ou tels polynômes avec les moindres carrés...


Pour la deuxième solution, à partir de tes deux dernières mesures, je calculerais la dernière tendance de ta courbe. Une approximation linéaire en attendant le dernier enregistrement.


A+


Thomas.










Marin Marais

Merci marinmarais 3

Avec quelques mots c'est encore mieux Ajouter un commentaire

Codes Sources a aidé 88 internautes ce mois-ci

Commenter la réponse de marinmarais
Meilleure réponse
Messages postés
106
Date d'inscription
lundi 11 avril 2005
Dernière intervention
16 juillet 2010
- 10 juil. 2009 à 12:15
3
Merci
OK Florent, j'ai pigé.


Dans ce cas là, avant tout, je me demande si tu ne devrais pas utiliser un autre logiciel de calcul. Tu sais utiliser MatLab ? Tu pourrais résoudre ton problème en 5 min avec ce soft si tu es relativement familier avec des calculs matriciel usuels.


Je vois toujours deux solutions :


- La FFT. Tu appliques une FFT avec toutes tes mesures. Si tu veux limiter lisser ton signal, et bien tu ne prends pas en compte les plus petites fréquences calculées. Si tu dois programmer ça en VB, tu vas passer des plombes. Commande MatLab : Y=FFT(X)... Ton résultat est dans l'ensemble des complexes mais tu as peu de boulot pour obtenir le résultat réel de la forme :


altitude(t)= a0 + [somme de 1 à n](a_i*sin(omega_i*t+Phase_i)


Tu pourras aisément dériver ton résultat pour étudier les forces en jeu et déduire ton coefficient de frottement.


- La méthode des moindres carrés : Tu ne connais pas, a priori, le degré optimal de ton polynôme. Je te conseille, avant d'implémenter une solution générale pour tous tes ballons, d'effectuer plusieurs ajustements en étudiant les résidus et surtout, en regardant à partir de quel degré de polynôme ceux-ci suivent la loi normale. C'est la méthode que j'ai adoptée pour optimiser le choix de polynômes de calibration de capteurs. Sans même le programmer, tu peux le faire par Excel, mais je maintiens qu'avec MatLab, c'est beaucoup moins laborieux, pour peu que tu sois familier de la méthode des moindres carrés.


Je vais essayer de te faire un exemple sur Excel avec les données que tu as envoyées.


Par contre, ton histoire de coefficient de friction qui varie, ça m'intéresse. Je fais des recherches en topo de précision et j'étudie de nouveaux systèmes de mesures dont un projet de système de mesure de différences d'altitude par un laché de bille. Je mesure le temps nécessaire au déplacement permettant, par une intégration et connaissant les forces en jeu, de déduire la dénivellée. C'est dans le cas d'un transfert de géométrie entre la surface et le fond d'un tunnel. Aurais-tu des références sur ce point ?


A+ et courage,


Thomas.








Marin Marais

Merci marinmarais 3

Avec quelques mots c'est encore mieux Ajouter un commentaire

Codes Sources a aidé 88 internautes ce mois-ci

Commenter la réponse de marinmarais
Messages postés
2117
Date d'inscription
lundi 11 avril 2005
Dernière intervention
14 mars 2016
- 7 juil. 2009 à 23:15
0
Merci
Bonsoir,

Bien compliqué tout cela...
Perso, j'arrive à :
y = 0,00023124*x^2 + 2,12946809*x + 342,20515405
pour la partie courbe, et à :
y=16286,1363636364
pour l'horizontale...

Reste à calculer le point de rencontre... pas très compliqué...

Je ne sais pas si c'est la réponse attendue, ou si tu veux qu'on reprenne un programme qui ne convient pas... mais dans ce cas, je dirais qu'il faut demander à l'auteur...

Ah, oui... C'est une utilisation de base d'Excel...

Amicalement,
Us.
Commenter la réponse de us_30
Messages postés
5
Date d'inscription
jeudi 2 juillet 2009
Dernière intervention
15 juillet 2009
- 10 juil. 2009 à 11:02
0
Merci
Salut Thomas,
tu as  bien compris pour ce qui est des ballons, en effet j'étudie leur déplacement.
En fait, j'ai déja les position en fonction du temps, je veux juste lisser un peu ça sans perdre les déplacemnts incertains du ballons, chose que j'étudie à vrai dire...
Pourquoi?
car je cherche à optimiser ces courbes pour calculer les meilleurs Cx(coeff de friccion de l'air sur le ballon) et Reynodls, nombre adimensionnel qui dépend du milieu et de la vitesse de déplacement du ballon dans celui ci. Ces renseignements me servent à étudier la mise au plafond des ballons qui parfois dépassent leur plafond car ils arrivent trop vite et éclatent... (le Cx diminue trop avant la mise au plafond)
La majeur partie du temps lorsqu'ils dépassent leur plafond, ils redescendent car l'enveloppe se comprime trop et fait redescendre le ballon qui va osciller autour de l'altitude du plafond pour s'y stabiliser.
C'est pas si facile à expliquer ^^
Pour le moment j'ai des résultats corrects. Je cherche à sortir le polynome calculé.
Le créateur du programme m'a conseillé de je cite :

Pour ta deuxième question, dans la sub main, tu peux
récupérer les coefficients du polynôme, après la
ligne

Call
regression_polynomiale(6, x(), y(), ai_p())

Les coefficients sont les valeurs du vecteur
ai_p()

 

Pour les récupérer, il suffit par exemple de remplacer
la dernière boucle par celle-ci-dessous :

For j = 1 To
(UBound(ai_p()))

   
Print #1, Format(ai_p(j), "####0.000")

Next
i

Je vais essayer de faire ce qu'il dit mais je ne comprend pas trop ce que Print #1, Format(ai_p(j), "####0.000") veut dire ^^

Merci de l'interêt que vous portez à mon sujet

Florent
Commenter la réponse de Nukix
Messages postés
106
Date d'inscription
lundi 11 avril 2005
Dernière intervention
16 juillet 2010
- 13 juil. 2009 à 15:58
0
Merci
Bonjour Florent,


J'ai jeté un coup d'oeil à tes données. Vu la rupture nette à 16100m, une régression polynomiale va te donner des erreurs importantes à cet endroit. Si l'altitude en fonction du temps avait été à peu près parabolique, ça aurait pu jouer mais là...


Bref, j'ai fait une FFT vite fait sur MatLab. Ca marche pas mal... Mon erreur max après modélisation est de 266m (écart-type sur les résidus de 82m). Bah ouais, la FFT ça marche bien sauf aux extrémités.


Si ça te tente, je peux t'envoyer tout ça par mail (fichier excel de résultat et le code MatLab de mon traitement. Dans ce cas, donne moi ton adresse mail par mp.


A+,


Thomas








Marin Marais
Commenter la réponse de marinmarais
Messages postés
2
Date d'inscription
lundi 10 novembre 2003
Dernière intervention
16 juillet 2009
- 13 juil. 2009 à 23:02
0
Merci
bonjour ,
tous roule impec pour moi avec ce programme.
le seul point un peu zarbi est l'appel a des valeurs dans un fichier texte dans la procédure de lancement main
ci dessous comment j'ai modifé la procédure main  pour
-charger les données a étudier dans une feuille excel
-choisir les degrés du polynome de regression

 Sub main()
Dim degpol As Integer
degpol = 2 ' degrés demandé du polynome de régression  , ici 2 ( si un on aura la une droite )
 'Appel de la procédure de régression polynômiale :
 'le premier argument est le degré du polynôme de tendance choisi, x() et y() sont les données à traiter
 'ai_p() est un argument de retour renvoyant les coef ai du polynôme
 nbp = 25 ' nombre de points couple x y
 ReDim x(nbp)
 ReDim y(nbp)
 ReDim ai_p(degpol)
 For i = 1 To nbp
 x(i) = Sheets("prob").Cells(i, 1) ' charge les données x dans une feuille appelée prob de la ligne 1 a 25 colonne 1 à adapter
 y(i) = Sheets("prob").Cells(i, 2) ' charge les données y dans une feuille appelée prob de la ligne 1 a 25 colonne 2 à adapter
  Next
 Call regression_polynomiale(degpol, x(), y(), ai_p()) ' passe les x , y et l'ordre du polynome au programme de calcul
' au retour les coefficients de regression  sont le tableau ai_p()  de taille degpol
 
 End Sub

sans trop de difficulté il doit être possibles de transformer ces programmes en fonction
si on pouvait avoir en rab le coeff de coreelation ce serait super
Commenter la réponse de hogsvoor
Messages postés
5
Date d'inscription
jeudi 2 juillet 2009
Dernière intervention
15 juillet 2009
- 15 juil. 2009 à 10:05
0
Merci
Ca y est Thomas je t'ai envoyé un message. Peux tu expliciter ce que tu cherches, en quoi je peux t'aider sur ton taff.

hogsvor, génial ce que tu as fait, tu peux préciser où mettre ça dans le programme?

à+
Commenter la réponse de Nukix
Messages postés
106
Date d'inscription
lundi 11 avril 2005
Dernière intervention
16 juillet 2010
- 15 juil. 2009 à 10:58
0
Merci
Salut !


En utilisant la FFT dans ton cas, il ne s'agit pas d'obtenir un spectre en fréquence. Fourier a inventé la transformée juste pour modéliser des fonctions, afin qu'elles soient facilement dérivables pour résoudre des équa diff...


Ici, le but est le même. Tu as l'altitude de ton ballon en fonction du temps. Tu cherches une fonction (pas trop compliquée si possible) modélisant au mieux le mouvement dans le temps de ton ballon.


Et bien tu obtiens comme modèle une somme de cosinus...


Le calcul que j'ai fait, sur MatLab, se déroule en 2 étapes.


- Je ré-échantillonne pour avoir un pas régulier de 60 s et un nombre pair de données


- Je lance la FFT sur les données ré-échantillonnées


Voila le code MatLab. Je suis pas un pro de ce langage, ça peut sans doute être amélioré, mais ça marche !


%**************************************************************************


%* Modelisation de l'altitude d'un ballon sonde par la FFT *


%**************************************************************************


% Thomas Touze


% CERN BE-ABP-SU


% 13 juillet 2009


% 1- Import des donnees


% Les donnees du ballon sonde (97 enregistrements) sont regroupes en deux


% vecteurs de 97 lignes nommes temps et altitude.


% 2- Reechantillonnage des donnees


% Le but est d'obtenir des enregistrements avec un pas regulier de 60 s


m=97;


T=temps(end,1);


n=T/60;


n=n+1;


alti=zeros(n,1);


seconde=zeros(n,1);


alti(1,1)=altitude(1,1);


alti(n,1)=altitude(m,1);





for i=1:n


seconde(i,1)=60*(i-1);


for j=2:97


if seconde(i,1)==temps(j,1)


alti(i,1)=altitude(j,1);


elseif temps(j-1,1)<seconde(i,1) && seconde(i,1)< temps(j,1)


alti(i,1)=altitude(j-1,1)+(seconde(i,1)-temps(j-1))*(altitude(j,1)-altitude(j-1,1))/(temps(j,1)-temps(j-1,1));


end


end


end





% clear m altitude temps i j





% Un nombre pair de donnees est plus pratique pour la FFT


dimension=size(alti);


dimension=[dimension(1,1)];





if mod(dimension,2)==1


seconde=[seconde(1:n-1,1)];


alti=[alti(1:n-1,1)];


n=n-1;


else


T=T+60;


end





% 3- FFT





N=n/2;


Z=fft(alti);


A=zeros(N-1,1);


P=zeros(N-1,1);





Ao=((conj(Z(1,1))*Z(1,1))^0.5)/(2*N);





for i=1:N-1


A(i,1)=((conj(Z(i+1,1))*Z(i+1,1))^0.5)/N;


P(i,1)=mod(-angle(Z(i+1,1))-2*pi*i/T,2*pi);


end





fourier=zeros(n,1);


residus=zeros(n,1);


% le vecteur fourier represente la modelisation de l'altitude. Le vecteur residu represente l'erreur de modelisation





for i=1:n


toto=Ao;


for j=1:N-2


toto=toto+A(j,1)*cos(2*pi*j*seconde(i,1)/T-P(j,1));


end


fourier(i,1)=toto;


residus(i,1)=alti(i,1)-fourier(i,1);


end





% 4 - Analyse residus modelisation par la fft





toto=0;


for i=1:n


toto=toto+residus(i,1);


end


Vm=toto/n;





toto=0;


for i=1:n


toto=toto+(residus(i,1)-Vm)^2;


end


Sm=(toto/(n-1))^0.5;


Voili voilou...



Thomas.

Marin Marais
Commenter la réponse de marinmarais
Messages postés
2
Date d'inscription
lundi 10 novembre 2003
Dernière intervention
16 juillet 2009
- 16 juil. 2009 à 00:41
0
Merci
pour nukix et les autres
j'ai fait un fichier exemple téléchargeable à l'adresse ci dessous
normalement çà marche ,
http://rapidshare.com/files/256258820/regresspolynom.xls.html



il n'y pas de limite au nombre de couples de données ( 32000 dans excel 2003 et infini dans excel 2007 sauf la longueur des calculs
ni au nombre de degrés du polynome , mais vu ma (très) longue expérience de modélisateur du vivant , quand un phénomène de régression a deux variables n'est pas toujours pas correctement modélisable par un polynome qui dépasse le 4 ième degré , il y a bien peu de chance que les deux variables soient corrélées ( les mouvements d'un couple de danseur de tangou tango font exception à cette règle)

Dans vba le ficher macro lancé par calcul est la procédure calcul()
j'ai essayé de documenter le truc

précision : n'ayant jamais eu le moindre cours d'excel et de visual basic j'utilise mon propre patois de programmation, désolé pour les pro et les puriste si j'ai un accent et une syntaxe pourris.
mais comme disait le vieux Deng Xiao Ping en 1978 , peu importe la couleur du chat l'important c'est qu'il chope les souris.

nb nikix t'a trouvé où les programmes de calcul parce que celui ou celle qui fait çà , il a fait l'essentiel du travail ce serait bien de lui faire savoir que ça sert à d'autres.
Commenter la réponse de hogsvoor
Messages postés
2
Date d'inscription
jeudi 22 avril 2010
Dernière intervention
25 mai 2010
- 25 mai 2010 à 09:09
0
Merci
Bonjour.
Je fais des maths sur les nuages de points, dont les nuages plans comme le tien.
Il faut que la fonction numérique f(x) qui représente ton nuage ait un sens physique. Il faut donc lui imposer des contraintes, par exemple sur les coordonnées d'au moins un point du nuage.
Ex:f(x)=7.4447664402*10^-15*x^5-0.000000000149*x^4+0.000000926419*x^3-0.00203323893*x^2+4.1391583705*x-39.4144584247 ne convient pas car pour x=0, ton ballon est enterré à une profondeur de 39.4144584247 m.(x est le temps).Somme des résidus yi-f(xi)=14986.23..
On peut choisir f(x)avec ln.exp.Arctan. ect...
Je te propose un autre polynôme qui satisfait aux contraintes des moindres carrés avec contrainte sur trois points du nuage:(0.200) (4860.16180) (6480.16299)ici.
f(x)=3.42288845778*10^-21*x^7-6.863274651*10^-17*x^6+5.25877769937*10^-13-0.000000001966*x^4+0.000003784328*x^3-0.003400684559*x^2+3.63474227891*x+200. Somme des résidus yi-f(xi)=15834.42.
J'ai pris trois points caractéristiques de ton nuage, dans cet exemple. Si tu as d'autres nuages en dim 2, dim 3..tu peux me les montrer si je peux te citer pour publication dans mon travail. Amicalement. Bernard.
Commenter la réponse de alusipdranreb

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.