Con lắc giảm chấn Python

Với mục đích của bài viết này, con lắc điều hòa là con lắc có điểm neo di chuyển theo thời gian $x_0[t] = A\cos\omega t$. Cũng như các bài viết trước, vị trí của con lắc dao động theo thời gian có thể được mô tả bằng cơ học Lagrange. Trong một hệ tọa độ với con lắc neo ban đầu ở $[0,0]$ và trục $y$ hướng lên trên, các thành phần của vị trí và vận tốc của quả lắc là một hàm của thời gian là

\begin{align*} x &= l\sin\theta + A\cos\omega t& \dot{x} &= l\dot{\theta}\cos\theta - A\omega\sin\omega t\\ y &= -l\cos\theta & \dot{y} &= l\dot{\theta}\sin\theta \end{align*}

Động năng, $T = \tfrac{1}{2}m[\dot{x}^2 + \dot{y}^2]$, và thế năng, $V = mgy$ rồi cho Lagrangian, $ \mathcal{L} = T - V$ như

\begin{align*} \mathcal{L} = T - V = \tfrac{1}{2}m[l^2\dot{\theta}^2 + A^2\omega^2\sin^2\ omega t - 2A\omega l \dot{\theta}\sin\omega t\cos\theta] + mgl\cos\theta, \end{align*}

\begin{align*} \ddot{\theta} = \frac{A\omega^2}{l}\cos\omega t - \frac{g}{l}\sin\theta. \end{align*}

Phương trình này không thể được giải một cách giải tích, nhưng trong giới hạn $\theta$ nhỏ, nó quy về phương trình vi phân không thuần nhất cấp hai

$$ \theta_\mathrm{approx} = \frac{A\omega^2}{l[\omega_0^2-\omega^2]}[\cos\omega t - \cos\omega_0 t], $$

đối với các điều kiện ban đầu đơn giản nhất, $\theta[0] = \dot{\theta}[0] = 0$, trong đó $\omega_0 = \sqrt{g/l}$ là tần số "tự nhiên" của con lắc

Đây là biểu đồ của cả $\theta[t]$ và $\theta_\mathrm{approx}[t]$ cho biên độ dao động nhỏ, $A$. Do cộng hưởng [khi $\omega \approx \omega_0$], ngay cả một $A$ nhỏ cũng không đảm bảo giữ $\theta$ đủ nhỏ để công thức gần đúng có thể giữ vô thời hạn

Con lắc đơn giản là một ví dụ về hệ dao động cổ điển. Chuyển động điều hòa cổ điển và tương tự lượng tử của nó đại diện cho một trong những mô hình vật lý cơ bản nhất. Mô hình dao động điều hòa có thể được sử dụng để mô tả và mô hình hóa các hiện tượng như nhiệt, liên kết phân tử và tinh thể, dao động mạng tinh thể, sóng điện từ, quang phổ dao động, sóng nước, giảm xóc, sóng âm, âm học, động đất, v.v.

Trong bài viết này, chúng tôi mô tả 3 phương pháp cơ bản có thể được sử dụng để giải phương trình ODE [phương trình vi phân thường] bậc hai cho một hệ dao động điều hòa đơn giản. Sau đó, chúng tôi triển khai 3 phương thức cơ bản bằng bộ giải python

Ở đây, chúng tôi nhận xét rằng 3 phương pháp được mô tả ở trên có thể được mở rộng để bao gồm các lực bên ngoài khác như lực giảm chấn hoặc lực ma sát

Nhập thư viện cần thiết

import numpy as np
import matplotlib.pyplot as plt

Bộ giải ODE

