Formes quadratiques et Maple

L3 Math : Algèbre et Géométrie  -  Feuille 7 exercice II (copie) -  2007-08

Instructions exécutées avec Maple 10 en ligne de commandes. Traduction pour Maxima ici.

Définition de B2, matrice de B2 relativement à (1,X,X^2,X^3), forme quadratique sur R^4 correspondante

> B2:=(f,g)-> -f(0)*g(0)+int(diff(f(x),x)*diff(g(x),x),x=0..1);
                                           1
                                          /
                                         |   /d      \ /d      \
           B2 := (f, g) -> -f(0) g(0) +  |   |-- f(x)| |-- g(x)| dx
                                         |   \dx     / \dx     /
                                        /
                                          0

> B2(x->x^3,x->x^3);
                                      9/5


> f:=(i,j)->B2(x->x^(i-1),x->x^(j-1));
                                        (i - 1)        (j - 1)
                f := (i, j) -> B2(x -> x       , x -> x       )

> M:=Matrix(4,f);
                            [-1    0     0      0 ]
                            [                     ]
                            [ 0    1     1      1 ]
                            [                     ]
                            [ 0    1    4/3    3/2]
                            [                     ]
                            [ 0    1    3/2    9/5]

> unassign('a','b','c','d'):q:=(Matrix(1,4,[a,b,c,d]).M.Matrix(4,1,[[a],[b],[c],[d]]))[1,1];
              2                   /    4 c   3 d\     /    3 c   9 d\
       q := -a  + (b + c + d) b + |b + --- + ---| c + |b + --- + ---| d
                                  \     3     2 /     \     2     5 /

> # On peut inversement obtenir la matrice de la forme bilineaire symetrique associee a q :
> with(linalg):hessian(q/2,[a,b,c,d]);
                            [-1    0     0      0 ]
                            [                     ]
                            [ 0    1     1      1 ]
                            [                     ]
                            [ 0    1    4/3    3/2]
                            [                     ]
                            [ 0    1    3/2    9/5]

> # B2 est elle degeneree ?
> Rank(M);
                                       4

>

Réduction de Gauss de q(a,b,c,d)

> var:=[a,b,c,d];
                              var := [a, b, c, d]

> q1:=q;var1:=var;liste:=[]:
               2                   /    4 c   3 d\     /    3 c   9 d\
       q1 := -a  + (b + c + d) b + |b + --- + ---| c + |b + --- + ---| d
                                   \     3     2 /     \     2     5 /

                             var1 := [a, b, c, d]

> printlevel:=2:
> while var1<>[] do
aa:=diff(q1,var1[1],var1[1])/2;l:=eval(diff(q1,var1[1]),var1[1]=0);
# Si q1 ne depend pas de var1[1], on supprime var1[1] de la liste var1
if aa=0 and l=0 then var1:=subsop(1=NULL,var1);next end if;
# Il y a un terme en (var1[1])^2
if aa<>0 then liste:=[op(liste),[aa,var1[1]+l/(2*aa)]];q1:=eval(q1,var1[1]=0)-l^2/(4*aa);var1:=subsop(1=NULL,var1);next end if;
i:=2;while eval(diff(l,var1[i]),var1[i]=0)=0 do i:=i+1 end do;
aa:=diff(q1,var1[i],var1[i])/2;
# Il y a un terme en var1[1]*var1[i] et en (var1[i])^2
if aa<>0 then liste:=[op(liste),[aa,var1[i]+l/(2*aa)]];q1:=eval(q1,var1[i]=0)-l^2/(4*aa);var1:=subsop(i=NULL,var1);next end if;
# Il y a un terme en var1[1]*var1[i] mais pas de terme en (var1[1])^2 ni en (var1[i])^2
aa:=diff(q1,var1[1],var1[i]);l1:=eval(diff(q1,var1[1]),[var1[1]=0,var1[i]=0]);li:=eval(diff(q1,var1[i]),[var1[1]=0,var1[i]=0]);
liste:=[op(liste),[aa/4,var1[1]+var1[i]+(l1+li)/aa],[-aa/4,var1[1]-var1[i]+(li-l1)/aa]];
q1:=eval(q1,[var1[1]=0,var1[i]=0])-l1*li/aa;
var1:=subsop(1=NULL,i=NULL,var1)
end do;
                                   aa := -1

                                    l := 0

                              liste := [[-1, a]]

                                /    4 c   3 d\     /    3 c   9 d\
          q1 := (b + c + d) b + |b + --- + ---| c + |b + --- + ---| d
                                \     3     2 /     \     2     5 /

                               var1 := [b, c, d]

                                    aa := 1

                                l := 2 c + 2 d

                      liste := [[-1, a], [1, b + c + d]]

                                                               2
                    /4 c   3 d\     /3 c   9 d\     (2 c + 2 d)
              q1 := |--- + ---| c + |--- + ---| d - ------------
                    \ 3     2 /     \ 2     5 /          4

                                var1 := [c, d]

                                   aa := 1/3

                                    l := d

                                                           3 d
              liste := [[-1, a], [1, b + c + d], [1/3, c + ---]]
                                                            2

                                          2
                                         d
                                  q1 := ----
                                         20

                                  var1 := [d]

                                  aa := 1/20

                                    l := 0

                                                      3 d
         liste := [[-1, a], [1, b + c + d], [1/3, c + ---], [1/20, d]]
                                                       2

                                    q1 := 0

                                  var1 := []



