chap2-exercice-valeurs-propres3D

Contenu

• 1. Calcul des valeurs propres
  • 1.1. Calcul des invariants
• 2. Intégration numérique
• 3. expo
  • 3.1. calcul avec boucle for

1. Calcul des valeurs propres¶

$$ \det(M - \lambda I) = 0 $$ Pour $n = 3$, on obtient un polynôme caractéristique de la forme : $$-\lambda^3 + I_1 \lambda^2 - I_2 \lambda + I_3 = 0$$

avec : $I_1$,, $I_2$ et $I_3$ sont appelés les invariants de la matrice

$$ \begin{align} I_1 &= \mathrm{Trace}(M) = M_{ii} \\ I_2 &= \frac{1}{2} (I_1^2-\mathrm{Trace}(M^2)) = \frac{1}{2} (M_{ii} M_{jj} - M_{ij} M_{ji}) \\ I_3 &= \mathrm{det}(M) = \varepsilon_{ijk} M_{1i} M_{2j} M_{3k} \end{align} $$

La résolution de l'équation cubique en $\lambda$ se fera avec la méthode de Newton-Raphson $$ \begin{align} f(\lambda) &= -\lambda^3 + I_1 \lambda^2 - I_2 \lambda + I_3 \\ f'(\lambda) &= -3\lambda^2 + 2I_1 \lambda - I_2 \end{align} $$

$$ \lambda_{k+1} = \lambda_k - f( \lambda_k ) / f'(\lambda_{k} ) $$

In [1]:

# Définition de la matrice M
M = [
    [8,  5,  3],    # O
    [5,  7,  6],    # 1
    [3,  6,  9]]    # 2
#    0   1   2    indices

M=[ [2, 1, 3],
    [1, 6, 1],
    [3, 1, 2]]
print(M)

[[2, 1, 3], [1, 6, 1], [3, 1, 2]]

In [2]:

# Test si n != 3 on soulève une erreur
n = len(M)
if n !=3:
    raise ValueError('n doit être égal à 3')

1.1. Calcul des invariants¶

In [3]:

# Calcul de I1 = M_ii
I1 = sum(M[i][i] for i in range(n))
print(f'{I1 = }')

# Calcul de I2 =  (M_{ii} M_{jj} - M_{ij} M_{ji}) / 2
I2 = sum(M[i][i] * M[j][j] - M[i][j] * M[j][i] for i in range(n) for j in range(n)) / 2
print(f'{I2 = }')

# Calcul de I3 =  eps_{ijk} M_{1i} M_{2j} M_{3k}
# on a besoin de epsilon de Levi-Civita
eps = [ [ [(i-j)*(j-k)*(k-i)//2 for k in range(3)] for j in range(3)] for i in range(3)]

I3 = sum(eps[i][j][k] * M[0][i] * M[1][j] * M[2][k]           # attention aux indices
          for i in range(n) for j in range(n) for k in range(n))
print(f'{I3 = }')

I1 = 10
I2 = 17.0
I3 = -28

Résolution avec Newton-Raphson¶

Remarques:

• Comme lambda est un mot clé de python, on nommera une variable x à la place
• On prendra $\frac{1}{3}I_1$ comme valeur initiale

En principe, ça ne demande que très peu d'itérations, donc elle est intéressant même pour un calcul manuel.

Par exemple pour une précision à 2 chiffres, on n'a que 3 itérations

In [4]:

# Méthode de Newton-Raphson pour l'extraction 
# d'une valeur propre d'une matrice 3x3

tol = 1e-6
it, maxit = 0, 100 

x, x1 = 1e30, I1/3    # initialisation de x et x1
while abs(x1-x) > tol and it < maxit:
    x = x1
    f  = x * ( x * (I1 - x) - I2 ) + I3  # Algo de Ruffini-Horner
    df = x * ( 2*I1 - 3*x ) - I2         # dérivée  
    x1 = x -  f/df
    it+= 1
    print(f'{it = }, {x1 = :.8f}')


lambda_1 = x1
print('itérations :', it)
print('valeur propre = ', lambda_1 )
#print(M)

it = 1, x1 = 3.98185941
it = 2, x1 = 3.99995712
it = 3, x1 = 4.00000000
it = 4, x1 = 4.00000000
itérations : 4
valeur propre =  4.0

Calcul des 2 autres valeurs¶

Une fois on obtient $\lambda_1$ comme solution avec Newton-Raphson, on réduit le polynôme en :

$$ f(\lambda) = (\lambda -\lambda_1) (a\lambda^2 + b \lambda + c)$$ on obtient : $$ a = -1 \quad; \quad b = I_1 - \lambda_1 \quad\quad c = b \lambda_1 - I_2 $$

Reste plus qu'à trouver les racines du polynôme d'ordre 2

In [5]:

# Calcul des 2 autres valeurs propres
a = -1 
b = I1 - lambda_1 
c = b*lambda_1 - I2

delta = b*b - 4*a*c
if delta < 0 :        #on soulève une erreur
    raise ValueError('Discriminant négatif')
lambda_2 = (- b + delta**0.5)/(2*a)
lambda_3 = (- b - delta**0.5)/(2*a) 

print(f'{lambda_1 = :.4f}')
print(f'{lambda_2 = :.4f}')
print(f'{lambda_3 = :.4f}')

lambda_1 = 4.0000
lambda_2 = -1.0000
lambda_3 = 7.0000

In [6]:

#tri des valeurs a < b < c
def tri(a,b,c):
    if a>b : a,b = b,a
    if a>c : a,c = c,a
    if b>c : b,c = c,b
    return a,b,c

sig1,sig2,sig3 = tri(lambda_1,lambda_2,lambda_3)
print(sig1,sig2,sig3)

-1.0 4.0 7.0