From d2f533dce9049236458172b1afbf3f7b9164dfdd Mon Sep 17 00:00:00 2001 From: Didictateur Date: Sun, 15 Feb 2026 17:15:25 +0100 Subject: [PATCH] working example --- example.py | 52 ++++++++++++++++++++--- legrandchien.py | 111 ++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 153 insertions(+), 10 deletions(-) diff --git a/example.py b/example.py index c276946..f013678 100644 --- a/example.py +++ b/example.py @@ -1,21 +1,59 @@ import legrandchien as lgc +import matplotlib.pyplot as plt +import matplotlib.animation as animation if __name__ == "__main__": - g = lgc.Const(-10) + g = lgc.Const(10) m = lgc.Const(1) x = lgc.Var("x") + y = lgc.Var("y") - T = 0.5 * m * x.diff() * x.diff() - V = m * g * x + T = 0.5 * m * (x.diff() * x.diff() + y.diff() * y.diff()) + V = m * g * y L = T - V - dL = L.diff() - dico = { "x" : 0, "d_x" : 1, - "d_d_x" : 1, + "y" : 0, + "d_y" : 1 } - print(dL.evaluate(dico)) \ No newline at end of file + 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 202792e..9ee6a6f 100644 --- a/legrandchien.py +++ b/legrandchien.py @@ -186,10 +186,20 @@ class Var: return Var(self.name, self.degree) def diff(self, n=1): - return Var("d_"+self.name, self.degree+n) + return Var("d_"*n+self.name, self.degree+n) def evaluate(self, dico): return dico[self.name] + + def __eq__(self, other): + if isinstance(other, str): + return False + if not (self.type == "Var" and other.type == "Var"): + return False + return self.name == other.name and self.degree == other.degree + + def __hash__(self): + return hash((self.name, self.degree)) def __add__(self, other): if other.type == "Const": @@ -350,9 +360,9 @@ class Equation: def copy(self): if self.op == Operation.CONST: - return Const(self.value.value) + return Equation(Operation.CONST, None, None, Const(self.value.value)) if self.op == Operation.VAR: - return Var(self.value.name, self.value.degree) + return Equation(Operation.VAR, None, None, Var(self.value.name, self.value.degree)) return Equation(self.op, self.left.copy(), self.right.copy(), None) def diff(self): @@ -402,6 +412,101 @@ class Equation: if self.op == Operation.DIV: return self.left.evaluate(dico) / self.right.evaluate(dico) + def solve(self, dico, tmax=10, dt=0.01): + variables = self.getAllVar() + + equations = [] + 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: + unknown = unknown.union(eq.getUnknown(dico)) + + t = 0 + values = {t : {var.name : dico[var.name] for var in variables}} + + for var in unknown: + dico[var.name] = 0. + + assert(len(unknown) == len(equations)) + n = len(unknown) + + variables = list(variables) + unknown = list(unknown) + equations = list(equations) + + while t <= tmax: + err = 1 + while err > 10**-6: + # jacobian + J = np.zeros((n, n)) + F = np.zeros(n) + for i in range(n): + F[i] = equations[i].evaluate(dico) + for j in range(n): + # print(equations[0]) + J[i][j] = equations[i].partial(unknown[j].name).evaluate(dico) + + # print(J, F) + + err = np.linalg.norm(F) + dX = -np.linalg.inv(J).dot(F) + + for i in range(n): + dico[unknown[i].name] += dX[i] + + values[t] = {var.name : dico[var.name] for var in variables} + t += dt + + # update all values + for i in range(n): + var = unknown[i] + + d = var.degree + name = var.name + while d > 0: + ivar_name = name[2:] + dico[ivar_name] += dt * dico[name] + + d -= 1 + name = ivar_name + + return values + + def getAllVar(self): + if self.op == Operation.CONST: + return set() + if self.op == Operation.VAR: + if self.value.degree == 0: + return {self.value.copy()} + else: + return {Var(self.value.name[2*self.value.degree:], 0)} + return self.left.getAllVar().union(self.right.getAllVar()) + + def getUnknown(self, dico): + if self.op == Operation.CONST: + return set() + if self.op == Operation.VAR: + if self.value.name not in dico: + return {self.value} + return set() + return self.left.getUnknown(dico).union(self.right.getUnknown(dico)) + + def __str__(self): + if self.op == Operation.CONST: + return str(self.value.value) + if self.op == Operation.VAR: + return self.value.name + if self.op == Operation.ADD: + return f"{str(self.left)} + {str(self.right)}" + if self.op == Operation.SUB: + return f"{str(self.left)} - {str(self.right)}" + if self.op == Operation.MULT: + return f"({str(self.left)})*({str(self.right)})" + if self.op == Operation.DIV: + return f"({str(self.left)})/({str(self.right)})" def __add__(self, other): if other.type == "Const":