Lagrangian Mechanics

By TITIPATA

ชอบฟังเพลงอินดี้และอัลเทอร์เนทีฟ เวลาว่างชอบเขียนโปรแกรมและเล่นเกม

Email Twitter Google+ Github

ในโพสต์นี้เราจะมาพูดถึงฟิสิกส์กันบ้างหลังจากโพสต์ก่อนหน้าเราพูดถึงโปรแกรมมิ่งไปซะเยอะ ทุกคนที่เรียนวิชาฟิสิกส์ต้องเคยผ่านการเรียนกฎของนิวตันกันมา ไม่ว่าจะรุ่นไหน เราก็เรียนกฎของนิวตัน $F = ma$ ท่องจำกันอย่างขึ้นใจตั้งแต่สมัย ม.2 หรือ ม.3 แถมแก้สมการการเคลื่อนที่จากกฎของนิวตันได้อย่างไม่ยากเย็น เราก็เช่นกัน

แต่เราพอโตขึ้นมา ได้เรียนหลายวิชามากขึ้น เราก็คนพบว่า เห้ย! จริงๆ แล้วไม่ใช่แค่กฎของนิวตันที่ใช้แก้สมการการเคลื่อนที่ได้นะ ยังมีสมการที่เรียกว่าสมการออยเลอร์-ลากรานจ์ (Euler-Lagrange equation) หรือเรียกสั้นๆว่าสมการลากรานจ์ ที่แก้สมการการเคลื่อนที่ของวัตถุได้เช่นกัน คือจริงๆสำหรับนักเรียนทั่วไปอาจจะไม่เจอสมการนี้กันมากเท่าไหร่ แต่สำหรับคนที่อ่านเปเปอร์ในสายวิศวกรรมควบคุมหรือ Control ทีไร มักจะเจอการพูดถึงสมการนี้โดยประจำ แถมไม่ค่อยจะแก้สมการมาให้อีกด้วย ชอบทิ้งสมการไว้ใช้แก้สมการการเคลื่อนที่เอง หรือ State equation

กลับมาที่สมการลากรานจ์ มันค่อนข้างจะมหัศจรรย์เล็กน้อย แต่ว่าหลักการของสมการนี้คือเราสามารถเขียนสมการการเคลื่อนที่ได้จากเพียงพลังงานศักย์ พลังงานจลน์ และความรู้แคลคูลัส (ดิฟ ดิฟ ดิฟ) เท่านั้น แน่นอนว่าเราจะไม่ลงรายละเอียดการพิสูจน์สมการเนื่องจากต้องใช้ความรู้ Calculus of Variation ถึงจะพิสูจน์สมการลากรานจ์ได้ อันนั้นไว้สำหรับให้คนที่เรียน Optimal Control อย่างนักศึกษาภาคไฟฟ้าหรือเครื่องกลเท่านั้น

ข้อดีขอการใช้สมการลากราจน์นี้คือเราสามารถเช็คสมการการเคลื่อนที่ได้จากเพียงพลังงานจลน์และพลังงานศักย์อย่างไม่ยากนัก และเราไม่ต้องแตกแรงและแก้หลายๆสมการอย่างการใช้กฎของนิวตันอีกด้วย ยกตัวอย่างสมการการเคลื่อนที่ ที่เขียนด้วยสมการลากรานจ์ได้ง่ายๆเลยก็เช่น single pendulum, double pendulum หรือแม้แต่ triple pendulum ก็ยังได้อยู่ หรือการเขียน inverted pendulum ก็ไม่ยากเพราะว่าเราเขียนสมการของพลังงานสำหรับ inverted pendulum ง่ายกว่าใช้กฎของนิวตันแน่ๆ

โดยถ้าเราเขียนสมการของพลังงานถูกต้อง (ซึ่งง่ายกว่าการเขียนสมการนิวตันเยอะมาก) ก็จะมีโอกาสเขียนสมการการเคลื่อนที่ได้ถูกต้องมากกว่าการใช้กฎของนิวตัน และเราสามารถเอาลากรานจ์ที่เขียนไปแก้แบบ Symbolic โดยใช้ solver ในโปรแกรมต่างๆ Mathematica, sympy บน Python หรือ Matlab ก็ได้ สำหรับ Mathematica การหาสมการการเคลื่อนที่โดยใช้โค้ดเขียนแค่ 3-4 บรรทัดเท่านั้น คือจริงๆไม่สั้นขนาดนั้นนะ แต่ว่าไม่ยาวมาก ในตัวอย่างด้านล่าง เราจะใช้ Mathematica เป็นตัวอย่างในการแก้สมการ