class ODESolver[object]:
    """Second-order ODE Solver.
    Parameters
    ------------
    omega_0 : float
            initial angular velocity
    theta_0 : float
            initial angular displacement
    eta : float
        time step size
    n_iter : int
           number of steps
        
    Attributes
    -----------
    time_ : 1d-array
        Stores time values for each time step.
    omega_ : 1d-array
        Stores angular velocity values for each time step.
    theta_ : 1d-arra
       Stores angular displacement values for each time step.
        
    Methods
    -----------
    euler[alpha]: Implements the Euler algorithm for the acceleration function alpha.
    
    midpoint[alpha]: Implements the Midpoint algorithm for the acceleration function alpha.
    
    verlet[alpha]: Implements the Verlet algorithm for the acceleration function alpha.
    """def __init__[self, omega_0 = 0, theta_0 = 10, eta=0.01, n_iter=10]:
        self.omega_0 = omega_0
        self.theta_0 = theta_0
        self.eta = eta
        self.n_iter = n_iter
        
    def euler[self,alpha]:
        """Implements Euler Method.
        
        Parameters
        ----------
        alpha : acceleration function
        
        Returns
        -------
        self : object
        """
        self.time_ = np.zeros[self.n_iter]
        self.omega_ = np.zeros[self.n_iter]
        self.theta_ = np.zeros[self.n_iter]
        self.omega_[0] = self.omega_0
        self.theta_[0] = self.theta_0*np.pi/180.0
        
        for i in range[self.n_iter-1]:
            self.time_[i+1] = self.time_[i] + self.eta
            self.omega_[i+1] = self.omega_[i] + self.eta*alpha[self.theta_[i]]
            self.theta_[i+1] = self.theta_[i] + self.eta*self.omega_[i]
        return self
    
    def midpoint[self,alpha]:
        """Implement Midpoint Method.
        
        Parameters
        ----------
        alpha : acceleration function
        Returns
        -------
        self : object
        """
        self.time_ = np.zeros[self.n_iter]
        self.omega_ = np.zeros[self.n_iter]
        self.theta_ = np.zeros[self.n_iter]
        self.omega_[0] = self.omega_0
        self.theta_[0] = self.theta_0*np.pi/180.0
        
        for i in range[self.n_iter-1]:
            self.time_[i+1] = self.time_[i] + self.eta
            self.omega_[i+1] = self.omega_[i] + self.eta*alpha[self.theta_[i]]
            self.theta_[i+1] = self.theta_[i] + 0.5*self.eta*[self.omega_[i]+self.omega_[i+1]]
        return self
    
    def verlet[self,alpha]:
        """Implement Verlet Method.
        
        Parameters
        ----------
        alpha : acceleration function
        Returns
        -------
        self : object
        """
        self.time_ = np.zeros[self.n_iter]
        self.theta_ = np.zeros[self.n_iter]
        self.theta_[0] = self.theta_0*np.pi/180.0
        self.time_[1]= self.eta
        self.theta_[1] = self.theta_[0]+self.omega_0*self.eta +0.5* [self.eta**2]*alpha[self.theta_[0]]
        
        for i in range[self.n_iter-2]:
            self.time_[i+2] = self.time_[i+1] + self.eta
            self.theta_[i+2] = 2.0*self.theta_[i+1] -self.theta_[i] + [self.eta**2]*alpha[self.theta_[i+1]]
        return self

Định nghĩa hàm gia tốc góc

def alpha[x]:
    return -np.sin[x]

ví dụ 1. Phương pháp Euler

time=ODESolver[omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300].euler[alpha].time_
theta=ODESolver[omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300].euler[alpha].theta_
plt.plot[time,theta*180/np.pi,lw=3,color='red']
plt.xlabel['time[s]',size=13]
plt.ylabel['angle [deg]',size=13]
plt.title['Euler Method',size=13]
plt.show[]

Chúng tôi quan sát thấy rằng với một bước thời gian là 0. 1, phương pháp Euler cho nghiệm không ổn định. Vấn đề này có thể được giải quyết bằng cách giảm bước thời gian xuống một giá trị nhỏ hơn, ví dụ 0. 001

ví dụ 2. Phương pháp trung điểm

time=ODESolver[omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300].midpoint[alpha].time_
theta=ODESolver[omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300].midpoint[alpha].theta_
plt.plot[time,theta*180/np.pi,lw=3,color='green']
plt.xlabel['time[s]',size=13]
plt.ylabel['angle [deg]',size=13]
plt.title['Midpoint Method',size=13]
plt.show[]

Chúng tôi quan sát thấy rằng với một bước thời gian là 0. 1, phương pháp Midpoint đưa ra một giải pháp không ổn định, nhưng tương đối tốt hơn khi so sánh với phương pháp Euler. Vấn đề này có thể được giải quyết bằng cách giảm bước thời gian xuống một giá trị nhỏ hơn, ví dụ 0. 001

ví dụ 3. Phương pháp Verlet

time=ODESolver[omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300].verlet[alpha].time_
theta=ODESolver[omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300].verlet[alpha].theta_
plt.plot[time,theta*180/np.pi,lw=3,color='blue']
plt.xlabel['time[s]',size=13]
plt.ylabel['angle [deg]',size=13]
plt.title['Verlet Method',size=13]
plt.show[]

Chúng tôi quan sát thấy rằng với một bước thời gian là 0. 1, phương pháp Verlet đưa ra một giải pháp hợp lý ổn định. Bởi vì phương pháp Verlet dựa trên đạo hàm trung tâm trong khi Euler và Midpoint sử dụng đạo hàm thuận, sai số trong phương pháp Verlet là khá nhỏ

Tóm lại, chúng tôi đã chỉ ra cách xây dựng một đối tượng python để triển khai 3 phương pháp cơ bản để giải ODE bậc hai. Chúng tôi cũng cung cấp một số đầu ra mẫu từ mã được phát triển. Dựa trên phân tích này, chúng tôi quan sát thấy rằng phương pháp Verlet là phương pháp tính toán hiệu quả nhất vì nó sử dụng đạo hàm trung tâm, đây là định nghĩa đối xứng hơn của đạo hàm so với đạo hàm kỳ hạn