Con lắc phi tuyến trăn

Một trong những điều thực sự khiến vật lý học "nhấp chuột" đối với tôi là học cách mô phỏng số các hệ thống mà tôi đã học trên lớp. Viết ra các phương trình chuyển động cho một bài toán nào đó là một chuyện - thực sự hình dung chúng trông như thế nào trên máy tính lại là một chuyện khác. Đó là mối liên kết còn thiếu giữa trừu tượng và thực tế

Trong bài đăng này, tôi sẽ giải bài toán con lắc đơn - một bài toán đồ chơi phổ biến trong vật lý - và chỉ ra cách mô phỏng nó trên máy tính. Để hiểu các khái niệm trong bài viết này, bạn cần có kiến ​​thức cơ bản về giải tích và làm quen với logic lập trình

Bước đầu tiên để mô phỏng bất kỳ loại vấn đề vật lý nào là viết ra các phương trình giải tích của chuyển động. Chúng ta sẽ sử dụng phương pháp Lagrangian, một cách đặc biệt hiệu quả để phân tích các hệ thống phức tạp. Sử dụng nó cho trường hợp của một con lắc đơn giản có lẽ là quá mức cần thiết, nhưng tôi thích sự sang trọng trong cách tiếp cận của nó. Bước đầu tiên của chúng tôi là xác định vấn đề của chúng tôi. Một sơ đồ của một con lắc đơn giản được trình bày dưới đây

Tuyên bố vấn đề là sau đó. Tìm phương trình chuyển động của con lắc có chiều dài l và khối lượng m. Cách tiếp cận Lagrangian để làm điều này như sau

1. Xác định tọa độ tổng quát của hệ thống,
   
   .
2. Viết Lagrange,
   
   trong đó
   
   là động năng và
   
   là thế năng, tính theo tọa độ tổng quát .
3. Các phương trình chuyển động của hệ khi đó sẽ được cho bởi phương trình Euler-Lagrange,
   
   .

Sẽ luôn có nhiều tọa độ tổng quát tương ứng với số bậc tự do trong hệ thống của bạn. Cách tôi muốn nghĩ về tọa độ tổng quát là tưởng tượng tôi đang nói chuyện điện thoại với ai đó và tôi cần mô tả trạng thái hệ thống của mình cho họ. Tập hợp thông tin tối thiểu mà tôi có thể nói với họ để họ có thể sao chép hoàn toàn hệ thống của tôi là gì? . Một nhỏ hơn hai, và do đó tọa độ tổng quát cho con lắc đơn của chúng ta đơn giản là

and positions of the mass, of course, but I could also just tell them the angle the pendulum makes with the vertical. One is fewer than two, and so the generalized coordinate for our simple pendulum is simply , tương ứng với một bậc tự do.

Bước tiếp theo là viết năng lượng của chúng ta theo góc . Động năng được viết dễ dàng bằng cách ghi nhớ hệ thức

đối với chuyển động tròn đều và lưu ý rằng trong trường hợp của chúng ta, bán kính chỉ đơn giản là chiều dài của con lắc.

Đối với thế năng, chúng ta định nghĩa "không" của mình là khi khối lượng hoàn toàn thẳng đứng, i. e. khi

. Khi đó "độ cao" phía trên điểm 0 sẽ chỉ là , do đó thế năng hấp dẫn là.

Bây giờ chúng ta có thể viết Lagrangian

và áp dụng phương trình Euler-Lagrange để có được phương trình chuyển động của chúng ta

Bây giờ chúng ta có phương trình chuyển động giải tích cuối cùng. Lưu ý rằng các số hạng khối lượng triệt tiêu, cho thấy chuyển động của con lắc không phụ thuộc vào khối lượng của nó. Ở giai đoạn này, nhiều khóa học vật lý nhập môn sẽ lấy xấp xỉ góc nhỏ

để thu được phương trình chuyển động điều hòa đơn giản, có thể giải bằng phương pháp giải tích. Không có phép tính gần đúng này, sẽ không có nghiệm giải tích cho con lắc đơn, nhưng điều đó sẽ không làm phiền chúng ta ở đây, vì chúng ta đang tìm cách giải nó bằng số. Ta có phương trình vi phân thường cấp hai, có thể viết thành hai phương trình vi phân thường cấp một.

