Algorithmique (Agrégation interne)Bernard.Parisse@ujf-grenoble.fr2010/11 |
N.B.: ce texte est disponible en ligne sur
www-fourier.ujf-grenoble.fr/~parisse/agregint.html
et en PDF sur
www-fourier.ujf-grenoble.fr/~parisse/agregint.pdf
Dans cette section, on passe en revue les notions de bases en algorithmique au programme (B.2), en commencant par travailler en ligne de commande (notion de variable, affectation), puis en complexifiant au fur et à mesure la ligne de commande (test, boucle) et on termine par l’écriture de fonctions. On pourra trouver d’autres tutoriels introduisant l’algorithmique via la géométrie ou l’arithmétique sur la page algorithmique seconde de Xcas
www-fourier.ujf-grenoble.fr/~parisse/algoxcas.html
La syntaxe sera donnée en Xcas, avec les éléments de syntaxe permettant d’adapter à d’autres languages utilisables à l’oral (Python, Scilab, Casio Classpad, TI Nspire CAS, Maple, Maxima), on ne reprendra par contre pas Free Pascal qui n’a pas d’interpréteur en mode interactif et se prête donc moins à la démarche proposée ici, ni Carmetal dont la version 3 permet d’écrire des programmes en javascript mais qui est trop spécialisé en géométrie par rapport aux thèmes abordés ici (sauf dans la section 1.2.2), ni bien sur les logiciels de géométrie ne permettant pas de programmer (en-dehors de macro-constructions). Quelques remarques sur ces logiciels et langages :
Pour stocker une donnée en mémoire, on utilise une variable. Une variable possède un nom, qui est composé d’au moins un caractère alphabétique (de A à Z ou de a à z) suivi ou non d’autres caractères alphanumériques, par exemple a, var, b12. Xcas est sensible à la différence entre majuscules et minuscules (c’est aussi le cas de Python, Maple, Maxima, Scilab mais ce n’est pas le cas chez TI ni Casio). Une variable peut contenir n’importe quel type de donnée.
Pour donner une valeur à une variable
:= suivi par
la valeur. Par exemple a := 1.2 (même syntaxe chez TI, Maple ;
avec interversion et le signe ⇒ chez Casio,
en utilisant = en Python et Scilab, avec : en Maxima).
saisir(a)
ou encore saisir(a,b) ou encore pour faciliter la saisiesaisir("Entrez votre nom",a,"Entrez votre age",b).a=input("Entrez une valeur"),
chez CasioInput "A=",A
Pour utiliser la valeur d’une variable
a:=2 puis a+1
on obtient 3. En Xcas, TI, Casio, Maple, Maxima
une variable non affectée n’est pas remplacée
et reste symbolique (calcul formel), en Python et Scilab, cela
déclenche une erreur.
afficher (en Maple, Maxima, Xcas
print("texte",variable),
en Python 2 print "texte",variable,
en Scilab disp(variable), chez
Casio Print ou PrintNatural).
Cela permet d’afficher un résultat intermédiaire dans une
zone située avant la ligne contenant le résultat du programme.
Par exemple afficher(a+1) affichera
3 en bleu dans cette zone.
Xcas peut manipuler différents types de données dont :
1, 2/3, 1.1
(on peut aussi utiliser la notation scientifique mantisse, exposant
comme dans 1e-7)
x, t, z,
x+1, x^2+y^2,
f:=x->x^2 ou f(x):=x^2,
s:="bonjour".
On accède à un caractère d’une chaine en donnant le nom de la chaine
indicié par un entier compris entre 0 et la taille de la chaine moins
1 (même convention qu’en langage C), par exemple s[0]
renvoie "b".
a:=[1,2,3].
On accède à un élément d’une liste en donnant le nom de la liste
indicié par un entier compris entre 0 et la taille de la liste moins
1 (même convention qu’en langage C), par exemple a[0] renvoie
1. On peut créer des listes de taille donnée avec une formule
de remplissage (par exemple l:=[0$10], l:=[j^2$j=1..10],
voir aussi seq, makelist, ranm)
Maple, Maxima et Python peuvent manipuler des entiers,
flottants, chaines de caractères,
et listes, avec les mêmes délimiteurs que Xcas. TI et Casio aussi,
avec des délimiteurs {} au lieu de [] pour
différencier les listes des vecteurs. Pour adresser
un élément d’une liste, on utilise le nom de variable
suivi de l’indice entre [], les indices
commencent à 0 en Xcas, Scilab, Python et à 1
en Maple, Maxima, TI, Casio. Attention, Scilab ne propose
pas de type entier (on peut représenter un entier par un flottant
sans erreur de représentation s’il est inférieur à 253).
Il faut
être bien conscient du type de données que l’on manipule,
par exemple en Python 7/2 renvoie 3 le quotient euclidien
des deux entiers, alors que 7./2.
renvoie 3.5 (le quotient des deux flottants)
Exercices :
saisir le montant en
euros, on convertir en francs et on affiche le résultat avec
afficher)
Une instruction valide dans Xcas doit suivre une syntaxe précise qui
ressemble à l’écriture d’une expression en mathématique. Les
différents éléments d’une instruction simple peuvent être des
nombres (entiers, réels, flottants), des noms de variables, des
opérateurs (par exemple +), des appels à des fonctions
(l’argument ou les arguments se trouvant entre les parenthèses
séparés par une virgule s’il y a plusieurs arguments) : par exemple
sqrt(3) pour désigner la racine carrée de 3 ou
irem(26,3) pour désigner le reste de
la division euclidienne de 26 par 3).
Il faut
prendre garde à la priorité entre les différentes opérations. Par
exemple dans 1+2*3 l’opération * est effectué
avant + comme en mathématiques. On peut utiliser
des parenthèses pour forcer l’ordre des opérations, par exemple
(1+2)*3.
Exercice : Calculer
5/2/3, 5/(2/3), 5/2*3, 5/(2*3),
5^2*3, 5^(2*3),5*2^3, (5*2)^3
Conclure sur les priorités relatives de la multiplication,
de la division et de la puissance.
Attention, la multiplication doit toujours
être indiquée explicitement contrairement aux mathématiques.
Ainsi l’écriture ab ne désigne pas le produit de 2 variables
a et b mais une variable dont le nom est ab.
Lorsqu’on veut exécuter plusieurs instructions à la suite, on
termine chaque instruction par le caractère ; (ou par
les caractères :; si on ne souhaite pas afficher le résultat
de l’instruction). En Python il suffit de
passer à la ligne, en Scilab, Maxima, Maple on utilise ;,
chez TI et Casio on utilise :.
Cette section est spécifique à Xcas, les autres systèmes n’ayant en général pas d’instruction géométrique interagissant comme du calcul numérique. Les habitués de Carmetal pourront noter une certaine similarité dans la programmation géométrique des deux systèmes (les instructions ont en général le même nom, avec une majuscule dans Carmetal et la syntaxe est très proche du C dans les 2 languages).
Toutes les instructions du menu Geo ont un résultat graphique
(à l’exception des instructions du menu Mesure), par exemple :
point(1,2) affiche le point de coordonnées 1 et 2,
droite(A,B) la droite passant par deux points A et
B définis auparavant.
Lorsqu’une ligne de commande contient une instruction
graphique, le résultat est affiché dans un repère 2-d ou 3-d selon
la nature de l’objet généré. On peut controler
le repère avec les boutons situés à droite du graphique,
par exemple orthonormaliser avec le bouton ⊥.
Si une ligne de commande contient des instructions graphiques et non
graphiques, c’est la nature de la dernière instruction qui décide
du type d’affichage.
On peut donner des attributs graphiques aux objets graphiques
en ajoutant à la fin de l’instruction graphique
l’argument affichage=..., argument dont la
saisie est facilitée par le menu Graphic->Attributs.
L’instruction A:=point(click()) permet de
définir une variable contenant un point du plan que l’on clique
avec la souris (click() renvoie le nombre complexe correspondant).
C’est un peu l’analogue géométrique
de l’instruction saisir.
Lorsqu’on veut voir sur un même graphique le résultat de plusieurs
lignes de commandes, on crée une figure avec le menu
Geo->Nouvelle figure->graph geo2d, la figure correspond alors
aux lignes de commande situées à gauche. Dans ce mode, on peut
créer les objets graphiques en ligne de commande ou pour certains
directement à la souris en sélectionnant un Mode et
en cliquant. On peut aussi déplacer un point en mode Pointeur
et observer les modifications des objets géométriques qui en
dépendent.
Exemples :
milieu.
A:=point(click()); B:=point(click()); xm:=(abscisse(A)+abscisse(B))/2; ym:=(ordonnee(A)+ordonnee(B))/2; C:=point(xm,ym)
A:=point(click()); D:=droite(x=abscisse(A))
Exercices (cf. le menu Geo pour les commandes de géométrie de Xcas) :
point, soit
en passant en mode point et en cliquant 3 points à la souris).
Afficher le centre du cercle
circonscrit du triangle formé par ces 3 points en
utilisant uniquement les instruction mediatrice et
inter_unique.
On peut tester l’égalité de 2 expressions en utilisant l’instruction
==, alors que != teste si 2 expressions ne sont pas
égales. On peut aussi tester l’ordre entre 2 expressions
avec <, <=, >, >=, il s’agit
de l’ordre habituel sur les réels pour des données numériques
ou de l’ordre lexicographique pour les chaines de caractères.
En Python et Scilab, les tests sont identiques, en Maple, Maxima on utilise
= pour l’égalité, chez TI et Casio
on teste l’égalité avec = ou ≠ et l’inégalité
avec ≤, ≥.
Un test renvoie 1 ou true s’il est vrai, 0 ou false s’il est faux.
On peut combiner
le résultat de deux tests au moyen des opérateurs logiques &&
(ou et ou and), || (ou ou ou or),
et on peut calculer la négation logique
d’un résultat de test avec ! (ou not).
En Python, on utilise les symboles &&, ||
ou !. En Maple, Maxima, Scilab, TI, Casio,
and, or et not.
On utilise ensuite souvent la valeur du test pour exécuter une instruction
conditionnelle comme l’instruction si :
si condition alors bloc_vrai sinon bloc_faux fsi;
(cela peut aussi servir pour arrêter une boucle en cours d’exécution).
if(test): passer à la ligne puis mettre le bloc
vrai indenté, et else: passer à la ligne puis
mettre le bloc faux indenté.
if condition then bloc_vrai else bloc_faux end.
if condition then bloc_vrai else bloc_faux fi;.
if condition then (bloc_vrai) else (bloc_faux):If condition Then : action1 : Else : action2: EndIf (Classpad IfEnd).
if (condition) { bloc_vrai } else { bloc_faux }.
Par exemple, on pourrait stocker la valeur absolue d’un réel x dans y par :
si x>0 alors y:=x; sinon y:=-x; fsi;
(on peut bien sur utiliser directement y:=abs(x)).
Exercices :
droite d’équation y=2x+1, puis le point
O origine, puis un point A quelconque, faire afficher si
A est du même côté de la droite que O. De même pour O
un point quelconque. Faire bouger le point A à la
souris (en passant en Mode pointeur) pour vérifier
les différentes possibilités.
distance
pour avoir la distance entre 2 sommets).
On peut exécuter des instructions plusieurs fois de suite en utilisant une boucle définie (le nombre d’exécutions est fixé au début) ou indéfinie (le nombre d’exécutions n’est pas connu). On utilise en général une variable de controle (indice de boucle ou variable de terminaison).
pour ... de ... jusque ... faire ... fpourf:=1; pour j de 1 jusque 10 faire f:=f*j; fpour;f:=1; pour j de 2 jusque 10 pas 2 faire f:=f*j; fpour;i comme indice
de boucle, Xcas ne pourra pas utiliser ailleurs sa valeur
prédéfinie comme étant √−1.
repeter ... jusqu_a ...repeter saisir("Entrer un nombre entre 1 et 10",a); jusqu_a a>=1 && a<=10;
tantque ... faire ... ftantquetantque b!=0 faire r:=irem(a,b); a:=b; b:=r; ftantque
Remarques :
for(init;condition;incrementation){ instructions }if (...) break;) dont
l’usage peut éviter l’utilisation de variables de contrôle
compliquées, et le passage immédiat à l’itération suivante
(mot-clef continue) qui évite des forêts d’if.
res:=NULL avant la boucle,
et on écrit dans la boucle res:=res,objet_graphique.
Notez que les objets graphiques intermédiaires sont de toutes facons affichés
dans la fenêtre DispG (menu Cfg->Montrer).
range,for var in range(start,stop):while condition: et le bloc
indenté.
for var from debut to fin do bloc; od;while condition do instructions od;.
for var:debut thru fin do (bloc)
for var=debut:fin; bloc; end;while condition do instructions end;.
:For I,A,B : action : EndFor (TI)For A⇒ I To B : action : Next (Classpad)
Exercices :
rotation).
Pour réaliser des taches un peu plus complexes, il devient
intéressant de
les subdiviser en plusieurs entités indépendantes
appelées fonctions, qui
permettent d’étendre les possibilités des instructions de Xcas.
Une fonction peut avoir 0, 1 ou plusieurs
arguments (comme la fonction racine carrée sqrt() en a un).
Elle calcule une valeur, appelée valeur de retour.
La définition de fonctions peut parfois se faire directement avec une formule (fonction définie algébriquement) mais nécessite parfois un algorithme plus complexe. C’est peut-être à cet endroit qu’il faut s’arrêter en classe de seconde (au moins en terme d’exigences), en réservant les fonctions plus complexes (et les variables locales) à la première scientifique et les fonctions récursives à la classe de terminale, conjointement à l’étude de la récurrence.
Il s’agit de fonctions définies par une formule algébrique.
Exemples :
f(x):=x^2+1f:=x->x^2+1,
TI et CasioDefine f(x)=x^2+1),
abs, on tape :absolu(x):= si x<0 alors -x sinon x fsi;
Exercices :
max2
calculant le maximum de 2 entiers (sans utiliser max)
puis de 3 entiers en appelant la fonction précédente ou sans
appeler la fonction
précédente. Observer l’intérêt d’utiliser la fonction max2.
Attention
Il faut faire la différence entre fonction et expression.
Par exemple a:=x^2-1 définit une expression alors que
b(x):=x^2-1 définit une fonction.
On peut donc écrire b(2)
qui renverra 4 mais l’analogue avec a serait subst(a,x=2).
On peut écrire factor(a) mais l’analogue avec b sera
factor(b(x)).
On peut composer des fonctions, par exemple la fonction c=b ∘ b
est c:=b@b et aussi construire une fonction à partir d’une
expression avec l’instruction unapply par exemple la
fonction b
définie par l’expression a est b:=unapply(a,x).
La plupart des fonctions ne peuvent pas être définies simplement
par une formule algébrique. On doit souvent calculer des données
intermédiaires, faire des tests et des boucles. Il faut alors
définir la fonction par une suite d’instructions,
délimitées par { ... }.
La valeur calculée par la fonction est alors la valeur calculée
par la dernière instruction ou peut être explicitée en utilisant le
mot-clef return suivi de la valeur à renvoyer, par exemple si
une fonction calcule plusieurs valeurs on peut les renvoyer dans
une liste ou une séquence.
Notez que l’exécution
de return met un terme à la fonction même s’il y a encore
des instructions après.
f:=proc(parametres) ...; end;f(parametres):=(bloc)
function nom_retour=f(parametres); ...; endfunctionnom_retour.
def f(parametres):, on saisit le bloc
définissant la fonction indenté où on renvoie la valeur de retour
par return.
Pour éviter que les données intermédiaires n’interfèrent
avec les variables de la session principale,
on utilise un type spécial de variables,
les variables locales, dont la valeur ne peut être
modifiée ou accédée
qu’à l’intérieur de la fonction. On utilise à cet effet le mot-clef
local suivi par les noms des variables locales séparés par
des virgules. En Python, toute variable intermédiaire est considérée
comme locale, sauf indication contraire (avec global).
Comme le texte définissant une telle fonction ne tient en général pas
sur une ou deux lignes, il est commode d’utiliser un éditeur
de programmes pour définir une fonction. Pour cela, on utilise
le menu Prg->Nouveau programme de Xcas. Le menu
Prg->Ajouter
facilite aussi la saisie des principales commandes
de programmation.
Exemple : encore le PGCD mais sous forme de fonction
pgcd(a,b):={
local r;
tantque b!= faire
r:=irem(a,b);
a:=b;
b:=r;
ftantque
return a;
}
On clique ensuite sur le bouton OK, si tout va bien, le programme
pgcd est défini et on peut le tester dans une ligne de commande
par exemple par pgcd(25,15).
On peut exécuter en mode pas à
pas un programme et visualiser l’évolution de la valeur
des variables avec l’instruction debug(),
ce qui peut servir à comprendre le déroulement d’un
algorithme, mais aussi à corriger un programme erroné
(un analogue existe en Scilab et Python, selon l’environnement de programmation).
Exercice :
On peut au cours du déroulement d’une fonction faire appel à cette même fonction avec un ou des arguments “plus simples”, de sorte qu’après un nombre fini d’appels récursifs, l’argument ou les arguments sont devenus triviaux, et on peut alors renvoyer le résultat. Par exemple, l’algorithme d’Euclide peut s’énoncer sous la forme
On peut le traduire en Xcas par
pgcdr(a,b):=si b==0 alors a sinon pgcdr(b,irem(a,b)); fsi;
Dans certains cas, la récursivité permet de simplifier grandement la conception d’un algorithme, par exemple pour le calcul du reste de an par un entier fixé m, en distinguant n pair et n impair et en utilisant pour n>0 pair :
| an (mod m ) = (an/2 (mod m ))2 (mod m ) |
et pour n>1 impair :
| an (mod m ) = (a × (a(n−1)/2 (mod m ))2) (mod m ) |
Exercice : implémenter le calcul de an (mod m ) de cette manière.
Attention à bien s’assurer que le programme récursif se termine. En particulier quand on teste un programme récursif, il est conseillé de commencer par tester le cas où la récursion doit se terminer. Notez aussi que Xcas limite par défaut à 25 le nombre de récursions.
Attention aussi, un algorithme récursif peut être très
inefficace, par exemple si on calcule les termes de la suite
de Fibonacci par la formule
F(n):=si n<2 alors 1; sinon F(n-1)+F(n-2); fsi;
alors le calcul de F5 nécessite le calcul de F4 et F3 mais
le calcul de F4 demande lui-même le calcul de F3 et F2,
on calcule donc deux fois F3, trois fois F2, cinq fois F1.
(exercice : le nombre total de calcul est lui-même une suite de Fibonacci!).
Il est plus efficace dans ce cas d’écrire une boucle.
Pour cette section, on pourra se référer au manuel de programmation de Xcas. Pour les instructions prédéfinies, utiliser le menu CAS, Arithmétique ou le menu Cmds, Entier et Cmds, Polynômes.
Exposés : 103, 104, 106, 157, 159. Exercices : 302, 303, 304, 305, 306, 309
Application :écriture d’un entier en base b.
Écrire un algorithme de PGCD itératif ou/et récursif. Même question pour Bézout. Algorithme itératif pour a et b :
Adapter pour calculer l’inverse d’un entier modulo un autre entier (on peut diminuer la longueur des listes de 1).
Test de primalité par division. Recherche des nombres premiers inférieurs à un entier donné (crible).
Recherche de a (mod ∏pi ) connaissant a (mod pi ). On peut commencer par 2 nombres premiers (et utiliser Bézout).
Rappel du principe de codage RSA :
| (ac (mod n ))s (mod n ) = (ac (mod n ))s (mod n ) = a (mod n ) (1) |
Exercice 1: Générer une paire de clefs
Génèrez deux nombres premiers p et q>256 au hasard, en utilisant par exemple les fonctions
nextprime et rand. Calculez n=p × q puis ϕ(n)=(p−1)(q−1)
puis choisissez une clef secrète s inversible modulo ϕ(n) et calculez son inverse c.
Vérifiez sur quelques entiers la propriété (1), on utilisera la fonction
powmod(a,c,n) pour calculer ac (mod n ).
Exercice 2: Codage et décodage d’un message
On transforme une chaine de caractère en une liste d’entiers et réciproquement
avec asc et char. Pour l’appliquer à une liste l, on peut utiliser
map(l,powmod,c,n).
En utilisant la paire de clefs de l’exercice 1,
codez un message puis décodez ce message pour vérifier.
Décodez le message authentifié situé à l’URL
http://www-fourier.ujf-grenoble.fr/~parisse/mat249/rsa1.
Exercice 3: Puissance modulaire rapide
Pour pouvoir crypter des messages de cette manière, il est nécessaire d’avoir une
fonction calculant ac (mod n ) rapidement. Comparer le temps de calcul de
a^c mod n et powmod(a,c,n) pour quelques valeurs de a, c, n (on
pourra utiliser time(instruction) pour connaitre le temps d’exécution d’une
instruction. Pour comprendre cela, programmez une fonction puimod calculant
ac (mod n ) en utilisant et en justifiant un des algorithmes
Exercice 4: Attaque simple
L’utilisation de 256 valeurs possibles pour a se prête à une attaque très simple :
la personne souhaitant décoder un message codé avec une clef publique
sans en connaitre la clef secrète calcule simplement
la liste des ac (mod n ) pour les 256 valeurs possibles de a et compare au message.
Décodez de cette manière le message situé à l’URL
http://www-fourier.ujf-grenoble.fr/~parisse/mat249/rsa2
Exercice 5: Groupement de lettres
Pour parer à cette attaque, on va augmenter le nombre de valeurs possibles de a pour que
le calcul de la liste de toutes les puissances des a possibles soit trop long. Pour cela, on
groupe par paquets de x caractères et on associe à un groupe de caractères
l’entier correspondant en base 256. Par exemple, si on prend des groupes de x=3 caractères,
"ABC" devient 65*256^2+66*256+67 car le code ASCII de A, B, C est respectivement
65, 66, 67.
Donner une condition reliant n et x pour que le décodage redonne le message original.
Choisissez une paire de clefs vérifiant cette condition pour x=8.
Ecrire un programme de codage et de décodage avec groupement (on commencera par compléter
le message original par des espaces pour qu’il soit un multiple de 8 caractères,
l’instruction size permet de connaitre la taille d’une chaine de caractères,
on pourra utiliser la fonction d’écriture en base de la feuille d’exercices 2).
Exercice 6: Sécurité du codage 1
Montrer que la connaissance
de ϕ(n) et de n permet de calculer p et q par résolution d’une
équation de degré 2.
La sécurité du codage repose donc sur la difficulté de factoriser n. Tester sur des entiers
de taille croissante le temps nécessaire au logiciel pour factoriser p et q.
Une valeur de n de taille 128 bits, 512 bits, 1024 bits parait-elle suffisante?
Exercice 7: Sécurité du codage 2
Le choix de c et de s est aussi important. Pour le comprendre, prenons p=11 et q=13.
Représentez pour différentes valeurs de c les points (a,ac (mod n )), plus
le dessin obtenu est aléatoire, plus il sera difficile à une personne mal intentionnée
de déchiffrer un message sans connaitre la clef.
On pourra utiliser les instructions seq pour générer une suite de terme général
exprimé en fonction d’une variable formelle, et scatterplot(l) qui représente
le nuage de points donné par une liste l de couples de coordonnées.
Observez en particulier les cas où c n’est pas premier avec ϕ(n) et c=3.
Conclusions ?
Exercice 8: Primalité
Pour créer une paire de clefs, il faut générer des nombres premiers et donc être
capable de déterminer si un nombre est premier ou non. Un test simple consiste à appliquer
la contrapposée du petit théorème de Fermat: si ap ≠ a (mod p ), alors p n’est
pas premier. Ecrire une fonction prenant en argument a et p et renvoyant 0 si
p n’est pas premier et 1 si ap =a (mod p ), puis une fonction prenant en argument
p et effectuant le test pour toutes les valeurs de a∈[2,p−1] jusqu’à ce
que le test échoue. Si tous le tests réussissent, p est peut-être premier mais ce n’est pas
certain: existe-t-il un nombre non premier pour lequel tous les tests réussissent ?
Ce test utilise le fait que si p est premier, ℤ/pℤ est un corps et pour a≠ 0, ap−1=1 (mod p ). On écrit p premier différent de 2 sous la forme 2t s avec s impair. Comme l’équation x2=1 n’admet que 1 et -1 comme racines modulo p, pour a fixé on peut avoir soit at=1 (mod p ), sinon en prenant t fois le carré modulo p on doit prendre la valeur -1. Si le test échoue, p n’est pas premier. Si le test réussit, p n’est pas forcément premier, mais on peut montrer qu’il y a au plus 1 entier sur 4 pour lequel il réussit. On reprend donc le test pour quelques autres valeurs de a jusqu‘à être raisonnablement sur que p est premier (il existe aussi des certifications de primalité).
Calcul du quotient et du reste de la division de 2 polynômes dont les coefficients sont donnés dans une liste. Exemple d’application : élimination de racines (pour trouver toutes les racines d’un polynôme)
Application : inverse de l’écriture en base b.
Programmation de la méthode de Horner (fonction horner de Xcas)
Il s’agit d’évaluer efficacement un polynôme
| P(X) = an Xn + ... + a0 |
en un point. On pose b0=P(α ) et on écrit :
| P(X)−b0=(X−α )Q(X) |
où :
| Q(X) = bn Xn−1 + ... +b2 X + b1 |
On calcule alors par ordre décroissant bn, bn−1, ..., b0.
horn effectuant ce calcul:
on donnera en arguments le polynôme sous forme de la
liste de ces coefficients (dans l’exemple [1,2,0,-1,5]) et la
valeur de α et le programme renverra P(α ).
(On pourra aussi renvoyer les coefficients de Q).
Se programme comme pour les entiers en utilisant la division
euclidienne des polynômes (quo, rem).
Applications :
Écrire une fonction qui détermine les racines rationnelles d’un polynôme P à coefficients entiers (elles sont de la forme p/q où q divise le coefficient dominant de P et ± p divise son coefficient de plus bas degré). Tester avec le polynôme P=12x5+10x4−6x3+11x2−x−6.
N.B.: ce n’est pas un algorithme efficace, cf. le manuel Algorithmes de calcul formel.
Interpolation de Lagrange (ne semble pas au programme, mais est
plus ou moins implicite dans les méthodes de quadrature), fonction
lagrange de Xcas. Référence : documentation de Xcas
et livre de Demailly.
Référence: cours mat249, livre de Demailly. Exposés : 201, 208, 217, 220, 250, 252, 253, 254, 255, 256 Exercices : 403, 406, 433, 443
Recherche de solutions de f(x)=0 sur [a,b] sachant que f(a)f(b)<0.
Recherche de solutions de f(x)=x par étude de la suite un+1=f(un). Exemple, résolution de l’équation du temps en mécanique céleste x − e sin(x) = y, e ∈ [0,1[.
Exercice 1
Écrire un programme iter prenant en argument
la fonction f, la valeur de u0, de N et de ε, et
qui s’arrête dès que l’une des conditions suivantes est satisfaite :
Dans le premier cas le programme renverra la valeur de un+1, dans
le second cas une séquence composée de uN et de N.
Tester votre programme avec f(x)=√2+x et f(x)=x2.
On suppose que la fonction f satisfait aux hypothèses du théorème du point fixe. On notera k<1 la constante de contractance. On peut alors trouver un encadrement de la limite l de la suite (un) en fonction de un, un−1 et k.
Exercice 2
Écrire un programme iter_k prenant en argument
la fonction f, la valeur de u0, la constante k et l’écart toléré ε,
et qui s’arrête dès que | un − l | ≤ ε.
Vérifier les hypothèses du théorème du point fixe pour f(x)=3cos(x/4)
sur [0,3] et expliciter une constante de contractance k.
Déterminer une valeur approchée de la limite de (un)
à 10−3 près en utilisant la fonction iter_k.
La convergence de ces suites est en général linéaire, le nombre de décimales exactes augmente de la même valeur à chaque itération. Par contre lorsqu’on est prêt d’une racine, la méthode de Newton permet en gros de multiplier par deux le nombre de décimales à chaque itération.
Exercice 3
| f(x)= |
|
iter, trouver un encadrement de
√5 à 10−6 près par les deux méthodes (on pourra
prendre une valeur initiale approchée puis entière exacte pour avoir
une valeur numérique approchée puis une fraction).
Combien d’itérations sont nécessaires?
Dans certains cas, la fonction f n’est pas contractante, mais on peut réécrire l’équation à résoudre sous une autre forme avec une fonction contractante, par exemple en utilisant une fonction réciproque.
Exercice 4
Donner un encadrement à 10−6 près d’une racine
de l’équation tan(x)=x sur l’intervalle ]5π/2,7π/2[ en
utilisant une méthode de point fixe.
| un+1=un− |
|
Programmation de la méthode de Newton. Utilisation de la convexité pour prouver la convergence. Exemple : racine carrée, recherche des racines d’un polynôme (illustration avec l’utilisation du calcul formel pour prouver la convexité), bassins d’attraction des racines.
Rectangles, trapèzes, Simpson. Programmation ou/et illustration. Majoration de l’erreur (illustration possible avec calcul formel pour calcul exact de l’intégrale)
Exercice 1 : Calculer une valeur approchée de
| ∫ |
|
|
par la méthode des rectangles, du point milieu et des trapèzes
en utilisant un pas de 1/10 et de 1/100 (vous pouvez la fonction
plotarea ou utiliser le
tableur ou écrire un programme effectuant ce calcul avec comme arguments
la fonction, la borne inférieure, la borne supérieure et le nombre
de subdivision). Observez numériquement
la différence entre les valeurs obtenues et la valeur exacte de
l’intégrale.
Exercice 2
Calculer le polynôme interpolateur P de Lagrange de f(x)=1/1+x2
aux points d’abscisse j/4 pour j variant de 0 à 4.
Donner un majorant de la différence entre P et f en un point
x ∈ [0,1]. Représenter graphiquement ce majorant.
Calculer une majoration de l’erreur entre l’intégrale de f
et l’intégrale de P sur [0,1].
En déduire un encadrement de π/4.
Exercice 3
On reprend le calcul de ∫01 dx/1+x mais en utilisant
un polynôme interpolateur de degré 4 sur N subdivisions de [0,1]
(de pas h=1/N). Déterminer une valeur de N telle que la valeur
approchée de l’intégrale ainsi calculée soit proche à 10−8 près
de ln(2). En déduire une valeur approchée à 10−8 de
ln(2).
Même question pour ∫01 dx/1+x2
et π/4 (pour majorer la dérivée n-ième de 1/1+x2,
on pourra utiliser une décomposition en éléments simples sur ℂ).
Méthode d’Euler (illustration avec interactive_odeplot
ou/et programmation).
Exercice 1.
Donner une majoration indépendante de x de l’erreur commise.
(A titre d’illustration, tracer la différence T8(x)−cos(x).)
Exercice 2.
On veut approcher sin(x) à 1e-6 près
en utilisant des développements en séries entières.
| T2k+1(x)= |
| (−1)j |
|
1e-6 de sin(x) sur [−100,100] en justifiant et en effectuant les
étapes suivantes:
Exercice 3
| α− |
| + |
| − |
| ≤ arctan(α) ≤ α− |
| + |
|
1e-8 près de π/4= arctan(1).
| = 4 arctan( |
| ) − arctan( |
| ) |
Exercice 4
| g(x)= |
|
| I= | ∫ |
| g(x) dx |
Méthode de Richardson lorsqu’on connait le développement asymptotique d’une suite convergeant vers la limite. Par exemple pour la constante d’Euler (cf. le manuel de programmation de Xcas)
Cas particulier : convergence de la méthode de trapèzes : méthode de Romberg pour approcher numériquement une intégrale.
Autres : méthode du Δ2 d’Aitken, séries alternées
(théorie voir le texte 585 de l’option
C de l’agrégation externe,
http://agreg.dnsalias.org/Textes/585.pdf,
cf. la session du menu Exemples, analyse de Xcas).
Référence: dans Xcas, menu Aide, puis manuel, puis Programmation (pour le pivot de Gauss) ou Algorithmes de calcul formel (Gauss-Bareiss, déterminant modulaire, Fadeev).
Exposés : 110, 151, 155, 160, 162 Exercices : 310, 313, 316, 319,
On rappelle qu’en mode Xcas, les indices commencent à 0 (en mode maple les indices commencent à 1). Etant donné une matrice M ayant L lignes et C colonnes, on demande de programmer l’algorithme du pivot de Gauss que l’on rappelle :
rowSwap(matrice,l1,l2))
mRowAdd(coeff,matrice,l1,l2))
Pour calculer l’inverse d’une matrice M carrée de taille n, on peut résoudre simultanément n systèmes linéaires du type Mxk=yk où yk représente (en colonne) les coordonnées du k-ième vecteur de la base canonique pour k variant de 1 à n. On écrit donc la matrice M puis les colonnes des coordonnées de ces n vecteurs, donc la matrice identité de taille n. On réduit complètement la matrice par l’algorithme du pivot de Gauss. Si M est inversible, les n premières colonnes après réduction doivent être la matrice identité de taille n, alors que les n colonnes qui suivent sont les coordonnées des xk donc ces n colonnes constituent M−1. Écrire un programme de calcul d’inverse de matrice par cet algorithme.
On présente ici deux méthodes, la première se généralise au cas des systèmes à coefficients entiers, la deuxième utilise un peu moins de mémoire (elle travaille sur une matrice 2 fois plus petite).
Première méthode Soir M la matrice dont on cherche le noyau. On ajoute à droite de la matrice transposée de M une matrice identité ayant le même nombre de lignes que Mt. On effectue une réduction sous-diagonale qui nous amène à une matrice composée de deux blocs
| ( Mt In ) → ( U L ) |
Attention, L n’est pas la matrice L de la décomposition LU de Mt, on a en fait
| L Mt = U |
donc
| M Lt = Ut |
Les colonnes de Lt correspondant aux colonnes nulles de Ut (ou si on préfère les lignes de L correspondant aux lignes nulles de U) sont donc dans le noyau de M et réciproquement si Mv=0 alors
| Ut (Lt)−1 v =0 |
donc, comme U est réduite, (Lt)−1 v est une combinaison linéaire des vecteurs de base d’indice les lignes nulles de U. Finalement, les lignes de L correspondant aux lignes nulles de U forment une base du noyau de M.
On peut faire le raisonnement ci-dessus à l’identique si M est une matrice à coefficients entiers, en effectuant des manipulations élémentaires réversibles dans ℤ, grâce à l’idendité de Bézout. Si a est le pivot en ligne i, b le coefficient en ligne j à annuler, et u, v, d les coefficients de l’identité de Bézout a u + b v =d on fait les changements :
| Li ← uLi +v Lj, Lj ← − |
| Li + |
| Lj |
qui est réversible dans ℤ car le déterminant de la sous-matrice élémentaire correspondante est
| ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ |
| ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ | = 1 |
Cette réduction (dite de Hermite) permet de trouver une base du noyau à coefficients entiers et telle que tout élément du noyau à coefficient entier s’écrit comme combinaison linéaire à coefficients entiers des éléments de la base.
Deuxième méthode On commence bien sûr par réduire la matrice (réduction complète en-dehors de la diagonale), et on divise chaque ligne par son premier coefficient non nul (appelé pivot). On insère alors des lignes de 0 pour que les pivots (non nuls) se trouvent sur la diagonale. Puis en fin de matrice, on ajoute ou on supprime des lignes de 0 pour avoir une matrice carrée de dimension le nombre de colonnes de la matrice de départ. On parcourt alors la matrice en diagonale. Si le i-ième coefficient est non nul, on passe au suivant. S’il est nul, alors tous les coefficients d’indice supérieur ou égal à i du i-ième vecteur colonne vi sont nuls (mais pas forcément pour les indices inférieurs à i). Si on remplace le i-ième coefficient de vi par -1, il est facile de se convaincre que c’est un vecteur du noyau, on le rajoute donc à la base du noyau. On voit facilement que tous les vecteurs de ce type forment une famille libre de la bonne taille, c’est donc bien une base du noyau.
Lorsque les coefficients de la matrice M=(mj,k)0 ≤ j <L, 0 ≤ k < C sont entiers, on peut souhaiter éviter de faire des calculs dans les rationnels, et préférer utiliser une combinaison linéaire de lignes ne faisant intervenir que des coefficients entiers. On peut par exemple effectuer l’opération (où l,c désignent la ligne et colonne du pivot)
| Lj ← ml,c Lj − mj,c Ll |
Tester la taille des coefficients obtenus pour une matrice carrée aléatoire de taille 5 puis 10. L’idée est-elle bonne ?
On peut montrer qu’on peut toujours diviser par le pivot utilisé pour réduire la colonne précédente (initialisé à 1 lors de la réduction de la première colonne)
| Lj ← |
| (ml,c Lj − mj,c Ll) |
Tester à nouveau sur des matrices carrées de taille 5, 10, vérifier que les calculs sont bien effectués dans ℤ. Comparer le dernier coefficient en bas à droite avec le déterminant de la matrice (si vous avez vu les propriétés des déterminants, démontrez ce résultat. Plus difficile: en déduire qu’on peut bien diviser par le pivot de la colonne précédente en considérant des déterminants de matrices extraites)
Exercice 1
Calcul de la somme de deux sous-espaces vectoriels.
On donne deux sous-espaces vectoriels E1 et E2
de ℝn par deux familles géneratrices (c’est-à-dire
une liste de vecteurs), il s’agit d’écrire un algorithme permettant
On pourra utiliser les fonctions rref ou/et ker de Xcas. Tester avec deux sous-espaces de dimension 2 de ℝ5. N’oubliez pas de rédiger une justification mathématique de la méthode mise en oeuvre.
Exercice 2
Calcul de l’intersection de deux sous-espaces vectoriels.
On donne deux sous-espaces vectoriels E1 et E2
de ℝn par deux familles géneratrices (c’est-à-dire
une liste de vecteurs), il s’agit d’écrire un algorithme permettant
de trouver une base de E1 ∩ E2.
On pourra utiliser les fonctions rref ou/et ker de Xcas.
Tester avec 2 sous-espace de dimension 3 de ℝ4.
N’oubliez pas de rédiger une justification mathématique de la
méthode mise en oeuvre.
Exercice 3
Mise en oeuvre du théorème de la base incomplète. On donne
une famille de vecteurs de ℝn (liste de vecteurs),
il s’agit d’écrire un algorithme permettant d’extraire
une base du sous-espace engendré,
puis de compléter par des vecteurs afin de former une base
de ℝn tout entier, et enfin d’écrire un vecteur de ℝn
comme combinaison linéaire des vecteurs de la base complétée.
On pourra utiliser les fonctions rref ou/et ker de Xcas.
Tester avec une famille de 3 vecteurs de ℝ5.
N’oubliez pas de rédiger une justification mathématique de la
méthode mise en oeuvre.
Calcul par réduction. Autres algorithmes (à la limite du programme) : développement (2n opérations), modulaire (sur ℤ).
On propose ici quelques algorithmes de calcul du polynôme minimal M et/ou caractéristique P d’une matrice carrée A de taille n. On peut bien sur calculer le polynôme caractéristique en calculant directement le déterminant (par exemple avec l’algorithme de Gauss-Bareiss pour éviter les fractions de polynômes), mais c’est souvent plus couteux que d’utiliser un des deux algorithmes ci-dessous.
On sait que le degré de P est n, il suffit donc de connaitre
la valeur de P en n+1 points distincts pour connaitre P.
Par exemple, on calcule P(k) pour k=0,...,n, et en utilisant
l’instruction lagrange de Xcas, on en déduit le
polynôme caractéristique. Programmer cet algorithme et testez
votre programme en comparant le résultat de votre programme
et de l’instruction charpoly de Xcas sur quelques
matrices.
Cet algorithme permet de calculer le polynôme caractéristique C dans presque tous les cas, en recherchant le polynôme minimal M de A (on note m le degré de M).
On sait que M(A)=0, donc pour tout vecteur v, M(A)v=0. Si M(x)=∑k=0m Mk xk, alors
| 0 = M(A)v= |
| Mk Ak v |
On va donc rechercher les relations linéaires qui existent entre les n+1 vecteurs v, ..., Anv. Cela revient à déterminer le noyau K de l’application linéaire dont la matrice a pour colonnes v,...,Anv. On sait que le vecteur (M0,...,Mm,0,...,0) ∈ ℝn+1 des coefficients de M (complété par des 0 si m<n) appartient à ce noyau K. Si le noyau K est de dimension 1, alors m=n, et les coefficients de M sont proportionnels aux coefficients du vecteur de la base du noyau calculé. On en déduit alors le polynôme minimal M et comme m=n, et aussi le polynôme caractéristique C=M.
Programmer cet algorithme en prenant un vecteur v aléatoire. Attention, pour calculer Akv on formera la suite récurrente vk=Avk−1, pourquoi?
Si on n’a pas de chances dans le choix de v, on trouvera un noyau de dimension plus grande que 1 bien que le polynome minimal soit de degré n. On peut alors recommencer avec un autre vecteur. On peut aussi chercher la relation de degré minimal pour un v donné (elle apparait automatiquement comme premier vecteur de la base dans l’algorithme de calcul du noyau donné au TP2), prendre le PPCM P des polynomes pour 2 ou 3 vecteurs aléatoires et conclure s’il est de degré n. Programmer cette variante.
Si le polynome minimal est de degré m<n, on peut tester si le PPCM P est le polynome minimal en calculant P(A) mais ce calcul est couteux. On peut aussi faire confiance au hasard, supposer que le polynome minimal est bien M=P et essayer de déterminer C/M par d’autres moyens, par exemple la trace si n=m+1. On dispose alors d’un moyen de vérification en calculant l’espace caractéristique correspondant à la valeur propre double. Programmer cette variante.
Algorithme de Fadeev (hors programme, référence, manuel algorithmes de calcul formel).
Si la matrice M possède une seule valeur propre de module maximal, alors la suite vn+1=Mvn se rapproche de l’espace propre correspondant. Écrire un programme mettant en oeuvre la méthode de la puissance dans le cas réel (penser à normaliser vn pour éviter les overflow). Utilisez ce programme pour trouver une valeur approchée de la valeur propre de norme maximale par la méthode de la puissance de B t B où B est une matrice aléatoire.
Représentation par la liste des images ou par une liste de cycles, chaque cycle étant la liste des éléments du cycle. Composition, inverse, décomposition en cycles.
Encadrement de π par un polygone régulier inscrit et par un polygone extérieur.
En prenant un couple de nombres aléatoires entre -1 et 1, on regarde si le couple est dans le disque de rayon 1, et on calcule la proportion.
Expérience : on tire n fois au hasard une boule dans une urne avec remise, la proportion de boules noires dans l’urne étant p (proche de 0.5), on compte le nombre de boules noires. On effectue N fois l’expérience. On compte le nombre de fois où la proportion de boules noires observée est dans un intervalle centré en p de taille proportionnelle à 1/√n (par exemple [p−1/√n,p+1/√n]). On peut aussi représenter graphiquement par des segments horizontaux les intervalles proportion de boules observées ± 1/√n et tracer verticalement la droite d’abscisse p.
Expérience : on ajoute n réels tirés au hasard pris entre 0 et 1 (loi uniforme, ou autre loi), par exemple pour n=30. On effectue N fois l’expérience (par exemple N=1000), et on trace un histogramme des résultats, et sur le même graphique on représente la loi normale de moyenne l’espérance de la loi et d’écart-type l’écart-type de la loi divisé par √n.
Par discrétisation.
Remarques :
www-fourier.ujf-grenoble.fr/~parisse/giac_fr.html
Numéros de leçons correspondants : 302, éventuellement 303. Exposés : 103.
Énoncé : : le but de l’exercice est de calculer rapidement an modulo p où a,n,p sont 3 entiers avec 0 ≤ a < p.
| n= |
| nj 2j, nj ∈ { 0,1 } |
| aj=a2j (mod p ) = (a2j−1 (mod p ))2 (mod p ) |
Solution :
base2(n):={
local l;
l:=NULL; // sequence vide
tantque n>0 faire
// ajoute n_0 puis n_1 puis ... en fin de l
l:=l,irem(n,2); // irem: reste euclidien
n:=iquo(n,2);
ftantque;
// renverse l'ordre des elements de l
return revlist(l);
}
N.B. :
l’instruction intégrée Xcas ou Maple permettant de faire la conversion
est convert(n,base,2), elle renvoie la liste des bits à l’envers.En syntaxe compatible Maple V, sans renverser l’ordre des éléments de
l :
base2:=proc(n0)
local l,n;
n:=n0;
l:=NULL;
while n>0 do
l:=l,irem(n,2);
n:=iquo(n,2);
od;
RETURN l;
end;
Attention, avec Maple il n’est pas possible de modifier les arguments d’une procédure et il n’existe pas de débugger interactif.
|
base2 précédente, on peut
aussi calculer simultanément les nj et an (mod p ). En syntaxe Xcas en
français:
puimod(a,n,p):={
local aj,an_modp;
aj:=irem(a,p); // precaution au cas ou a>=p
an_modp:=1;
tantque n>0 faire
si irem(n,2)==1 alors
an_modp:=irem(an_modp*aj,p);
fsi;
n:=iquo(n,2);
aj:=irem(aj*aj,p);
ftantque;
return an_modp;
}
En syntaxe compatible Maple V :
puimod:=proc(a,n0,p)
local aj,an_modp,n;
n:=n0;
aj:=irem(a,p);
an_modp:=1;
while n>0 do
if irem(n,2)=1 then
an_modp:=irem(an_modp*aj,p);
fi;
n:=iquo(n,2);
aj:=irem(aj*aj,p);
od;
RETURN an_modp;
end;
Commentaires
Numéros de leçons correspondants : 305, 302. Exposés : 104, 103.
Énoncé :
isprime pour tester si un entier
est premier.
Solution :
| (a+1)P = |
| ⎛ ⎝ |
| ⎞ ⎠ | aj = 1+ap + |
| ⎛ ⎝ |
| ⎞ ⎠ |
|
test(p,nmax):={
local a;
nmax := min(nmax,p-1);
pour a de 2 jusque nmax faire
si powmod(a,p,p)!=a alors retourne a; fsi;
fpour
retourne 0;
}:;
On utilise ici la fonction powmod (puissance modulaire rapide).
Le paramètre nmax sert à limiter le nombre d’essais, ici on essaie les
entiers de 2 à nmax, mais on pourrait aussi tirer nmax entiers
au hasard. La boucle sur a se poursuit jusqu’à ce que le test de Fermat
échoue, return provoque en effet l’arrêt prématuré
de la fonction (donc de la boucle) ou va à son terme (test réussi, valeur
de retour 0).
L’équivalent en syntaxe maple V en utilisant la fonction
puimod ci-dessus (aussi utilisable dans Xcas) :
test:=proc(p,nmax0)
local a,nmax;
nmax := min(nmax0,p-1);
for a from 2 to nmax do
if puimod(a,p,p)<>a then RETURN a; fi;
od;
RETURN 0;
end
Lorsque la réponse du programme est non nulle, on est assuré que le nombre n’est pas
premier, la valeur renvoyée est un certificat de non-primalité (un entier
tel que ap ≠ a (mod p )). Par contre, si la valeur renvoyée est
0, on ne peut pas en déduire que p est premier (il y a seulement
des présomptions).
carmi(nmax):={
local a,p;
pour p de 3 jusque nmax faire
pour a de 2 jusque p-1 faire
si powmod(a,p,p)!=a alors break; fsi;
fpour
// a la sortie de la boucle,
// a vaut p si le test a reussi
// ou est inferieur a p (arret par break)
si a==p et !isprime(p) alors return a; fsi;
fpour;
return "non trouve";
}
En syntaxe Maple V :
carmi:=proc(nmax)
local a,p;
for p from 3 to nmax do
for a from 2 to p-1 do
if puimod(a,p,p)<>a then break; fi;
od;
if a=p and not isprime(p) then RETURN a; fi;
od;
RETURN "non trouve";
end
On exécute par exemple carmi(100) qui renvoie “non trouvé”
(il n’y a pas de nombre de Carmichael inférieur à 100) alors
que carmi(1000) renvoie après quelques secondes 561, le premier
nombre de Carmichael.
Commentaires :
pari("isprime",9856989898997789789,1).
Numéros de leçons correspondants : 306, Exposés: 159, 157
Énoncé : Soient a et b 2 entiers. On pose r0=a,r1=b, et pour n≥ 0 tel que rn+1 ≠ 0, on note qn+2 et rn+2 le quotient et le reste de la division euclidienne de rn par rn+1. On rappelle qu’il existe un (premier) entier N tel que rN= 0, et que pgcd(a,b)=rN−1 (dernier reste non nul).
| u0=1, u1=0, un+2=un−qn+2 un+1, v0=1, v1=0, vn+2=vn−qn+2 vn+1 |
| aun+bvn=rn |
| au+bv=auN−1+bvN−1=pgcd(a,b) |
| un rn+1− un+1 rn |
Solution
|
l0 contiendra les valeurs de rn,un,vn, l1 les
valeurs de rn+1,un+1,vn+1 et l2 les valeurs
de rn+2,un+2,vn+2, en fin de boucle l’incrémentation de n
se traduit par la recopie de l1 dans l0 et de
l2 dans l1.
bez(a,b):={
local l0,l1,l2,q;
l0:=[a,1,0];
l1:=[b,0,1];
tantque l1[0]!=0 faire
q:=iquo(l0[0],l1[0]);
l2:=l0-q*l1;
l0:=l1;
l1:=l2;
ftantque;
// l0 contient les coefficients de Bezout
// l1 contient les cofacteur du PPCM
return l0;
}
On peut observer le déroulement de l’algorithme pas à pas en tapant par exemple
bez:=proc(a,b)
local u0,u1,u2,v0,v1,v2,r0,r1,r2,q;
r0:=a; u0:=1; v0:=0;
r1:=b; u1:=0; v1:=1;
while r1<>0 do
q:=iquo(r0,r1);
r2:=r0-q*r1;
u2:=u0-q*u1;
v2:=v0-q*v1;
r0:=r1; u0:=u1; v0:=v1;
r1:=r2; u1:=u2; v1:=v2;
od;
RETURN [r0,u0,v0];
end
| au+bv=1 |
|
| unrn+1 −un+1 rn = (−1)n (u0r1 −u1 r0 ) =(−1)n b |
| uN a =(−1)N |
| = (−1)N ppcm(a,b) |
| auN+bvN=0 |
Commentaires
Numéros de leçons correspondants : 310, 319. Exposés: 110.
Énoncé : Soit A une matrice carrée d’ordre n à coefficients complexes. On se propose de déterminer le polynôme caractéristique de A en recherchant des relations linéaires entre un vecteur v0, et ses images successives par A: vk+1=Avk
Solution
| pj vj=0 |
ker
de calcul du noyau.
polmin(A):={
local k,n,vk,B,noyau;
n:=size(A); // nombre de lignes
vk:=ranm(n); // initialise v0 aleatoire (reel)
B:=[vk]; // initialise B (tranposee)
pour k de 1 jusque n faire
vk:=A*vk;
B[k]:=vk;
fpour;
B:=tran(revlist(B)); // vn ... v0
noyau:=ker(B);
si size(noyau)>1 alors
return "echec de l'algorithme";
fsi;
return noyau[0];
}
On peut le tester avec une matrice aléatoire, par exemple
A:=ranm(3,3). Le même algorithme est plus
délicat à écrire en syntaxe Maple selon qu’on utilise
with(linalg); ou les nouvelles fonctions d’algébre
linéaire sur les versions plus récentes, une matrice n’est pas
une liste de listes avec linalg, le produit matrice-vecteur ou
matrice-matrice n’est pas simplement *, etc.
Commentaires
Numéros de leçons correspondants : 443. Exposés: 250, 201, 208
Énoncé :
Soit e∈[0,1[ et t∈[−π,π]
On se propose de résoudre l’équation
| x=t + e sin(x) |
par dichotomie et par la méthode du point fixe et comparer ces deux méthodes. (N.B. e ne désigne pas la base des logarithmes ici, mais l’excentricité d’une ellipse dont on cherche à résoudre l’équation du temps).
1e-8 près.
Combien d’étapes sont-elles nécessaires ?
1e-8 près. Combien d’étapes sont-elles nécessaires ?
Solution
dicho(f,a,b,eps):={
local m,fa,fb,fm;
// assure une resolution numerique
a:=evalf(a); b:=evalf(b);
fa:=f(a); fb:=f(b);
// teste si a ou b est solution
si fa=0 alors return a; fsi;
si fb=0 alors return b; fsi;
// si fa*fb>=0 pas forcement de solution
si fa*fb>=0 alors return "echec"; fsi;
// debut de la dichotomie
tantque b-a>eps faire
m:=(a+b)/2;
fm:=f(m);
si fm==0 alors return m fsi; // solution trouvee
si fa*fm<0 alors b:=m; sinon a:=m; fsi;
ftantque;
return [a,b]; // intervalle contenant la solution
}
Puis on essaie par exemple pour t=1 et e=0.1.
PI
au lieu de pi avec Maple) :
dicho:=proc(f,a0,b0,eps)
local m,fa,fb,fm,a,b;
a:=evalf(a0); b:=evalf(b0);
fa:=f(a); fb:=f(b);
if fa=0 then RETURN a; fi;
if fb=0 then RETURN b; fi;
if fa*fb>=0 then RETURN "echec"; fi;
while b-a>eps do
m:=(a+b)/2;
fm:=f(m);
if fm=0 then RETURN m fi;
if fa*fm<0 then b:=m; else a:=m; fi;
od;
RETURN [a,b];
end;
Le nombre d’étapes est le plus petit entier n tel que
| ≤ ε |
| n ≥ |
|
| |un−l| ≤ 2(π+e) en |
| |un−l| ≤ |
|
fixe(f,u0,k,eps):={
local u1;
u0:=evalf(u0);
u1:=f(u0);
eps:=eps*(1-k);
tantque abs(u1-u0)>eps faire
u0:=u1;
u1:=f(u0);
ftantque;
return u0;
}
On teste ensuite avec par exemple u0=0.0 :
fixe:=proc(f,u,k,eps0)
local u0,u1,eps;
u0:=evalf(u);
u1:=f(u0);
eps:=eps0*(1-k);
while abs(u1-u0)>eps do
u0:=u1;
u1:=f(u0);
od;
RETURN u0;
end;
g:=x->1+0.1*sin(x); fixe(g,0.0,0.1,1e-8);
Le nombre d’étapes se majore avec l’estimation à priori
| n ≥ |
|
Commentaires
Dans la documentation de Xcas (donc disponible le jour de l’oral 2), aller dans le menu Aide, sous-menu Manuels, puis Programmation. Vous y trouverez une grande partie des algorithmes ci-dessus. Le manuel Exercices et le manuel Amusement peuvent aussi contenir des exercices à vocation algorithmique ou illustrable par Xcas. Enfin, le manuel tableur, statistiques peut servir pour une illustration d’exercices de proba-stats. Vous pouvez aussi lancer une recherche avec un mot clef (menu Aide, Rechercher un mot, par exemple Bezout).
En livres :
Sur le web :
http://www-fourier.ujf-grenoble.fr/~parisse/giac_fr.html,
la page de Xcas, avec en particulier la page algorithmiquehttp://www-fourier.ujf-grenoble.fr/~parisse/algoxcas.html
http://tehessin.tuxfamily.org/, site de
Guillaume Connan, contient un livre téléchargeable librement
couvrant pratiquement tout le programme
http://www-fourier.ujf-grenoble.fr/~parisse/#mat249,
B. Parisse, mon cours Math Assistés par ordinateur
de licence 2ìeme année
http://www.inria.fr/rrrt/rr-5105.html
http://www.loria.fr/~zimmerma/mca/pub226.html
| Instructions Xcas en français | |
| affectation | a:=2; |
| entrée expression | saisir("a=",a); |
| entrée chaine | saisir_chaine("a=",a); |
| sortie | afficher("a=",a); |
| valeur retournée | retourne(a); |
| arrêt dans boucle | break; |
| déclaration fonction | f(parametres):={ ... } |
| alternative | si <condition> alors <inst> fsi; |
si <condition> alors <inst1> sinon <inst2> fsi; | |
| boucle pour | pour j de a jusque b faire <inst> fpour; |
pour j de a jusque b pas p faire <inst> fpour; | |
| boucle répéter | repeter <inst> jusqua <condition>; |
| boucle tantque | tantque <condition> faire <inst> ftantque; |
| Instructions Xcas "C-like" | |
| affectation | a:=2; |
| entrée expression | input("a=",a); |
| entrée chaine | textinput("a=",a); |
| sortie | print("a=",a); |
| valeur retournée | return(a); |
| arrêt dans boucle | break; |
| déclaration fonction | f(parametres):={ ... } |
| alternative | if (<condition>) {<inst>}; |
if (<condition>) {<inst1>} else {<inst2>}; | |
| boucle pour | for (j:= a;j<=b;j++) {<inst>} |
for (j:= a;j<=b;j:=j+p) {<inst>} | |
| boucle répéter | repeat <inst> until <condition>; |
| boucle tantque | while (<condition>) {<inst>}; |
| Instructions Maxima | |
| affectation | a:2; |
| sortie | print("a=",a) |
| valeur retournée | return(a) |
| déclaration fonction | f(parametres):=( ... ) |
| alternative | if (<condition>) then <inst> |
if (<condition>) then <inst1> else <inst2> | |
| boucle pour | for j:1 thru 26 do (inst1,inst2) |
for j:2 thru 20 step 2 do (inst,...,inst) | |
| boucle tantque | for i:1 while <condition> do (inst,...,inst) |
| Instructions Maple V ou Xcas en mode Maple | |
| affectation | a:=2; |
| sortie | print("a=",a); |
| valeur retournée | RETURN(a); |
| arrêt dans boucle | break; |
| déclaration fonction | f:=proc(parametre) ... end; |
| alternative | if <condition> then <inst> fi; |
if <condition> then <inst1> else <inst2> fi; | |
| boucle pour | for j from a to b do <inst> od; |
for j from a to b step p do <inst> od; | |
| boucle tantque | while <condition> do <inst> od; |
Avec des versions plus récentes de Maple, on utilisera return
et non RETURN.
| Instructions Python | |
| affectation | a=2 |
| entrée expression | a=input("a=") |
| sortie | print "a=",a |
| déclaration fonction | def f(parametres): |
instructions | |
| valeur retournée | return a; |
| arrêt dans boucle | break; |
| alternative | if <condition>: |
instructions1 | |
else: | |
instructions2 | |
| boucle pour | for j in range(debut,fin,pas): |
instructions | |
| boucle tant que | while condition: |
instructions | |
| Instructions Scilab | |
| affectation | a=2 |
| entrée expression | a=input("a=") |
| sortie | disp("a=",a) |
| déclaration fonction | function nom_var=f(parametres) ... endfunction |
| valeur retournée | nom_var=...; |
| arrêt dans boucle | break; |
| alternative | if <condition> then <inst> end; |
if <condition> then <inst1> else <inst2> end; | |
| boucle pour | for j=debut:fin do <inst> end; |
| boucle tantque | while <condition> do <inst> end; |
Ce document a été traduit de LATEX par HEVEA