{"nbformat":4,"nbformat_minor":0,"metadata":{"anaconda-cloud":{},"colab":{"name":"mac0209-movimento-2D-stateSpace.ipynb","provenance":[],"collapsed_sections":[]},"kernelspec":{"display_name":"Python [default]","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.5.2"},"toc":{"base_numbering":1,"nav_menu":{},"number_sections":true,"sideBar":true,"skip_h1_title":false,"title_cell":"Table of Contents","title_sidebar":"Contents","toc_cell":false,"toc_position":{},"toc_section_display":true,"toc_window_display":true}},"cells":[{"cell_type":"markdown","metadata":{"id":"JAcSMzHxSxR5","colab_type":"text"},"source":[" MAC0209 - Modelagem e Simulação \n","\n","Roberto M. Cesar Jr. (IME-USP)\n","\n","Roberto Hirata Jr. (IME-USP)\n","***\n"," Movimento 2D usando Sistemas dinâmicos, State Vectors e State Space\n","***"]},{"cell_type":"markdown","metadata":{"id":"uEcmOVX9SxR7","colab_type":"text"},"source":["# Introdução: modelagem por sistemas dinâmicos e vetores de estado"]},{"cell_type":"markdown","metadata":{"id":"MhfLFFY0SxR8","colab_type":"text"},"source":["Equações diferenciais do movimento:\n","\n","$$\\frac{d^2\\vec{x}}{dt^2} = \\frac{d\\vec{v}}{dt} = \\vec{a}(t)$$\n","\n","$$\\frac{d\\vec{x}}{dt} = \\vec{v}(t)$$\n","\n","Euler:\n","\n","$$\\vec{x}(t+\\Delta t) = \\vec{x}(t) + \\vec{v}(t) \\; \\Delta t $$\n","\n","$$\\vec{v}(t+\\Delta t) = \\vec{v}(t) + \\vec{a}(t) \\; \\Delta t $$\n","\n","Assim, o movimento 2D da partícula pode ser representado pelo vetor de estados \n","\n","$\\vec{s} = (\\vec{x},\\vec{v}) = ([x_1, x_2], [v_1, v_2])$\n"," \n"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"zqAj0xhX4BVu","colab":{}},"source":["# funcoes genericas que podem ser re-usadas em outros problemas\n","\n","import math\n","import matplotlib.pyplot as pyplot\n","import numpy as np\n","import matplotlib as mpl\n","from mpl_toolkits.mplot3d import Axes3D\n","\n","\n","# funcoes base para implementar o Euler. \n","# Deve-se implementar a funcao rates, que depende de cada modelo.\n","\n","def initStateVector(s):\n"," return(np.array(s))\n","\n","def updateStateVectorEuler(s,dt):\n"," return(s + rates(s,dt))\n","\n","# State Vector Trajectories store state space evolution. Uses list to init empty.\n","\n","def initSVTrajectory():\n"," return([])\n","\n","# append s a svt\n","def updateSVTrajectory(svt,s):\n"," svt.append(s)\n"," return(svt)\n","\n","def extractSVTrajectory(svt,i): # returns the trajectory as numpy array\n"," foo = np.array(svt)\n"," return(foo[:,i])\n"," \n","def plotCompareAnalyticalEuler(vxa, vva, vxe, vve, vtime):\n","\n"," fig, ax = pyplot.subplots()\n"," pyplot.plot(vtime, vxe, label='Euler',linestyle='',marker='o') \n"," pyplot.plot(vtime, vxa, label='Analytical') \n"," pyplot.title('Posição')\n"," ax.set_xlabel('Tempo (segundos)')\n"," ax.set_ylabel('Posição (metros)')\n"," pyplot.show(block=False)\n"," \n"," fig, ax = pyplot.subplots()\n"," pyplot.plot(vtime, vve,label='Euler',linestyle='',marker='o')\n"," pyplot.plot(vtime, vva,label='Analytical')\n"," pyplot.title('Velocidade')\n"," ax.set_xlabel('Tempo (segundos)')\n"," ax.set_ylabel('Velocidade (metros / segundo)')\n"," pyplot.show()\n","\n"," fig, ax = pyplot.subplots()\n"," pyplot.plot(vxe, vve,label='Euler',linestyle='',marker='o')\n"," pyplot.plot(vxa, vva,label='Analytical')\n"," pyplot.title('Dynamical System Trajectory')\n"," ax.set_xlabel('Posição (metros)')\n"," ax.set_ylabel('Velocidade (metros / segundo)')\n"," pyplot.show()\n","\n","\n","def erroTrajetorias(v1,v2,tipoErro):\n"," if (tipoErro == 0): # erro com sinal\n"," return(np.array(v1) - np.array(v2))\n"," elif (tipoErro == 1): # erro quadratico\n"," return((np.array(v1) - np.array(v2))**2)\n"," elif (tipoErro == 2): # erro em modulo\n"," return(fabs((np.array(v1) - np.array(v2))))\n"," \n","\n","def easyPlot(v,title):\n"," pyplot.figure()\n"," pyplot.plot(v)\n"," pyplot.title(title)\n"," pyplot.show()\n","\n","def easyPlot2D(x,y,title):\n"," pyplot.figure()\n"," pyplot.plot(x,y,'*')\n"," pyplot.title(title)\n"," pyplot.show()\n"," \n","def easyPlot3D(x,y,z,title,xl,yl,zl):\n"," mpl.rcParams['legend.fontsize'] = 10\n"," fig = pyplot.figure()\n"," ax = fig.gca(projection='3d')\n"," ax.plot(x, y, z, label=title)\n"," ax.set_xlabel(xl)\n"," ax.set_ylabel(yl)\n"," ax.set_zlabel(zl)\n"," ax.legend()\n"," pyplot.show()\n"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"HNEmHeOp5lab","colab_type":"code","colab":{}},"source":[""],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"ZzsRkUuw4BV1"},"source":["## Exercício: $\\frac{d^2\\vec{x}}{dt^2} = \\vec{a}(t)$\n","\n","Escreva a solução para as EDOs:\n","\n","$\\frac{d^2x_0}{dt^2} = \\sin (k_0 \\pi t)$\n","\n","$\\frac{d^2x_1}{dt^2} = \\cos (k_1 \\pi t)$\n","\n","\n"]},{"cell_type":"markdown","metadata":{"id":"H5aEXF0MSxSc","colab_type":"text"},"source":["### Exercício: Resolva na célula abaixo antes de olhar a solução"]},{"cell_type":"markdown","metadata":{"id":"NOsdjhKiSxSd","colab_type":"text"},"source":["O programa abaixo implementa a solucão desse problema **com** a modelagem por sistemas dinâmicos e vetores de estado."]},{"cell_type":"code","metadata":{"code_folding":[],"id":"Odk1OQRiSxSd","colab_type":"code","colab":{}},"source":["# Escreva sua solucao aqui."],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"n4iKEtZDSxSf","colab_type":"text"},"source":["### Resolução: Compare sua solução"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"xDtL5OZS4BV6","colab":{}},"source":["# Implementa o exercício da integração de Euler de \n","# compara com o resultado analítico\n","\n","# Euler:\n","def rates(s,dt):\n"," k0 = 2\n"," k1 = 3\n"," r0 = s[2] * dt\n"," r1 = s[3] * dt\n"," r2 = math.sin(k0 * math.pi * s[4]) * dt\n"," r3 = math.cos(k1 * math.pi * s[4]) * dt\n"," r4 = dt\n"," return(np.array([r0,r1,r2,r3,r4]))\n","\n","def main():\n"," t=0\n"," tf = 4\n"," dt=0.05\n"," x0=0\n"," x1=0\n"," v0=0\n"," v1=0\n"," \n"," # state vector: [position x, velocity v, time t]\n"," stateVectorEuler = initStateVector([x0,x1,v0,v1,t])\n"," \n"," svtEuler = initSVTrajectory() \n"," \n"," while (stateVectorEuler[4]