Les vecteurs se notent entre [ ]
.
Attention, les vecteurs sont par notés en ligne, et les
matrices sont simplement des vecteurs de lignes.
Pour entrer une matrice sous xcas, plusieurs options :
donner une fonction de coefficients:
matrix(6,6,(k,l)->1/(1 + k + l));
entrer les lignes à la main ou avec des séquences, en mettant bien les crochets autour:
M:=[[1,2,3,4],[5,6,7,8]];
M:= [ seq( [ seq( 1/(1 + k + l), l = 0..5) ], k = 0..5) ];
(à éviter si possible) créer une matrice nulle et la remplir:
M := matrix(6,6);
for (k:=0; k < 6; k++) {
for (l:=0; l < 6; l++) {
M[k,l] = 1/(1+k+l);
}
}
M;
En sage, on utilisera par exemple:
M = matrix(3, lambda i,j: 1/(1+i+j))
Consulter la documentation xcas de la fonction déterminant. Quelle est l’option qui correspond au pivot de Gauss usuel ? au développement selon les lignes ou les colonnes ?
Pour faire des essais, définir une variable pour la dimension n:=10
et créer une matrice A
de taille n
à coefficients tirés uniformément
dans \([-20,20]\).
Elle servira pour toute la suite. On pourra faire varier n
ou considérer
la matrice Ap = A%1009
pour se placer sur un corps fini.
n := 10;
A := matrix(n,n,(j,k)->rand(41)-20);
ou bien:
A := ranm(n,n)
À l’aide de l’instruction time
, étudier (rapidement) la complexité
des différents algorithmes disponibles en xcas pour des matrices
Dans quels cadres ci-dessus peut-on modéliser le coût d’une multiplication par une constante ?
Par exemple:
1 2 3 4 5 6 7 8 9 10 11 12 13 | complexite(n) := {
local A;
A := ranm(n,n);
[ n, time(det(A,linsolve))[0],
(n<=300)? time(det(A,bareiss))[0] : 0,
(n<=100)? time(det(A,rational_det))[0] : 0,
(n<= 20)? time(det(A,minor_det))[0] : 0
];
}
M := [ seq(complexite(5*2^k),k=0..6) ];
scatterplot(M);
|
Le coût de multiplication est constant pour les flottants et les corps finis. Attention, dans le premier cas la méthode utilisée a une influence énorme sur les propagations d’erreurs lors du calcul (stabilité numérique).
Poser à présent n:=4
, entrer purge(a)
pour être sûr d’avoir une
variable formelle et définir A:=matrix(n,n,(k,l)->a[k,l])
.
Comparer det(A,rational_det)
à det(A, bareiss)
. Comparer également
les façons de développer : simplify
, normal
et expand
.
en sage on pourra écrire matrix(n,n,lambda i,j:var("a_%s_%s"%(i,j)))
Sur cet exemple avec \(n=5\) ou \(6\), comparer le pivot, le pivot de Bareiss et la méthode des mineurs. Conclusion ?
A présent, le résultat est la somme sur toutes les permutations, et l’algorithme le plus efficace est le développement des mineurs ! Le pivot de Gauss donne un résultat particulièrement atroce qui illustre bien la progression du calcul.
Calculer de manière formelle le déterminant de Vandermonde de taille 6.
Calculer la liste des déterminants des matrices \(kI_n-A\) pour \(0\leq k\leq n\).
ly := seq( det(k-A), k, 0, n );
En déduire le polynôme caractéristique de \(A\) en utilisant
la fonction d’interpolation lagrange
.
lx := seq( k, k, 0, n );
normal(lagrange(lx,ly));
Remarque: xacs fait pareil si on demande:
det(x-A,lagrange)
Soient \(x_1,\dots x_n\) des indéterminées.
Il existe plusieurs familles d’expressions invariantes par permutation des \(x_i\), en particulier
les fonction symétriques élémentaires
qui tirent leur nom d”«élémentaires» du fait que tout polynôme \(f(x_1,\dots x_n)\) invariant par permutation des \(x_i\) s’exprime en fonction des \(\sigma_k\) ; (un théorème de Gauss)
les sommes de Newton
\[S_k(x_1,\dots x_n) = \sum_{i=1}^n x_i^k.\]
Ces familles ont une expression élégante en terme de série génératrice : si l’on considère le polynôme
les fonction symétriques élémentaires sont (à un signe près) les coefficients du polynôme développé
tandis que les sommes de Newton sont données par le développement en série
Vérifier les identités ci-dessus sur l’exemple du polynôme
P:=(1-a*T)(1-b*T)(1-c*T)(1-d*T)
.
On pose \(P(x)=1-6x-7x^2+10x^4-9x^5\). Déterminer la liste des complexes \(x_i\) tels que \(P(T)=\prod(1-x_iT)\).
Ce sont les inverses des racines de \(P\), ou les racines du polynôme réciproque
P:=1-6x-7x^2+10x^4-9x^5
R:=col(roots(polrecip(P)),0)
Calculer les 15 premières sommes de Newton de \(P\) à l’aide de ces nombres complexes. Que constate-t-on sur la précision, pourquoi ?
SP:=[ seq( sum(R[j]^k,j,0,4) ,k,0,15) ]
On voit que la précision absolue diminue au fur et à mesure : au début on est proche de la valeur entière, mais les sommes sont de plus en plus grandes et à la fin on ne peut en connaître qu’une valeur approchée. Cela vient du fait que \(P\) a une «grande» racine qui absorbe le comportement des autres.
Calculer les quinze premières sommes de Newton de \(P\) à l’aide d’un développement en série. Qu’en est-il de la stabilité numérique ?
series(diff(P,x)/P,x,0,15)
À présent les sommes sont calculées comme des entiers exacts.
La série génératrice des sommes de Newton sert aussi dans l’autre sens pour retrouver \(P\) : en effet l’équation différentielle se résout en
(on sait que \(P(0)=1\))
Retrouver \(P\) à l’aide de ses 6 premières sommes de Newton:
SN := 6+50x+342x^2+2362x^3+16551x^4+115394x^5;
On utilisera
int( ,x)
pour intégrer des polynômes, series( ,x)
pour obtenir un
développement en serie, et convert( ,polynom)
pour passer d’une série à
un polynôme.
ISN := int(SN,x)
SP := series(exp(-ISN),x)
P := convert(SP,polynom)
A.<x> = QQ[['x']]
SN = 6+50*x+342*x^2+2362*x^3+16551*x^4+115394*x^5
exp(integral(SN+O(x^6)))
Soit \(A\in M_n(\k)\), et \(\chi_A(x)\) son polynôme caractéristique.
Montrer que les sommes de Newton de \(\chi_A\) sont données, pour tout \(k\geq 0\), par
Calculer les dix premières sommes de Newton de \(\chi_A\).
En déduire un algorithme de calcul du polynôme caractéristique. Quelle est sa complexité, lors de quelle étape est-elle atteinte ?
On peut écrire un algorithme en deux étapes :
de sorte que c’est l’étape de calcul des sommes qui est la plus lourde et donne la complexité de l’algorithme.
charpoly_newton(A) = {
my(d = matsize(A)[1]);
my(B = matid(d));
my(S = x*Ser(vector(d,i,B*=A;-trace(B)/i)));
polrecip(Pol(exp(S)));
};
Remarque: il est possible de calculer les traces des \(A^k\) plus rapidement à l’aide d’une stratégie «pas de bébé-pas de géant» : en effet, si l’on décompose \(A^k = A^{aq+r}\), où \(a=\floor{\sqrt{n}}\) et que l’on précalcule les puissances \(A^{aq}\) et \(A^r\), la trace de \(A^k\) nécessite \(O(n^2)\) opérations (seuls les coefficients de la diagonale nous intéressent).
Programmer cet algorithme, vérifier en pratique sa complexité.
Comparer avec l’algorithme ci-dessus pour des matrices à coefficients
dans \(\Z/p\Z\), avec \(p>n\) (la matrice Ap
).
Soit \(A\in M_n(\k)\), on définit \(M=x I_n -A\in M_n(\k[x])\) de sorte que \(\det(M)=\chi_A(x)\) soit le polynôme caractéristique de \(A\).
On dispose de deux équations :
d’une part l’expression du produit matrice-comatrice :
de l’autre la dérivée du déterminant, obtenue en utilisant la multilinéarité (on somme les dérivées partielles selon chaque colonne) :
Si l’on décompose ces deux équations dans \(M_n(\k)[x]\) en introduisant les coefficients \(c_k\in\k\) et \(A_k\in M_n(\k)\)
on obtient les formules de récurrence suivantes
en partant de \(C_0=I_n\).
Utiliser cette méthode pour calculer \(\chi_A(x)\).
Puisque les expressions sont imbriquées (\(C_k\) dépend de \(c_k\) qui dépend de \(C_k\)), on résout la première ligne en \(c_k=-\trace(AC_{k-1})/k\) ce qui invite à décomposer l’itération de la manière suivante
1 2 3 4 5 6 7 8 9 10 11 12 13 | Faddeev(A,x) := {
local n, c, C, P, k;
n := size(A);
C := identity(n);
P := x^n;
for(k:=1;k<n;k++) {
C := A*C;
c := -trace(C)/k;
C := C+c;
P += c*x^(n-k);
}
return(P);
}
|
Note
exercice
La technique de Faddeev et celle de Le Verrier sont en réalité rigoureusement identiques si l’on retrouve itérativement les fonctions symétriques à partir des sommes de Newton.
La vision «Le Verrier» des choses est plus globale et se prête mieux à l’adoption d’algorithmes rapides pour la calcul simultané des traces de \(A^k\) et l’obtention de \(\chi_A\) à partir de ses sommes de Newton.
De son côté, la présentation «Faddeev» fournit gratuitement l’inverse de \(A-XI_n\).
On considère la matrice
A := [[1,1,-1,2,-1],[2,0,1,-4,-1],[0,1,1,1,1],[0,1,2,0,1],[0,0,-3,3,-1]];