{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pêndulo e caos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 1: Pêndulo simples (regime linear)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O pêndulo simples é constituído por uma haste indeformável de comprimento $\\ell$ e de massa desprezível presa a um ponto em torno do qual pode rodar livremente em um plano fixo e com uma massa pontual $m$ na outra extremidade. A massa é sujeita à ação de um campo gravitacional uniforme com aceleração $g$. O grau de liberdade do pêndulo pode ser representado por um ângulo $\\theta$ em relação à vertical, como na figura abaixo\n",
    "\n",
    "![Pêndulo](./pendulum-small.png)\n",
    "\n",
    "Como a haste é indeformável, a tensão na haste compensa a componente da força da gravidade paralela à haste. A força resultante na massa será então a componente perpendicular à haste:\n",
    "\n",
    "$$ F = - m g \\sin \\theta. $$\n",
    "\n",
    "(O sinal negativo é devido ao fato de que a força age na direção oposta ao sentido positivo do ângulo $\\theta$). Por outro lado,\n",
    "\n",
    "$$F = m \\frac{d^2 x}{dt^2},$$\n",
    "\n",
    "onde $d x$ é o deslocamento infinitesimal linear da massa. Este pode ser calculado como $dx = \\ell d\\theta$, o que nos leva a\n",
    "\n",
    "$$ m \\ell \\frac{d^2\\theta}{dt^2} = - m g \\sin \\theta,$$\n",
    "\n",
    "resultando na conhecida equação diferencial de segunda ordem:\n",
    "\n",
    "$$ \\frac{d^2\\theta}{dt^2} = -\\frac{g}{\\ell}\\sin\\theta.$$\n",
    "\n",
    "Em princípio, como dispomos de equação diferencial que descreve o nosso modelo, seríamos capazes de, dadas as condicões inciais, calcular o comportamento futuro do pêndulo. Infelizmente, a não-linearidade introduzida pelo $\\sin\\theta$ do lado direito da equação impede que sejam usados os métodos conhecidos de solução de equações diferenciais."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linearização"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A saída tradicional para esse dilema é usar um método perturbativo e nos preocuparmos apenas com variações pequenas em torno de um ponto conhecido. No caso expandimos em torno do ponto de equilíbrio $\\theta=0$ e consideramos que o ângulo $\\theta$ é pequeno, portanto podemos aproximar\n",
    "\n",
    "$$\\sin\\theta \\approx \\theta,$$\n",
    "\n",
    "o que pode ser entendido lembrando a expansão de Taylor do seno em torno do ângulo 0:\n",
    "\n",
    "$$\\sin\\theta = \\theta - \\frac{\\theta^3}{3!} +  \\frac{\\theta^5}{5!} + \\cdots,$$\n",
    "\n",
    "onde vemos que para $\\theta$ pequeno (próximo de zero) os termos não lineares podem ser desprezados.\n",
    "\n",
    "Outra forma de ver isto é através de um gráfico."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "θ = np.linspace(0, 20, 100+1) # Variando theta de 0 a 20 graus, a cada 0.2 graus\n",
    "θ_rad =  np.deg2rad(θ) # Precisamos dos valores em radianos\n",
    "pl_linear, = plt.plot(θ, θ_rad, label=r'$\\theta$') # Plota θ no eixo x e θ_rad no eixo y.\n",
    "                                                   # label dá o texto que descreve a curva (aqui usando LaTeX)\n",
    "pl_sin, = plt.plot(θ, np.sin(θ_rad), label=r'$\\sin\\theta$') # Plota θ no eixo x e sin(θ) no y\n",
    "plt.xlabel(r'$\\theta$') # xlabel indica o nome que aparece no eixo x. ylabel seria para o eixo y.\n",
    "plt.legend(handles=[pl_linear, pl_sin]) # Solicita a criação de uma legenda para as duas curvas.\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos ver que para ângulos pequenos (menores que 10°) há pouca diferença entre $\\sin\\theta$ e $\\theta$.\n",
    "\n",
    "A equação aproximada fica então\n",
    "\n",
    "$$ \\frac{d^2\\theta}{dt^2} = -\\frac{g}{\\ell}\\theta$$\n",
    "\n",
    "que é uma equação linear e pode ser resolvida por métodos tradicionais, resultando em\n",
    "\n",
    "$$\\theta(t) = \\theta_0\\cos(\\omega_0 t - \\phi),$$\n",
    "\n",
    "onde $\\theta_0$ é o ângulo inicial, $\\omega_0 = \\sqrt{g/\\ell}$ é a denominada frequência angular natural do pêndulo e $\\phi$ depende da velocidade inicial do pêndulo (pode ser encontrado de $v_0 = \\theta_0\\sqrt{g\\ell}\\sin\\phi$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixando alguns parâmetros do pêndulo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A partir de agora vamos fixar um pêndulo com $\\ell=9.8$ usando $g = 9.8m/s^2$ (de modo que $\\omega_0=1$), sempre considerando a massa sendo largada do repouso (que resulta em $\\phi=0$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ℓ = 9.8 # ℓ é o comprimento\n",
    "g = 9.8 # g a aceleração da gravidade\n",
    "ω0 = np.sqrt(g/ℓ)\n",
    "T0 = 2 * np.pi / ω0 # Período da oscilação natural"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discretização do tempo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ao realizar cálculos numéricos ou visualização no computador, precisamos discretizar todas as funções contínuas. Ao discretizar no tempo precisamos fixar a definição temporal (intervalo entre dois instantes de tempo computados), o tempo total e o número total de instantes de tempo avaliados. Vamos definir uma função para ajudar a calcular fatores de temporização. Como estamos lidando com funções oscilantes, nos baseamos no tempo de um ciclo de oscilação. Então a função recebe o número de ciclos desejados, o tempo um ciclo e, opcionalmente, o número de intervalos de tempo em um ciclo; caso este último argumento não seja fornecido, vamos assumir o valor 100. Com esses valores a função calcula e retorna (na ordem) o tempo total de simulação, o número total de pontos simulados e o intervalo de tempo entre dois pontos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def timing(ncycles, T, n_per_cycle = 100):\n",
    "    return (ncycles * T, ncycles * n_per_cycle, T / n_per_cycle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizando a solução analítica para alguns valores de $\\theta_0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculamos e mostramos então o solução para 5 ciclos de oscilação:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax, nintervals, _ = timing(5, T0) # 5 ciclos. Uso de _ diz que estamos descartando este valor (sem interesse aqui).\n",
    "t = np.linspace(0, tmax, nintervals)\n",
    "hndl = []\n",
    "some_θs = [3, 5, 7]\n",
    "for θ0 in some_θs:\n",
    "    h, = plt.plot(t, θ0*np.cos(ω0 * t), label=r'$\\theta_0 = ' + str(θ0) + '$')\n",
    "    hndl.append(h)\n",
    "plt.xlabel(r'$t$')\n",
    "plt.ylabel(r'$\\theta$')\n",
    "plt.legend(handles=hndl)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solucionando a equação diferencial numericamente"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uma atividade interessante é solucionar computacionalmente essa equação e comparar o resultado com o esperado pela solução apresentada acima.\n",
    "\n",
    "Para isso, precisamos reescrever a equação diferencial como um sistema de equações de primeiro grau (pois isso é requerido pela função para solução numérica de equações diferenciais que vamos usar), da seguinte forma:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "  \\frac{dx_1}{dt} & = & f_1(x_1, x_2, \\ldots, x_n, t)\\\\\n",
    "  \\frac{dx_2}{dt} & = & f_2(x_1, x_2, \\ldots, x_n, t)\\\\\n",
    "  \\vdots \\\\\n",
    "  \\frac{dx_n}{dt} & = & f_n(x_1, x_2, \\ldots, x_n, t)\\\\\n",
    "\\end{eqnarray}\n",
    "\n",
    "O sistema é então definido pelo conjunto de funções $f_i$, pois dado o estado (valores dos diversos $x_i$) e tempo atuais, podemos calcular o valor de todas as derivadas usando as expressões das funções $f_i$, possibilitando o conhecimento da trajetória do sistema.\n",
    "\n",
    "No nosso caso, temos uma equação de segundo grau, então precisamos escrevê-la como duas equações de primeiro grau. A forma natural de fazer isto neste caso é definirmos as duas variáveis:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "  x & = & \\theta \\\\\n",
    "  y & = & \\frac{d\\theta}{dt}\n",
    "\\end{eqnarray}\n",
    "\n",
    "Com isso, ficamos com o sistema de equações:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "  \\frac{dx}{dt} & = & y \\\\\n",
    "  \\frac{dy}{dt} & = & - \\frac{g}{\\ell} x\n",
    "\\end{eqnarray}\n",
    "\n",
    "Com isso feito, podemos usar a função `odeint` do SciPy para fazer a integração das equações diferenciais, lembrando que $x(t=0) = \\theta_0$ e $y(t=0) = 0$ (pois partimos do repouso).\n",
    "\n",
    "Primeiro precisamos definir uma função para o cálculo das funções $f_i$ dado um tempo e os valores correntes de x e y. Conforme requerido pela `odeint`,  o primeiro parâmetro para essa função é tal que, se o parâmetro se chama `xy` então `xy[i]` deve dar o valor atual da correspondente variável. Qual variável corresponde a qual posição deve ser escolhido pelo programador, e mantido consistente por todo o código.  Aqui, escolhemos que `xy[0]` corresponde à variável $x(t)$ enquanto `xy[1]` corresponde à variável $y(t)$. O segundo parâmetro é o tempo corrente. Demais parâmetros são opcionais, e dependem do que estamos modelando: eles podem ser usados para passar parâmetros próprios do sistema. Em nosso caso, o valor da aceleração da gravidade e comprimento da haste.\n",
    "\n",
    "A função deve retornar, em ordem, os valores das $f_i$ calculadas para os valores passados das variáveis e do tempo. No nosso exemplo atual, o sistema é independente do tempo, então o tempo não será usado."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simple_linear_pendulum(xy, t, g, ℓ):\n",
    "    x, y = xy # Separamos os componentes x e y do estado\n",
    "    return y, - g / ℓ * x # retornamos funções das derivadas em x e derivada em y nas equações"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Agora basta usar essa função para resolver o sistema:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Definimos uma função para fazer o cálculo deste sistema, dados os parâmetros.\n",
    "# Isto evita repetição de código abaixo.\n",
    "def simulate_simple_linear_pendulum(θ0, tmax, nintervalos, g, ℓ):\n",
    "    θ0_rad = np.deg2rad(θ0)\n",
    "    t = np.linspace(0, tmax, nintervalos + 1) # Instantes de tempo para os quais queremos resultados\n",
    "    xy0 = [θ0_rad, 0.0] # Estado inicial\n",
    "    xy = odeint(simple_linear_pendulum, xy0, t, args=(g, ℓ)) # Resolve sistema de equações definido por deriv_simples,\n",
    "                                                             # Partindo do estado inicial xy0 para os instante t\n",
    "    return t, np.rad2deg(xy) # Retorna um par: instantes e correspondentes estados (em graus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "θ0 = 5.\n",
    "ncycles = 5 # número de ciclos a simular\n",
    "tmax, nintervals, _ = timing(ncycles, T0)\n",
    "t, xy = simulate_simple_linear_pendulum(θ0, tmax, nintervals, g, ℓ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Agora basta plotar o resultado juntamente com a solução exata para comparação:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exact = θ0 * np.cos(np.sqrt(g/ℓ)*t) # Valor da solução analítica\n",
    "pn, = plt.plot(t, xy[:, 0], label='Numérico')\n",
    "pe, = plt.plot(t, exact, label='Exato')\n",
    "plt.xlabel(r'$t$')\n",
    "plt.ylabel(r'$\\theta$')\n",
    "plt.legend(handles=[pn, pe])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Não é possível distinguir visualmente entre as duas soluções, o que é um bom sinal. Vejamos a diferença entre os resultados diretamente:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t, xy[:, 0] - exact)\n",
    "plt.xlabel(r'$t$')\n",
    "plt.ylabel(r'Diferença em $\\theta$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note que, como a amplitude da oscilação é 5.0, a amplitude do erro é várias ordens de grandeza menor. No entanto, notamos uma tendência a crescimento do erro com o tempo, o que indica que a solução numérica não é confiável para intervalos de tempo muito grandes.\n",
    "\n",
    "Vamos primeiramente tentar aumentar o número de intervalos simulados no mesmo tempo total (aumentando a resulução temporal) para ver se isso tem algum efeito."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax, nintervals, _ = timing(ncycles, T0, n_per_cycle=1000)\n",
    "t, xy = simulate_simple_linear_pendulum(θ0, tmax, nintervals, g, ℓ) # Simlando para intervalos 10 vezes menores\n",
    "exact = θ0*np.cos(np.sqrt(g/ℓ)*t)\n",
    "plt.plot(t, xy[:,0] - exact)\n",
    "plt.xlabel(r'$t$')\n",
    "plt.ylabel(r'$\\theta$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Não se nota diferença significativa, o que mostra que a definição em tempo da simulação anterior já era adequada. Vamos agora realizar a simulação por mais tempo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax, nintervals, _ = timing(50, T0)\n",
    "t, xy = simulate_simple_linear_pendulum(θ0, tmax, nintervals, g, ℓ) # Simula por 50 ciclos\n",
    "exact = θ0*np.cos(np.sqrt(g/ℓ)*t)\n",
    "plt.plot(t, xy[:, 0] - exact)\n",
    "plt.xlabel(r'$t$')\n",
    "plt.ylabel(r'$\\theta$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Veja que há um crescimento quase linear da amplitude do erro com o tempo. Portanto se desejarmos podemos avaliar, dada uma precisão aceitável, até que ponto a simulação pode ser levada adiante."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Espaço de fase (espaço de estado)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Um outro gráfico interessante é o chamado *gráfico de espaço de fase* (ou *espaço de estado*), que consiste em plotar posição contra velocidade. No nosso caso, ângulo contra velocidade angular. Como sabemos, se comenhecemos as equações diferenciais que descrevem o sistema, com a posição e a velocidade em um dado instante todas as posições e velocidades posteriores são determinadas, formando uma trajetória no espaço de fase. Isto é, a trajetória no espaço de fase diz tudo o que desejamos saber sobre o movimento do sistema."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "plt.plot(xy[:,0], xy[:,1])\n",
    "plt.xlabel(r'$\\theta$')\n",
    "plt.ylabel(r'$\\dot{\\theta}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No nosso caso, vemos que a trajetória é fechada em um **ciclo**, isto é, o sistema repete continuamente seu estado."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 2: Pêndulo simples (regime não-linear)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A questão agora é a seguinte: O que acontece se o ângulo $\\theta$ não é necessariamente pequeno? Neste caso, não temos uma solução exata, mas podemos usar o mesmo método para a solução numérica. Voltando à esquação original,\n",
    "\n",
    "$$\\frac{d^2\\theta}{dt^2} = - \\frac{g}{\\ell}\\sin\\theta,$$\n",
    "\n",
    "fazemos a transformação\n",
    "\n",
    "\\begin{eqnarray}\n",
    "  x & = & \\theta \\\\\n",
    "  y & = & \\frac{d\\theta}{dt}\n",
    "\\end{eqnarray}\n",
    "\n",
    "Chegando a:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "  \\frac{dx}{dt} & = & y \\\\\n",
    "  \\frac{dy}{dt} & = & - \\frac{g}{\\ell} \\sin x\n",
    "\\end{eqnarray}\n",
    "\n",
    "Com o qual podemos montar a solução usando SciPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simple_nonlinear_pendulum(xy, t, g, ℓ):\n",
    "    x, y = xy\n",
    "    return y, -g/ℓ * np.sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simulate_simple_nonlinear_pendulum(θ0, tmax, nintervals, g, ℓ):\n",
    "    θ0_rad = np.deg2rad(θ0)\n",
    "    t = np.linspace(0, tmax, nintervals + 1)\n",
    "    xy0 = [θ0_rad, 0.]\n",
    "    xy = odeint(simple_nonlinear_pendulum, xy0, t, args=(g, ℓ))\n",
    "    return t, np.rad2deg(xy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Uma função para fazer comparação da solução linear com a não-linear, dados os parâmetros.\n",
    "def compare_nonlinearity(θ0, tmax, nintervals, g, ℓ):\n",
    "    t, xy = simulate_simple_nonlinear_pendulum(θ0, tmax, nintervals, g, ℓ)\n",
    "    linear = θ0*np.cos(np.sqrt(g/ℓ) * t)\n",
    "    pn, = plt.plot(t, xy[:,0], label='Não-linear')\n",
    "    pl, = plt.plot(t, linear, label='Linear')\n",
    "    plt.xlabel(r'$t$')\n",
    "    plt.ylabel(r'$\\theta$')\n",
    "    plt.legend(handles=[pn, pl])\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fixando o intervalo de tempo e número de intervalos:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncycles = 5\n",
    "tmax, nintervals, _ = timing(ncycles, T0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_nonlinearity(5.0, tmax, nintervals, g, ℓ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note como para $\\theta_0$ pequeno não se vê diferença entre o sistema não-linear e o linear. Vejamos agora com um ângulo inicial maior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_nonlinearity(25.0, tmax, nintervals, g, ℓ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aqui já começamos a ver que o sistema não linearizado apresenta uma frequência menor de operação. Vejamos para um ângulo inicial ainda maior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_nonlinearity(45., tmax, nintervals, g, ℓ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mantendo a tendência de diminuir a frequência com o aumento de $\\theta_0$. Mais dois ângulos maiores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_nonlinearity(90., tmax, nintervals, g, ℓ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_nonlinearity(160., tmax, nintervals, g, ℓ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note como neste último caso fica também claro que o formato da onda, antes aparentemente senoidal, é dependente do ângulo inicial.\n",
    "\n",
    "Vemos que o modelo linear, como esperado, é completamente inadequado para descrever o sistema real, que apresenta não-linearidade, quando fora do regime perturbativo.\n",
    "\n",
    "Para concluir esta parte, vejamos o espaço de fase do sistema não linear para este último valor de $\\theta_0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t, xy00 = simulate_simple_nonlinear_pendulum(160., tmax, nintervals, g, ℓ)\n",
    "plt.figure(figsize=(5,5))\n",
    "plt.plot(xy00[:,0], xy00[:,1])\n",
    "plt.xlabel(r'$\\theta$')\n",
    "plt.ylabel(r'$\\dot{\\theta}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note como a figura não tem mais um formato circular, porém continua sendo uma figura fechada, característica de um sistema periódico (volta regularmente ao mesmo estado inicial, de onde a trajetória se repete)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 3: Atrito"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Já que não precisamos mais nos prender a estudar o que pode ser resolvido analiticamente, vamos agora introduzir um novo elemento no sistema. Até agora, as forças atuando nos pêndulos eram conservativas. Num pêndulo real, haverá atuação de outras forças não-conservativas, as mais importantes se referindo a atrito na junção onde o pêndulo se prende ao teto e resistência do ar. Felizmente para nosso modelamento, nas baixas velocidades típicas de um pêndulo, essas duas forças são aproximadamente proporcionais à velocidade. (Para altas velocidades a resistência do ar é proporcional ao quadrado da velocidade). Vamos então introduzir uma força da forma $-q\\frac{d\\theta}{dt}$, onde $q>0$ é um coeficiente efetivo que engloba as duas contribuições e o sinal negativo indica que a força é contrária à direção de deslocamento do pêndulo.\n",
    "\n",
    "Ficamos então com o seguinte sistema (usando nossas coordenadas $x$ e $y$ como anteriormente):\n",
    "\n",
    "\\begin{eqnarray}\n",
    "  \\frac{dx}{dt} & = & y \\\\\n",
    "  \\frac{dy}{dt} & = & -\\frac{g}{\\ell}\\sin x - q y.\n",
    "\\end{eqnarray}\n",
    "\n",
    "Esse sistema pode ser simulado como abaixo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dissipative_nonlinear_pendulum(xy, t, g, ℓ, q):\n",
    "    x, y = xy\n",
    "    return y, -g/ℓ*np.sin(x)-q*y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simulate_dissipative(θ0, tmax, nintervals, g, ℓ, q):\n",
    "    θ0_rad = np.deg2rad(θ0)\n",
    "    t = np.linspace(0, tmax, nintervals + 1)\n",
    "    xy0 = [θ0_rad, 0.]\n",
    "    xy = odeint(dissipative_nonlinear_pendulum, xy0, t, args=(g, ℓ, q))\n",
    "    return t, np.rad2deg(xy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixando um valor de $\\theta_0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos a partir de agora fixar $\\theta_0 = 12$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "θ0 = 12.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparando alguns valores de $q$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncycles = 4\n",
    "tmax, ninvervals, _ = timing(ncycles, T0)\n",
    "some_qs = [0.0, 0.5, 1.0, 2.0] # q=0 é o mesmo que sem atrito"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hndl = []\n",
    "xys = []\n",
    "for q in some_qs:\n",
    "    t, xy = simulate_dissipative(θ0, tmax, nintervals, g, ℓ, q)\n",
    "    h, = plt.plot(t, xy[:, 0], label=r'$q=' + str(q) + '$')\n",
    "    hndl.append(h)\n",
    "    xys.append(xy)\n",
    "\n",
    "plt.xlabel(r'$t$')\n",
    "plt.ylabel(r'$\\theta$')\n",
    "plt.legend(handles=hndl)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vemos que as oscilações decaem rapidamente, e mais rapidamente quanto maior o coeficiente de atrito, conforme esperado. Vejamos esses mesmos resultados no espaço de fase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "hndl = []\n",
    "for i, q in enumerate(some_qs):\n",
    "    h, = plt.plot(xys[i][:,0], xys[i][:,1], label=r'$q='+str(q)+'$')\n",
    "    hndl.append(h)\n",
    "plt.xlabel(r'$\\theta$')\n",
    "plt.ylabel(r'$\\dot\\theta$')\n",
    "plt.legend(handles=hndl)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As trajetórias no espaço de fase tendem rapidamente para o ponto $(0,0)$, que é o ponto de energia nula, devido à perda de energia por atrito. Um ponto para o qual as trajetórias tendem é denominado um **ponto fixo**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 4: Força externa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos agora acrescentar um novo elemento: uma força externa atuando sobre o sistema. Para termos um comportamento mais interessante, essa força externa terá um formato senoidal com frequência angular $\\omega_F$ e amplitude $f$ dadas, isto é, terá a forma $ f\\sin(\\omega_F t)$.\n",
    "\n",
    "Nosso sistema fica então:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "  \\frac{dx}{dt} & = & y \\\\\n",
    "  \\frac{dy}{dt} & = & -\\frac{g}{\\ell}\\sin x - q y + f \\sin \\omega_F t.\n",
    "\\end{eqnarray}\n",
    "\n",
    "Para simular este sistema seguimos o mesmo procedimento anterior, mas iremos tomar um cuidado adicional: Como existe uma força externa além da gravidade, é possível que o pêndulo passe acima de seu ponto de origem, podendo inclusive \"rodar\" para o outro lado. Para manter os ângulos dentro de uma faixa específica (para facilitar os gráficos) vamos então continuamente manter os ângulos na faixa $[-180º, 180º)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalize_angle(angle):\n",
    "    return (angle + 180.) % 360. - 180."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def forced_dissipative_nonlinear_pendulum(xy, t, g, ℓ, q, f, ωf):\n",
    "    x, y = xy\n",
    "    return y, -g/ℓ*np.sin(x)-q*y+f*np.sin(ωf*t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simulate_forced(θ0, tmax, nintervals, g, ℓ, q, f, ωf):\n",
    "    θ0_rad = np.deg2rad(θ0)\n",
    "    t = np.linspace(0, tmax, nintervals + 1)\n",
    "    xy0 = [θ0_rad, 0.]\n",
    "    xy = odeint(forced_dissipative_nonlinear_pendulum, xy0, t, args=(g, ℓ, q, f, ωf))\n",
    "    return t, np.rad2deg(xy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos fixar $q=0.5$ e $\\omega_F = 2/3$, e usar $f$ como parâmetro e ajustar o intervalo de $t$. Também, agora calculamos a temporização com base na frequência forçante:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 0.5\n",
    "ωf = 2./3.\n",
    "Tf = 2 * np.pi / ωf\n",
    "ncycles = 10\n",
    "tmax, nintervals, _ = timing(ncycles, Tf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Primeiro vejamos um exemplo para um valor pequeno de $f$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = 0.01\n",
    "t, xy = simulate_forced(θ0, tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "q01 = normalize_angle(xy[:,0])\n",
    "p01, = plt.plot(t, xy[:, 0])\n",
    "plt.xlabel(r'$t$')\n",
    "plt.ylabel(r'$\\theta$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note como o sistema começa oscilando de acordo com a sua frequência natural, mas com o amortecimento que desgasta a energia inicialmente fornecida, passa a ser dominado pela frequência da força externa.\n",
    "\n",
    "Vamos comparar alguns valores diferentes para $f$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uma função auxiliar para plotar resultados para um f específico\n",
    "def subplot_f(f, subindex):\n",
    "    return p, xy\n",
    "\n",
    "fs = [0.0, 0.5, 1.085, 1.2]\n",
    "xys = []\n",
    "for i, f in enumerate(fs):\n",
    "    t, xy = simulate_forced(θ0, tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "    xy[:, 0] = normalize_angle(xy[:,0])\n",
    "    plt.subplot(2, 2, i+1)\n",
    "    plt.title(r'$f='+str(f)+r'$')\n",
    "    plt.plot(t, xy[:, 0])\n",
    "    plt.xlabel(r'$t$')\n",
    "    plt.ylabel(r'$\\theta$')\n",
    "    xys.append(xy)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "É bastante interessante ver estes resultados no espaço de fase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, f in enumerate(fs):\n",
    "    plt.subplot(2, 2, i+1)\n",
    "    plt.title(r'$f=' + str(fs[i])+'$')\n",
    "    plt.plot(xys[i][:, 0], xys[i][:,1], 'o', markersize=1)\n",
    "    plt.xlabel(r'$t$')\n",
    "    plt.ylabel(r'$\\theta$')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aqui plotamos apenas pontos, ao invés de linhas, para evitar linhas inexistentes nos dois gráficos de baixo. O resultado à esquerda acima, usando $f=0$ se reduz aos resultados anteriores, com um ponto fixo na origem. No caso de $f=0.5$ já vemos um comportamento mais interessante: ao invés de um ponto-fixo, a trajetória tende ao que é denominado um ciclo-limite (tende ao ciclo de oscilação determinado pela força externa). Mas para forças maiores, o comportamento é ainda mais interessante. No caso $f=1.085$ o sistema percorre uma maior região do espaço de estado antes de se aproximar do ciclo limite. No caso de $f=1.2$ não vemos tendência a um ponto fixo ou a um ciclo limite. Vamos avaliar o sistema por um tempo maior, para verificar o que ocorre."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncycles = 200\n",
    "tmax, nintervals, _ = timing(ncycles, Tf)\n",
    "t, xy200 = simulate_forced(θ0, tmax, nintervals, g, ℓ, q, 1.2, ωf)\n",
    "f200 = normalize_angle(xy200[:,0])\n",
    "plt.plot(f200, xy200[:,1], '.', markersize=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Não vemos uma repetição da trajetória, ao contrário, uma ampla região do espaço de estado é explorada. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 5: Sensibilidade a condições iniciais"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos agora avaliar agora um outro aspecto do sistema: Vejamos o que acontece quando o sistema se desenvolve a partir de condições iniciais ligeiramente diferentes (continuamos com velocidade inicial zero, mas variamos ligeiramente o ângulo inicial). Avaliaremos então a diferença (absoluta) entre as duas trajetórias com a passagem do tempo.\n",
    "\n",
    "Primeiramente, para $f=0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "θ0 = 12. # Ângulo inicial\n",
    "Δθ0 = 0.2\n",
    "θ0prime = θ0 + Δθ0 # Ângulo inicial ligeiramente modificado\n",
    "\n",
    "ncycles = 10\n",
    "tmax, nintervals, _ = timing(ncycles, Tf)\n",
    "\n",
    "f = 0.1\n",
    "\n",
    "t, xy_0 = simulate_forced(θ0,      tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "t, xy_1 = simulate_forced(θ0prime, tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "\n",
    "plt.plot(t, np.abs(normalize_angle(xy_0[:, 0] - xy_1[:, 0])))\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$t$')\n",
    "plt.ylabel(r'Diferença')\n",
    "plt.title(r'$f=' + str(f) + r'$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note como a diferença entre as trajetórias decai rapidamente. Na verdade, usamos a escala logaritmica no eixo das diferenças para mostrar que a queda tem tendência exponencial. Isto é,\n",
    "\n",
    "$$ | \\Delta\\theta | \\approx e^{\\lambda t}, $$\n",
    "\n",
    "onde $\\lambda$ é negativo neste caso.\n",
    "\n",
    "Agora vejamos o que acontece para $f=1.2$. (Para mostrar melhor o comportamento, começamos com um valor bem menor de diferença nos ângulos iniciais.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "θ0 = 12.\n",
    "Δθ0 = 0.00000001\n",
    "θ0prime = θ0 + Δθ0\n",
    "\n",
    "ncycles = 30\n",
    "tmax, nintervals, _ = timing(ncycles, Tf)\n",
    "\n",
    "f = 1.2\n",
    "\n",
    "t, xy_0 = simulate_forced(θ0,      tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "t, xy_1 = simulate_forced(θ0prime, tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "\n",
    "plt.plot(t, np.abs(normalize_angle(xy_0[:, 0] - xy_1[:, 0])))\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$t$')\n",
    "plt.ylabel(r'Diferença')\n",
    "plt.title(r'$f=' + str(f) + r'$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vemos que neste caso o comportamento é oposto: a diferença entre as trajetórias continua seguindo uma tendência da forma\n",
    "\n",
    "$$ |\\Delta\\theta| \\approx e^{\\lambda t}, $$\n",
    "\n",
    "mas agora o coeficiente $\\lambda$ é **positivo**! Isto significa que **qualquer pequena diferença nas condições inciais será amplificada rapidamente**, resultando em que o comportamento do sistema é na prática impossível de prever, visto que não é possível uma precisão perfeita na determinação das condições iniciais. Esta é uma característica dos sistemas denominados **caóticos**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 6: Seção de Poincaré"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uma técnica útil na análise de sistemas dinâmicos é a denominada *seção de Poincaré*, ou *mapa de Poincaré*, que é essencialmente a interseção de uma trajetória no espaço de fase com um espaço de menor dimensão (um subespaço) específico.\n",
    "\n",
    "O subespaço a escolher depende do sistema e da análise desejada. No nosso caso, vamos plotar apenas os pontos que correspondem a instantes que satisfazem a relação $\\omega_F t = 2n\\pi$ onde $n$ é um inteiro. Isto é, estamos plotando apenas os pontos onde a trajetória está em fase com a força externa aplicada. Devido às aproximações numéricas, devemos na verdade plotar um ponto se ele satisfizer $|t - 2n\\pi/\\omega_F| < \\Delta t/2$ para algum valor de $\\Delta t$ apropriado. No nosso caso, usaremos para $\\Delta t$ o valor do intervalo entre amostras no sistema simulado, visto que não há precisão temporal maior do que essa em nossos resultados.\n",
    "\n",
    "Vamos fazer uma função que dado um array de valores de tempo, o tempo de ciclo e a precisão de tempo $\\Delta t$ seleciona todos os índices (retornando um vetor de booleanos) que estão sincronizados com o ciclo de acordo com a expressão acima:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def select_synchronized(t, Tf, Δt):\n",
    "    return np.abs(np.mod(t+Tf/2, Tf) - Tf/2) < Δt/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 0.5\n",
    "ncycles = 200\n",
    "tmax, nintervals, Δt = timing(ncycles, Tf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = 0.1\n",
    "\n",
    "t, xy = simulate_forced(θ0, tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "traj = normalize_angle(xy[:, 0])\n",
    "\n",
    "selected = select_synchronized(t, Tf, Δt)\n",
    "\n",
    "plt.plot(traj[selected], xy[selected, 1], 'o', markersize=3, alpha=0.3) # Plota apenas esses pontos\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note como, com exceção de alguns pontos devidos ao transiente inicial, os pontos agora se acumulam em um mesmo local. A razão é simples: depois do transiente inicial, o sistema assume uma trajetória periódica de acordo com a frequência da força externa, e portanto, se olhamos a intervalos determinados pelo período dessa oscilação, encontraremos o pêndulo sempre no mesmo ponto.\n",
    "\n",
    "Vejamos agora o que acontece no caso caótico."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = 1.2\n",
    "\n",
    "t, xy = simulate_forced(θ0, tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "traj = normalize_angle(xy[:, 0])\n",
    "\n",
    "selected = select_synchronized(t, Tf, Δt)\n",
    "\n",
    "plt.plot(traj[selected], xy[selected, 1], '.', markersize=3, alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O que acontece se deixamos o sistema executar por mais tempo?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncycles = 2000\n",
    "tmax, nintervals, Δt = timing(ncycles, Tf)\n",
    "\n",
    "t, xy = simulate_forced(θ0, tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "\n",
    "traj = normalize_angle(xy[:, 0])\n",
    "selected = select_synchronized(t, Tf, Δt)\n",
    "\n",
    "plt.plot(traj[selected], xy[selected, 1], '.', markersize=3, alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note como o formato ficou o mesmo. Quer dizer, apesar de toda a característica caótica do sistema, ele é atraído a uma região específica do espaço de estado, denominado um **atrator**. A estrutura dos atratores no regime caótico é bastante complexa, tendo dimensionalidades fracionárias (o número de dimensões do espaço do atrator não é um número inteiro), e portanto são denominados **fractais**; os atratores fractais são denominados de **atratores estranhos**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 7: Duplicação de período"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos analisar o regime estacionário para diferentes parâmetros de amplitude da força externa ($f$). Para isso, simulamos o sistema para um grande número de ciclos e descartamos todos os resultados com exceção dos últimos ciclos, que mostrarão o funcionamento do sistema após o término do transiente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncycles = 212 # Número total de ciclos a simular\n",
    "samples_per_cycle = 100 # Número de pontos amostrados por ciclo\n",
    "tmax, nintervals, _ = timing(ncycles, Tf, samples_per_cycle)\n",
    "ncycles_discard = 200 # Número de ciclos a descartar (evitando o transiente)\n",
    "ndiscard = ncycles_discard * samples_per_cycle # Número de amostras descartadas\n",
    "\n",
    "plt.figure(figsize=(9,9))\n",
    "fs = [1.065, 1.075, 1.082, 1.0826] # Alguns valores de f\n",
    "for i, f in enumerate(fs):\n",
    "    t, xy = simulate_forced(θ0, tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "    traj = normalize_angle(xy[:, 0])\n",
    "    \n",
    "     # Seleciona os valores não descartados\n",
    "    to_keep = np.arange(nintervals+1) > ndiscard # Calcula os indices a manter\n",
    "    traj_steady = traj[to_keep]\n",
    "    t_steady = t[to_keep]\n",
    "    \n",
    "    # plota os gráficos\n",
    "    plt.subplot(len(fs), 1, i+1)\n",
    "    plt.plot(t_steady, traj_steady)\n",
    "    plt.ylabel(r'$\\theta$')\n",
    "    plt.title(r'$f='+str(f)+r'$')\n",
    "plt.xlabel(r'$t$')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note que para $f=1.065$, no regime estacionário a oscilação tem um período igual ao da força externa (no caso dos parâmetros usados, $T_F = 3\\pi$).\n",
    "\n",
    "Já para $f=1.075$, os picos oscilam alternadamente de altura, o que significa que o período de repetição é dobrado. Para $f=1.082$ temos 4 picos de valores diferentes antes da repetição, o que implica um período quadruplicado em relação ao original. Para $f=1.0826$ temos 8 picos distintos, apesar de isso ser muito difícil de enxergar na figura. Esse processo de duplicação se repete indefinidamente, com intevalos entre os valores de $f$ para a duplicação cada vez menores, até que chegamos no regime caótico. Essa **duplicação de período** é típica em vários sistemas caóticos."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Diagrama de bifurcação"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uma forma tradicional de analisar o processo de duplicação de período na transição para o caos é através da construção de um **diagrama de bifurcação** para analisarmos os pontos onde a duplicação de período ocorre.\n",
    "\n",
    "Para o nosso sistema usamos a mesma técnica da seção de Poincaré para extrair os pontos em fase com força externa e plotamos, para cada valor de $f$, todos os valores encontrados nesses pontos dentro do regime estacionário.\n",
    "\n",
    "*OBS:* A execução do código abaixo é demorada, pois realiza um grande número de simulações."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Faixa de valores de f a simular\n",
    "fs = np.linspace(1.065,  1.087, 200+1)\n",
    "\n",
    "# Parâmetros da simulação\n",
    "ncycles = 250\n",
    "ncycles_discard = 200\n",
    "samples_per_cycle = 100\n",
    "tmax, nintervals, Δt = timing(ncycles, Tf, samples_per_cycle)\n",
    "ndiscard = ncycles_discard * samples_per_cycle\n",
    "\n",
    "plt.figure(figsize=(9,9))\n",
    "for f in fs:\n",
    "    t, xy = simulate_forced(θ0, tmax, nintervals, g, ℓ, q, f, ωf)\n",
    "    traj = normalize_angle(xy[:, 0])\n",
    "    selected = select_synchronized(t, Tf, Δt) # Selecionados pela seção de Poincaré\n",
    "    traj_steady = traj[selected & (np.arange(nintervals+1) > ndiscard)] # Pega apenas os não-descartados\n",
    "    # Plota todos os pontos no mesmo valor de f\n",
    "    plt.plot(f * np.ones_like(traj_steady), traj_steady, 'k.', markersize=1)\n",
    "plt.minorticks_on()\n",
    "plt.grid(which='both', color='r', alpha=0.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note como até um pouco acima de $f=1.066$ a trajetória atinge sempre o mesmo ponto (para um dado $f$). Um pouco acima disso já são atingidos dois pontos distintos. Isso é uma marca de que o período de oscilação é o dobro do período de amostragem. Por volta de $f=1.079$ há uma nova duplicação (cada um dos pontos anteriores é dividido em dois). Próximo a $f=1.082$ há outra duplicação. Esse processo em princípio se repete indefinidamente, com o intervalo entre um ponto de duplicação e outro diminuindo exponencialmente, até que se atinge o ponto de caos.\n",
    "\n",
    "Vale a pena você fazer mudanças na faixa de $f$ usada para entender melhor o sistema."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}