?DoublePendulum类的equations()用于计算各个未知函数的导数，输入参数w数组中的变量依次为th1、th2、v1和v2，它们分别表示上球角度、下球角度、上球角速度和下球角速度。

返回值为每个变量的导数dth1、dth2、dv1和dv2，它们分别表示上球角速度、下球角速度、上球角加速度和下球角加速度。其中，dth1和dth2很容易计算，它们直接等于传入的角速度变量。

为了计算dv1和dv2，需要将微分方程组变形为如下格式：

由于两个微分方程对于和来说是两个如下形式的一次方程：

因此可以求出其中的系数a、b、c、d、e、f，?然后调用linalg.solve()解出和。

在double_pendulum_odeint()中，?调用odeint()对双摆的方程组求数值解。将odeint()得到的最终的小球状态保存到pendulum.init_status中，作为下一次调用odeint()的初始值，因此多次调用double_pendulum_odeint()可以生成连续的运动轨迹。函数返回4个数组，分别是两个小球的X-Y轴的坐标。

最后是主程序部分，我们使用小角度和大角度的初始值分别计算双摆的摆动轨迹。小初始角度时小球的运动轨迹如图18-5所示(见文前彩插)。大初始角度时小球的运动轨迹如图18-6所示(见封三彩插)。可以看出，当初始角度很大时，摆动出现混沌现象。