ไม่ใช่แค่เพียงเท่านี้ ยังมีอีกหลายอย่างที่สมการลากรานจ์สามารถทำได้ เช่น การใส่โมเมนตัม การใส่ constraint เข้าไปในวิถีที่วัตถเคลื่อนที่เช่นวัตถุวิ่งอยู่บนราง หรือแม้แต่ใส่แรงเข้าไปในสมการก็ได้เช่นกัน ซึ่งเราจะไม่พูดถึงในโพสต์นี้

โอเค จะให้อธิบายสมการออยเลอร์-ลากรานจ์ ไปมากกว่านี้ก็ยังได้ แต่เนื้อหาโดยละเอียดน่าจะมีตามหนังสือวิชาการหรือแม้แต่บน Wikipedia เราจะมาลองเขียนสมการกันจริงๆเลยในโพสต์นี้ เราจะใช้ single pendulum และ double pendulum เป็นตัวอย่าง พร้อมตัวอย่างโค้ด snippet ให้ไปลองเขียนกันบน Mathematica กัน

Euler-Lagrange equation basic and terminology

ก่อนจะเริ่มต้นกันด้วยการเขียนสมการลากรานจ์ เราต้องรู้จักการเขียนตัวแปรที่เรียกกันว่าลากรานเจียน (Lagrangian) หรือใช้สัญลักษณ์ว่า $L$ กันก่อน ซึ่งคือพลังงานจลน์ของระบบลบกับพลังงานศักย์

เมื่อ $KE$ คือพลังงานจลน์ ส่วน $PE$ คือพลังงานศักย์

ส่วนสมการออยเลอร์-ลากรานจ์ (Euler-Lagrange Equation) เขียนได้ดังนี้

เมื่อ $q$ คือเวกเตอร์ของตัวแปรในระบบ สมมติเราแค่เขียนสมการของ single pendulum ตัวแปรก็คือแค่ $\theta$ นั่นเอง ถ้า double pendulum ก็จะเป็น $\theta_1, \theta_2$ เป็นต้น

เมื่อเราเขียนและแก้สมการลากรานจ์ทั้งหมดแล้ว ก็จะได้มาการการเคลื่อนที่ของวัตถุมานั่นเอง เช่น สมมติว่า $q = [x, y]^T$ ก็จะมีสมการลากรานจ์ 2 สมการ ซึ่งแก้ได้ทั้งหมด 2 differential equation

ตัวอย่างที่ 1: ก้อนมวลตกลงบนพื้น

จุดมวล $m$ ตกลงบนพื้นด้วยแรงโน้มถ่วง $g$ โดยเรากำหนดให้ $x$ เป็นความสูงของวัตถุเทียบกับพื้น ในที่นี้ตัวแปร $q$ คือ $x$ นั่นเอง เขียนพลังงานจลน์ได้ว่า $KE = \frac{1}{2}m\dot{x}^2$ ส่วนพลังงานศักย์ $PE = mgx$ ในที่ ดังนั้นเราเขียน Lagrange ได้ดังต่อไปนี้

ด้านซ้ายของสมการลากรานจ์เขียนได้ดังนี้ $\cfrac{d}{dt}\cfrac{\partial L}{\partial \dot{x}} = \cfrac{d}{dt}m\dot{x} = m\ddot{x}$

ส่วนด้านขวาของสมการคือ $\cfrac{\partial L}{\partial x} = -mg$

ดังนั้นเมื่อใส่เข้าไปในสมการลากรานจ์ จะได้ $m\ddot{x} = -mg, \ddot{x} = -g$ เป็นสมการการเคลื่อนที่ของก้อนมวลนี้ ซึ่งได้เหมือนกฎของนิวตันเป๊ะเลย เจ๋งมาก!

ตัวอย่างที่ 2: ตุ้มนาฬิกาหรือ single pendulum

สำหรับ single pendulum $q$ ในที่นี้คือ $\theta$ นั่นเอง โดยเราเขียนพลังงานจลน์ได้ว่า $KE = \frac{1}{2} mv^2 = \frac{1}{2}m(L \dot{\theta})^2$ ส่วนพลังงานศักย์เมื่อให้จุดต่ำสุดที่ pendulum ไปถึงเป็นจุดเทียบ เขียนได้ว่า $PE = mgh = mgL(1 - \cos(\theta))$ ดังนั้นเราเขียนลากรานจ์ได้ดังต่อไปนี้

ด้านซ้ายของสมการลากรานจ์ $\cfrac{d}{dt}\cfrac{\partial L}{\partial \dot{\theta}} = \cfrac{d}{dt}(mL^2\dot{\theta}) = mL^2 \ddot{\theta}$

