EDIT: le code est updaté dans son entier dans le 5ème message pour les intéressés.

Bonjour,

je m'excuse d'abord de poster dans cette partie du forum, je ne suis pas sur s'il s'agit plus d'un problème numérique ou mathématique mais comme la compréhension du code en Python est importante pour résoudre mon problème, je me suis dis qu'il valait mieux poster ici.

Alors voila, mon objectif est:

* créer une matrice tridiagonale carrée dont la diagonale principale a une valeur, et dont les deux diagonales l'entourant ont une autre valeur

* en tirer les vecteurs propres

* modéliser les vecteurs propres de façon à avoir: Amplitude de l'élément une fonction de la position de cet élément dans le vecteur propre

Ce que je dois obtenir:

(a) montre la modélisation des vecteurs propres 1, 2, 3 et 4, (b) montre la modélisation du vecteur propre 200, pour une matrice tridiagonale 200x200.

ce que j'obtiens:

La première image représentant le premier vecteur propre, la deuxième du deuxième vecteur propre etc. Et la 5ème image, du 200ème vecteur propre.

On voit immédiatement un déphasage sur les images 2, 3, qui disparaît pour l'image 4.

On voit également que l'enveloppe du signal est rectangulaire au lieu d'être sinusoïdale comme dans l'image du résultat attendu.

Mon code pour cette partie est le suivant:

from matplotlib import pyplot as plt # Import pyplot as plt from the library matplotlib
import numpy as np # Import the library numpy under the name np
from IPython.display import Image # Import Image from IPython's module Display to insert images in the notebook
from scipy.sparse import diags # Here for tridiagonal matrices
import random # Will be used to generate random numbers
%matplotlib inline

# Allows us to plot graphs in the notebook

n = int(input("Size of the matrix: ")) # Size of the matrix
I0 = float(input("Main diagonal's value: ")) # Set value for the main diagonal
I1 = float(input("Secondary diagonals' value: ")) # Set value for the second and third diagonals
print("\n")

A = np.zeros((n,n)) # Matrix without disorder
np.fill_diagonal(A, I0) # Changing value of main diagonal for input value
i = np.arange((n-1)) # Changing value of off-diagonal for input value
A[i,i+1] = I1
A[i+1,i] = I1

print("Ordered matrix: \n", A) # Shows the ordered matrix

eigenvalues, eigenvectors = np.linalg.eig(A) # Get the eigenvectors and eigenvalues of the ordered matrix
print("\n")
print("Matrix of Eigenvectors: \n", eigenvectors) # Matrix of ordered eigenvectors
print("\n")
print("Matrix of eigenvalues: \n", eigenvalues) # Array of ordered eigenvalues
print("\n")

eigenvector = [] # Array of one eigenvector from the eigenvectors matrix

j = int(input("Eigenvector's index: \n")) # Chooses which eigenvector will be plotted
print("\n")

j = j-1 # Set the index to start at 0

for i in range(n): # Append every element of the jth eigenvector to its array for plotting 
    eigenvector.append(eigenvectors[i,j])

print("Array of a single eigenvector: \n", eigenvector) # To check the eigenvector array

plt.plot(eigenvector) # Plot the ordered eigenvector
plt.grid() # Add a grid to the plot

Mes questions est donc:

* Est-ce que je prend les mauvais éléments pour isoler un unique vecteur propre?

* Pourquoi le déphasage dans ce qui est supposé être trivial et retourner le même graph?

* Pourquoi l'enveloppe du graph n'est-elle pas sinusoïdale?

J'espère pouvoir trouver la solution ici, je vous remercie pour avoir jeté un œil au problème.

Merci,

Timothy.

-
Edité par TymeToTry 13 juin 2017 à 20:54:30

Bonjour,

Pour simplifier un petit bout de code tu peux extraire les colonnes de tes vecteurs propres comme ci dessous :

for j in range(0, n):
    eigenvector = eigenvectors[:,j]
    print("Array of a single eigenvector: \n", eigenvector) # To check the eigenvector array
    plt.plot(eigenvector) # Plot the ordered eigenvector

