114 lines
3.3 KiB
Python
Executable File
114 lines
3.3 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
|
|
"""
|
|
auteur: Dylan Voisin - M1 Informatique Luminy
|
|
MNI: TP 4
|
|
"""
|
|
|
|
import numpy as np
|
|
from numpy.polynomial import polynomial as P
|
|
from copy import deepcopy as dp
|
|
from math import isclose
|
|
|
|
def vect_propre(A, µ=1, nb_iter=15):
|
|
I = np.identity(A.shape[0])
|
|
v = np.ones(A.shape[0])
|
|
inv = np.linalg.inv(A - µ*I)
|
|
for i in range(nb_iter):
|
|
v = inv @ v / np.linalg.norm(inv@v)
|
|
return v
|
|
|
|
def vect_propre2(A, eps = 0.1):
|
|
v = np.random.random(A.shape[0])
|
|
v /= np.linalg.norm(v)
|
|
lambda_ = v.transpose() @ A @ v
|
|
eps = 1e-13
|
|
écart = float('inf')
|
|
while écart > eps:
|
|
w = lambda_ * A @ v
|
|
v = w / np.linalg.norm(w)
|
|
old = lambda_
|
|
lambda_ = (v.transpose() @ A @ v) / (v.transpose() @ v)
|
|
écart = abs(old - lambda_)
|
|
return v
|
|
|
|
def vect_propre_généralisé(A, eps = 1-9):
|
|
A_ = A @ A.transpose()
|
|
v = vect_propre2(A_, eps)
|
|
v = (v.transpose() @ A @ v) / (v.transpose() * v)
|
|
return v
|
|
|
|
def eval_poly(p, x):
|
|
"""un polynome sera représenté par une liste pleine de ses coefficients
|
|
f = a_0 + a_1 * x ** 1 + … + a_n * x ** b_n
|
|
[a_0, a_1, a_2, …, a_n]
|
|
"""
|
|
return sum(p[i] * x ** i for i in range(len(p)))
|
|
|
|
def dichotomie(size, step):
|
|
l = [0, size]
|
|
lc = dp(l)
|
|
for i in range(step):
|
|
for j in range(len(l)-1):
|
|
m = (l[j+1] + l[j]) / 2
|
|
lc.insert(2*(j+1)-1, m)
|
|
yield m
|
|
l = dp(lc)
|
|
|
|
def newton(p, epsilon_f = 1e-13, epsilon_c = 1e-9):
|
|
"""peut faire une boucle infinie à cause de l'imprécision
|
|
à régler avec les paramètres epsilon, mais l'imprécision va
|
|
se répercuter sur les valeurs des racines"""
|
|
x = np.random.random() * 100
|
|
pp = np.polyder(p[::-1])[::-1]
|
|
roots = []
|
|
while True:
|
|
# toujours de l'imprécision à gérer
|
|
if abs(eval_poly(p, x)) < epsilon_f:
|
|
roots.append(x)
|
|
p = P.polydiv(p, (-x, 1))[0]
|
|
# on purge l'imprécision
|
|
p = list(map(lambda x: x if abs(x) > epsilon_c else 0, p))
|
|
pp = np.polyder(p[::-1])[::-1]
|
|
if len(p) == 1:
|
|
return roots # f = c
|
|
fpx = eval_poly(pp, x)
|
|
if fpx == 0:
|
|
print(x)
|
|
# extremum (local ou non)
|
|
old = x
|
|
step = 5
|
|
maxrange = 1000
|
|
eps = 1*10**(-8)
|
|
sg = np.sign(eval_poly(pp, x-eps))
|
|
sd = np.sign(eval_poly(pp, x+eps))
|
|
for diff in dichotomie(maxrange, step):
|
|
if np.sign(eval_poly(pp, x-diff)) != sg:
|
|
x = x - diff
|
|
break
|
|
elif np.sign(eval_poly(pp, x+diff)) != sd:
|
|
x = x + diff
|
|
break
|
|
if x == old:
|
|
return roots # on considère qu'on a toutes les racines
|
|
x -= eval_poly(p, x) / eval_poly(pp, x)
|
|
return roots
|
|
|
|
if __name__ == '__main__':
|
|
m = np.array(
|
|
[
|
|
[2, 3],
|
|
[3, 2]
|
|
]
|
|
)
|
|
v = vect_propre(m)
|
|
v /= np.linalg.norm(v)
|
|
print(v)
|
|
v = vect_propre2(m)
|
|
print(v)
|
|
v = vect_propre_généralisé(m)
|
|
v /= np.linalg.norm(v)
|
|
print(v)
|
|
# p = (0, 1, -2) # -2x²+x
|
|
p = (30, -19, -15, 3, 1) # (x+5)(x-3)(x-1)(x+2)
|
|
print(newton(p)) |