Cherche coup de main

Question

Bonjour à tous,
Je cherche un coup de main pour m'aider à programmer un algorithme en VB6. 
Il s'agit du mds de Torgerson. On peut le trouver dans cet article
http://acoustics.aau.dk/~fw/mds03.pdf 

A la page 10 l'auteur présente les 4 étapes de l'algorithme et à la page 11 il donne un exemple.
J'ai fait les étapes 1 et 2. Je retrouve bien la matrice B de la page 11, mais je coince sur la troisième étape lorsqu'il s'agit de trouver les 2 plus grandes valeurs propres positives.

Quelqu'un pourrait-il m'aider ?
Merci

Flachy Joe · Answer

Salut, voila ce qu'on obtient avec Maple :
_______________________________________________________________________
Initialisation :

> reset:
> with (linalg):

_______________________________________________________________________
Matrice P :

> P:= linalg[matrix](4,4,[
> 0   ,93 ,82  ,133 ,
> 93  ,0  ,52  ,60  ,
> 82  ,52 ,0   ,111 ,
> 133 ,60 ,111 ,0   
> ]);

                         [  0    93     82    133]
                         [                       ]
                         [ 93     0     52     60]
                    P := [                       ]
                         [ 82    52      0    111]
                         [                       ]
                         [133    60    111      0]

_______________________________________________________________________
Matrice P2 :
> P2:=map(x->x^2,P);

                     [    0    8649     6724    17689]
                     [                               ]
                     [ 8649       0     2704     3600]
               P2 := [                               ]
                     [ 6724    2704        0    12321]
                     [                               ]
                     [17689    3600    12321        0]

_______________________________________________________________________
Matrice J :

> J:=linalg[matrix](4,4):
> for i from 1 to 4 do
>  for j from 1 to 4 do
>   J[i,j] := -0.25 + `if`(i=j,1,0):
>  od:
> od:
> J:=eval(J);

                    [0.75     -0.25    -0.25    -0.25]
                    [                                ]
                    [-0.25    0.75     -0.25    -0.25]
               J := [                                ]
                    [-0.25    -0.25    0.75     -0.25]
                    [                                ]
                    [-0.25    -0.25    -0.25    0.75 ]

_______________________________________________________________________
Matrice B :

> B:= map(x->-x/2,multiply(J,P2,J));

  B :=

        [5035.062500 , -1553.062500 , 258.9375000 , -3740.937500]

        [-1553.062500 , 507.8125000 , 5.312500000 , 1039.937500]

        [258.9375000 , 5.312500000 , 2206.812500 , -2471.062500]

        [-3740.937500 , 1039.937500 , -2471.062500 , 5172.062500]

_______________________________________________________________________
                      Equation caractéristique :

> I4:=Matrix(4,4,shape=identity):
> det(B-lambda*I4)=0;

                                  10                         8       2                     3         4
  6000. - 0.112490345 10   lambda + 0.3120950838 10  lambda  - 12921.75000 lambda  + lambda  = 0

_______________________________________________________________________
    Valeurs propres, multiplicités et vecteurs propres associés :

> eigenvectors(B);

  [9724.167602, 1,{[-0.6371596899, 0.1866206893, -0.2531171389, 0.7036561391]}]

  [36.59655948, 1,{[0.006882387287, -0.8181761912, 0.4324934135, 0.3788004753]}

                  -5
  [0.2199994167 10  , 1,{[0.4999999997, 0.5000000345, 0.4999999817, 0.4999999843]}]

  [3160.985840, 1,{[-0.5864982205, 0.2139168147, 0.7063152123, -0.3337338072]}]

Je suis pas sur que ça s'affiche bien, il faudrat sans doute que tu décale les lignes des exposants pour les ramener sur les bons termes...

Il n'y a donc pas de probleme, juste un polynome de degre 4 à résoudre...

 ;)  Flachy Joe  ;)

Flachy Joe · Answer

En effet, ça c'est vraiment affiché n'importe comment : l'equation caractéristique est 
 6000. - 0.112490345*10^10*lambda + 0.3120950838*10^8 * lambda^2 - 12921.75000*lambda^3 + lambda^4 = 0

 ;)  Flachy Joe  ;)

kaputzan · Answer

