Shared2.sagewsOpen in CoCalc

Thème 2 : Algèbre matricielle avec Sagemath, systèmes d'équations linéaires

Sujet de devoir : Application à la théorie des jeux : Ecrire un algorithme prenant en entrée une matrice M (la matrice de paiement du premier joueur J1 d'un jeu fini à somme nulle), rendant les sommets (au sens de la convexité) de l'ensemble des stratégies mixtes prudentes de J1. Dessiner le découpage des stratégies mixtes prudentes de J1, le paiement minimal sur chacun des morceaux, les stratégies mixtes prudentes. Commencer par le cas où J1 a deux stratégies pures (la matrice M a deux lignes).

Document de cours sur la recherche des équilibre de l'extension mixte d'un jeu (fini à deux joueur à somme nulle), exemple de recherche algorithmique des solutions ; voir la page web du cours.

Autre sujet de devoir (difficile) : Application aux Chaînes de Markov : Ecrire un algorithme prenant en entrée la matrice M des probabilités de transition d'une chaîne de Markov (homogène, à temps discret, ayant un ensemble fini d'états), rendant la ou les formes limites de Mn pour n grand. Simuler la chaîne de Markov de matrice M pour un état initial donné ; estimer les probabilités de présence en un état après un temps long, les temps d'atteinte d'un état donné.

Documents de cours sur cette page.


TP2 : Expérimenter l'algèbre matricielle avec Sagemath

Tutoriel sur l'algèbre matricielle avec sagemath

D'autres ressources sur internet (google "sagemath matrices")

Soit M la matrice

\left(\begin{array}{rr} \frac{1}{3} & \frac{1}{2} \\ \frac{2}{3} & \frac{1}{2} \end{array}\right)
et I la matrice unité. Exhiber avec Sagemath une base du noyau de M-I (l'espace propre de M pour la valeur propre 1 ). Quelles sont les méthodes disponibles ? Comment utiliser la méthode (M_I).kernel() ?

Quel est le système d'équations linéaires correspondant ? (Appliquer la méthode solve() ). Illustration graphique.

Calculer M^k pour k=2,16,128. Calculer une valeur numérique approchée de M^k avec la méthode (M^k).n()

Mêmes questions avec la matrice

N=\left(\begin{array}{rrr} \frac{1}{2} & 0 & 0 \\ \frac{1}{4} & 0 & 1 \\ \frac{1}{4} & 1 & 0 \end{array}\right)


Au passage listes, vecteurs, matrices dans Python, listes imbriquées, génération de listes

Corrigé du TP

M=matrix([[1/3,1/2],[2/3,1/2]]);I=matrix([[1,0],[0,1]])
show("M=",M," , I=",I) #latex(M)
M= \displaystyle \left(\begin{array}{rr} \frac{1}{3} & \frac{1}{2} \\ \frac{2}{3} & \frac{1}{2} \end{array}\right) , I= \displaystyle \left(\begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array}\right)
M-I
[-2/3 1/2] [ 2/3 -1/2]
K=(M-I).right_kernel();K
X=K.basis();X;
M*X[0]
Vector space of degree 2 and dimension 1 over Rational Field Basis matrix: [ 1 4/3] [ (1, 4/3) ] (1, 4/3)
M*K.basis_matrix().transpose()
[ 1] [4/3]
M.eigenvalues()
[1, -1/6]
M.eigenvectors_right()
[(1, [ (1, 4/3) ], 1), (-1/6, [ (1, -1) ], 1)]
M.eigenvectors_right? #mieux que M.eigenvectors_right()?
File: /projects/sage/sage-7.5/src/sage/matrix/matrix2.pyx
Signature : M.eigenvectors_right(self, extend=True)
Docstring :
Compute the right eigenvectors of a matrix.

For each distinct eigenvalue, returns a list of the form (e,V,n)
where e is the eigenvalue, V is a list of eigenvectors forming a
basis for the corresponding right eigenspace, and n is the
algebraic multiplicity of the eigenvalue. If "extend = True" (the
default), this will return eigenspaces over the algebraic closure
of the base field where this is implemented; otherwise it will
restrict to eigenvalues in the base field.

EXAMPLES: We compute the right eigenvectors of a 3 x 3 rational
matrix.

   sage: A = matrix(QQ,3,3,range(9)); A
   [0 1 2]
   [3 4 5]
   [6 7 8]
   sage: es = A.eigenvectors_right(); es
   [(0, [
   (1, -2, 1)
   ], 1),
   (-1.348469228349535?, [(1, 0.1303061543300932?, -0.7393876913398137?)], 1),
   (13.34846922834954?, [(1, 3.069693845669907?, 5.139387691339814?)], 1)]
   sage: A.eigenvectors_right(extend=False)
   [(0, [
   (1, -2, 1)
   ], 1)]
   sage: eval, [evec], mult = es[0]
   sage: delta = eval*evec - A*evec
   sage: abs(abs(delta)) < 1e-10
   True

