exemple pendule à ressort

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()

example.py

@@ -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()

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())

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()