working example
This commit is contained in:
Didictateur
parent
431570ce73
commit
d2f533dce9
2 changed files with 153 additions and 10 deletions
52
example.py
52
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))
|
||||
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()
|
||||
111
legrandchien.py
111
legrandchien.py
|
|
@ -186,11 +186,21 @@ 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":
|
||||
return Equation(
|
||||
|
|
@ -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":
|
||||
|
|
|
|||
Loading…
Reference in a new issue