单摆问题可表示为

$$\ddot{\theta }+\frac{g}{L}\sin \theta =0$$

其中，$g$是重力加速度；$L$是摆长；$\theta$是摆动角度，这个方程是没有解析解的。

当角度比较小时，有$\sin (\theta )\approx \theta$，则钟摆问题可简化为线性形式

$$\ddot{\theta }+\frac{g}{L}\theta =0$$

当角度较大时，可通过数值求解。首先，将其拆分成一阶的微分方程组

$$\begin{array}{rl}& \dot{\theta }=\omega \\ & \dot{\omega }=-\frac{g}{L}\sin \theta \end{array}$$

using OrdinaryDiffEq, Plots
plotlyjs()

const g = 9.81  #g是常量
L = 1.0

u₀ = [0,π/2]
tspan = (0.0,6.0)

#定义问题
function simplependulum(du,u,t)
    θ = u[1]
    ω = u[2]        #\omega + tab = ω
    du[1] = ω
    du[2] = -(g/L)*sin(θ)
end

prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob,Tsit5())

# plotlyjs绘图面板中有个相机的图标，点击可保存图片
plot!(sol, title ="钟摆问题", xaxis = "时间", label = ["θ" "ω"])

如图所示

其中Tsit5是一阶非线性微分方程组的求解器，之前在求解衰减曲线时曾经用过，但并未做出说明。Tsit5是SciML优先推荐的求解器，是5阶Tsitouras方法，是一种经过改进的龙格库塔法。

在物理学中，除了建立物理量随时间的变化关系之外，有时还需要建立物理量之间的关系，尤其是速度与坐标的关系，二者构成了系统可能存在的所有状态，俗称相空间。

在三维空间中，速度和动量分别有三个方向，所以三维空间中的相空间一般有6个维度。但对于摆动问题来说，可以将速度和位置简化至一维，所以相空间有两个维度。

由于我们的绘图空间也有两个维度，所以只能画线，也就是说，只能从相空间中抽出一些线来绘制。

p = plot(sol,vars = (1,2), xlims = (-15,15), title = "相空间", xaxis = "速度", yaxis = "位置", leg=false)
# 用于绘图的函数
function phase_plot(func, u0, p)
    prob = ODEProblem(func,u0,(0.0,2pi))
    sol = solve(prob,Vern9()) # Vern9精度更高
    plot!(p,sol,vars = (1,2))
end
for i in -4pi:pi/2:4π
    for j in -4pi:pi/2:4π
        phase_plot(simplependulum, [j,i], p)
    end
end
plot(p,xlims = (-15,15))

结果如下