Méthode solve, exploitation de la réponse

var('x,y')
show(1/3*x+1/2*y==x);show(2/3*x+1/2*y==y)
S=solve([1/3*x+1/2*y==x,2/3*x+1/2*y==y],[x,y]);S
(x, y)
\displaystyle \frac{1}{3} \, x + \frac{1}{2} \, y = x
\displaystyle \frac{2}{3} \, x + \frac{1}{2} \, y = y
[[x == 3/4*r1, y == r1]]
#le paramètre r* est créé par Sage et change à chaque exécution. Peut on lui assigner une valeur ?
#solution après visualisation du résultat
S[0][1].rhs() #rhs pour "right hand side"
vector([x,y]).subs(S).subs(S[0][1].rhs()==1) #[x,y].subs(S).subs(S[0][1].rhs()==1) donne une erreur
r1 (3/4, 1)
#alternative avec solution_dict
var('x y');T=solve([1/3*x+1/2*y==x,2/3*x+1/2*y==y],[x,y],solution_dict=True);T
[T[0][x],T[0][y]]
vector([T[0][x],T[0][y]]).subs(T[0][y]==1)
(x, y) [{y: r2, x: 3/4*r2}] [3/4*r2, r2] (3/4, 1)
#Une alternative sans substitution explicite du paramètre mais avec un choix d'une des inconnues
var('y t');s1=solve([1/3*x+1/2*y==x,2/3*x+1/2*y==y,y==t],[x,y]);s1
(y, t) [[x == 3/4*t, y == t]]
f(t)=s1[0];f
t |--> (x == 3/4*t, y == t)
f(1)
(x == (3/4), y == 1)
f(t)=[s1[0][0].rhs(),s1[0][1].rhs()];f
t |--> (3/4*t, t)
#Liste des paramètres dans la réponse de solve :
#maxima_calculus("%rnum_list") cf https://ask.sagemath.org/question/26205/questions-about-the-parameters-in-the-output-of-solve/
var('y');solve([1/3*x+1/2*y==x,2/3*x+1/2*y==y],[x,y])
a=maxima_calculus("%rnum_list")
print a, a[0].str()[1:]
y [[x == 3/4*r12, y == r12]] [%r12] r12
#application :
var('x y');S=solve([1/3*x+1/2*y==x,2/3*x+1/2*y==y],[x,y]);S
param=[i.str()[1:] for i in maxima_calculus("%rnum_list")]
t=[var('t'+str(i)) for i in range(len(param))];t
r=[var(i) for i in param];r
v=vector([x,y])
for i in range(len(param)):v=v.subs(S).subs(r[i]==t[i])
print "v=", v
(x, y) [[x == 3/4*r21, y == r21]] [t0] [r21] v= (3/4*t0, t0)
g(x)=[i.subs(t0==x) for i in v];g #mais g(x)=v.subs(t0==x) donne une erreur
x |--> (3/4*x, x)
#Autre méthode
g(x)=v(t0=x).list();g
x |--> (3/4*x, x)
g(1);g(1).parent()
(3/4, 1) Vector space of dimension 2 over Symbolic Ring
diff(g,x)
x |--> (3/4, 1)

Dessin

show(implicit_plot(1/3*x+1/2*y==x,(x,-2,2),(y,-2,2))+implicit_plot(2/3*x+1/2*y==y,(x,-2,2),(y,-2,2))+point(v,size=20))
%md
Puissances de $M$

Puissances de M

for n in [2,16,128]:show(M^n)
\displaystyle \left(\begin{array}{rr} \frac{4}{9} & \frac{5}{12} \\ \frac{5}{9} & \frac{7}{12} \end{array}\right)
\displaystyle \left(\begin{array}{rr} \frac{302261775799}{705277476864} & \frac{403015701065}{940369969152} \\ \frac{403015701065}{705277476864} & \frac{537354268087}{940369969152} \end{array}\right)
\displaystyle \left(\begin{array}{rr} \frac{429856276558674694292684291539547551359685980829479881352399816713084541777734705200945845241605559}{1002997978636907620016263346925610953172600621935453056488932905663863930814714312135540305563746304} & \frac{573141702078232925723579055386063401812914641105973175136533088950779389036979606934594460322140745}{1337330638182543493355017795900814604230134162580604075318577207551818574419619082847387074084995072} \\ \frac{573141702078232925723579055386063401812914641105973175136533088950779389036979606934594460322140745}{1002997978636907620016263346925610953172600621935453056488932905663863930814714312135540305563746304} & \frac{764188936104310567631438740514751202417219521474630900182044118601039185382639475912792613762854327}{1337330638182543493355017795900814604230134162580604075318577207551818574419619082847387074084995072} \end{array}\right)
for n in [2,16,128]:
    print
    show("M^"+str(n),"=",(M^n).n(digits=3))
