{ "cells": [ { "cell_type": "markdown", "id": "1620df10-4f75-4a4b-adee-01330b7d750e", "metadata": { "tags": [] }, "source": [ "# Kapitza's Inverted Pendulum" ] }, { "cell_type": "markdown", "id": "3e8a58d2-77d2-4ebf-81ca-0b3a088d3511", "metadata": {}, "source": [ "An inverted pendulum --- one where the rod points upward so that the bob is above the pivot --- is usually unstable. However, if the pivot is driven to oscillate up and down at a high frequency, the pendulum can be stabilized in the upward position. This system, called Kapitza's pendulum, is shown in the Figure below. We will solve its equation of motion and illustrate the stability of the pendulum. \n", "![Kapitza](source/Kapitza.png)" ] }, { "cell_type": "markdown", "id": "6fd94eb6-36b6-492e-acfa-730ada348f83", "metadata": {}, "source": [ "The Lagrangian of the system is given by (ignoring terms that do not depend on $\\theta$ or $\\dot{\\theta}$):\n", "\\begin{equation*}\n", "L = \\frac{1}{2} m l^2 \\dot{\\theta}^2 - m l \\sin\\theta \\, \\dot{\\theta} \\, \\dot{z} - m g l \\cos\\theta\n", "\\end{equation*}\n", "Let us assume that the pivot is driven to move as $z(t) = A \\cos \\omega t$. The Euler-Lagrange equation for $\\theta$ leads to:\n", "\\begin{equation*}\n", "\\ddot{\\theta} = (\\omega_0^2 - \\alpha \\omega^2 \\cos \\omega t) \\sin\\theta\n", "\\end{equation*}\n", "where $\\omega_0^2 \\equiv g/l$ and $\\alpha \\equiv A/l$." ] }, { "cell_type": "markdown", "id": "594bbe16-bfd3-4a1b-a431-1fe036ef9676", "metadata": {}, "source": [ "To solve this equation numerically, we first cast it into a set of first order ordinary differential equations (ODE). Let us define an auxilliary variable $\\eta$ to be the time derivative of $\\theta$, i.e., $\\eta \\equiv \\dot{\\theta}$, then we have:\n", "\\begin{align*}\n", "\\dot{\\theta} &= \\eta \\\\\n", "\\dot{\\eta} &= (\\omega_0^2 - \\alpha \\omega^2 \\cos \\omega t) \\sin\\theta\n", "\\end{align*}\n", "We will solve these equations using the `odeint` function from the `scipy.integrate` module, as illustrated below." ] }, { "cell_type": "code", "execution_count": 1, "id": "6bcaed1e-42f3-48f5-a53f-2c52db5c14a5", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.integrate as intgr" ] }, { "cell_type": "markdown", "id": "cb44f561-534a-4f75-be2f-89b0b6a60b37", "metadata": {}, "source": [ "We first have to define a function to calculate the time derivatives of the variables. All variables are stacked into an array `X`, which is the input to the function, besides time and other parameters." ] }, { "cell_type": "code", "execution_count": 2, "id": "6f154eb4-84b4-4c23-9353-08ba62849e79", "metadata": {}, "outputs": [], "source": [ "def equations(X, t, alpha, omega): # calculate time derivatives\n", " theta, eta = X # parse variables\n", " d_theta = eta\n", " d_eta = (1 - alpha * omega**2 * np.cos(omega*t)) * np.sin(theta)\n", " return [d_theta, d_eta] # assemble derivatives into an array" ] }, { "cell_type": "markdown", "id": "4ee1e88d-1d68-41e1-a4bb-70a6dd5d7ec1", "metadata": {}, "source": [ "Now we choose some values for the parameters. We have set $\\omega_0 = 1$ by rescaling time. In order for Kapitza's pendulum to be stable, the driving frequency $\\omega$ must be large compared to $\\omega_0$ (you can test this later by changing the value of $\\omega$)." ] }, { "cell_type": "code", "execution_count": 3, "id": "ab95d350-7385-4a56-b841-96de8da5c569", "metadata": {}, "outputs": [], "source": [ "alpha = 0.1 # amplitude of pivot << 1\n", "omega = 30. # omega >> 1\n", "\n", "theta0 = 0.1 # initial position\n", "eta0 = 0.0 # initial speed\n", "\n", "T = 10. # total time to solve for\n", "time = np.arange(0, T, 0.01) # time points to evaluate solution at\n", "\n", "sol = intgr.odeint(equations, [theta0, eta0], time, args=(alpha, omega)) # solve equations\n", "theta = sol[:,0] # theta is the first component of solution" ] }, { "cell_type": "markdown", "id": "1afbcaaa-18fc-4372-955a-0204f50150a5", "metadata": {}, "source": [ "We have solved the equations using the `odeint` function. Let us plot our result." ] }, { "cell_type": "code", "execution_count": 4, "id": "b773e153-a7a5-465d-93fd-80b6607f7020", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEICAYAAACNn4koAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABMvUlEQVR4nO29d5hkZ3Xn/zlV1V3VsapzDhM1QZpRGISEkAxIAoGxhJdgsTJIXmyt7WUdWB4bHMDG5vmB7V3W68VeywrINhhhCSNhgqyACALEjKTR5NATOudQHau7q+v8/rj31tS0umd6pqvqVng/z1PPVN16761z5+37ft9w3nNEVTEYDAaDIRV43DbAYDAYDLmLERmDwWAwpAwjMgaDwWBIGUZkDAaDwZAyjMgYDAaDIWX43DYgk6iurtb29na3zTAYDIas4uWXXx5R1ZqVvjMik0B7ezv79u1z2wyDwWDIKkSkc7XvzHSZwWAwGFKGERmDwWAwpAwjMgaDwWBIGUZkDAaDwZAyjMgYDAaDIWVkhMiIyB0iclxEOkTkEyt8f4uIvCIiURF537Lv7hWRk/br3oTj14nIQfua/0dEJB33YjAYDIZzuC4yIuIFvgi8E9gBfFBEdiwr1gXcB3xl2bmVwKeBNwLXA58WkQr7678Dfg3YYr/uSNEtGAwGg2EVXBcZLHHoUNXTqroAfBW4K7GAqp5V1QNAbNm57wCeUdUxVR0HngHuEJEGoFxVf6pWLoN/BN6TqhvoGJri8989hkmbkFlMzC7wb6/2sBQz9ZJJRBaX+PbBfhaXlj/OhlwkE0SmCehO+NxjH1vPuU32+4teU0TuF5F9IrJveHh4zUYn8sLxYf7uhVN87/jQZZ1vSA1//ORhfvex1/j6Kz0XL2xIG3/z/El+88uv8KUXz7ptiiENZILIuIqqPqCqe1R1T03NilERLsqHb2wnUODhBydGkmydYT387MwoAD84aeolk3Cek++fuLxOnSG7yASR6QVaEj4328fWc26v/f5yrnnJFPo8bKsv5+TQVKp+wnCJhOcWGZycB+DkoKmXTCEWU44PWPXRMTTtsjWGdJAJIrMX2CIiG0SkELgbeGqN5z4NvF1EKuwF/7cDT6tqPzApIjfYXmUfBp5MhfEObVXFdI7OpvInDJdA38QcAPXlAbrGZs16WYYwPD3PwlKM2jI/A5MR5haW3DbJkGJcFxlVjQIfxRKMo8DXVPWwiHxGRO4EEJE3iEgP8H7g70XksH3uGPBnWEK1F/iMfQzgN4EHgQ7gFPCdVN5HW2UxfRNzLETNYmYm0DtuicybNlUxu7DE8PS8yxYZAHrserlpczUAXWOmY5bruC4yAKr6bVXdqqqbVPWz9rFPqepT9vu9qtqsqiWqWqWqOxPOfVhVN9uvRxKO71PVK+1rflRT3JVtrSohptBr96AN7uLUw42bqgDoMqPMtLO4FOPsyMx5x5wR5pvsejmz7HtD6nm1a5xf/6eXGZqKpOX3MkJkcoG2qmIAOkfPf2hm5qNumJNX/ODEMJ/55hGiCS6xvRNzFPo8XNNqbZsyPeb08/tPHOAtf/UCL3eOx48tF/+ecVMv6eb/+84xvnt4gH/6yarR+ZOKEZkk0RgqAqA/fK538LW93Vz5J0/z3NFBt8zKeVSV+x75GQ+/eIZnjpz7f+4dn6MpVETTCvUyFVnkwR+eZmJ2Ie325guxmPLN1/oA+NaB/vjx3vE5gkUFNIWKCBR4zqsXQ+pRVQ71hgH46enRtPymEZkkUVfmxyPnpgMA/mVvF6rwhNmnkTJ6J+Zw9lr++NToecebQkUUFXqpKC44r14e+MFp/vxbR/nLp4+n29y8YWR6nsUlq2IO9k7Ejzv1IiI0BovoD5+rl9PD07zjCz/gPw4PpNvcvGF8dpFZ29nixOB0WhxijMgkCZ/XQ115gL4Jq2e2uBTjYI/VYzjWb1xoU8WJBPdkxzUWrMasMRQAoCFYdF6P+YXj1v6Mn5xKT08uH+mz/79ry/ycHDrXmPWOz8VH/Q2hc88LwD//tIvjg1P8r2dOpN/gPMHpbF2/oZLw3CIj06kfzRuRSSKNoaJ4JXaOzhKNKY1By4U2akJopIQzI9ac/tt31NExbO27iCwuMTw1T1PIWidrDAXi9RJdinFsYNI6d3TGuNCmiAF7hHLzlhomZhcZnVlAVemdmKO5whaZZSOZn521RL9jaJr5qKmXVDBgi/8tWyzvvnTsVTIik0QaQ0X02Q+NU3m3bq8jGlPjdZYihiYjFPo8XL+hkrGZBUan5+Ojlia7MUsU/66xWRaXlJ/bWoMqdJuF55Tg1MEtW63G7OTgNJNzUabno/F1ssZggKGpeRaXYizFlBOD0zQEA0Rjyulh43WWCvonrXp58xYruknHsDXKfLFjhLGZ1IxqjMgkkcZggP5whFhMOWX3qt+2vRYwrpqpYmAyQl25ny11ZQCcHJqO75FxGrOGYBGTkSgz89G4+N+y1XrIjPinhv6wJf572isBOD0yHRd0R/wbQkWowtDUPD3jsyxEY7x1m/W8OHVoSC4D4Tm8HuHKxnKKCrycHp5meGqeex58iW8d6EvJbxqRSSKNoSIWojFGZxY4NWz1yrbXlwPnNqEZkstAOEJ9eYDNtaWANYJ0Ri3xHrO9NtMfnotPqd1sTxf0GZFJCf3hCA3BAA3lAYoLvZwamom7K7dWWtOYDUG7XibmODlo1cvP2eLfFzb1kgr6wxHqyvz4vB42VJdwZmSGTtu9v8Wul2RjRCaJnHNjnuPU0DSbakqpKfPj9Uh8LtSQXIam5qkrD9AYtBqzjqFpesZnEYF6uxFz6qVvIkLH0DR15X421ZTi9YgRmRQxEJ6jIRjA4xE2VJdwang6vlfJaczi9RKOxMX/hg1VFPo8ZiSTIvonIjTY/+8ba0o4PTwTD4fVVlWSkt80IpNEnJ5Zz/gcHUPTbK61GjInTpMhuagqA+EIdeUBRIRNNaWcGp7m1MgMzRVFFPqsP2+nXvomzom/1yPUl5/v3WRIHn0TERqCVmPm1Ev3mLVHJlhUAJw/kukYmqamzE+wuIDGYMBMY6aIgclIvPO1saaUnvFZOoam8ci5kX+yMSKTRJxK+tmZMWYWlthWb60T1JUHzEgmBUxGoswtLlFfbj00m2tL6Ria5tTQNJtrSuPl6soDeMQS/+ODU1xh10tjyDRmqSAWUwYnI3ER2VRTSu/EHCcGp2ipPNeQlQUKKPP76JuY4+TgFFvsKc+miiIzwkwBqkp/eI4G+3nZWG2FwvpRxzCNoXOdsmRjRCaJhIqtXtq3Dlo7nJ3GrCEYOM9V05AcBu3RYV3wnMj0hyMcG5iKOwIAFHg9tFWV8PThASKLMXY0WOtkdeUBhqdM4MxkMzI9TzSm50SmtgRVeOnMGJsSxB+sKZsj/ZMcHZhiZ+O5enHSNBiSR3hukchiLD6S2VJn1cWh3sn4mmYqMCKTRESEbfVlDE/NIwJb7YauPmgemlTgiIwzknHEA3jdQ7OrOchJ27Nse0NiYxYxaQCSjOO+XG9Pl13dEop/t6s5dF7ZnU1B9p4dZyEa48qmIGDVy9CUqZdk49SLsxZ2RUJH7JqWipT9rhGZJPOmTZbX0s7Gckr8PsBqBKfno0xFFt00LedwpiDryv2AtYvZ4ZYt52c5TWzcHPGvK/czu7DEtAlimlScxswZyTRXFHNdWwU+j/COnXXnlb2yMRh/v9N+X1fmZ3FJGZ81z0syGYiLv1UvPq+Hj751M1Ulhfyna9ea8f7S8aXsynnKh29s4+zoDB+6sS1+zKnUgXCEskCBW6blHPHpMnskU+L38cRv3Mh8wpSAw89f1cBDPzzNB97QEp97ds4bnJw39ZJEnKnhhoQ6eORX3sDs/NLr6uXmLdV4xPJs2lRjeTfVxuslQmVJYZqszn36VqiXj7/jCj7+jitS+rtGZJJMRUkhX/ilq8875njZDExGzlsrMKyPwcl5QsUFBAq88WPXtVWuWLY+GODFT7wNK1GqRW2Z9bANTUZSOiedbwzYGzETBaI8UED5CkLeUlnMf/zuLVSV+ON144xMBycj8alNw/oZCEfwCNSU+tP6uxkxXSYid4jIcRHpEJFPrPC9X0Qes79/SUTa7eP3iMj+hFdMRK62v3vBvqbzXW167+occVdN42GWVAYmI/H1mLWQKDCQ0JilKXlTvuBsxFz+/70am2vLqEgQpHPib9Yxk0m/7e7v86a32XddZETEC3wReCewA/igiOxYVuwjwLiqbga+AHweQFW/rKpXq+rVwIeAM6q6P+G8e5zvVXUoxbeyKrV2Y2bcmJPL4GQkPrVyOdQmTJcZkkd/eO6SxH85tQkjGUPyGAhHXjddmQ5cFxngeqBDVU+r6gLwVeCuZWXuAh613z8O3Cqv7yZ90D434/D7vFSVFJoNmUnGCilz+UP/Ur+PUr/PNGZJxtqIefmNmd9n5QAyI8z1MRCO8PVXeliyEy71jM/GPcvSSSaITBPQnfC5xz62YhlVjQJhoGpZmV8C/mXZsUfsqbI/XkGUABCR+0Vkn4jsGx4evtx7uCj1QbMhM5lEFpcYnj4Xzv9yqS33G5FZJ91js3znYD+qykI0Rn94Lh6f7HIxe2XWz0e/8gof+9prfHVvF9GlGD3jc7SlKD7ZhcgEkVk3IvJGYFZVDyUcvkdVrwJutl8fWulcVX1AVfeo6p6ampqViiSF+vKAWZNZJ12js7zYMQJY0ZNVOW8H+eVQbxqzdaGqfOihl/iNL7/Cc0eH6LMzlbauMw5WbXmAISP+l014dpGXu8YB+PorvfRNRIjGlPYUxSe7EJkgMr1AS8LnZvvYimVExAcEgcS0hnezbBSjqr32v1PAV7Cm5VzD2pBpPTQL0RgvHB9i0SQyWzNLMeUX//ZF7nnwJV7pGqd77PyIvpeLsyHTcHl0jc1y1g6w+MQrPfEgmOuulzI/QyYaw2VzuC+MKuxpq+CVrnFe7bYEp7UqP0cye4EtIrJBRAqxBOOpZWWeAu61378PeF7t7cAi4gE+QMJ6jIj4RKTafl8AvBs4hIs0BAOMzSwQWVziC8+e4L5H9vIFk2Z2zRzomWDUTqr0r/t64iKz3vDkDfY05uJSDFXl1a5xI/6XwAE7xfietgqePzYUT4edDPEfmpo3GWUvk1N2/qrfeMsmVOHBH54BYEN1Ho5k7DWWjwJPA0eBr6nqYRH5jIjcaRd7CKgSkQ7gY0Cim/MtQLeqnk445geeFpEDwH6skdA/pPZOLoyz8a9vYo4nXu4B4JspShKUiziN2fUbKvnOoX5ODc8QKPCs2+d/Y00p0ZjSPTbLoz8+yy/+7Y/57LeOJsPkvMBJxvdff24T89EYj7x4lrKAj9qy9dVLa2UxSzGN52HqHpuNL2AbLs5AeA6fR3jLFbXUlwc42Bumvjyw7nq5HFwXGQBV/baqblXVTar6WfvYp1T1Kft9RFXfr6qbVfX6REFR1RdU9YZl15tR1etUdZeq7lTV31ZVV5OGOz2IL7/UxdDUPFe3hOgem2N02kwJrIWzozOUFHq598Z2JmYX+dKPz7KjoRyPZ217MVYjMdnZP7/UBcC/vdpLzDRoa6I/HKG6tJCbt1RTUuild2KOa1sr1l0vG+3d/6dHpnn68AA3/8X3+KNvuDoZkVX0T1h7Yrwe4V1XNQBw+466Ne9dSiYZITL5wPaGcjwCD/3oDCWFXn7jLZsAOG3SMq8JJ9nSz11Rg9duwG7YuNzB8NJxGrMn9/fRMTTN7pYQ4blFzoyaelkL/eE56oMBAgVe7rupHYD//MbWdV93ox2t+fTwDP/wA6tP+fjL3cxHXe0rZg39CXtiPv6OrfzP9+/mk+/a5ootRmTSRInfFw/g+PO7GuJBGp2sdIYL02/nJyn1+/iTX9jBruYgv3xD28VPvAjlgQJaKov41sF+Cr0e/sftWwE4ZUdsNlwYa6+S5eH38bdfwb4/uo137Kxf93UrSwqpKC7g6cMD7Osc59rWEItLyqHeyXVfOx9ITE5WXOjjvdc1U1zoThQxIzJp5M/uupL73tTO792xjaZQER4h7o1juDD9E3M02jHgPnRjO0999M1J21jmTCe8e3cDu5qtSMCmXtZG38QcjSGrMRMRqpMYF+va1gr2nrW8ov7gXdsBODVsxP9iOMnJGl3Y3b8SJkBmGtlSV8af3Lkz/rkhWESXmZa5KAvRGMPT8ykLifG7t21lZ2OQ27fXUVToJVhUwFlTLxdlZj7KZCSasnq5+/pWnj8+xPuubebqlhA+j3DWTC9flHPJydK/u38ljMi4SGtlMd3jJmPmxbASWBHvMSebQIGXO3c3xj+3VxWbacw1sDxvTLK5fUcdL33yVmrKrAjNLZWmXtZCquvlUjHTZS7SEDKhZtbCuYcmPT2zlspiM122BgbSUC+15eeiObdVFcddpg2rszw5mdsYkXGRhmCAgcmI8f+/CH0Tr0+2lEoaQ0X0h03634uxUnKyVNJaWUzPuBH/i2FGMoY4DcEilmLKiNkrc0HiD02aIsg2BgMsRGPxCAOGlemPp79OT2PWECxiMhI16bIvQn94zpXkZKthRMZFTDKztWGlrbbC8qcDR8ycEZRhZfrDEapKCs/LTJpKnDW5flMvF6Q/HKG2LP3JyVYjM6zIU5w504GweWguRN/EXFqH/k1xkTHifyGcjZjpwnFZ7zOdsgsyEI7QkCInmcvBiIyLOAumpjG7MFY63/S5YzqCZkYyF2YgzfXSaEaYa6I/nN5O2cUwIuMiFcUF+H0ekzHzIvSHIylzX16JypJC/D5PfGHbsDKW+KevXurK/HjETJddCGsj5rkoDJmAERkXEREagiaZ2YWYjy4xMj2f1odGRGgMFZkR5gWYXYgSnltM63SZz+uhrjxAr6mXVZmajzK7sGRGMoZz1AcDpmd2AQbDludduueYG0MB+sxIZlWcjlE6R5iA3Skz9bIa/bYA1xmRMTg0BIvMSOYCOA1KY5pDZDQEi8zc/wVwGrN0T8tYI0xTL6vh7CNqrjDTZQabBjsts8lfsjLOaCLdu5cbQ0UMTc2bLJmr0O1SY9YYKqLPbJRdlXjG2Ir0p1lejYwQGRG5Q0SOi0iHiHxihe/9IvKY/f1LItJuH28XkTkR2W+//l/COdeJyEH7nP8jbmTrWQMNwQDRmDIyYzZkrkTPmCUyaW/MggFUMWF/VqF7bBavR9I+9282yl6Y7vE5igq8VJcWum1KHNdFRkS8wBeBdwI7gA+KyI5lxT4CjKvqZuALwOcTvjulqlfbr19POP53wK8BW+zXHam6h/XgRErtN4uZK9IzPkdNmT9tG/4cHHdZM5W5Mj3jVoj/dG/4czbKmudlZbrGZmmpLHIlA+ZquC4ywPVAh6qeVtUF4KvAXcvK3AU8ar9/HLj1QiMTEWkAylX1p2qNq/8ReE/SLU8CZtf/heken6XFhfllZ0HbzP+vjFUv6Z+ScTbK9pp6WZGzIzO0VmbOVBlkhsg0Ad0Jn3vsYyuWUdUoEAac3LsbRORVEfm+iNycUL7nItcEQETuF5F9IrJveHh4fXdyGTSYXf8XpGd8jmYXGrP4RllTL69DVTkzMkNblRv14nTKTL0sZ3YhyqnhaXY0Bt025TwyQWTWQz/QqqrXAB8DviIi5ZdyAVV9QFX3qOqempqalBh5ISpLCin0esxIZgUii0v0jM/S7kJjVuL3ESwqMCOZFeidmGNidtGVxszZKGvq5fUc6ZskprCryYjMcnqBloTPzfaxFcuIiA8IAqOqOq+qowCq+jJwCthql2++yDUzAhGhvbrYpJVdgeMDU8QU13pmjaEiM/e/Aof7JgHY2XhJ/bmkEN8oazplr+NnZ8cA4inEM4VMEJm9wBYR2SAihcDdwFPLyjwF3Gu/fx/wvKqqiNTYjgOIyEasBf7TqtoPTIrIDfbazYeBJ9NxM5fDtvpyjvZPuW1GxnGk373GDCxPJjP3/3pe7hynwCtsr3epXkIBM5JZxo9OjvDwj85ybWuI2jSlXlgrrouMvcbyUeBp4CjwNVU9LCKfEZE77WIPAVUi0oE1Lea4Od8CHBCR/VgOAb+uqmP2d78JPAh0YI1wvpOO+7kctjWU0Tsxx2Rk0W1TMopXOscJFhW4trGsMVRE7/ic2ZORwEA4wjde7eXNm6spKkyvx59DQ9CMMBP5t1d7+OWHXiJQ4OFP7tzptjmvIz0JOi6Cqn4b+PayY59KeB8B3r/CeU8AT6xyzX3Alcm1NDXsagoBsPfMGLdur3PXmAxhaCrC04cHeOu2WtfcMXc2lvNPP+3kzMgMG2tKXbEhk9jfPcGHHnqJ+WiM37p1i2t2NIWKGJyKEFlcSrtreybyN893sLs5yGP/9caM/P9wfSRjgOs3VFIe8PHN1/rcNiUj6J2Y471/92MWl5Rfu3mja3Zc11YBwCtdE67ZkEn88TcOUeb38d3fvplrWitcs2NnYzmqcLgv7JoNmULvxBynh2e48+qmjBQYMCKTERT6PPyna5t58rU+ntyfkf4JaeVPnzrM2PQCX73/Bq500VNmU00plSWFPHtk0DUbMoX+8BwHe8Pcd1O766O6q1tDALxqxJ+fnBoF4KbNVRcp6R5GZDKE379jG29or+R3H9vPa90TbpvjGkNTEf7jyCAfefMGdreEXLXF4xHufkML/3FkgFe6xl21xW1e7LAaszdvTr+b/3JqywK0Vhbz/RPp39eWaRwfmMTv87C1tsxtU1bFiEyGUFTo5eH73kB5UQGPvHjGbXNc48d2Y3b7jnqXLbH4tZs30lxRzP3/uC+vPc0O9kxQ6vexrT4zGrP3XtvMD0+OcHIwv70yO0dnaaksxuPJnDAyyzEik0GU+n289YpafnhyJG+jMv/s7BjlAZ9rbsvLqSgp5OH73sBUJMrffq/DbXNco3NslraqzGnM7rmhlWBRAb//xAGW8vRZAStWWVuGhZFZjhGZDOPmLdWMzixwdGDSbVNc4ezIDJtqSzOmMQPYXFvKbdvrePboYN66M3eNzroSRmY1qkv9/OmdO3mlayJvR/6qStfYLK0ZVC8rYUQmw7jaXoc40pefItM5Okt7VYnbZryOW7ZWMzg5z4nB/IvMsBRTusdnaa3MrHq56+pGbtpcxZd+fDYvxX90ZoHZhaWMC4i5HCMyGUZbVQl+n4fjA/k31zwfXaIvPJeRD80NGy3vnVfz0AGgb2KOxSXNqJEMWCFm7riygZ7xOc6MzLhtTtrpHLUSlGVavSzHiEyG4fUIW+pKOZ6HC5rdY3OoQnt15j00LRXFFBd6OZaH4t9lZ1vMxLn/mzdXA/DT02MXKZl7OFkwM7FTlogRmQxkS20Zp4byb1qma8zqjbZl4HSZxyNsqSvLyxGm02POxLn/tqpiSv0+juXhGqZTL26kwrgUjMhkIK2VxfRPRpiPLrltSlo5O5K5PWaAbXVlnMjDEWbn2AwFXonn2MkkRIRt9WUcy8MAs51jM9SXBzJ2p7+DEZkMpKWyGFXoy7MggF1js5T5fVSWZE5+8kTaqosZnVlgej7qtilppWvUyoLpzSCPv0S2NZRxdGAy7xb/u7PAswyMyGQkzhyrMxeeL5wdnaG1qjij8pMn0mZ7V3WN5le9dI5mdmO2sbqUqUiUsZkFt01JK52jsxm/HgNGZDKSfBWZrgx1X3bIx3px9mJk6hQmnPOu6syjeplbWGJoaj6j68XBiEwGUlvmp9DniXuP5APRpZi1FyODe8yOyORTvYzZ04OtGSz+jsjk0wizezxznTGWkxEiIyJ3iMhxEekQkU+s8L1fRB6zv39JRNrt47eLyMsictD+920J57xgX3O//apN4y2tC49HaK4oyqvGrD8cYXFJac/ghyZYXECwqIDOsfzZk+GMDjK5Xhzvqs48EhlHULNhusz1pGV2+uQvArcDPcBeEXlKVY8kFPsIMK6qm0XkbuDzwC8BI8AvqGqfiFyJlV2zKeG8e+zkZVlHa2VxXk3LxN1kM2xX+XKsesmfQJldWbDhL1Dgpb48kJfinw0ikwkjmeuBDlU9raoLwFeBu5aVuQt41H7/OHCriIiqvqqqTqavw0CRiPjTYnWKaa0spmt0Nm88Zs6OWg1EJm7ETKS1qjivRpido7OIZP5ejNaq4ryaLusanaE0gz0xE8kEkWkCuhM+93D+aOS8MqoaBcLA8iw97wVeUdX5hGOP2FNlfyyruCyJyP0isk9E9g0PZ05+itbKYqbmo0zMLrptSlroGpul0OehrizgtikXpLWymJ7x2byJ/JstezHaKovzauH/zOgs7dWZ64mZSCaIzLoRkZ1YU2j/NeHwPap6FXCz/frQSueq6gOqukdV99TUuJ+QycEZBufLg3N2ZIa2DM+LAVa9LC4pfXmSWyZb3GTbqooZnppndiE/9jCdGZlmQ7W7GUrXSiaITC/QkvC52T62YhkR8QFBYNT+3Az8G/BhVT3lnKCqvfa/U8BXsKblsoYN1dbaxNk8CfzXNZZZoeRXw3GxzpdF5s4MC/G/Go73Wz6sY85Hl+gdn2NDFtQLZIbI7AW2iMgGESkE7gaeWlbmKeBe+/37gOdVVUUkBHwL+ISqvugUFhGfiFTb7wuAdwOHUnsbyaWlshiRc2sVuYyq2o1ZZi/6w7k1o3yol8nIIiPT87RXZ0G92A2uE5ool+kemyWmsKEm8+sFMkBk7DWWj2J5hh0Fvqaqh0XkMyJyp13sIaBKRDqAjwGOm/NHgc3Ap5a5KvuBp0XkALAfayT0D2m7qSQQKPDSGCzKi5FMfzjC3OJSVjRmdWUBAgWevKgXJ0jr5prMn5Zx/nbyQfxPD1v3mC3TZa67MAOo6reBby879qmE9xHg/Suc9+fAn69y2euSaaMbtFUVczYPpmVODWdPY+bxCG2VJXlRLx2OyNRmfr2UBwqoKinMC/E/bd/jhiwY+UMGjGQMq9NeXUJnHvTMnMZsU212PDTt1cV50WPuGJ6m0OvJioV/cDpluV8vh3rDNAQDBIsL3DZlTRiRyWDaq4oZn10knONuzB1D05QHfNSUZscWp/aqErpGc9+N+dTQNO3Vxfi82dFMtFeX5MWazP7uiXia9mwgO/568hTHk+lMjvfOTg1Ps7m2NCt8/sFqzBaWYvSHc9uN+eTQdFZMlTlsqCphYDLC3ELu5mEamZ6nZ3zOiIwhOWyyH/CTOZ4oq2Nohk1ZsB7j4Ih/LveaJ2YX6BydZWdj0G1T1kybvfify+FlXuueADAiY0gO7VUlFBV4OdKfu6llB8IRRqbn2dFY7rYpa2aT7Tp6cih3xX+/3Zhdk0WN2Yaq3N9btr97Aq9HuKo5e8TfiEwG4/UI2xrKONKXuyLjNGa7s6gxqynzU1lSyNEcFv/93ROIkFWNmbNvxHEkyUX2d0+wta6M4sKMcAxeE0ZkMpwdDeUc6c/d1LIHeibweYQdDdkzkhERtjeUcTSH88rv755gS20pZYHs8GACKPX7aKks4thAbtZLLKZZt+gPRmQynh2N5UxFovSM5+Yi82s9E2xrKMv4AIzL2V5fzvHBKaJLMbdNSTqqymtZ2JgBbKsvz9kR5umRGaYi0ayawgQjMhmP08PPxXWZWEw50BNmV3PIbVMume0N5SxEY5zJwfn/rrFZxmcXubqlwm1TLpntDeWcGZkhsph7HmbO1PLVrSFX7bhUjMhkONvqy/EIHO4Nu21K0jkzavXMrs5SkYHcFP9z62TZsx7jsL2+jJjCycHcW5fZ3z1Oqd+XVZ6YYEQm4ykq9LKltowDOSgyB3omANiVhY3Z5tpSCrySkyLzatcERQVerqgrc9uUS2abLf5HB3KvXvZ3T7CrOYg3w9NhLMeITBawqznIgZ5wzi3+v9YdpqjAmxUxy5ZT6POwta6Mw72515gd6JngyqbyrNnpn0hrZTFFBd6cW5eJLC5xrH8qq7wwHbLvrygP2dUSYmxmIecW/1/rmeCqpmBWNmYAVzUFOdibW+IfXYpxuG8yK9fJwHL7396Qe+J/fGCKaEzZ1ZR9o/7sfLrzjN32XoUDPbkzZba4FONI3yS7smgfxnKubAoSnlvMKfE/MTjNfDSW1fWyqznEob5wTsWWO2hPl19pRMaQCrbVl1Po9cTXMHKB4wNTVmOWhcN/h105KP4HeycAsnYkA1a9zC4sxVNI5AIHe8KEigtorihy25RLxohMFlDo87C9oYzXckhknJ7Z7izuMV9RX0aBV+L3kgsc6AlTFvDRliXh/VciN8U/zFVNwawJIptIRoiMiNwhIsdFpENEPrHC934Recz+/iURaU/47pP28eMi8o61XjPb2NUc4lDvJLEcmQI40DNBsKgga3KVrITf52VrXRmHckhknMbMk2UeTIlsqC6lpNDLwRzplEUWlzgxOMVVWThVBhkgMiLiBb4IvBPYAXxQRHYsK/YRYFxVNwNfAD5vn7sDuBvYCdwB/K2IeNd4zaxiV3OQ6fkop0dyYwrA2oSZnT2zRHJp8X8+usTR/smsile2El6PsLMpmDNu/86if86KjIi0i8hfisjXReRBEfmoiLQl0YbrgQ5VPa2qC8BXgbuWlbkLeNR+/zhwq1it013AV1V1XlXPAB329dZyzazCcV18rTv7H5zI4hLHB7K3Z5bIVc25s/h/YmCaxSVlV1PIbVPWze7mIEf6JlnMgbA/2bzoD2sbyTwJHMMaGdwO7AZ+ICJfFJFkpDJsAroTPvfYx1Yso6pRIAxUXeDctVwTABG5X0T2ici+4eHhddxGatlUU0pxoTcnFv+P9k9a7phZ3mMG4kKZC/P/zppfTtRLc4j5aIwTOZCL6VBv9i76w9pExquqD6nqc8CYqv4asAk4CzyQSuPSgao+oKp7VHVPTU2N2+asitcjXNkU5LUcaMycnlk2ezA55NLi/8GeMBVZ3Jgl4uwnOZgjz0u2LvrD2kTmWRH5qP1ewRpNqOpfAjcmwYZeoCXhc7N9bMUyIuIDgsDoBc5dyzWzjt3NQY70T7IQze4pgNe6w1SXFtIQDLhtyrrJpcX/A71hrmoOZW1jlkhbVTFlAV/Wr8s4U8vZOlUGaxOZjwFBEdkHNNrTS78sIl/EaujXy15gi4hsEJFCrIX8p5aVeQq4137/PuB5tVZanwLutr3PNgBbgJ+t8ZpZx+6WEAs5MAVwsHeCXTnSmEFuLP47HkzZuKN8JUTEDsc04bYp6yKbd/o7XFRkVDWmqp8FbgHuB+qB64BDWN5b68JeY/ko8DRwFPiaqh4Wkc+IyJ12sYeAKhHpwBK9T9jnHga+BhwBvgv8N1VdWu2a67XVbXbb00vZvF9mZj5Kx9B0Tiz6O+TC4v+R/kmWYpr1nmWJ7GoOcXxgKqvD/mf7oj/AmnN4quos1mgg6SMCVf028O1lxz6V8D4CvH+Vcz8LfHYt18x2miuKqCgu4EB3mHve6LY1l8fhvklimhuLyw6Ji/8tWbrvx1m72J0D62QOu5qCLC4pxweyM7AkZP+iP2TAPhnD2rGmAEJZPZJxpi9yqcecC4v/r/VMUFPmp648GQ6jmYHzN5bN6zLZvugPRmSyjt3NQU4MTjG7EHXblMviQE+YhmCA2rLsX/R3yIXF/4M9YXZleWO2nKZQEZUlhRywk7BlG876azZPlYERmaxjV3OImFrTTtmI0zPLNXY1Z+/i/+xClI7h6axvzJbjLP5n6wjz9Ii1OXZbffYlj0vEiEyW4WSRfC0Le2fhuUXOjMxk7fz4hcjmsP/HB6ZQhR2N5W6bknR2NVkj/5n57Bv5Hx+wvEi31Wd3vRiRyTJqywI0BANZucPcmU7KxZFMNu/8dxqz7VnemK3ENa0VxDQ76+XYwBQFXmFDdYnbpqwLIzJZSLb6/+eyyGTz4v+xgSmKC71Z7cG0Glfbo+ZXu8fdNeQyOD4wxcbqUgp92d1MZ7f1ecqu5hBnR2cJzy66bcolcXxgirpyPxUlhW6bknSyefH/2MAkV9SXZXV4/9WoKClkY3UJr3ZNuG3KJXN8YIorsnw9BozIZCXOXoYDdhbDbOHYwBRX5OCUjEM2Lv6rWvtIsn1x+UJc3Rri1a7xrKqXqcgivRNzRmQM7nBlk9VQH+3PHg+z6FKMjuHpnG7MsnHxf2hqnvHZxaxfXL4Q17RWMDK9kFX1cnLIyht1RV32Py9GZLKQUHEhNWV+Tg5mTwKzs6MzLERjOfHQrEY2Lv47cfC25nC9XNsaAuCVruxZl+kcnQGgPcsX/cGITNaypbY03tvJBk7YgpgLw//VyMbF/7MjVmO2sSb7G7PVuKKuDL/Pk1Vh/ztHZxGBlsrsd8YwIpOlbKktpWNoOmvmmc+M5E7PbDX8Pi9X1GfX4n/n6CyBAg+1ZbkTTmY5Pq+HzbWlnMiiTlnn6CyNwSL8Pq/bpqwbIzJZyua6Mqbno/SHI26bsia6x2apLi2k1L/mmKxZSbaF/e8cm6W1sjinwsmsxNa6Mk5mUYqMztEZWrM02OpyjMhkKVtqSwGyZsqsc3Q2Zx6aC5Fti/9do7O0Vubu6NJhS10p/eEIk5HscPvvHJ2lvTo3nhcjMllKe5XVMHSNzbpsydroGssPkXEW/7NhykxV6Rybob0q9+tla621FpgNzjLT81FGZxZyRvyNyGQptWV+Cn0eurNAZOajS/SF52ityo2H5kJsrLFGmGdHM79ehqbmiSzGaMsHkalzRCbzp8zinmU5Ui+uioyIVIrIMyJy0v63YpVy99plTorIvfaxYhH5logcE5HDIvK5hPL3iciwiOy3X7+arntKFx6P0FJRRFcWNGY943OoQlsejGRK/T4qSwqzYoTZaf/t5IP4N1UU4fMInVlVL7nxvLg9kvkE8JyqbgGesz+fh4hUAp8G3ghcD3w6QYz+SlW3AdcAN4lIYjrox1T1avv1YErvwiVaKovpHs/8h8YRwnzoMYNVLz1ZUC9OjzkfxN/rEZoqirJi5N8Zf15yQ/zdFpm7gEft948C71mhzDuAZ1R1TFXHgWeAO1R1VlW/B6CqC8ArQHPqTc4cWiuL6RqdzXhPJqcxy5We2cVoqSjKmpGM0/jmAy0VxXRngUNG5+hMTnliui0ydarab78fAOpWKNMEdCd87rGPxRGREPALWKMhh/eKyAEReVxEWlYzQETuF5F9IrJveHj4cu7BNVori5majxKey2yPmc6xWYoLvdSU5u5ejERaK4vpHZ9jKZbh4j82S1OoiAKv281AemipLM6KkczZHHJfhjSIjIg8KyKHVnjdlVhOre74JT+VIuID/gX4P6p62j78TaBdVXdhjXweXe18VX1AVfeo6p6amppL/XlXabH/EDO919ydJ3sxHFori4nGlP5wZveau0Zn8mYKE6zd82MzCxmfwOzsyGxObVpOucio6m2qeuUKryeBQRFpALD/HVrhEr1A4kik2T7m8ABwUlX/d8JvjqrqvP3xQeC6JN5SxtBSYTUQ3WOZ3Zjlyx4ZB0f8M71ezo7O5pfIOM9LBq+XzS0sMTAZYUOOrMeA+9NlTwH32u/vBZ5coczTwNtFpMJe8H+7fQwR+XMgCPxO4gmOcNncCRxNrtmZgRPXqHNsxmVLVicWUzrHcqtndjFa4yKTuY3ZxOwC4bnF+H6rfMCpl0z2yDybQ4ExHdwWmc8Bt4vISeA2+zMiskdEHgRQ1THgz4C99uszqjomIs3AHwI7gFeWuSr/lu3W/BrwW8B96bypdFEWKKCypDCjG7PBqQgL0fzYi+HQEAzg9UhGT2M6+3jySWTiI8wMXvx3ApZme8rlRFx1X1DVUeDWFY7vA3414fPDwMPLyvQAK07yq+ongU8m1dgMpbWyOO7ymImcHbHdMXNk9/Ja8Hk9NIYCmS0y8YCl+SP+FcUFlBR6M7pTdsaMZAyZRltVZotMfC9GHo1kwBLVjBaZ0RlEoLkif+pFRDLew+zM8AzVpf6ccV8GIzJZT1tlMf3hORaiMbdNWZGzo7MUeIXGUH7sxXDI9MbMCSUfKMj+UPKXQmtlcUaL/5H+yZzLHmtEJstprSohpmTsDvOusRlaKorxevLDfdmhtbKY0ZkFpjPUXfb4wFROJypbDSdKRiZuYJ6PLnFicIor7SCruYIRmSzHmYbK1JhMp4Zmcmp+ea1ksidTrjZma6G1spjIYozh6fmLF04zx/qnWFzSeCTvXMGITJbjiIyzkJtJRBaX6BieZkdDudumpB2nXroy0L38xMA00ZhyZWNuNWZrIZPdy793fAgRuH5DpdumJBUjMllOTamfypJCjvVnXgjzE4NTLMWUnY35JzKOu2wmOmX8+NQIAFe3htw1xAWcvWWZtlF2PrrEk/v7eENbJTU5lgrbiEyWIyLsaCjnSP+k26a8jle7JgDYmYc95mBRATVlfo5nYP6S7xwa4MqmcpryzBkDznnTZZr4/69nTnBmZIbfeMsmt01JOkZkcoAdjeUcH5xiPrrktinn8fThATbVlORN9OXlbG8o52iGjTC/8Wov+7sneM/VTRcvnIMECrw0hYo4NZw5GTI7hqZ44Aen+eD1Lbx1W63b5iQdIzI5wA0bK1mIxvjJqVG3TYnzatc4Pz09yruuarh44Rxle0MZHUNTGeNefmp4mt97/ABv3FDJh25sc9sc19hWX8axgcwZ+T/84lkKvR4+/vYr3DYlJRiRyQFu2lxNqd/HY3u7L144DZwenua+R/bSVFHEf7lpg9vmuMY1LRUsLikvd467bQoAX3y+A59X+OI91+L35df+mES2NZRxangmY0b+3z8+zNu21VKVo6kwjMjkAH6fl/tv2ch3Dg3wqScPub4H4AvPniSmypc/cgMVJYWu2uImt2ytxu/z8K2DfW6bwlJMefboIO/e1UB1jjZma2VnY5ClmHKoN+y2KfSMz9I7Mccbc8yjLBEjMjnCf3vrZn7lpnb+8Sed/MeRQdfsmIws8vThAd57bXPersU4FBf6uHN3I195qYsvv9Tpqvgf6ZtkMhLlTZuqXbMhU7hpczU+j/DMkZUyi6SX/d0TAOxpNyJjyHC8HuEP37WdhmCAf93X45odr3VPsBCNcdv2lZKc5h+fvnMnN22u5g//7RDfPNB/8RNSxP6eCQDekMM95rUSLCrgzVuq+ereLl447q7QnB629lFtqil11Y5UYkQmh/B5Pfzc1hr2dY651ms+0GNNQeTaruXLpdTv40u/cj2bakr4p5+cdc2OU0PTlBR6aQwGXLMhk/jjd++gsriQ+x7Zy9dfca9Tdnp4mqZQEUWFubtGZkQmx9jVHGJidtG1IICHesO0VRUTLC5w5fczEa9HuH1HPfu7J4gsurPYfGp4mo01pXmTAvtibKop5Tu/czNXNpXz/75/yjU7To/M5HwMOVdFRkQqReQZETlp/1uxSrl77TInReTehOMviMhxO2HZfhGptY/7ReQxEekQkZdEpD1Nt+Q6u1usEcRrPe4sap4ZmWFzDg/9L5drWkMsLimH+9xxnT09PMOmHG/MLhW/z8uduxs5MTjN8JQ7sczOjszkfOI4t0cynwCeU9UtwHP25/MQkUrg08AbgeuBTy8To3tU9Wr75UywfgQYV9XNwBeAz6fyJjIJZ273zHD6Y2apKj3jc/GQKoZzXNMSAqw1q3QzH12id2IuLwOVXozr2qym5JWu9LuZT89HmYxEcz4NhtsicxfwqP3+UeA9K5R5B/CMqo6p6jjwDHDHJVz3ceBWyZN5gkCBl/ryAJ0uBGacnIsyPR+luSK3H5rLoabMT1nAF8/hnk6GJq1eemPQ1MtyrmwK4vOIK+I/EI4AVrruXMZtkalTVcflZgBYySWpCUjcZdhjH3N4xJ4q++MEIYmfo6pRIAxUJdXyDKa1yp2EWd12Tpt8jIl1MUTEtSymfRNWMMj6HG/MLge/z+taCvN8EZmU5/gUkWeB+hW++sPED6qqInKpLlH3qGqviJQBTwAfAv7xEu27H7gfoLW19RJ/PjNpqyzm+yeG0/67PeNWY5ZPKX0vhbbKEg73pX+tbGDSaswaQ7ndmF0ubVXFnHEhVUZf2HpeGnJ8hJnykYyq3qaqV67wehIYFJEGAPvflZzWe4GWhM/N9jFU1fl3CvgK1prNeeeIiA8IAisG9lLVB1R1j6ruqampWe/tZgQtlcUMTc2n3ZPJyc5ppstWprWqmJ7xOaJL6Y1l1jdhiUx9jjdml0tbVQmdozNpd/t3RjK15bkdgcHt6bKnAMdb7F7gyRXKPA28XUQq7AX/twNPi4hPRKoBRKQAeDdwaIXrvg94Xt2OtZJGnIVE5484XfSMz1FS6CVk3JdXpKWimGhMGUqzJ9NAeI6ygI9Sf8onLrKS9qpiZhaWGJleSOvv9ocjVJUUEijI3T0y4L7IfA64XUROArfZnxGRPSLyIICqjgF/Buy1X5+xj/mxxOYAsB9r9PIP9nUfAqpEpAP4GCt4reUyzrSIMxefLvom5mgMFZm9GKvgWr2EIzk/778emuzp3XTXS394joY8mMJ0tWujqqPArSsc3wf8asLnh4GHl5WZAa5b5boR4P1JNTaLcBbee9P80AxORszi8gVwRph9aR5hDoQjOT/vvx4cAe6bmGO37WqeDgbCkbyYWnZ7JGNIAU5D35/uxmwyQn25EZnVSGzM0km/GclcELfEvz9PxN+ITA7i93mpLvWntTGLLsUYnpo3I5kLUBYooCzgoz+N9bIQjTEyberlQlQUF+D3edJaL7MLUcJzi3lRL0ZkcpSmUCCt02Uj0wvEFOrMSOaCNAaL6J1IX495cDI/9mKsBxGhMVSU1pF/vuyRASMyOUtjqCitI5n+uM9/7j8066ExFIj/X6UDR2SM+F+YxlAgvm8lHTgiY0YyhqylIWj1zNLluW0as7XRkOYec3+8x5z7c//roSFYRH8aR5j5VC9GZHKUxlCA2YUlwnOLafm9fOqZrYemUBFjMwvMLaRno6wj/qZeLkxjMMDgVITFNG2UHcijaUwjMjlKut2YBybnKfAKlcWFafm9bKUh7vmXnnrpD0coKvBSHjAbMS9EQ6gI1XOinGr6w3NUFBfk/EZMMCKTszTYIpOuKYDByQi1ZQE8HrMR80I40yN9aaqXAdt92WyQvTANaXb775+I5E2YHyMyOUp8d3maeswDZi/GmmiK78lI10hmzqyTrYH4XpkUjfxVlR+dHGF8xgpdc3Z0hpY82IgJRmRyluoSP4VeT9qmywYnI9QZkbkodUErGGKqRphTkUX+5rmT8VQPnaOztFWZqNgXI9UjmW8d7OeXH3qJDz/8M+ajS5wdnWVLXX5kkDUik6N4PEJ9MJCyxuxAzwT/5Ut7OTYwyeJSjO7xWdpMRsyL4vd5qSlL3UbZ//t8B//zmRP8/hMHGJ9ZYHRmgc21+dGYrYeyQAFl/tRtlH1yfx8AB3vD/Ou+HpZiypbaspT8VqZhRCaHaQwF4o3ZXz19nDd//vmkjWw+951jPH9siM995xjdY7MsLikba0xjthYag+f2ZESXYhztn0zatZ8/ZmXL+PGpUZ45OgicS8ltuDCNoaKUhZbZ3z3Bu66qpzzg44++YQWLzxfxNyKTwzgbMiOLS/z9D07RMz7HV3/Wte7rLsU0nhP9ByeG+clpK1XPphqTQ34tJG6U/Z3H9vPOv/4h3z7Yf5GzLk5kcYnTIzN8YE8zHoG/+O4xIH8as/XSkNApG56a54vf6yA8u/4tAFORRYan5tndHOK91zUDUFTgZVu9GckYspzGYBGDU/N899AAi0vWpsyfnRlb93XPjMwQWYzxG2/ZREzhr589CWBGMmvE2Sg7NrPAvx+wxOUbr/au+7odQ9MsxZRbttbwpk3VjEwvUF8eyItIv8nAqReAP/rGQf7y6eN84dkT675uPIRMqIiPvHkDVzaV88l3bcPnzY/mNz/uMk9pqihiKab89XMnaQoV8f7rmjk1vP40s04K4Tt3N7KtvoyhqXm21JYSLDLJytaCs1H266/0ALC1rpRXuyfWfd2TQ1MAXFFXxu/ctoXq0kJ+862bjPvyGmkMBhibWWBkep7vHbPSlz93bHDd1+1PiFPWXFHMv//3m/nwje3rvm62YEQmh7nazo1xZmSG917bxMaaUkam55mKrG8KoHPU8lzaWFPCb9+6BZ9H+JWbNqzX3LzBcWP++x+cpilUxF1XNzE8Nc/0fHRd13X23jRXFLOnvZJ9f3R7XjVm68XZW/av+3pYWIrx5s3VdI/NMbPOeolHw8hTV3JXRUZEKkXkGRE5af9bsUq5e+0yJ0XkXvtYmYjsT3iNiMj/tr+7T0SGE7771ZWum+tsqy/j2tYQ1aWF3HNDGxuqLe+vsyOz67ruwKSVNtbv8/LOqxo49md38J/f2JoMk/OCa1qtP/PhqXnevrMuvpZ1Zp2jzMHJCOUBH0WFub+LPBU404oP/vA0waICPvCGFsDqpK0HZySTr/uV3B7JfAJ4TlW3AM+xQppkEakEPg28Ebge+LSIVKjqlKpe7byATuDrCac+lvD9gym/kwxERPjq/Tfyw997G3XlAdqr7cZsdJ2NWThCbcIDky9zy8miPhjg5i3VFBd6+eD1rWyottayTo9Mr+u6g5ORvG3IksGVTUEARmcWuHVbbXxh/tTw+uplYHKO6tJCCn35+Zy4HdDoLuAt9vtHgReA319W5h3AM6o6BiAizwB3AP/iFBCRrUAt8MPUmpt9JP5hNycpl7mVAdO/rmvkOw/eu4fIYoxgUQGRRStYpjMNebkMTJrkZOuh1O/jF69p4t8P9HH39a20VRXjEda9jtkfzu+05G5La52qOr6bA0DdCmWagO6Ezz32sUTuxhq5JMa1f6+IHBCRx0WkZTUDROR+EdknIvuGh4cv4xayh1K/jzK/Lz5HfLkMTub3Q5MM/D5v3FEiUOClurRw3UEzB8NmJLNe/uf7d/PTT97K9Rsq4xlmB9ZZLwPhCPXl+evhl3KREZFnReTQCq+7EsvZAnG5yU/uJmFkA3wTaFfVXcAzWKOkFVHVB1R1j6ruqampucyfzx7qg4F1icziUoyR6QVqy0xjlkwagkXrCpq5FFOGp+epMyPMdeHxCFWl5/4P64MBBifn13XNgcn8juuX8ukyVb1tte9EZFBEGlS1X0QagKEVivVybkoNoBlrWs25xm7Ap6ovJ/zmaEL5B4G/uDzrc4/6YID+dYQzH5qaj1/HkDwaQwFOr2NaZnR6nqWY5q0HU6qoLQvQM37505hzC0tMzC7m9fPi9nTZU8C99vt7gSdXKPM08HYRqbC9z95uH3P4IOePYrAFy+FO4GjSLM5y6ssD6xr+57s7ZqqwRjJzl53J1Oltm+my5FJX7o93rNbCxOzCeVsE8ik52Wq4LTKfA24XkZPAbfZnRGSPiDwIYC/4/xmw1359xnECsPkAy0QG+C0ROSwirwG/BdyX0rvIIhqCAYan5okmZADsHpulY+j1HjTz0SX+9oUODvaE48dMmuXU0BQqYmZhicnIuT0ZsZjGnQISmYoscs+DP+VvX+iIHxsw9ZIS6sqtDZrz0YtnMp2YXeCWv/geH/j7n8Y7C846mxnJuISqjqrqraq6RVVvc8RDVfep6q8mlHtYVTfbr0eWXWOjqh5bduyTqrpTVXer6luXf5/P1AUDxBSGp63e2fR8lLf+1Qvc+X9/9LrNgE+83MtffPc4H//X1+LHBuI+/2buP5mslM/kvi/t5S1/+cLrNs9+5+AAL3aM8hffPR4XIZNmOTU4I/ahhHWZR148wwf+/ievq5fvHR9iMhLlaP8kh/usoKfxkDJ5kqBsJdweyRjSjDNsd/74X+wYIRpTZheWeOn06Hllf3bG+nxiaCouQINTEQq9HipLTJrlZNIQOj8t89mRGX5wYpiByQjfO36+1+NrPRPx9/vtcDSDkxE8AtWlRvyTSa3dmRqasp6X+egSf/rNI/zszFg87pzDa93nRvxOmKB+M71sRCbfcKZTHJHZnxAz60DCtBjAsQErFpYqHO61vrM2YvpNPKwk44Sa6bU9zPZ1jse/e7Vr/LyyJwenabVz9xy362ggHKGmzI/XpL9OKs7z4qx5He2fin+39+z5wWZPDk2xuzlIWcDH8YFzI5lQcUFeR2EwIpNnOMN2p4d1cnCarXWlNIWK6EyIBLC4FOPU8DTvs0OTOwJkbcTM315Zqqgu9ePzSDxp1vGBSfw+D9vqy84La6KqnBia4qbNVQSLCjg+aDV6g1Pzpl5SwPJO2SG7s9VWVfy6cDN9ExGaK4rZVl/GMVuM+sPmeTEik2dUFBdQ6PPE5/BPj0yzqaaU9upizibsOO+yE5HduLGKmjI/J5zGbHLeLC6nAK+dydRZkzk9PMOG6hI21ZSeFwlgbGaBidlFNteWcUVdWXwkszzUjyE5vO55GZ6hqMDLTZurXyf+fRNzNIYCbKsv5/jAFKrKwORc3q+TGZHJM0SE+vIA/eEIi0sxukZn2VhTQmtlCV1j5xqz/nhE3yK21JZycmjaemjMrvKUYSUzc8R/ho01JbRVFdM9Nhv3BnRGoE2hIrbUlXJi0GrM+sJzee0mmypEhLpyf1xkzo7O0G6L/8TsIuMzCwC2B1qMxlARW+vLmJqPMjAZoWd8Lu/z+RiRyUOcXf+do7NEY2qNZKqKGZtZYNL2mHHSAzcELZHpGJpmdGaBucUlmvL8oUkVTlrmhWiMrrFZNlSX0F5dQjSm9Ixb9eGMdBpDAbbWlTEViXJicJqpSDS+TmNILnVlgbiL+JmRGTZWl8Qjmp+2RzNO56AhWMRWOxPpvrPjTMwu5n29GJHJQ+rLrYfmtB1ddmNNaTx4Zq/dmDkjmbqgn811ZUzPR9lnL3S2GJFJCY2hIlv8Z1iKKRurS+O9YEdcBhJclbfUWY3Zs0etxFpOHRqSS50dWmZx6Zz4N4XODzbrdMqaQkVsrbOiNzv1YkTGkHc02COZk0OOyJS8zoXWCk/ux+/zssXumT1/zIr605LnD02qaAgVEY0pL9kpsjfWlCR4ndniH45Q4BWqS/zxxsypl3xvzFJFfXmAwckI3WOzLMWUDdUlNNrPS1xk7H8bQgEqSgqpKfPz7BFHZErcMTxDMCKTh9QHAywsxfjJqVGaQkWUBwpe50LbN3EuqJ8jMs5+DSMyqaHJbri+f8L6f95YUxpfNHamY/on5qgrD+DxCNWlfqpKCnnZdnduqTQjzFTQELTSZb/aNQFAe3UJZYECygK++BpZfzhCoc9Dlb1/7Iq6MmYWlhCB1qr8fl6MyOQhjkvljzpG2N5g9YZrSv0UeCXeI+tPWEiuKvVTWVLI8NQ81aWFlPrdTkOUmzi7/r93bIimUBHBogL8Pi81Zf74CLM/fH5E3yvsxFpNoSLKAgXpNzoP2FRjdbK+fdDafOlMUzaFiuIjzN6JORqDgfj+sT3tVvbT9qqSvH9ejMjkIVvsaRaAbfXlgBXivK48EN+n0T8RiTd6ADsarHK7mkPpMzTPcPYwRWMaF3+wHAJ6E9ZkEkOUvPWKWgBu2lyVRkvzi832SP65Y0M0V1gjf3C8Ac9NlyU+Lx+8vpXr2yv55Du3pd/gDCO/JTZP2Vh9bo74jRsr4+8dF9qpyCJT89Hz/Pt/+YY2Xjozyi+9YdX8b4Z1Uh449zhut0UdrHpxXJX7wxHu2HmuXj78pjZKAz7u2FmfVlvziaZQEX6fh/loLN4pA8vDz4nG0D8R4abN1fHv6soDfO3Xb0y7rZmIGcnkIR6P8OCH9/Dxt2/lpk3nHgxn+O/0mhP9+++4sp7Df3oH7zCNWcoQEd5zdSMi8O5djfHjjvgPT82zEI2d50Lu93n54PWtVJhYcinD45H4FNl1bRXx442hIsZnFwnPLjI4Fcn7/TCrYUYyecptO+q4bcf52a4bgpYXjbPDfLlLbKHP9ElSzefeu4s/eNf283bvN4aKmFtc4lCfFdLENGbp5w/euZ3H9nXzgT3N8WOOs8y+zjFUTb2shhEZQ5xG24XW8aIxD036CRR4CRScH0zR8Tr72Rlrasbsh0k/b9pczZsSpsPgnKOG43Ju6mVlTNfUEMfpme09O0ag4Jw7psFdnMbMSb1gxD8ziIuMnSLDuJCvjKsiIyKVIvKMiJy0/61Ypdx3RWRCRP592fENIvKSiHSIyGMiUmgf99ufO+zv29NwO1mP03i93DlOc0WxCeefIbTYPeRXuiaoKimkuNBMQGQC9eUBCn0eXusJ4/NI3kdbXg23RzKfAJ5T1S3Ac/bnlfhL4EMrHP888AVV3QyMAx+xj38EGLePf8EuZ7gI7dUlFHqtP4lt9WUXKW1IFxUlhfFMpNsaTL1kCl6PxDcqb6krw+d1uznNTNz+X7kLeNR+/yjwnpUKqepzwFTiMbG62W8DHl/h/MTrPg7cKqZbflEKvB621lsPzRvaKy9S2pBOrm4JAbCnzdRLJvHGDdb+pOvbV5yEMeD+wn+dqjo5TAeAugsVXkYVMKGqTmL6HqDJft8EdAOoalREwnb5keUXEZH7gfsBWltbL/kGco0/+YWdPPFKD794bdPFCxvSxh+8azvVpX4+dGOb26YYEvjNt27C64GPvHmj26ZkLCkXGRF5Flhpc8UfJn5QVRURTbU9y1HVB4AHAPbs2ZP238809rRXsseMYjKOtqoSPvuLV7lthmEZ1aV+/vDnd7htRkaTcpFR1dtW+05EBkWkQVX7RaQBGLqES48CIRHx2aOZZqDX/q4XaAF6RMQHBO3yBoPBYEgjbq/JPAXca7+/F3hyrSeqqgLfA963wvmJ130f8Lxd3mAwGAxpxG2R+Rxwu4icBG6zPyMie0TkQaeQiPwQ+FesBfweEXmH/dXvAx8TkQ6sNZeH7OMPAVX28Y+xuteawWAwGFKImA7+Ofbs2aP79u1z2wyDwWDIKkTkZVXds9J3bo9kDAaDwZDDGJExGAwGQ8owImMwGAyGlGFExmAwGAwpwyz8JyAiw0DnZZ5ezQoRBXIcc8/5gbnn/GA999ymqjUrfWFEJkmIyL7VvCtyFXPP+YG55/wgVfdspssMBoPBkDKMyBgMBoMhZRiRSR4PuG2AC5h7zg/MPecHKblnsyZjMBgMhpRhRjIGg8FgSBlGZAwGg8GQMozIJAERuUNEjotIh4jkfMRnEWkRke+JyBEROSwiv+22TelARLwi8qqI/LvbtqQLEQmJyOMickxEjorIjW7blGpE5Hftv+tDIvIvIhJw26ZkIyIPi8iQiBxKOFYpIs+IyEn736TklDYis05ExAt8EXgnsAP4oIjkeqq8KPA/VHUHcAPw3/LgngF+GzjqthFp5q+B76rqNmA3OX7/ItIE/BawR1WvBLzA3e5alRK+BNyx7NgngOdUdQvwHElKkWJEZv1cD3So6mlVXQC+Ctzlsk0pRVX7VfUV+/0UVsPT5K5VqUVEmoGfBx68WNlcQUSCwC3YeZpUdUFVJ1w1Kj34gCI7q24x0OeyPUlHVX8AjC07fBfwqP3+UeA9yfgtIzLrpwnoTvjcQ443uImISDtwDfCSy6akmv8N/B4Qc9mOdLIBGAYesacJHxSREreNSiWq2gv8FdAF9ANhVf0Pd61KG3Wq2m+/HwDqknFRIzKGy0ZESoEngN9R1Um37UkVIvJuYEhVX3bbljTjA64F/k5VrwFmyPEss/Y6xF1YAtsIlIjIL7trVfqx09UnZX+LEZn10wu0JHxuto/lNCJSgCUwX1bVr7ttT4q5CbhTRM5iTYe+TUT+2V2T0kIP0KOqzij1cSzRyWVuA86o6rCqLgJfB97ksk3pYlBEGgDsf4eScVEjMutnL7BFRDaISCHWIuFTLtuUUkREsObpj6rq/3LbnlSjqp9U1WZVbceq3+dVNed7t6o6AHSLyBX2oVuBIy6alA66gBtEpNj+O7+VHHd2SOAp4F77/b3Ak8m4qC8ZF8lnVDUqIh8FnsbyRHlYVQ+7bFaquQn4EHBQRPbbx/5AVb/tnkmGFPHfgS/bHajTwK+4bE9KUdWXRORx4BUsL8pXycEQMyLyL8BbgGoR6QE+DXwO+JqIfAQr5ckHkvJbJqyMwWAwGFKFmS4zGAwGQ8owImMwGAyGlGFExmAwGAwpw4iMwWAwGFKGERmDwWAwpAwjMgaDwWBIGUZkDAaDwZAyjMgYDBmOiDSLyC+5bYfBcDkYkTEYMp9byf2YYYYcxez4NxgyGBF5M1YMqQlgCvhPqnraVaMMhkvAiIzBkOGIyHeBj6vqoYsWNhgyDDNdZjBkPlcAx9w2wmC4HIzIGAwZjIhUY2VnjLpti8FwORiRMRgym3ZyMMe8IX8wImMwZDbHsHJ+HBKRfMnQaMghzMK/wWAwGFKGGckYDAaDIWUYkTEYDAZDyjAiYzAYDIaUYUTGYDAYDCnDiIzBYDAYUoYRGYPBYDCkDCMyBoPBYEgZ/z+56p6DSZHKSAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(time, theta) # plot angle as a function of time\n", "plt.xlabel(r'$t$')\n", "plt.ylabel(r'$\\theta$')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "3cc80c42-9be9-4a7c-b15e-53b5d5c8c7b4", "metadata": {}, "source": [ "The plot shows that the pendulum stays near the upright position, and is therefore stable! Notice that the motion apparently consists of a fast component (with the same frequency as the driving) and a slow component that oscillates around the upright position. If not for the high-frequency component that generates an \"effective restoring force\", the system would not have been stable. You can change the parameter $\\omega$ and repeat the steps above to see when the system stops being stable." ] }, { "cell_type": "markdown", "id": "3b556d39-9036-4130-b4bf-0b81d8cd3805", "metadata": {}, "source": [ "Just for fun, we can make an animation of the pendulum using the solution we got." ] }, { "cell_type": "code", "execution_count": 5, "id": "4c361047-02ed-4a24-b645-611b76530b09", "metadata": {}, "outputs": [], "source": [ "import matplotlib.animation as anim\n", "\n", "plt.rcParams[\"animation.html\"] = \"jshtml\"\n", "fig, ax = plt.subplots(figsize=(4,4))\n", "ax.set_xlim(-1, 1)\n", "ax.set_ylim(-0.5, 1.5)\n", "ax.axis('off')\n", "p0, = ax.plot([], [], '-')\n", "p1, = ax.plot([], [], 'o')\n", "\n", "def animate(t):\n", " Zt = alpha * np.cos(omega * time[t])\n", " Xt = np.sin(theta[t])\n", " Yt = np.cos(theta[t]) + Zt\n", " p0.set_data([Xt, 0], [Yt, Zt])\n", " p1.set_data([Xt], [Yt])\n", "\n", "mov = anim.FuncAnimation(fig, animate, frames=len(time), interval=10)\n", "plt.close()" ] }, { "cell_type": "code", "execution_count": 6, "id": "bf6849e7-3ad0-4e9a-b63e-a48d1216d856", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mov" ] }, { "cell_type": "markdown", "id": "eb3244b5-0eca-44e9-9c23-28f86a45e4b5", "metadata": {}, "source": [ "```{admonition} Exercise\n", ":class: tip\n", "**Solving the problem using Hamilton's equations** \n", "From the Lagrangian, you can derive the Hamiltonian for the inverted pendulum, then write down the equations of motion using Hamilton's equations. These equations will already be first order ODE's. Solve them numerically and compare your results with the above.\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "8e5bb857-818a-493e-ba3e-3d36506d0ef3", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.8.10" } }, "nbformat": 4, "nbformat_minor": 5 }