Merci Flachy Joe, je suis carrément impressionné ! Que dis-je ...époustouflé
Mais aussi .. . pas plus avancé qu'avant la bataille. Car je cherche à résoudre ce type de problème pour n'importe quelle matrice en entrée et surtout . . . tout ça en VB . . . 
Si tu as une petite idée, aussi époustouflante, que ta demo . . . je suis preneur.

Flachy Joe · Answer

Salut !

Ce n'est que des produits matriciel, qui ne devrait pas poser de probleme si tu as fait un chouillat de math post bac. (dans vb, il faut traiter ça comme tableau a 2d et programmer avec les formules indicielles).
Le seul (et gros) probleme c'est de trouver les racines du polynome du quatrieme ordre qui donne les valeurs propres, il faut que tu cherche un algo capable de te le faire, perso je ne sais pas. Tu peux toujours utiliser une méthode numérique pour obtenir des valeurs approchées...
(algo de Newton, etc) Si les valeurs approchées te conviennent, je pense pouvoir programmer ça, faut voir...

 ;)  Flachy Joe  ;)

Flachy Joe · Answer

Voila, j'ai fait quelque recherches, je n'ai pas codé la recherche de toutes les racines, mais si ça te pose problème, n'hesite pas...

Algo de Newton :
_____________________________________________________________________

'Ces fonction utilisent des polynomes dont les coéfficient sont stockés dans un
'tableau. Le coefficient de x^i doit être stocké dans la ieme valeur du tableau.
' par exemple :
' p(x) = 1-x+2*x^2 + 5*x^6
'devient :
' dim p(6) as double : p(0)=1 : p(1)=-1 : p(2)=1 : p(6)=5

'Ne pas utiliser de tableaux commençant à l'indice 1 :
Option Base 0

'évaluer la valeur du polynome au point x
Public Function CalculP(Polynome() As Double, ByVal x As Double) As Double
    Dim Result As Double
    For i = 0 To UBound(Polynome)
        Result = Result + Polynome(i) * x ^ i
    Next
    CalculP = Result
End Function

'créer le polynome dérivé du polynome initial
Public Sub Derive(Polynome() As Double, Result() As Double)
    ReDim Result(UBound(Polynome) - 1)
    For i = 1 To UBound(Polynome)
        Result(i - 1) = i * Polynome(i)
    Next
End Sub

'chercher une racine de p(x) a partir de x0
Public Function Newton(x0 As Double, Polynome() As Double, Precision As Double) As Double
    Dim d() As Double
    Derive Polynome, d
    If CalculP(d, x0) = 0 Then _
        Err.Raise 99999, "Algo de Nawton", _
        "Impossible d'utiliser ce point comme initialisation"
    Dim xk As Double
    xk = x0
    k = 0
    While Abs(CalculP(Polynome, xk)) > Precision
        xk = xk - CalculP(Polynome, xk) / CalculP(d, xk)
    Wend
    Newton = xk
End Function
_____________________________________________________________________

Pour trouver toutes les racines d'un polynomes, voila la marche à suivre :

On considere un polynome de degré N
On cherche ses N-2 premieres derivés, la dérivé (N-2)ieme est un polynome d'ordre 2 dont on trouve facilement les 2 racines exactes.

En choisissant des points à droite du 1er puis entre le premier et le 2d puis a gauche du 2d comme valeur initiale pour la méthode de Newton (on prendra par exemple x1+m, m, x2-m avec m=(x1-x2)/2), on obtient les 3 racines de la derivé (N-3)ieme.

On continu jusqu'a arriver à la dervié 0ieme, le polynome initial, qui à N racines (racines eventuellement multiples, dans ce cas 2 points d'une des derivées feront convergé l'algo vers le même zero).

L'algo de Newton n'est pas trés rapide (encore moins comme je l'ai programmé, avec des tableaux) mais ç'est suffisant pour des polynomes de degré de l'ordre de 100.

 ;)  Flachy Joe  ;)

kaputzan · Answer

Merci pour ton aide,
Super ! Je vais regarder cela d'un peu plus près
En fait j'étais parti sur la méthode QR pour trouver les valeurs propres mais je trouve des valeurs approximatives qui ne sont pas assez prècises. 
J'essaye ta méthode

Cherche coup de main

6 réponses

Discussions similaires