Snippet for Lagrangian Mechanics

Douple Pendulum Snippet (Thai)

(*ตัวแปรสำหรับสมการลากรานจ์*)
q = {{theta1[t]}, {theta2[t]}};
dq = D[q, t];
ddq = D[dq, t];

(*เริ่มด้วยการเขียน Lagrangian ก่อนซึ่งคือพลังงานจลน์กับพลังงานศักย์*)
PE = -m1 g R1 Cos[theta1[t]] + m2 g ( -R1 Cos[theta1[t]] - R2 Cos[theta1[t] + theta2[t]]); (*พลังงานศักย์*)
x1[t]  = R1 Sin[theta1[t]];
y1[t] = -R1 Cos[theta1[t]];
x2[t] = R1 Sin[theta1[t]] + R2 Sin[theta1[t] + theta2[t]];
y2[t] = -R1 Cos[theta1[t]] - R2 Cos[theta1[t] + theta2[t]];
KE = (1/2)*m1*(D[x1[t], t]^2 + D[y1[t], t]^2) + (1/2)*m2*(D[x2[t], t]^2 + D[y2[t], t]^2 ); (*พลังงานจลน์*)

L = KE - PE // Simplify;

(*Euler-Lagrange Equation โดย ELtemp คือคำตอบของสมการแล้ว*)
D[D[L, Transpose[dq]], t] - D[L, Transpose[q]] == 0;
Eq = (Thread[D[D[L, Transpose[dq]], t] - D[L, Transpose[q]] == 0]);
ELtemp = Solve[Eq[[1]] && Eq[[2]], Flatten[ddq]];
EL = {theta1''[t] == ELtemp[[1, 1, 2]],
   theta2''[t] == ELtemp[[1, 2, 2]]};

(*ด้านล่างนี้เรากำหนดจุดเริ่มต้นให้กับ double pendulum และให้  Mathematica แก้สมการถึงเวลา t=50*)

(*พารามิเตอร์สำหรับ simulation*)
m1 = 1;
m2 = 1;
R1 = 1;
R2 = 1;
g = 9.8; (*gravitational constant*)

(*เวลาแก้สมการ second-order เราจำเป็นจะต้องมีจุดเริ่มต้นหรือ initial condition*)
InitCon = {theta1[0] == 0.1, theta1'[0] == 0, theta2'[0] == 0, theta2[0] == Pi};

(*แก้สมการโดยใช้ Mathematica ND Solve*)
sol = NDSolve[Join[EL, InitCon], {theta1, theta2}, {t, 0, 50}];

(*Assign Solution and Plot*)
theta1sol = theta1[t] /. sol[[1, 1]];
theta2sol = theta2[t] /. sol[[1, 2]];

(*create the function to solve position of pendulums*)
x1[theta1_] := R1 Sin[theta1];
y1[theta1_] := -R1 Cos[theta1];
x2[theta1_, theta2_] := R1 Sin[theta1] + R2 Sin[theta2];
y2[theta1_, theta2_] := -R1 Cos[theta1] - R2 Cos[theta2];

(*Animate Part!*)
X1[tau_] := R1 Sin[theta1sol] /. t -> tau;
Y1[tau_] := -R1 Cos[theta1sol] /. t -> tau;
X2[tau_] :=  (R1 Sin[theta1sol] + R2 Sin[theta1sol + theta2sol]) /. t -> tau;
Y2[tau_] := (-R1 Cos[theta1sol] - R2 Cos[theta1sol + theta2sol]) /. t -> tau;

(*มาถึงบรรทัดนี้เราจะวาดแอนิเมชั่นของ double pendulum*)
prop1 = Line[{{-4, 0}, {4, 0}}];
prop2 = Line[{{0, -4}, {0, 4}}];
prop3 = {PointSize[Large], Red, Point[{0, 0}]};
anim = Animate[Graphics[{{AbsoluteThickness[2], Blue,
      Line[{{0, 0}, {X1[tau], Y1[tau]}}]}, {AbsoluteThickness[2], Blue,
      Line[{{X1[tau], Y1[tau]}, {X2[tau], Y2[tau]}}]},
      {PointSize[Large], Red, Point[{X1[tau], Y1[tau]}]},
      {PointSize[Large], Red, Point[{X2[tau], Y2[tau]}]}, prop1, prop2, prop3}],
      {tau, 0, 50},
   AnimationRate -> 0.8]

ลองก้อปโค้ดไปใส่กันที่ Mathematica Online ที่นี่ https://lab.open.wolframcloud.com/app/ > Create a New Notebook

แล้วก็มีน้องส่งเมล์มาถามว่า ถ้าอยากจะพล้อตเส้นทางเดินของ double pendulum ต้องทำยังไง ข้างล่างนี้เป็น snippet เพื่อใช้พล้อตทางเดินของ pendulum ก้อนล่างจากเวลา 0 ถึง 10 วินาทีฮะ

ParametricPlot[{x2[theta1sol, theta2sol], y2[theta1sol, theta2sol]}, {t, 0, 10}, AxesLabel -> {x2, y2}]

titipata/lagrangian.md

Douple Pendulum Snippet (Thai)