ด้านขวาของสมการลากรานจ์เป็นประมาณนี้ (โน้ตไว้ว่า $\cfrac{d}{d\theta}\cos(\theta) = -\sin(\theta)$ อันนี้คนเรียนแคลคูลัสน่าจะจำได้ ขนาดเมาก็ยังต้องดิฟได้ ถือเป็นการเช็ค) $\cfrac{\partial L}{\partial \theta} = -mgL\sin(\theta)$

ดังนั้นสมการการเคลื่อนที่จึงเขียนได้ดังนี้ $mL^2 \ddot{\theta} = -mgL\sin(\theta)$ หรือว่า $\ddot{\theta} = -\cfrac{g}{L}sin(\theta)$

ตัวอย่างที่ 3: ตุ้มนาฬิกาสองอัน หรือ double pendulum

วิธีการก็เหมือนกันกับ single pendulum เป๊ะ แค่เราต้องเขียนพลังงานจลน์และพลังงานศักย์ของมวลอีกก้อนเพิ่ม diagram ของ double pendulum ดูได้จากข้างบน สำหรับสมการนั้นเราจะไปแก้กันโดยใช้ Mathematica เลย เราจะไม่เขียนลงไปในโพสต์ แต่เชื่อว่าผู้ที่อ่านมาถึงจุดนี้สามารถลองเช็คคำตอบของตัวเองได้กับ Mathematica โดยตรงเลย

ในที่นี้ $q = [\theta_1, \theta_2]$ นั่นเองดังนั้นหลังจากเขียน $L$ เราจะได้สมการลากรานจ์ทั้งหมดสองสมการ

และ

จากสองสมการนี้เราสามารถแก้หา $\ddot{\theta}_1, \ddot{\theta}_2$ ได้โดยใช้คำสั่ง Solve ของ Mathematica

ก่อนจะไปดูโค้ดเพื่อแก้สมการ เราลองมาเขียนกันลากรานจ์ดูก่อน จะได้อ่านโค้ดเข้าใจมากขึ้น

โดยที่

และ

ส่วนพลังงานศักย์ เขียนได้ดังนี้

จากนั้นเมื่อเราแทนทุกอย่างลงไปในสมการจะได้ลากรานจ์ $L$ ตามในโค้ดที่ให้ด้านล่าง

(*ตัวแปรสำหรับสมการลากรานจ์*)
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]

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

จริงๆส่วนของการแก้สมการจบแค่ช่วงแรกเท่านั้น แต่เนื่องจากเราต้องการเอาค่าของ $\theta_1, \theta_2$ จากเวลา 0 ถึง 50 มาพล้อตเป็นแบบแอนิเมชั่นด้วย จึงต้องใช้ NDSolve ของ Mathematica เพื่อแก้สมการ

สรุป

มาถึงจุดนี้คนเขียนโพสต์ก็ได้สลบลงหน้าคอมพิวเตอร์ (พูดเล่นนะ)

ถ้าใครลองเอาโค้ดไปรันบน Mathematica ก็ลองเปลี่ยนพารามิเตอร์ได้ตามชอบเลยตั้งแต่ขนาดของมวล รวมถึง gravitational constant เป็นต้น เพื่อดูว่าส่งผลต่อ double pendulum อย่างไรบ้าง ใครที่สนใจลองเขียนสมการอื่นๆ เราแนะนำให้ลองเขียนลากรานจ์ของระบบที่เราสนใจ อาจจะเป็น inverted pendulum ก็ได้ ($x(t), \theta(t)$) แล้วก็ลองใส่และแก้สมการใน Mathematica ดูว่าสมการที่เราเขียนนั้นถูกหรือไม่

สำหรับคนที่ชอบฟิสิกส์ก็สามารถเอาเทคนิคนี้ไปเขียนสมการการเคลื่อนที่ในระบบที่แปลกๆได้อย่างไม่ยากนัก ส่วนคนที่ชอบ simulation ก็สามารถลองเขียน simulation ฟิสิกส์ง่ายๆได้ผ่าน Mathematica และสำหรับคนที่เรียน Control โดยตรงอย่างเช่นนักศึกษาภาคไฟฟ้าหรือภาคเครื่องกล คราวนี้ก็ไม่ต้องกลัวเวลาเจอเปเปอร์ Optimal Control แล้วมีสมการหน้าตาประมาณ​ $\cfrac{d}{dt}\cfrac{\partial L}{\partial \dot{q}} = \cfrac{\partial L}{\partial q}$ นี้แล้ว :)