> liste[];
                                                 3 d
              [-1, a], [1, b + c + d], [1/3, c + ---], [1/20, d]
                                                  2

> q1:=0;for i from 1 by 1 to nops(liste) do q1:=q1+liste[i,1]*liste[i,2]^2 end do:q1;
                                         /    3 d\2
                                         |c + ---|      2
                      2              2   \     2 /     d
                    -a  + (b + c + d)  + ---------- + ----
                                             3         20


> # Verification :
> simplify(q1);
                                              2              2
                  2    2                   4 c            9 d
                -a  + b  + 2 b c + 2 b d + ---- + 3 c d + ----
                                            3              5

> simplify(q);
                                              2              2
                  2    2                   4 c            9 d
                -a  + b  + 2 b c + 2 b d + ---- + 3 c d + ----
                                            3              5

Avec q1:=a*b+b*c+a*c et var1:=[a,b,c] on obtiendrait à l'exécution :

                            q1 := a b + b c + a c

                               var1 := [a, b, c]

                                    aa := 0

                                  l := b + c

                                    i := 2

                                    aa := 0

                                    aa := 1

                                    l1 := c

                                    li := c

                 liste := [[1/4, a + b + 2 c], [-1/4, a - b]]

                                           2
                                   q1 := -c

                                  var1 := [c]

                                   aa := -1

                                    l := 0

             liste := [[1/4, a + b + 2 c], [-1/4, a - b], [-1, c]]

                                    q1 := 0

                                  var1 := []



On peut faire de cette suite d'instructions une fonction qui prend comme argument la forme quadratique et ses variables dans l'ordre choisi et qui rend la liste des formes lineaires et leur coefficients :

Gauss:=proc(quad,variable) local q,var,liste,i,a,l;
q:=quad;var:=variable;liste:=[];
while var<>[] do
 a:=diff(q,var[1],var[1])/2;l:=eval(diff(q,var[1]),var[1]=0);
 if a=0 and l=0 then var:=subsop(1=NULL,var);next end if;
 if a<>0 then liste:=[op(liste),[a,var[1]+l/(2*a)]];q:=eval(q,var[1]=0)-l^2/(4*a);var:=subsop(1=NULL,var);next end if;
 i:=2;while eval(diff(l,var[i]),var[i]=0)=0 do i:=i+1 end do;
 a:=diff(q,var[i],var[i])/2;
 if a<>0 then liste:=[op(liste),[a,var[i]+l/(2*a)]];q:=eval(q,var[i]=0)-l^2/(4*a);var:=subsop(i=NULL,var);next end if;
 a:=diff(q,var[1],var[i]);l1:=eval(diff(q,var[1]),[var[1]=0,var[i]=0]);li:=eval(diff(q,var[i]),[var[1]=0,var[i]=0]);
 liste:=[op(liste),[a/4,var[1]+var[i]+(l1+li)/a],[-a/4,var[1]-var[i]+(li-l1)/a]];
 q:=eval(q,[var[1]=0,var[i]=0])-l1*li/a;
 var:=subsop(1=NULL,i=NULL,var)
end do;
q:=add(liste[i,1]*liste[i,2]^2,i=1..nops(liste));
[liste,q] end proc;

> # Test (les calculs a l'interieur de la procedure Gauss n'apparaissent plus a l'ecran) :
> unassign('x','y','z');
> r:=Gauss(x*y+y*z+x*z,[x,y,z]):print('r[1]'=r[1]);print(x*y+y*z+x*z=r[2]);
              r[1] = [[1/4, x + y + 2 z], [-1/4, x - y], [-1, z]]

                                              2          2
                                 (x + y + 2 z)    (x - y)     2
               x y + y z + x z = -------------- - -------- - z
                                       4             4

> # Verification
> is(x*y+y*z+x*z=r[2]);
                                     true

>

Valeurs propres de la matrice de B2 : on sait que cette matrice est diagonalisable sur R

> with(LinearAlgebra):P:=CharacteristicPolynomial(M,X);
                           4   47  3   57  2   19
                     P := X  - -- X  - -- X  + -- X - 1/60
                               15      20      15

> # Maple ne trouve que -1 comme racine rationnelle de P, de multiplicite 1 :
> roots(P,X);
                                   [[-1, 1]]

> #localisation exacte des racines reelles de P (cf l'aide Maple '? realroots')
> realroot(P,1/10);
                                                61
                [[0, 1/16], [5/16, 3/8], [15/4, --], [-1, -1]]
                                                16

