教程1：单摆

介绍

本教程的目的是建模和解决尚未古典但不那么简单的钟摆问题。下面给出了一个陈述（来源： Wikipedia ）

从机械的角度来看， \(M\) 质量 \(m\) 应该集中在刚性臂的下端。记录臂的长度 \(l\) .框架 \(R_0\) 假定为惯性。注意臂与垂直方向之间的角度 \(\theta\) .使用动力学进行简单建模会导致：

\[m\vec a（m/0）=\vec p+\vec t\]

在哪里？

\(\vec A(M/0)\) 质量加速度，
\(\vec P\) 如果质量的重量，
\(\vec T\) 如果手臂的反作用力。

该方程沿垂直于手臂的方向的投影给出了一个更简单的方程：

\[\ ddot\theta=-\dfrac g l \sin\theta\]

这个方程是一个二阶非线性方程。只有当方程通过假设线性化时，才能知道闭式解。 \(\theta\) 小到可以写那个 \(\sin \theta \approx \theta\) .在本教程中，我们将使用不需要这种简化的数值方法来解决这个问题。

# Setup
%matplotlib nbagg
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint

###数值

l = 1.   # m
g = 9.81 # m/s**2

第1部分：问题的重新表述

这个问题可以重新制定以符合标准的配方。 \(\dot X = f(X, t)\) ：

\[\begin{split}X=\begin{bmatrix}\theta\\\dot\theta\end{bmatrix}\end{split}\]

编写函数 \(f\) 在python中：

def f(X, t):
    """
    The derivative function
    """
    # To be completed
    return

第2部分：小角度数值解

用欧拉积分、RK4积分和欧丁积分求解，并与闭式解进行比较。首先假设摆锤是在没有速度的情况下释放的。 (\(\dot \theta = 0 ^o/s\) ）在 \(\theta = 1 ^o\) .时间离散化如下：

持续时间：10秒，
时间步长：0.01秒。

def Euler(func, X0, t):
    """
    Euler integrator.
    """
    dt = t[1] - t[0]
    nt = len(t)
    X  = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt-1):
        X[i+1] = X[i] + func(X[i], t[i]) * dt
    return X

def RK4(func, X0, t):
    """
    Runge and Kutta 4 integrator.
    """
    dt = t[1] - t[0]
    nt = len(t)
    X  = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt-1):
        k1 = func(X[i], t[i])
        k2 = func(X[i] + dt/2. * k1, t[i] + dt/2.)
        k3 = func(X[i] + dt/2. * k2, t[i] + dt/2.)
        k4 = func(X[i] + dt    * k3, t[i] + dt)
        X[i+1] = X[i] + dt / 6. * (k1 + 2. * k2 + 2. * k3 + k4)
    return X

# ODEint is preloaded.

# Define the time vector t and the initial conditions X0

第3部分：大角度解决方案

用欧拉积分、RK4积分和欧丁积分求解，并与闭式解进行比较。首先假设摆锤是在没有速度的情况下释放的。 (\(\dot \theta = 0 ^o /s\) ）在 \(\theta \in [30 ^o, 90 ^o, 179 ^o]\) .

第四部分：初角速度下的大角度解

使用Euler、RK4和ODEint积分器求解问题，并将结果与闭合形式解进行比较。首先假设摆以非零的初始角速度释放，并且 \(\dot \theta = 179 ^o/s\) 。