from IPython.display import display, Latex
from IPython.core.display import HTML
%reset -f
%matplotlib inline
%autosave 300
from matplotlib.pylab import *
En important la librairie matplotlib, la librairie numpy est également chargée. Elle permet la manipulation de tableaux et contient des opérations empruntées aux opérations matricielles mathématiques.
On définit des tableaux comme suit, soit par remplissage de tous les coefficients, soit par des tableaux de 1 (ones) ou 0 (zeros).
a=array([[1,2],[3,4]])
b=5*ones((2,2))
print('a=',a)
print('b=',b)
print('a+b=',a+b)
On profite du tableau "zeros" pour intialiser un tableau à la taille souhaitée et modifier les coefficients en position ligne colonne de notre choix. Attention, les indices démarrent à 0 et non 1.
c=zeros((3,3))
c[2,1]=2
print(c)
On peut également affecter un sous-tableau d'un tableau.
d=zeros((3,3))
d[0:2,0:2]=a
print(d)
On peut effectuer des produit composantes par composantes:
print(d*d)
La commande "arange" permet de créer sous forme de vecteur numpy l'équivalent de l'itérateur "range". Il est possible de modifier la forme du tableau pour le transformer (reshape) en un tableau du même nombre d'éléments mais avec un nombre de ligne et colonne différent. Notez bien l'ordre dans lequel sont rangés les éléments.
vec=arange(3,21,2)
print(vec)
mat3x3=vec.reshape(3,3)
print("matrie 3x3 : ")
print(mat3x3)
print("mis en une unique colonne : ")
print(mat3x3.reshape(9,1))
print("mis en une unique ligne : ")
print(mat3x3.reshape(1,9))
print("reshape après transposition, l'ordre est bousculé :")
print(mat3x3.T.reshape(1,9))
On peut utiliser la commande "diag" afin de positionner en diagonale, un vecteur ou un sous-tableau dans un tableau dont on précise la taille. Les autres coefficients non précisés etant nuls. On choisit de combien on décale la diagonale par rapport à la diagonale principale.
vec=arange(1,10)
M=diag(vec,1)+diag(ones(9),-1)
print(M)
La multiplication matricielle utilise la fonction dot:
print(dot(d,d))
print("une variante d'écriture pour le même résultat:")
print(d.dot(d))
On peut également utiliser l'opérateur "@" qui signifie produit matriciel entre deux tableaux de taille compatible:
print(d@d)
Le produit scalaire de deux vecteurs n'est autre que le produit matriciel à une transposition près (produit d'un vecteur ligne par un vecteur colonne).
Le produit matrice vecteur n'est rien d'autre qu'un produit matriciel. Exemples:
vec=array([[1],[2],[3]])
print("vecteur colonne ou tableau 3 lignes 1 colonne")
print(vec)
vectrans=vec.T
print("son transposé:",vectrans)
print("2 produits possibles:")
print(vec@vectrans)
print("le deuxième")
print(vectrans@vec)
print("produit matrice vecteur:")
print(d@vec)
Attention aux tableaux qui ne sont qu'une conversion d'une liste de nombres et non une liste de liste de nombres. Dans ce cas, la notion de ligne ou colonne est ambigue. La transposée n'agit pas, le produit dot ou @ a un sens de produit scalaire.
a=array([1,2,3])
print(a.dot(a),a@a)
print("a=",a)
print("a.T=",a.T)
Coder et afficher les calculs suivants: $$ \left (\matrix{1&2&1\\ 2&1&1\\ 1&1&1}\right)\left (\matrix{1\\ 2\\ 1}\right) $$ $$ \left (\matrix{ 1&0&0&0&0&0&0\\ 0&1&0&0&0&0&0\\ 0&0&2&2&1&0&0\\ 0&0&2&2&1&0&0\\ 0&0&1&1&2&0&0\\ 0&0&0&0&0&1&0\\ 0&0&0&0&0&0&1\\ }\right)+ \left (\matrix{ 1&0&0&0&0&0&0\\ 0&1&0&0&0&0&0\\ 0&0&1&0&0&0&0\\ 0&0&0&1&0&0&0\\ 0&0&0&0&1&0&0\\ 0&0&0&0&0&1&0\\ 0&0&0&0&0&0&1\\ }\right)+ \left (\matrix{ 0&0&0&0&0&0&1\\ 0&0&0&0&0&1&0\\ 0&0&0&0&1&0&0\\ 0&0&0&1&0&0&0\\ 0&0&1&0&0&0&0\\ 0&1&0&0&0&0&0\\ 1&0&0&0&0&0&0\\ }\right) $$
Commençons par le produit des deux premières matrices.
a=array([[1,2,1],[2,1,1],[1,1,1]])
b=array([[1],[2],[1]])
print("le produit des deux matrices est :")
print(dot(a,b))
Poursuivons avec l'addition des trois matrices ci-dessus. Commençons par définir les trois matrices intervenant dans la somme, notée $M1$, $M2$, $M3$ dans l'ordre d'apparition.
M1=diag(ones(7))
M2=M1
M1[2:5,2:5]=M1[2:5,2:5]+a
M3=zeros((7,7))
for i in range(7):
M3[i,6-i]=1
print("Somme =")
print(M1+M2+M3)
Construire et afficher la matrice de taille 9x9 qui ne contient que des zéros sauf des $-4$ sur la deuxième sous-diagonale.
vec=-4*ones(7)
M=diag(vec,-2)
print(M)
La fonction diag a permis de placer le vecteur "vec" sur la deuxième sous-diagonale de la matrice. On remarque que "vec" est de longueur 7 quand la matrice est de taille 9x9.on a utiliser la fonction "ones" pour créer le vecteur constant composé de 1, multiplié par les nombres souhaités.
Calculer la somme des carrés des 60 premiers entiers non nuls à l'aide d'un vecteur et de la commande "dot". Calculer la somme de ces entiers toujours à l'aide de vecteurs et de "dot" ou "@".
vec=arange(1,61,1)
a=dot(vec,vec)# ou a=vec@vec
print("La somme des carrés des 60 premiers entiers non nuls est égale à :",a)
b=ones(60)
c=dot(vec,b)
print("La somme des 60 premiers entiers non nuls est égale à :",c)
DOT permet de faire le produit matriciel lorsqu'il s'agît de matrice (tableaux à 2 indices). Néanmoins, lorsqu'il s'agît de vecteurs (tableaux à 1 indice), il calcule le produit scalaire.
On peut factoriser, sous quelques hypothèses, une matrice $A$ sous la forme $A=PLU$ où $L$ (lower) est une matrice triangulaire inférieure avec des $1$ sur la diagonale, $U$ (upper) une matrice triangulaire supérieure et $P$ une matrice de permutation. Cette factorisation correspond aux opérations de transformation de la matrice par la méthode de Gauss avec pivot. Le but étant de résoudre alors un système linéaire dont $A$ est la matrice qu'on supposera alors carrée. L'hypothèse $A$ inversible suffit pour assurer cette factorisation. La méthode de Gauss (sans pivot, $P=Id$) s'applique dès lors que toutes les sous-matrices de $A$ sont inversibles.
from scipy import linalg
#U=triu([[1,2,3],[4,5,6],[7,8,9]])
#L=tril([[1,2,3],[4,5,6],[7,8,9]])
A=array([[1,2,3],[4,5,6],[7,8,10]])
P,L,U=linalg.lu(A)
print(P)
print()
print(L)
print()
print(U)
print()
print("vérification :")
print(P.dot(L.dot(U)))
En n'utilisant que la transposition et un algorithme de résolution de système triangulaire (inférieur ou supérieur) que vous coderez avec une seule boucle, résoudre le système linéaire $Ax=b$ à partir de $P$, $L$, $U$ de la factorisation $A=PLU$. Pour vérifier votre calcul, cet outil est déjà codé dans scipy avec la commande suivante:
#Forme compactee de la factorisation PLU
lu=linalg.lu_factor(A)
#appel de la resolution du système depuis la factorisation lu
b=ones((3,1))
x=linalg.lu_solve(lu,b)
print("x est le vecteur : ")
print(x)
print("on vérifie Ax=b : ")
print(A@x)
# ECRIRE VOTRE PROPRE ALGO QUI TROUVE x depuis P,L,U,b
#sans savoir ce que scipy choisit comme solveur,
#il existe la commande de résolution de systeme linéaire
print(linalg.solve(A,b))
On cherche à resoudre le système $Ax=b$.
Or $A=PLU$, on obtient l'équation $PLUx=b$ ce qui équivaut à $LUx=P^tb=z$ car la matrice de permutation est orthogonale.
On pose alors $y=Ux$ ; $y$ satisfait $Ly=z$.
Il faudra ensuite résoudre $Ly=z$ puis $Ux=y$
#Première étape, LUx=z
pt=P.T
z=pt@b
print("la transposée de P est ;")
print(z)
#Deuxième étape, résoudre Ly=z, avec y=Ux
n=len(b)
y=zeros((n,1))
for i in range (n):
y[i]=z[i]-L[i,0:i].dot(y[0:i,0])
print(y)
print()
print("vérifions si Ly=z :")
print(L@y)
print("l'égalité est vérifiée")
#troisième étape, résoudre Ux=y
n=len(b)
x=zeros((n,1))
for i in range(n-1,-1,-1):
x[i,0]=(y[i,0]-U[i,i+1:n].dot(x[i+1:n,0]))/U[i,i]
print(x)
Lorsque l'on obtient un résultat à $10^{-16}$, on considère que ça vaut 0. On obtient alors le bon résultat ici.
La commande "eig" permet de disposer des valeurs propres et des vecteurs propres associés. La commande "inv" permet d'inverser la matrice. On ne s'en sert pas pour résoudre un système linéaire, c'est beaucoup trop cher en calcul!
eig(A)
inv(A)
Calculer les valeurs propres de la matrice $A$ de taille $(20,20)$ constituée de $3$ sur la diagonale principale, de $-1$ sur la diagonale supérieure et inférieure.
vec=3*ones(20)
M=diag(vec)
vect=-1*ones(19)
N=diag(vect,1)
O=diag(vect,-1)
P=M+N+O
eig(P)
Pourquoi la matrice est inversible? Les valeurs propres sont elles réelles et strictement positives?
Supposons la matrice A(20.20) diagonalisable elle possède donc 20 valeurs propres. Prouvons que ces valeurs propres sont réelles. On définit les valeurs propres a priori complexes et on montrera qu'elles sont réelles. On sait que $A x_{\lambda} = \lambda x_{\lambda}$ où $x_{\lambda}$ est le vecteur propre associé à la dite valeur propre $\lambda$. Soit $<x,y>$ le produit scalaire dans $\mathbb{C}$ (produit hermitien) tel que $$<x,y> = \sum x_{i} \overline{y_{i}}$$ avec $\overline{y_{i}}$ le conjugué de $y_{i}$. On a alors que $$<A x_{\lambda},x_{\lambda}> = \sum_{i=1}^{20}\lambda x_{\lambda_i}\overline{x_{\lambda_i}}.$$ Or $\lambda$ ne dépendant pas de la somme nous pouvons l'en sortir, i.e : $$<A x_{\lambda},x_{\lambda}> = \lambda \sum_{i=1}^{20}x_{\lambda_i}\overline{x_{\lambda_i}}=\lambda\|x_{\lambda}\|^2.$$ D'où $$<A x_{\lambda},x_{\lambda}> = \lambda \|x_{\lambda}\|^2.$$ Démonstration similaire dans l'autre sens : on a $$<x_{\lambda},A x_{\lambda}> = \sum_{i=1}^{20}x_{\lambda_i} \overline{\lambda x_{\lambda_i}}.$$ Donc $$<x_{\lambda},A x_{\lambda}> = \sum_{i=1}^{20}x_{\lambda_i}\overline{\lambda}\overline{x_{\lambda_i}}$$ On a alors que : $$<x_{\lambda},A x_{\lambda}> = \overline{\lambda} \sum_{i=1}^{20}x_{\lambda_i}\overline{x_{\lambda_i}}.$$ D'où $$<x_{\lambda},A x_{\lambda}> = \overline{\lambda} \|x_{\lambda}\|^2.$$ Comme $A$ est symétrique à coefficient réel, $$ <A x_{\lambda},x_{\lambda}>=< x_{\lambda},A x_{\lambda}>. $$ Donc on a l'égalité $\lambda\|x_{\lambda}\|^2 = \overline{\lambda}\|x_{\lambda}\|^2$ qui impose que $\lambda = \overline{\lambda}$.
Les valeurs propres sont donc forcément réelles. CQFD
Prouvons maintenant que ces valeurs propres sont strictement positives, on pose $A x_{\lambda}=\lambda x_{\lambda}$ avec $$x_{\lambda} \ne 0_{\mathbb{R^{20}}}$$ par définition d'un vecteur propre. On note $$x_{l} = max_{i\in[0;20]_{\mathbb{N}}} x_{i}.$$ Quitte à changer $x_{\lambda}$ en $-x_{\lambda}$, on peut supposer que $x_{l}>0$.
d'où pour la ligne $l$ qui n'est ni à la dernière ligne ni à la première,
on a l'égalité suivante, $$3 x_{l}-x_{l-1}-x_{l+1} = \lambda x_l.$$ Or, $2 x_{l}-x_{l-1}-x_{l+1}\ge 0$, ainsi $0<x_l\le \lambda x_l$. On conclut donc que $ x_l(\lambda -1)>0$ et donc que $\lambda >1>0$.
On note que si $l=0$ ou $l=20$, on a {\sl a fortiuri} $0<x_l\le \lambda x_l$ et on conclut de même.
On retiendra qu'une matrice symétrique à diagonale strictement dominante a nécessairement des valeurs propres réelles strictement positives.
Dans les exercices qui suivent, la solution numérique proposée ne doit faire appel à aucune boucle, mais s'appuyer sur des produits matriciels ou scalaires entre tableaux de type array (du module numpy)
1- Calculer $$ \sum_{i=1}^{49}(-u_{i-1}+2u_i-u_{i+1})u_i-f(u_i)u_i, $$ où $U=(u_0,u_1,\cdots,u_{50})$ est donné par $$ u_i=\sin(i)\frac {i^2}{1+i^2}, $$ où $f$ est la fonction définie par $$ f(x)=\sin(x)\cos(x^2). $$
def u(x):
return sin(x)*(x**2/(1+x**2))
def f(x):
return sin(x)*cos(x**2)
a=arange(0,49)
b=arange(1,50)
c=arange(2,51)
d=(-u(a)+2*u(b)-u(c))*u(b)
e=f(u(b))*u(b)
print(sum(d-e))
2- On se donne $U=(u_0,u_1,\cdots,u_{100})$, $V=(v,v_1,\cdots,v_{100})$ avec, $$ u_i=\sin(i)\frac {i^2}{1+i^2},~v_i=cos(i)i. $$ Calculer
1) $$ \sum_{i=0,~i\text{ pair}}^{100}u_iv_{100-i}, $$
2) $$ \sum_{i=0,~i\text{ impair}}^{97}u_iv_{i+3}. $$
def v(x):
return cos(x)*x
a1=arange(0,101,2)
a2=100-a1
a3=u(a1)@v(a2)
print("1)")
print(a3)
b1=arange(1,98,2)
b2=b1+3
b3=u(b1)@v(b2)
print("2)")
print(b3)