Để giải số một hệ, bước đầu tiên là rời rạc hóa nó. Cách thông thường để làm điều này là viết ra chuỗi Taylor cho một hàm liên tục và cắt ngắn nó ở một số hạng. Trong trường hợp của chúng tôi, chúng tôi muốn tìm một giá trị gần đúng của

với điều kiện ban đầu đã biết . Mở rộng về , mang lại cho đơn hàng đầu tiên.

Viết lại, ta có

và cắm cái này vào hệ phương trình vi phân của chúng tôi, chúng tôi nhận được

Tuy nhiên, có một câu hỏi – đối với vế phải, chúng ta nên tính giá trị

và tại thời điểm nào? . Di chuyển tất cả các "đã biết" của chúng tôi sang phía bên tay phải, chúng tôi có. , else we end up with multiple unknowns to solve for. Moving all our “knowns” to the right-hand side, we have:

Bây giờ chúng ta có một phương pháp để xác định và

với các điều kiện ban đầu và where . Sử dụng và cho và và và cho và , chúng ta có thể dịch điều này sang logic lập trình như sau.

omega = omega_old - (g/l)*sin(theta_old)*tau
theta = theta_old + omega_old*tau

trong đó chúng tôi đã xác định

là một số không đổi. Đây là phương trình cập nhật của chúng tôi, mà chúng tôi sẽ lặp lại mỗi khi chúng tôi muốn "tiến tới" kịp thời để có được giá trị trong tương lai cho . Điều duy nhất chúng ta cần xác định trước, ngoài các hằng số, là các điều kiện ban đầu. Mã có liên quan, được viết bằng Python, trông như thế này.

import numpy as np
from math import *

def simplePendulumSimulation(theta0,omega0,tau,m,g,l,
                             numSteps,plotFlag):

    # initialize vectors

    time_vec = [0]*numSteps
    theta_vec = [0]*numSteps
    omega_vec = [0]*numSteps

    # set initial conditions

    theta = theta0
    omega = omega0
    time = 0

    # begin time-stepping

    for i in range(0,numSteps):
        omega_old = omega
        theta_old = theta
        # update the values
        omega = omega_old - (g/l)*sin(theta_old)*tau
        theta = theta_old + omega_old*tau
        # record the values
        time_vec[i] = tau*i
        theta_vec[i] = theta
        omega_vec[i] = omega

Chạy mã này cho các điều kiện ban đầu

và , chúng tôi ngay lập tức thấy có vấn đề.

Tổng năng lượng đang tăng lên không giới hạn. Nói cách khác, năng lượng của hệ không được bảo toàn. Bạn có thể thấy điều này bằng sự dịch chuyển ngày càng tăng của con lắc của chúng ta, nó dao động ngày càng cao một cách nghịch lý với mỗi vòng quay. Vấn đề nằm ở phép tính xấp xỉ bậc nhất thô sơ mà chúng tôi đã thực hiện trước đó khi rời rạc hóa giải pháp của mình. Rốt cuộc, phương trình chuyển động của chúng ta là phi tuyến tính, làm cho các xấp xỉ tuyến tính đặc biệt không hiệu quả trong việc giải nó. Thay vì sử dụng các phương pháp bậc cao, chúng ta chỉ cần áp dụng một biến đổi đơn giản gọi là phương pháp Euler-Cromer, đảm bảo tiết kiệm năng lượng. Chúng tôi chỉ cần sử dụng giá trị cập nhật của vận tốc góc khi nó có sẵn

omega = omega_old - (g/l)*sin(theta_old)*tau
theta = theta_old + omega*tau

Bạn có thể thấy rằng ở dòng thứ hai đã được thay đổi đơn giản thành . Chạy mô phỏng lại ta được kết quả.

Lúc này con lắc trở về đúng độ dời ban đầu sau mỗi chu kỳ. Chúng ta thấy rằng mặc dù tổng năng lượng vẫn dao động một chút (do phép tính gần đúng bậc nhất của chúng ta), nhưng hiện tại nó đã được giới hạn đúng. Để giảm dao động năng lượng, chúng ta chỉ cần giảm bước thời gian. Tập lệnh Python để thực hiện mô phỏng bằng phương pháp Euler-Cromer và tạo ảnh được đưa ra tại đây. quả lắc. py