> # Avec des intervalles plus petits :
> realroot(P,1/100);
                                41   21    485  243
               [[1/128, 1/64], [---, --], [---, ---], [-1, -1]]
                                128  64    128  64

> # Resolution numerique approchee :
> fsolve(P);
                 -1., 0.01357893318, 0.3232993539, 3.796455046

> # Comparer avec l'instruction 'eigenvals' : (I designe le nombre complexe i)
> evalf(eigenvals(M),5);
                      -5
  -1., 3.7965 - 0.7 10   I, 0.01357 - 0.000082605 I, 0.32329 + 0.000090605 I

Dans certains cas, on connaît les valeurs propres exactes de la matrice de la forme quadratique. On peut alors trouver une base orthonormée pour le produit scalaire canonique et orthogonale pour q. Exemple avec q(x,y,z)=xy+yz+xz :

> q:=x*y+y*z+x*z;
                             q := x y + y z + x z

> with(VectorCalculus):M:=Hessian(q/2,[x,y,z]);
                                [ 0     1/2    1/2]
                                [                 ]
                           M := [1/2     0     1/2]
                                [                 ]
                                [1/2    1/2     0 ]

> eigenvals(M);
                                 1, -1/2, -1/2

> v:=eigenvectors(M);
                       v := [1, 1, {[1, 1, 1]}], [-1/2, 2, {[-1, 1, 0], [-1, 0, 1]}]

> # orthonormalisation des vecteurs propres dans chaque espace propre :
> w:=op(v[1,3] union v[2,3]);
                                  w := [1, 1, 1], [-1, 1, 0], [-1, 0, 1]

> n1:=GramSchmidt(v[1,3],normalized);
                                               [ 1/2   1/2   1/2]
                                               [3     3     3   ]
                                        n1 := {[----, ----, ----]}
                                               [ 3     3     3  ]

> n2:=GramSchmidt(v[2,3],normalized);
                            [   1/2  1/2     1/2  1/2   1/2  1/2]  [   1/2   1/2   ]
                            [  3    2       3    2     3    2   ]  [  2     2      ]
                     n2 := {[- ---------, - ---------, ---------], [- ----, ----, 0]}
                            [      6            6          3    ]  [   2     2     ]


> # Matrice des coordonnees des vecteurs de la nouvelle base n1 union n2 :
> P:=Matrix(3,(i,j)->(n1 union n2)[j][i]);
                                        [ 1/2       1/2  1/2       1/2]
                                        [3         3    2         2   ]
                                        [----    - ---------    - ----]
                                        [ 3            6           2  ]
                                        [                             ]
                                        [ 1/2       1/2  1/2      1/2 ]
                                   P := [3         3    2        2    ]
                                        [----    - ---------     ---- ]
                                        [ 3            6          2   ]
                                        [                             ]
                                        [ 1/2      1/2  1/2           ]
                                        [3        3    2              ]
                                        [----     ---------       0   ]
                                        [ 3           3               ]


> # P est orthogonale et diagonalise M
> P^%T.P;
                                               [1    0    0]
                                               [           ]
                                               [0    1    0]
                                               [           ]
                                               [0    0    1]

> Diag:=P^%T.M.P;
                                                [1     0       0  ]
                                                [                 ]
                                        Diag := [0    -1/2     0  ]
                                                [                 ]
                                                [0     0      -1/2]

> # Coordonnees de (x,y,z) dans la nouvelle base en utilisant P^%T=P^(-1) :
> X:=P^%T.Matrix(<x,y,z>);
                                  [         1/2      1/2      1/2           ]
                                  [        3    x   3    y   3    z         ]
                                  [        ------ + ------ + ------         ]
                                  [          3        3        3            ]
                                  [                                         ]
                                  [   1/2  1/2      1/2  1/2      1/2  1/2  ]
                             X := [  3    2    x   3    2    y   3    2    z]
                                  [- ----------- - ----------- + -----------]
                                  [       6             6             3     ]
                                  [                                         ]
                                  [               1/2      1/2              ]
                                  [              2    x   2    y            ]
                                  [            - ------ + ------            ]
                                  [                2        2               ]

> # q est la somme des Diag[i,i]*X[i,1]^2 (comparer avec X^%T.Diag.X) :
> q1:=add(Diag[k,k]*X[k,1]^2,k=1..3);
                                      /   1/2  1/2      1/2  1/2      1/2  1/2  \2   /   1/2      1/2  \2
                                      |  3    2    x   3    2    y   3    2    z|    |  2    x   2    y|
        / 1/2      1/2      1/2  \2   |- ----------- - ----------- + -----------|    |- ------ + ------|
        |3    x   3    y   3    z|    \       6             6             3     /    \    2        2   /
  q1 := |------ + ------ + ------|  - -------------------------------------------- - --------------------
        \  3        3        3   /                         2                                  2

> # Verification
> is(q=q1);
                                                   true



F-X Dehon, 5 dec. 2007, http://math.unice.fr/~dehon