M^2 = \displaystyle \left(\begin{array}{rr} 0.444 & 0.417 \\ 0.556 & 0.583 \end{array}\right)
M^16 = \displaystyle \left(\begin{array}{rr} 0.429 & 0.429 \\ 0.571 & 0.571 \end{array}\right)
M^128 = \displaystyle \left(\begin{array}{rr} 0.429 & 0.429 \\ 0.571 & 0.571 \end{array}\right)
show((v/sum(v)).column().n()) # la limite de chacune des colonnes de M^n qd n tend vers l'infini (telle que prédite par la théorie des chaînes de Markov)
\displaystyle \left(\begin{array}{r} 0.428571428571429 \\ 0.571428571428571 \end{array}\right)

Avec la matrice N :

N=matrix([[1/2,0,0],[1/4,0,1],[1/4,1,0]]);N;latex(N)
[1/2 0 0] [1/4 0 1] [1/4 1 0] \left(\begin{array}{rrr} \frac{1}{2} & 0 & 0 \\ \frac{1}{4} & 0 & 1 \\ \frac{1}{4} & 1 & 0 \end{array}\right)
#points fixes de la matrice (ou espace propre pour la valeur propre 1)
I=matrix.identity(3)
(N-I).right_kernel().basis()
[ (0, 1, 1) ]
var('x y z')
c=['gray','darkkhaki','aquamarine']
show(sum(implicit_plot3d((N-I)[k]*vector([x,y,z])==0,(x,-1,1),(y,-1,1),(z,-1,1),color=c[k]) for k in range(3)))
(x, y, z)
3D rendering not yet implemented
#Les limites de N^(2n) et N^(2n+1) quand n tend vers l'infini existent et diffèrent. Approximation :
(N^16).n(digits=3)
print
(N^17).n(digits=3)
[0.0000153 0.000 0.000] [ 0.500 1.00 0.000] [ 0.500 0.000 1.00] [7.63e-6 0.000 0.000] [ 0.500 0.000 1.00] [ 0.500 1.00 0.000]

Complément : Système d'équations avec un nombre n d'inconnues, où n est un paramètre

#Un exemple
reset()
n=3
x=var(['x'+str(i) for i in range(n)]);x # fait une erreur si n n'a pas de valeur ou si n=1 !
x=[var('x'+str(i)) for i in range(n)] # ok si n=1
A=matrix(3,range(9));A
for k in range(n):show(A[k]*vector(x)==0)
c=['gray','darkkhaki','aquamarine']
sum(implicit_plot3d(A[k]*vector(x)==0,*([(x[i],-1,1) for i in range(n)]),color=c[k]) for k in range(n))
#cf http://stackoverflow.com/questions/3394835/args-and-kwargs pour l'expansion des arguments *arg
#for k in range(n-1):
#    for l in [k+1..n-1]:
#        eqs=eqs+[A[k]*vector(x)==A[l]*vector(x)]
#eqs
#show(sum(implicit_plot3d(*([eq]+[(x[i],-1,1) for i in range(n)])) for eq in eqs))
(x0, x1, x2) [0 1 2] [3 4 5] [6 7 8]
\displaystyle x_{1} + 2 \, x_{2} = 0
\displaystyle 3 \, x_{0} + 4 \, x_{1} + 5 \, x_{2} = 0
\displaystyle 6 \, x_{0} + 7 \, x_{1} + 8 \, x_{2} = 0
3D rendering not yet implemented
sol=solve(*([[A[k]*vector(x)==0 for k in range(2)]]+[inc for inc in x]),solution_dict=True);sol
sol[0][x0]
param=[i.str()[1:] for i in maxima_calculus("%rnum_list")]
t=[var('t'+str(i)) for i in range(len(param))];t
r=[var(i) for i in param];r
v=vector([sol[0][inc] for inc in x])
for i in range(len(param)):v=v.subs(r[i]==t[i])
print "v=", v
v(t0=1)
[{x0: r26, x1: -2*r26, x2: r26}] r26 [t0] [r26] v= (t0, -2*t0, t0) (1, -2, 1)
#dessin
sol[0]
parametric_plot3d([sol[0][x[i]] for i in range(n)],(-1/2,1/2))
3D rendering not yet implemented
# n comme argument d'une fonction
def f(n):
    y=[var('x'+str(i)) for i in range(n)]
f(10)
show_identifiers()
['x9', 'x2', 'x3', 'x0', 'x6', 'x7', 'x4', 'x5', 'x8', 'x1', 'f']
x=[var('x'+str(i)) for i in range(2)];x
x[1].subs(x1==1)
def f(*x):return(sum(x))
print f,',',f(*[1,2]),',',f(*x)
[x0, x1] 1 <function f at 0x7f10c19ba488> , 3 , x0 + x1