170 lines
4.0 KiB
Python
170 lines
4.0 KiB
Python
from dataclasses import dataclass
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
|
|
|
|
@dataclass
|
|
class Point:
|
|
x: float
|
|
y: float
|
|
|
|
|
|
def interpolation_lagrange(points: list[Point], x: float, show_polynome: bool=False) -> float:
|
|
"""
|
|
|
|
:param points: Les points de la courbe représentant une fonction
|
|
:param x: le x en entrée de la fonction où on veut interpoler
|
|
:return: L'interpolation
|
|
"""
|
|
|
|
interpolation: float = 0.
|
|
n = len(points) # nombre de points
|
|
for i in range(n):
|
|
term: float = 1. # 0.0 will always give 0, it's a product ! So lets put 1
|
|
# Lagrange L(x)
|
|
for k in range(n):
|
|
if k != i:
|
|
term *= ((x - points[k].x) / (points[i].x - points[k].x))
|
|
|
|
# erreur ici, voir le chat
|
|
if show_polynome:
|
|
print(f'P(x) = (({x} - {points[k].x}) / ({points[i].x} - {points[k].x}))')
|
|
|
|
# interpolation with a piece of Lagrange (y_i)
|
|
interpolation += term * points[i].y
|
|
|
|
return interpolation
|
|
|
|
|
|
def calc_coef(x: list[float],y: list[float]) -> float:
|
|
"""
|
|
Renvois f[x0,x1,...,xk] (diff divisé de Newton)
|
|
avec la relation récursive :
|
|
f[x0...xk] = (f[x1...xk] - f[x0..x{k-1}] / (xk - x0)
|
|
(calcul alpha)
|
|
:param x:
|
|
:param y:
|
|
:return:
|
|
"""
|
|
if len(x) == 1:
|
|
return y[0]
|
|
elif len(x) == 2:
|
|
return (y[1]-y[0]) / (x[1] - x[0])
|
|
return (calc_coef(x[1:], y[1:]) - calc_coef(x[:-1], y[:-1])) / (x[-1] - x[0])
|
|
|
|
|
|
def P_newton(x_i: list[float], y_i: list[float], x: float) -> float:
|
|
"""
|
|
Calcul le polynôme de Newton
|
|
|
|
P_n(x) = y0 + sum_{k=1..n} a_k * prod_{j=0..k-1}(x - x_j),
|
|
avec a_k = f[x0,...,xk] calculé via calc_coef()
|
|
:param x_i: ensemble des x
|
|
:param y_i: ensemble des y
|
|
:param x: là ou on veut interpolé
|
|
:return:
|
|
"""
|
|
p = y_i[0]
|
|
e = 1.
|
|
|
|
for k in range(1, len(x_i)):
|
|
e *= (x - x_i[k-1])
|
|
a_k = calc_coef(x_i[:k+1], y_i[:k+1])
|
|
p += a_k * e
|
|
|
|
return p
|
|
|
|
|
|
def exercice1() -> None:
|
|
points = [Point(-1,0), Point(0,-1), Point(1,0), Point(3,70)]
|
|
|
|
inter_x2 = interpolation_lagrange(points, 2, True)
|
|
inter_x3 = interpolation_lagrange(points, 3, True)
|
|
print(f"X2={inter_x2}, X3={inter_x3}")
|
|
|
|
x = np.linspace(0, 4, 300)
|
|
y_lagrange = [interpolation_lagrange(points, x[i]) for i in range(len(x))]
|
|
|
|
y_newton = []
|
|
|
|
for i in range(4):
|
|
y_newton.append([P_newton([point.x for point in points], [point.y for point in points], i)])
|
|
|
|
print(y_lagrange)
|
|
print(y_newton)
|
|
|
|
# extrapolation
|
|
plt.plot(x, y_lagrange)
|
|
plt.title('fonction f(x) ex 1')
|
|
plt.show()
|
|
|
|
|
|
def exercice2() -> None:
|
|
x_i = [0, np.pi / 4, np.pi / 2, (3 * np.pi) / 4, np.pi]
|
|
|
|
x = [i for i in range(len(x_i))]
|
|
|
|
f = lambda input_of_f : np.sin(input_of_f)
|
|
|
|
results = []
|
|
real_results = []
|
|
|
|
for i in x:
|
|
res = P_newton(x_i, [f(x_i[i]) for i in range(len(x_i))], i)
|
|
print(res)
|
|
real_results.append(f(i))
|
|
results.append(res)
|
|
|
|
print(real_results)
|
|
|
|
plt.plot(results)
|
|
plt.plot(real_results, '-.')
|
|
plt.show()
|
|
|
|
|
|
def moindre_carre(points: list[Point]) -> np.ndarray[np.float64]:
|
|
"""
|
|
|
|
:param points:
|
|
:return: Le vecteur contenant les coefficients du polynome
|
|
"""
|
|
|
|
# Construction de la matrice A et du vecteur y
|
|
T = np.array([p.x for p in points])
|
|
Y = np.array([p.y for p in points])
|
|
|
|
A = np.column_stack((T**2, T, np.ones(len(T))))
|
|
|
|
# Équations normales
|
|
X = np.linalg.inv(A.T @ A) @ A.T @ Y
|
|
|
|
return X
|
|
|
|
|
|
def exercice3() -> None:
|
|
t = [1.5, 3.5, 4.5, 11, 12, 13]
|
|
h = [80, 160, 184, 184, 149, 116]
|
|
|
|
points = [Point(t[i], h[i]) for i in range(len(t))]
|
|
|
|
# fonction de calcul de h
|
|
f_h = lambda a, b, c, t: a * (t ** 2) + b * t + c
|
|
|
|
a, b, c = moindre_carre(points)
|
|
print(a,b,c)
|
|
|
|
results = []
|
|
|
|
for i in range(1, 7):
|
|
results.append(f_h(a,b,c,t[i-1]))
|
|
|
|
plt.plot(results)
|
|
plt.title('trajectoire du missile')
|
|
plt.show()
|
|
|
|
|
|
if __name__ == '__main__':
|
|
#exercice1()
|
|
#exercice2()
|
|
exercice3()
|