{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# The double pendulum problem\n\nThis animation illustrates the double pendulum problem.\n\nDouble pendulum formula translated from the C code at\nhttp://www.physics.usyd.edu.au/~wheat/dpend_html/solve_dpend.c\n\nOutput generated via `matplotlib.animation.Animation.to_jshtml`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nfrom numpy import cos, sin\n\nimport matplotlib.animation as animation\n\nG = 9.8  # acceleration due to gravity, in m/s^2\nL1 = 1.0  # length of pendulum 1 in m\nL2 = 1.0  # length of pendulum 2 in m\nL = L1 + L2  # maximal length of the combined pendulum\nM1 = 1.0  # mass of pendulum 1 in kg\nM2 = 1.0  # mass of pendulum 2 in kg\nt_stop = 2.5  # how many seconds to simulate\nhistory_len = 500  # how many trajectory points to display\n\n\ndef derivs(t, state):\n    dydx = np.zeros_like(state)\n\n    dydx[0] = state[1]\n\n    delta = state[2] - state[0]\n    den1 = (M1+M2) * L1 - M2 * L1 * cos(delta) * cos(delta)\n    dydx[1] = ((M2 * L1 * state[1] * state[1] * sin(delta) * cos(delta)\n                + M2 * G * sin(state[2]) * cos(delta)\n                + M2 * L2 * state[3] * state[3] * sin(delta)\n                - (M1+M2) * G * sin(state[0]))\n               / den1)\n\n    dydx[2] = state[3]\n\n    den2 = (L2/L1) * den1\n    dydx[3] = ((- M2 * L2 * state[3] * state[3] * sin(delta) * cos(delta)\n                + (M1+M2) * G * sin(state[0]) * cos(delta)\n                - (M1+M2) * L1 * state[1] * state[1] * sin(delta)\n                - (M1+M2) * G * sin(state[2]))\n               / den2)\n\n    return dydx\n\n# create a time array from 0..t_stop sampled at 0.02 second steps\ndt = 0.01\nt = np.arange(0, t_stop, dt)\n\n# th1 and th2 are the initial angles (degrees)\n# w10 and w20 are the initial angular velocities (degrees per second)\nth1 = 120.0\nw1 = 0.0\nth2 = -10.0\nw2 = 0.0\n\n# initial state\nstate = np.radians([th1, w1, th2, w2])\n\n# integrate the ODE using Euler's method\ny = np.empty((len(t), 4))\ny[0] = state\nfor i in range(1, len(t)):\n    y[i] = y[i - 1] + derivs(t[i - 1], y[i - 1]) * dt\n\n# A more accurate estimate could be obtained e.g. using scipy:\n#\n#   y = scipy.integrate.solve_ivp(derivs, t[[0, -1]], state, t_eval=t).y.T\n\nx1 = L1*sin(y[:, 0])\ny1 = -L1*cos(y[:, 0])\n\nx2 = L2*sin(y[:, 2]) + x1\ny2 = -L2*cos(y[:, 2]) + y1\n\nfig = plt.figure(figsize=(5, 4))\nax = fig.add_subplot(autoscale_on=False, xlim=(-L, L), ylim=(-L, 1.))\nax.set_aspect('equal')\nax.grid()\n\nline, = ax.plot([], [], 'o-', lw=2)\ntrace, = ax.plot([], [], '.-', lw=1, ms=2)\ntime_template = 'time = %.1fs'\ntime_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)\n\n\ndef animate(i):\n    thisx = [0, x1[i], x2[i]]\n    thisy = [0, y1[i], y2[i]]\n\n    history_x = x2[:i]\n    history_y = y2[:i]\n\n    line.set_data(thisx, thisy)\n    trace.set_data(history_x, history_y)\n    time_text.set_text(time_template % (i*dt))\n    return line, trace, time_text\n\n\nani = animation.FuncAnimation(\n    fig, animate, len(y), interval=dt*1000, blit=True)\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.13.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}