From 6abe888da4598a57962f5d5b9660823ec422e432 Mon Sep 17 00:00:00 2001 From: Didictateur Date: Mon, 16 Feb 2026 11:33:29 +0100 Subject: [PATCH] =?UTF-8?q?exemple=20pendule=20=C3=A0=20ressort?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- chute_libre.py | 58 ++++++++++++++++++ example.py | 59 ------------------ legrandchien.py | 138 ++++++++++++++++++++++++++++++++++++++++++- pendule_a_ressort.py | 56 ++++++++++++++++++ 4 files changed, 249 insertions(+), 62 deletions(-) create mode 100644 chute_libre.py delete mode 100644 example.py create mode 100644 pendule_a_ressort.py diff --git a/chute_libre.py b/chute_libre.py new file mode 100644 index 0000000..0303c25 --- /dev/null +++ b/chute_libre.py @@ -0,0 +1,58 @@ +import legrandchien as lgc +import matplotlib.pyplot as plt +import matplotlib.animation as animation + +g = lgc.Const(10) +m = lgc.Const(1) + +x = lgc.Var("x") +y = lgc.Var("y") + +T = 0.5 * m * (x.diff() * x.diff() + y.diff() * y.diff()) +V = m * g * y +L = T - V + +dico = { + "x" : 0, + "d_x" : 1, + "y" : 0, + "d_y" : 1 +} + +values = L.solve(dico, 10, 0.1) + +# Extraire les données pour x et y +times = sorted(values.keys()) +positions_x = [values[t]["x"] for t in times] +positions_y = [values[t]["y"] for t in times] + +# Créer la figure +fig, ax = plt.subplots(figsize=(8, 8)) +ax.set_xlim(min(positions_x) - 5, max(positions_x) + 5) +ax.set_ylim(min(positions_y) - 5, max(positions_y) + 5) +ax.set_xlabel("X (m)") +ax.set_ylabel("Y (m)") +ax.grid(True) + +# Point animé +point, = ax.plot([], [], 'ro', markersize=15, label="Particule") +# Trajectoire +line, = ax.plot([], [], 'b-', alpha=0.3, linewidth=2, label="Trajectoire") +# Texte temps +time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12) + +ax.legend() + +def animate(frame): + # Point courant + point.set_data([positions_x[frame]], [positions_y[frame]]) + # Trajectoire jusqu'au point courant + line.set_data(positions_x[:frame + 1], positions_y[:frame + 1]) + # Texte + time_text.set_text(f"t = {times[frame]:.2f}s\nx = {positions_x[frame]:.2f}m, y = {positions_y[frame]:.2f}m") + return point, line, time_text + +ani = animation.FuncAnimation(fig, animate, frames=len(times), + interval=50, blit=True, repeat=True) + +plt.show() \ No newline at end of file diff --git a/example.py b/example.py deleted file mode 100644 index f013678..0000000 --- a/example.py +++ /dev/null @@ -1,59 +0,0 @@ -import legrandchien as lgc -import matplotlib.pyplot as plt -import matplotlib.animation as animation - -if __name__ == "__main__": - g = lgc.Const(10) - m = lgc.Const(1) - - x = lgc.Var("x") - y = lgc.Var("y") - - T = 0.5 * m * (x.diff() * x.diff() + y.diff() * y.diff()) - V = m * g * y - L = T - V - - dico = { - "x" : 0, - "d_x" : 1, - "y" : 0, - "d_y" : 1 - } - - values = L.solve(dico, 10, 0.1) - - # Extraire les données pour x et y - times = sorted(values.keys()) - positions_x = [values[t]["x"] for t in times] - positions_y = [values[t]["y"] for t in times] - - # Créer la figure - fig, ax = plt.subplots(figsize=(8, 8)) - ax.set_xlim(min(positions_x) - 5, max(positions_x) + 5) - ax.set_ylim(min(positions_y) - 5, max(positions_y) + 5) - ax.set_xlabel("X (m)") - ax.set_ylabel("Y (m)") - ax.grid(True) - - # Point animé - point, = ax.plot([], [], 'ro', markersize=15, label="Particule") - # Trajectoire - line, = ax.plot([], [], 'b-', alpha=0.3, linewidth=2, label="Trajectoire") - # Texte temps - time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12) - - ax.legend() - - def animate(frame): - # Point courant - point.set_data([positions_x[frame]], [positions_y[frame]]) - # Trajectoire jusqu'au point courant - line.set_data(positions_x[:frame + 1], positions_y[:frame + 1]) - # Texte - time_text.set_text(f"t = {times[frame]:.2f}s\nx = {positions_x[frame]:.2f}m, y = {positions_y[frame]:.2f}m") - return point, line, time_text - - ani = animation.FuncAnimation(fig, animate, frames=len(times), - interval=50, blit=True, repeat=True) - - plt.show() \ No newline at end of file diff --git a/legrandchien.py b/legrandchien.py index 9ee6a6f..5b7544b 100644 --- a/legrandchien.py +++ b/legrandchien.py @@ -1,5 +1,6 @@ from enum import Enum import numpy as np +from tqdm import tqdm class Operation(Enum): VAR = 1 @@ -363,6 +364,8 @@ class Equation: return Equation(Operation.CONST, None, None, Const(self.value.value)) if self.op == Operation.VAR: return Equation(Operation.VAR, None, None, Var(self.value.name, self.value.degree)) + if self.op in {Operation.SIN, Operation.COS}: + return Equation(self.op, None, None, self.value.copy()) return Equation(self.op, self.left.copy(), self.right.copy(), None) def diff(self): @@ -379,6 +382,10 @@ class Equation: if self.op == Operation.DIV: return (self.left.diff()*self.right.copy() - self.left.copy()*self.right.diff()) \ / (self.right.copy()*self.right.copy()) + if self.op == Operation.SIN: + return self.value.diff() * Equation(Operation.COS, None, None, self.value.copy()) + if self.op == Operation.COS: + return Const(-1) * self.value.diff() * Equation(Operation.SIN, None, None, self.value.copy()) def partial(self, var_name): if self.op == Operation.CONST: @@ -397,6 +404,10 @@ class Equation: if self.op == Operation.DIV: return (self.left.partial(var_name)*self.right.copy() - self.left.copy()*self.right.partial(var_name)) \ / (self.right.copy()*self.right.copy()) + if self.op == Operation.SIN: + return self.value.partial(var_name) * Equation(Operation.COS, None, None, self.value.copy()) + if self.op == Operation.COS: + return Const(-1) * self.value.partial(var_name) * Equation(Operation.SIN, None, None, self.value.copy()) def evaluate(self, dico): if self.op == Operation.CONST: @@ -411,17 +422,23 @@ class Equation: return self.left.evaluate(dico) * self.right.evaluate(dico) if self.op == Operation.DIV: return self.left.evaluate(dico) / self.right.evaluate(dico) + if self.op == Operation.SIN: + return np.sin(self.value.evaluate(dico)) + if self.op == Operation.COS: + return np.cos(self.value.evaluate(dico)) def solve(self, dico, tmax=10, dt=0.01): + self.simplify() variables = self.getAllVar() equations = [] - print([v.name for v in variables]) + # print([v.name for v in variables]) for var in variables: equations.append(self.partial(var.diff().name).diff() - self.partial(var.name)) unknown = self.getUnknown(dico) for eq in equations: + eq.simplify() unknown = unknown.union(eq.getUnknown(dico)) t = 0 @@ -430,6 +447,8 @@ class Equation: for var in unknown: dico[var.name] = 0. + # print([var.name for var in unknown]) + # print(equations[0]) assert(len(unknown) == len(equations)) n = len(unknown) @@ -437,9 +456,13 @@ class Equation: unknown = list(unknown) equations = list(equations) - while t <= tmax: + for t in tqdm([n * dt for n in range(int(tmax/dt))]): err = 1 - while err > 10**-6: + iterations = 0 + max_iteration = 50 + + while err > 10**-6 and iterations < max_iteration: + iterations += 1 # jacobian J = np.zeros((n, n)) F = np.zeros(n) @@ -483,17 +506,102 @@ class Equation: return {self.value.copy()} else: return {Var(self.value.name[2*self.value.degree:], 0)} + if self.op in {Operation.SIN, Operation.COS}: + return self.value.getAllVar() return self.left.getAllVar().union(self.right.getAllVar()) def getUnknown(self, dico): + # print(self) if self.op == Operation.CONST: return set() if self.op == Operation.VAR: if self.value.name not in dico: return {self.value} return set() + if self.op in {Operation.SIN, Operation.COS}: + return self.value.getUnknown(dico) return self.left.getUnknown(dico).union(self.right.getUnknown(dico)) + def simplify(self): + if self.left != None: + self.left.simplify() + if self.right != None: + self.right.simplify() + if self.op in {Operation.SIN, Operation.COS}: + self.value.simplify() + + if self.op == Operation.CONST: + pass + if self.op == Operation.VAR: + pass + + if self.op == Operation.ADD: + if self.left.op == Operation.CONST and self.right.op == Operation.CONST: + self.op = Operation.CONST + self.value = Const(self.left.value.value + self.right.value.value) + self.left = None + self.right = None + elif self.left.op == Operation.CONST and self.left.value.value == 0: + self.op = self.right.op + self.left = self.right.left + self.value = self.right.value + self.right = self.right.right + elif self.right.op == Operation.CONST and self.right.value.value == 0: + self.op = self.left.op + self.right = self.left.right + self.value = self.left.value + self.left = self.left.left + + if self.op == Operation.SUB: + if self.left.op == Operation.CONST and self.right.op == Operation.CONST: + self.op = Operation.CONST + self.value = Const(self.left.value.value - self.right.value.value) + self.left = None + self.right = None + elif self.right.op == Operation.CONST and self.right.value.value == 0: + self = self.left + + if self.op == Operation.MULT: + if self.left.op == Operation.CONST and self.right.op == Operation.CONST: + self.op = Operation.CONST + self.value = Const(self.left.value.value * self.right.value.value) + self.left = None + self.right = None + elif self.left.op == Operation.CONST and self.left.value.value == 0: + self.op = Operation.CONST + self.left = None + self.right = None + self.value = Const(0) + elif self.right.op == Operation.CONST and self.right.value.value == 0: + self.op = Operation.CONST + self.left = None + self.right = None + self.value = Const(0) + elif self.left.op == Operation.CONST and self.left.value.value == 1: + self.op = self.right.op + self.left = self.right.left + self.value = self.right.value + self.right = self.right.right + elif self.right.op == Operation.CONST and self.right.value.value == 1: + self.op = self.left.op + self.right = self.left.right + self.value = self.left.value + self.left = self.left.left + + if self.op == Operation.SIN: + if self.value.op == Operation.CONST: + self.op = Operation.CONST + self.left = None + self.right = None + self.value = Const(np.sin(self.value.value)) + + if self.op == Operation.COS: + if self.value.op == Operation.CONST: + self.op = Operation.CONST + self.left = None + self.right = None + self.value = Const(np.cos(self.value.value)) + def __str__(self): if self.op == Operation.CONST: return str(self.value.value) @@ -507,6 +615,10 @@ class Equation: return f"({str(self.left)})*({str(self.right)})" if self.op == Operation.DIV: return f"({str(self.left)})/({str(self.right)})" + if self.op == Operation.SIN: + return f"sin({str(self.value)})" + if self.op == Operation.COS: + return f"cos({str(self.value)})" def __add__(self, other): if other.type == "Const": @@ -595,3 +707,23 @@ class Equation: def __rtruediv__(self, other): if isinstance(other, (int, float)): return Const(other) / self + +def cos(value): + if isinstance(value, (float, int)): + return Equation(Operation.CONST, None, None, np.cos(value)) + if value.type == "Const": + return Equation(Operation.CONST, None, None, np.cos(value.value)) + if value.type == "Var": + return Equation(Operation.COS, None, None, Equation(Operation.VAR, None, None, value.copy())) + if value.type == "Equation": + return Equation(Operation.COS, None, None, value.copy()) + +def sin(value): + if isinstance(value, (float, int)): + return Equation(Operation.CONST, None, None, np.sin(value)) + if value.type == "Const": + return Equation(Operation.CONST, None, None, np.sin(value.value)) + if value.type == "Var": + return Equation(Operation.SIN, None, None, Equation(Operation.VAR, None, None, value.copy())) + if value.type == "Equation": + return Equation(Operation.SIN, None, None, value.copy()) \ No newline at end of file diff --git a/pendule_a_ressort.py b/pendule_a_ressort.py new file mode 100644 index 0000000..463a85e --- /dev/null +++ b/pendule_a_ressort.py @@ -0,0 +1,56 @@ +import legrandchien as lgc +import matplotlib.pyplot as plt +import matplotlib.animation as animation +import numpy as np + +g = lgc.Const(10) +m = lgc.Const(1) +k = lgc.Const(10) +r0 = lgc.Const(1) + +theta = lgc.Var("theta") +r = lgc.Var("r") + +x = r * lgc.cos(theta) +y = r * lgc.sin(theta) + +T = 0.5 * m * (x.diff()*x.diff() + y.diff()*y.diff()) +V = m * g * r * lgc.cos(theta) + k * (r - r0) * (r - r0) + +L = T - V + +dico = { + "theta" : -3.14/2, + "d_theta" : 0, + "r" : 1.2, + "d_r" : 0, +} + +data = L.solve(dico, 10, 0.01) + +# Extraire les positions +times = sorted(data.keys()) +positions = [data[t] for t in times] + +x_positions = [p["r"] * np.sin(p["theta"]) for p in positions] +y_positions = [p["r"] * np.cos(p["theta"]) for p in positions] + +# Créer la figure et l'axe +fig, ax = plt.subplots(figsize=(8, 8)) +ax.set_xlim(-3, 3) +ax.set_ylim(-3, 3) +ax.set_xlabel("X (m)") +ax.set_ylabel("Y (m)") +ax.grid(True) + +line, = ax.plot([], [], 'o-', lw=2, markersize=8) +trace, = ax.plot([], [], 'r-', alpha=0.3, lw=1) + +def animate(frame): + line.set_data([0, x_positions[frame]], [0, y_positions[frame]]) + trace.set_data(x_positions[:frame+1], y_positions[:frame+1]) + return line, trace + +anim = animation.FuncAnimation(fig, animate, frames=len(x_positions), interval=50, blit=True, repeat=True) + +plt.show() \ No newline at end of file