Apres j'ai tester le code, bizarrement j'ai bien certaines fréquences qui oscillent dans une enveloppe sinusoïdale (peut être que c'est lié au valeurs de ta diagonale ou le vecteur que tu affiches. Perso j'ai tester d'afficher les 100 vect. propres d'une matrice de taille n. C'est pas beau mais j'ai bien l'enveloppe sur certaines freq.).

Pour le déphasage par contre... Bonne question

Bonjour,

je te remercie de la réponse, j'ai repris le code en prenant en compte ce que tu as écris, pour donner plus de choix à l'utilisateur.

J'ai remplacé tout à partir de la ligne "eigenvector = []" :

eigenvector = [] # Array of one eigenvector from the eigenvectors matrix

k = int(input("# of eigenvectors graphed: \n")) # Get the number of eigenvectors to graph on a same plot, starting with the 1st eigenvecto
print("\n")

y = [] # Creates an empty list that'll be used to store which eigenvectors will be graphed

for o in range(0, k): # Creates the command to the user that asks which eigenvectors will be graphed
    x = int(input("Which eigenvector do you want to graph?: \n")) # Asks the index of each eigenvectors the user wants to graph
    y.append(x-1) # Append the index of the eigenvector to the list
    print("\n")

for j in y: # For j taking each value in the list of input eigenvectors index
    h = j+1 # Used to store the actual index number of the eigenvector that is being considered
    eigenvector = eigenvectors[:,j] # Creates an array with the value of each element constituting the eigenvector that is being considered
    print("Array of eigenvector #%s: \n" %(h), eigenvector) # To check the eigenvector array
    print("\n")
    plt.plot(eigenvector) # Plot the ordered eigenvector
plt.title("Graph of the matrix' eigenvectors as function of their elements' position") # Graph's title
plt.xlabel("Element position") # x axis' label
plt.ylabel("Arbitrary unit") # y axis' label
plt.yticks([]) # Removes y axis' values since they are arbitrary

J'ai regardé ce que tu disais, à propos d'obtenir l'enveloppe sinusoïdale pour une taille de matrix de 100, sur certaines fréquences seulement. J'ai en effet obtenu ce résultat pour une matrice 100x100, avec le 1er vecteur propre.

Après, l'algèbre et moi ça fait deux à l'heure actuelle, alors je peux pas l'expliquer, ni savoir si je dois réellement obtenir ce résultat. ^.^ Il reste que je n'obtiens pas les résultats escomptés et j'ai toujours ce problème de déphasage.

Merci encore, surtout pour l'astuce au niveau du code.

-
Edité par TymeToTry 13 juin 2017 à 20:53:03

Pour les intéressés, voici le code complet actuel:

from matplotlib import pyplot as plt # Import pyplot as plt from the library matplotlib
import numpy as np # Import the library numpy under the name np
%matplotlib inline
 
# Allows us to plot graphs in the notebook

n = int(input("Size of the matrix: \n")) # Size of the matrix
print("\n")
I0 = float(input("Main diagonal's value: \n")) # Set value for the main diagonal
print("\n")
I1 = float(input("Secondary diagonals' value: \n")) # Set value for the second and third diagonals
print("\n")
 
A = np.eye(n, n, k=-1)*I1 + np.eye(n, n)*I0 + np.eye(n, n, k=1)*I1  # Tridiagonal matrix without disorder

#print("Ordered matrix: \n", A) # Shows the ordered matrix
 
eigenvalues, eigenvectors = np.linalg.eigh(A) # Get the eigenvectors and eigenvalues of the ordered matrix
#print("\n")
#print("Matrix of eigenvectors: \n", eigenvectors) # Matrix of ordered eigenvectors
#print("\n")
#print("Matrix of eigenvalues: \n", eigenvalues) # Array of ordered eigenvalues
#print("\n")

K = int(input("# of eigenvectors graphed: \n")) # Get the number of eigenvectors to graph on a same plot, starting with the 1st eigenvecto
print("\n")

y = [] # Creates an empty list that'll be used to store which ordered eigenvectors will be graphed

for o in range(0, K): # Creates the command to the user that asks which eigenvectors will be graphed
    x = int(input("Which eigenvector do you want to graph?: \n")) # Asks the index of each eigenvectors the user wants to graph
    y.append(x-1) # Append the index of the eigenvector to the list
    print("\n")

h = float(input("Max randomness: \n")) # Coefficient setting the maximum randomness
print("\n")

isvalid = False # Will test that the following user input is correct
while not isvalid: # not will test if isvalid is false, if so, it'll run the following code till isvalid is true
    I2bis = (input("Disorder in the main diagonal? (y/n) \n")) # Will decide what coefficient I2 will be (0 or 1)
    I2bis = I2bis.upper() # Ensures uniformity in the way the command is written
    if I2bis == 'YES' or I2bis == '1' or I2bis == 'Y': # Checks for positive input of the argument
        I2 = 1 # Will ensure the main diagonal of the matrix 't' is made of ones
        isvalid = True # Stops the loop
    elif I2bis == 'NO' or I2bis == '0' or I2bis == 'N': # Checks for negative input of the argument
        I2 = 0 # Will ensure the main diagonal of the matrix 't' is made of zeros
        isvalid = True # Stops the loop
    else: print("Invalid entry, please type a valid entry (y, yes, 1, n, no, 0).") # If input is incorrect, rerun the loop
    print("\n")

isvalid = False # Will test that the following user input is correct
while not isvalid: # not will test if isvalid is false, if so, it'll run the following code till isvalid is true
    I3bis = (input("Disorder in the secondary diagonals? (y/n) \n")) # Will decide what coefficient I2 will be (0 or 1)
    I3bis = I3bis.upper() # Ensures uniformity in the way the command is written
    if I3bis == 'YES' or I3bis == '1' or I3bis == 'Y': # Checks for positive input of the argument
        I3 = 1 # Will ensure the secondary diagonals of the matrix 't' is made of ones
        isvalid = True # Stops the loop
    elif I3bis == 'NO' or I3bis == '0' or I3bis == 'N': # Checks for negative input of the argument
        I3 = 0 # Will ensure the secondary diagonals of the matrix 't' is made of zeros
        isvalid = True # Stops the loop
    else: print("Invalid entry, please type a valid entry (y, yes, 1, n, no, 0).") # If input is incorrect, rerun the loop
    print("\n")

d = np.random.rand(n,n)*h # Generates a matrix of random numbers in range of 0 to 1, multiplied by the actual max randomness
#print("Random matrix: \n", d) # Shows the array of random numbers
#print("\n")

t = np.eye(n, n, k=-1)*I3 + np.eye(n, n)*I2 + np.eye(n, n, k=1)*I3  # Unitary matrix that'll multiply the random matrix to get a disordered tridiagonal matrix
#print("Multiplicative (t-)matrix to multiply the random matrix: \n", t)
#print("\n")

w = d*t # (Tri)diagonal random matrix
#print("Diagonal random matrix: \n", w) # Checks the (tri)diagonal random matrix
#print("\n")

C = A + w # Add the random matrix and the ordered matrix together
#print("Localized matrix: \n", C) # Show the localized matrix
#print("\n")

p, q = np.linalg.eigh(C) # Get the eigenvectors and eigenvalues of the localized matrix

eigenvector = [] # Array of one eigenvector from the ordered eigenvectors matrix
Q = [] # Array of one eigenvector from the random eigenvectors matrix

idx = eigenvalues.argsort()[::-1]  # The 2 following blocks are here to return the order
eigenvalues = eigenvalues[idx] # in which the eigenvectors are indexed so that the graph will show
eigenvectors = eigenvectors[:,idx] # eigenvector #1 as eigenvector #1 and not eigenvector #200 etc.
 
idx = p.argsort()[::-1] # 2nd block. This line sorts the indexes in reverse
p = p[idx] # This line sorts the eigenvalues in reverse
q = q[:,idx] # This line sorts the eigenvectors in reverse

for j in y: # For J taking each value in the list of input eigenvectors index
    h = j+1 # Used to store the actual index number of the eigenvector that is being considered
    eigenvector = eigenvectors[:,j] # Creates an array with the value of each element constituting the eigenvector that is being considered
#    print("Array of ordered eigenvector #%s: \n" %(h), eigenvector) # To check the eigenvector array
#    print("\n")
    plt.plot(eigenvector, label="Ordered eigenvector #%s" %(h))  # Plot the ordered eigenvector
    Q = q[:,j] # Creates an array with the value of each element constituting the eigenvector that is being considered
#    print("Array of disordered eigenvector #%s: \n" %(h), Q) # To check the eigenvector array
#    print("\n")
    plt.plot(Q, label="Disordered eigenvector #%s" %(h)) # Plot the ordered eigenvector

if K != 0: # So that it doesn't return an empty graph if K = 0
    #plt.grid() # Add a grid to the graph
    plt.title("Graph of the matrix' eigenvectors as function of their elements' position") # Graph's title
    plt.xlabel("Element position") # x axis' label
    plt.ylabel("Arbitrary unit") # y axis' label
    plt.yticks([]) # Removes y axis' values since they are arbitrary
    plt.legend(loc="upper left", bbox_to_anchor=(1,1)) # Legend is outside the plot
else: # If K = 0, do nothing
    pass

EDIT: j'ai rajouté la partie de la matrice random etc. etc. que j'ai pu terminer grâce a coni63. Il me reste cependant à comprendre pourquoi mes vecteurs propres qui sont normalement #1, #2, #3... deviennent #200, #199, #198... et rajouter la partie largeur de localisation, etc. etc. qui constitue la suite du stage. Du coup si jamais, je ré-update le post, ou je reposterai à nouveau selon combien de temps ça fait.

-
Edité par TymeToTry 19 juin 2017 à 11:51:50

Bonjour,

Je viens de tester avec une matrice avec la diag. haute différente de la diag. basse. Dans ce cas, ça donne des vecteurs propres amorties... Je ne sais pas pourquoi mais bon peut etre que ca peut t'eclairer

D'ailleurs on peut aussi simplifier la façon de générer la matrice tridiagonale :

n = 100
I0 = 5
I1 = 2
I2 = 7
 
A = np.eye(n, n, k=-1)*I1 + np.eye(n, n)*I0 + np.eye(n, n, k=1)*I2

EDIT 1 : Ah c'est bon j'ai trouvé

Comme ta matrice est symmetrique, il faut utiliser eigh et non eig:

Tu peux tester le code ci dessous :

from matplotlib import pyplot as plt # Import pyplot as plt from the library matplotlib
import scipy
import numpy as np # Import the library numpy under the name np
import random # Will be used to generate random numbers

n = 150
I0 = 5
I1 = 2
I2 = 2
sym = True

A = np.eye(n, n, k=-1)*I1 + np.eye(n, n)*I0 + np.eye(n, n, k=1)*I2

print("Ordered matrix: \n", A) # Shows the ordered matrix

if sym:
    eigenvalues, eigenvectors = np.linalg.eig(A) # Get the eigenvectors and eigenvalues of the ordered matrix
    print("\n")
    print("Matrix of Eigenvectors: \n", eigenvectors) # Matrix of ordered eigenvectors
    print("\n")
    print("Matrix of eigenvalues: \n", eigenvalues) # Array of ordered eigenvalues
    print("\n")
else:
    eigenvalues, eigenvectors = np.linalg.eigh(A) # Get the eigenvectors and eigenvalues of the ordered matrix
    print("\n")
    print("Matrix of Eigenvectors: \n", eigenvectors) # Matrix of ordered eigenvectors
    print("\n")
    print("Matrix of eigenvalues: \n", eigenvalues) # Array of ordered eigenvalues
    print("\n")
  
# for j in range(0, 1, 1):
#     eigenvector = eigenvectors[:,j]
#     print("Array of a single eigenvector: \n", eigenvector) # To check the eigenvector array
#     plt.plot(eigenvector) # Plot the ordered eigenvector

f, axarr = plt.subplots(2, 2)
x = np.linspace(0, n, n)

axarr[0, 0].plot(x, eigenvectors[:,0])
axarr[0, 0].set_title('Subplot 1: 1er Vect')

axarr[0, 1].plot(x, eigenvectors[:,1])
axarr[0, 1].set_title('Subplot 2: 2nd Vect')

axarr[1, 0].plot(x, eigenvectors[:,4])
axarr[1, 0].set_title('Subplot 3: 5eme Vect')

axarr[1, 1].plot(x, eigenvectors[:,100])
axarr[1, 1].set_title('Subplot 4: 101eme Vect')

#plt.grid() # Add a grid to the plot
plt.show()

Et ca donne ca !

https://img15.hostingpics.net/pics/46567120170614171942Figure1.png

-
Edité par coni63 14 juin 2017 à 17:20:56

Bonjour,

je te remercie pour tes infos, je regarde ca ce soir et je te tiens informé! j'apprécie particulièrement la manière dont tu codes la matrice tridiagonale.

Bonjour,

Suite a ton update de post, j'ai regardé pourquoi il y a un mix de vecteurs propres... Comme un vecteur propre "ne peut pas être trié" mais est lié uniquement a sa valeur propre. Je pense que si l'ordre des valeurs propre n'est pas la même alors les vecteurs propres ne le seront pas. 

Voila la représentation de mes valeurs propres pour A et C en fonction de la colonne

Si tu tries les vecteurs propres par rapport a leur valeurs propres décroissante, ca marche. Pour cela tu peux utiliser : 

idx = eigenValues.argsort()[::-1]   
eigenValues = eigenValues[idx]
eigenVectors = eigenVectors[:,idx]

 Apres tri, ca donne ca :

Quand au resultats je pense que l'on a ce que tu recherches car le 1er et dernier vect propres ont la même enveloppe mais une "fréquence incluse" ca doit etre pareil pour tous les vecteur propres "opposés" (a verifier...)

Autres petites améliorations que tu peux apporter :

- C n'est pas symétrique a cause du random donc utilise eig et non eigh qui est fait uniquement pour les matrices symétriques.

- Tu peux remplacer chaque bloc :

isvalid = False # Will test that the following user input is correct
while not isvalid: # not will test if isvalid is false, if so, it'll run the following code till isvalid is true
    I3bis = (input("Disorder in the secondary diagonals? (y/n) \n")) # Will decide what coefficient I2 will be (0 or 1)
    I3bis = I3bis.upper() # Ensures uniformity in the way the command is written
    if I3bis == 'YES' or I3bis == '1' or I3bis == 'Y': # Checks for positive input of the argument
        I3 = 1 # Will ensure the secondary diagonals of the matrix 't' is made of ones
        isvalid = True # Stops the loop
    elif I3bis == 'NO' or I3bis == '0' or I3bis == 'N': # Checks for negative input of the argument
        I3 = 0 # Will ensure the secondary diagonals of the matrix 't' is made of zeros
        isvalid = True # Stops the loop
    else: print("Invalid entry, please type a valid entry (y, yes, 1, n, no, 0).") # If input is incorrect, rerun the loop
    print("\n")

par :

I3bis = "Yes" # Will decide what coefficient I2 will be (0 or 1)
if I3bis.upper() in ['YES', '1', 'Y']: 
    I3 = 1 # Will ensure the secondary diagonals of the matrix 't' is made of ones
else: 
    I3 = 0 # Will ensure the secondary diagonals of the matrix 't' is made of zeros

Si l'utilisateur ne mets pas oui, alors c'est non et le in est plus joli

Voila le code si besoin (désolé mais je suis faignant donc j'ai supprimé quelques inputs etc

from matplotlib import pyplot as plt # Import pyplot as plt from the library matplotlib
import scipy
import numpy as np # Import the library numpy under the name np
import numpy.random # Necessaire mais je comprends pas pourquoi ^_^
 
n = 200
I0 = 1
I1 = 0.1
K = 2 # Get the number of eigenvectors to graph on a same plot, starting with the 1st eigenvecto
y = [0, 199] # Creates an empty list that'll be used to store which ordered eigenvectors will be graphed
h = 0.01 # Coefficient setting the maximum randomness

I2bis = "Yes" # Will decide what coefficient I2 will be (0 or 1)
if I2bis in ['YES', '1', 'Y']: 
    I2 = 1 # Will ensure the main diagonal of the matrix 't' is made of ones
else:
    I2 = 0 # Will ensure the main diagonal of the matrix 't' is made of zeros

I3bis = "Yes" # Will decide what coefficient I2 will be (0 or 1)
if I3bis.upper() in ['YES', '1', 'Y']: 
    I3 = 1 # Will ensure the secondary diagonals of the matrix 't' is made of ones
else: 
    I3 = 0 # Will ensure the secondary diagonals of the matrix 't' is made of zeros

A = np.eye(n, n, k=-1)*I1 + np.eye(n, n)*I0 + np.eye(n, n, k=1)*I1  # Tridiagonal matrix without disorder
eigenvalues, eigenvectors = np.linalg.eigh(A) # Get the eigenvectors and eigenvalues of the ordered matrix
print(eigenvalues)

d = np.random.rand(n, n)*h # Generates a matrix of random numbers in range of 0 to 1, multiplied by the actual max randomness
t = np.eye(n, n, k=-1)*I3 + np.eye(n, n)*I2 + np.eye(n, n, k=1)*I3  # Unitary matrix that'll multiply the random matrix to get a disordered tridiagonal matrix
 
w = d*t # (Tri)diagonal random matrix
 
C = A + w # Add the random matrix and the ordered matrix together
p, q = np.linalg.eig(C) # Get the eigenvectors and eigenvalues of the localized matrix
print(p)

#tri des vecteurs propres en fct des valeurs propres
idx = eigenvalues.argsort()[::-1]   
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:,idx]

idx = p.argsort()[::-1]   
p = p[idx]
q = q[:,idx]

# plt.plot(eigenvalues)
# plt.plot(p)

for j in y: # For J taking each value in the list of input eigenvectors index
    plt.plot(eigenvectors[:,j], label="Ordered eigenvector #%s" %(j+1))  # Plot the ordered eigenvector
    plt.plot(q[:,j], label="Disordered eigenvector #%s" %(j+1)) # Plot the ordered eigenvector


plt.title("Graph of the matrix' eigenvectors as function of their elements' position") # Graph's title
plt.xlabel("Element position") # x axis' label
plt.ylabel("Arbitrary unit") # y axis' label
plt.yticks([]) # Removes y axis' values since they are arbitrary
plt.legend(loc="upper left", bbox_to_anchor=(1,1)) # Legend is outside the plot

plt.show()

-
Edité par coni63 16 juin 2017 à 14:29:01

Petit point aussi qui peux te servir pour les questions que tu pose a l'utilisateur, tu peux utiliser distutils.util.strtobool comme suit:

try:
    question=input(" foo? ")
    i = not strtobool(question)
    break
except ValueError:
    print(" Please only enter yes or no. ")

strtobool retourne True pour ['y', 'yes', 't', 'true', 'on', '1'] et False pour ['n', 'no', 'f', 'false', 'off', '0']

Salut!

Désolé j'avais plus le temps de m'occuper de ça ces derniers jours. Je regarde tes dernières updates pour optimiser le code.

Par rapport a eig vs eigh quand il y a du random dans la matrice, en fait pour obtenir l'effet désiré il faut utiliser eigh. Eigh annule la composante complexe, et à priori ca donne ce que l'on est supposé trouvé, me demande pas comment ca se fait par contre, je sais pas ^.^

Pour réarranger la liste des valeurs propres qui donne une implication sur les vecteurs propres, j'y avais vraiment pas pensé, j'allais passer par un moyen détourné pour renverser la matrice vecteur propre ou quelque chose du style. ^.^ Merci beaucoup pour ton input.

Par rapport à la boucle if/elif/else et celle avec in, j'ai préféré faire comme ca car je souhaite que ca redemande à l'utilisateur d'envoyer une donnée correcte si l'input est pas bon, ce qui est pas possible avec I3bis = "Yes"...

Je test tout ça dès que j'arrive au bureau!

Pour le strtobool, ca me retourne une erreur quand je test ce code:

  File "<ipython-input-1-16dbf690c887>", line 4
    break
    ^
SyntaxError: 'break' outside loop

P.S. : T'inquiètes pas, les prints étaient là pour vérifier que j'avais les bonnes matrices et que je créais les bonnes listes de vecteurs propres. Ca n'a aucune utilité sinon.

-
Edité par TymeToTry 19 juin 2017 à 11:53:14