{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\tHiroshi TAKEMOTO\n", "\t(take.pwave@gmail.com)\n", "\t\n", "\t

微分方程式

\n", "\t

\n", "\t\tSageを使って微分方程式を解く方法を紹介します。微分方程式の性質から解析的な解が求まらない場合等に備え、\n", "\t\t数値的に微分方程式を方法についても取り上げます。\n", "\t

\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

desolveを使った微分方程式の解法

\n", "\t

\n", "\t\t最初にSageで微分方程式を解くときに一般的に使用されるdesolve関数の使い方について説明します。\t\t\n", "\t

\n", "\t

\n", "\t\t例として以下のような2次の微分方程式をSageで表してみます。\n", "$$\n", "\t\t\\ddot z + \\dot z + z = u(t)\n", "$$\n", "\t\tここで外部からの変位u(t)を、$u(t) = 1, t \\ge 0$のステップ関数とすると上記の微分方程式は\n", "\t\t以下のように定義されます。\n", "\t

\n", "\t

\n", "\t\tdesolveでは、変数tをvar関数で、変数tの関数zをfunction関数で定義し、\n", "\t\t微分方程式deに$\\ddot z + \\dot z + a = 1$の関係式をセットし、\n", "\t\t初期条件icsに$[ t_0, z(t_0), \\dot z(t_0)]$を与えて解きます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/usr/lib/sagemath/local/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2885: DeprecationWarning: Calling function('f',x) is deprecated. Use function('f')(x) instead.\n", "See http://trac.sagemath.org/17447 for details.\n", " exec(code_obj, self.user_global_ns, self.user_ns)\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "-1/3*(sqrt(3)*sin(1/2*sqrt(3)*t) + 3*cos(1/2*sqrt(3)*t))*e^(-1/2*t) + 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "t = var('t')\n", "# deprecateのエラーがでますが、このままで実行します。\n", "z(t) = function('z', t)\n", "de = diff(z, t, 2) + diff(z, t) + z == 1\n", "sol = desolve(de, z, ics=[0,0,0])\n", "show(sol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t求まった解$z(t)$とその微分$\\dot z(t)$をプロットしてみます。\n", "\t

\n", "\t

\n", "\t\t$z(t)$は、最初は振動していますが、t=10になると定常状態に収まって1となっています。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFnCAYAAACRo/HLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl8VNX9//H3hAQSQMIedgkoi6KAbBGpBUERlFiLYqIU\ngm2FgkWLisX2p1JtpYvlqxVrtcWVzaJFVBA1giCyiYE2EBKEKItmJCJJCEhCcn9/HEMykGUmmZl7\nZ/J6Ph7zGHK5c+8n63vOueee47IsyxIAAAi6CLsLAACgviKEAQCwCSEMAIBNCGEAAGxCCAMAYBNC\nGAAAmxDCAADYhBAGAMAmjgxhy7KUn58v5hEBAIQzR4ZwQUGBYmNjVVBQYHcpAAAEjCNDGACA+oAQ\nBgDAJoQwAAA2IYQBALBJpN0FoG4sSzp6VDpyRIqNlVq0kKKj7a4KAOANQjgElZRI69dLixdLr79u\nQriiSy+Vxo0zj0GDpAj6OwDAkVyWA2/Gzc/PV2xsrPLy8tSsWTO7y3GUtDQpJUX673+lrl2l5GRp\n4ECpTRspP19yu6W1a6VVq0w49+kjzZ0r/ehHhDEAOA0t4RBRVCQ9+qj0hz+YYF2/Xho2THK5zt33\n9tul06fNPvPmSePHS/37SwsWSJdfHvzaAQCVo20UAo4fl8aOlR57TPp//0/aulX6wQ8qD+AykZHS\nVVdJ775rwjgyUrriCunuu6VTp4JXOwCgao7ujh4zZowiIyOVnJys5ORku8uyxfHj0tVXS7t2SW+9\nJV15Ze2OU1IiPfGENGeO1Lev9O9/S+ef799aAQC+cXQI1/drwiUl5lruhx9KH3xgrv3W1SefSDfd\nJBUUSK+8Io0ZU/djAgBqh+5oB7vvPjPA6tVX/RPAkjnOp59KCQnSddeZ1jEAwB6EsEMtXy7Nn29C\n8tpr/Xvsli2lN9+U7r3XXCOePdvcbwwACC5GRztQTo70859LN98szZgRmHNEREh/+pPUsaMJ4pMn\nTeBzGxMABA8h7ECzZ0sNGkh//3v1I6D94a67pJgYaepUM2r6mWcIYgAIFkLYYdavl15+WfrnP6VW\nrYJzzjvukBo2NPcXFxVJ//qXeRMAAAgsQthBioul6dPNoKkpU4J77pQUE8STJpkW8UsvSVFRwa0B\nAOobQthB/vY3KSPD3EZkR5fwrbdKjRpJSUmmG/zll2kRA0AgcfXPIQoKzLSUd9xhppi0y/jx0pIl\n0rJl0p13MmoaAAKJEHaIZ54xs2M98IDdlZjJPJ57ztTkhHoAIFzRHe0AJ09Kjz8uTZ4sde5sdzXG\n7bdLeXnSrFlmneJf/9ruigAg/BDCDrBwoXTkiHT//XZX4ulXv5KOHTPzTTdvLk2bZndFABBeHB3C\nSUlJYb+AQ3GxmTQjKUm64AK7qznXww+bIJ4+XWrWzAzeAgD4h6NDeOnSpWG/gMOiRdKBA6a16UQu\nl5k+My/PdJc3b26WVQQA1B2rKNnIsqR+/cySgitX2l1N9U6fNgO21qyR3ntPGjbM7ooAIPQxOtpG\nH38s/fe/5lYgp4uMlJYuNROJXH+9tHOn3RUBQOgjhG309NPmOvCoUXZX4p3oaOmNN0zNo0dLn31m\nd0UAENoIYZscOya99ppZLSmUFkxo1kxavdpcG776aunwYbsrAoDQFUJ//sPL8uVmZPTEiXZX4rs2\nbcx14ZIS0yI+etTuigAgNBHCNnn5ZWnkSKlDB7srqZ3OnU0Qu91mtPTx43ZXBAChhxC2weefmyUL\nJ02yu5K66dlTeucdafdu6cc/NqsvAQC8Rwjb4JVXpCZNpBtvtLuSuhswwNxetWGDWfyBIAYA7xHC\nQWZZpit6/HgTxOFg+HBpxQrp/felCROkoiK7KwKA0EAIB9nWrVJWlvSTn9hdiX+NHi29/rrpnk5K\nMoPOAADVI4SD7JVXpI4dpREj7K7E/8aONbddvfWWmWOaIAaA6jk6hJOSkpSYmKglS5bYXYpflJaa\n1uLNN0sNGthdTWBcf73073+b7umJE810lwCAyjF3dBBt3SoNGSJ9+KF05ZV2VxNYr78u3XKLlJgo\nLV4sNWpkd0UA4DyObgmHmxUrpNatpaFD7a4k8H78YxPEb78tXXedlJ9vd0UA4DyEcBCtWCGNG2cW\nQ6gPxo0zqy5t22ZGULvddlcEAM5CCAdJZqaUkSH96Ed2VxJcP/yhuYc4J8f0ALDoAwCUI4SDZMUK\nqXFjs+hBfXPppWbZxshIE8Qff2x3RQDgDIRwkKxYYe6ljYmxuxJ7dO0qbdwo9ehhuqb//nczcQkA\n1GdBCeENGzYoMTFRHTt2VEREhFauXBmM0zrGV19JmzeHxzSVddG6tfTBB9K0adL06dLtt0snT9pd\nFQDYJyghXFhYqH79+unpp5+Wy+UKxikd5Z13JJdLGjPG7krs17Ch9OST0ksvSUuXSsOGSV98YXdV\nAGCPoIzTvfbaa3XttddKkhx4W3LArVkjDRxoWoIwfvITqU8fcyvTgAHSwoXmnmIAqE+4JhxgJSVm\n3d3Ro+2uxHn695e2bzeDtW64wcywdfSo3VUBQPDUkztW7bN9uwkWQrhyLVtKb7xh5tSeOdOsxPSP\nf5hQRu2dOiV9+6352fv223P/XVhorsefOOH5+O47M9Xo6dPmDWRlD5fLTLvaoIEUEVH+78o+btjQ\nzJZW8VHZNm+3N2xY8yNcp4RFeCKEA2zNGqlZMzNdJSrncpnu6VGjpKlTzb3UN9wgzZsn9epld3XO\nUlpqBvodPCh9+WX54/Bhz4+PHav89TExUosWUtOm5t+NG5tHTIzUqpUUHS1FRXkGaWSk58eW5RnK\npaXnBnVpqQnyoiLzhqCw0LwJOHWqfNvZj7LtJSV1+xpFRHgX1mWfZ8U3DmX/Pvu5pm0VnX3FzV8f\nW1b1//Z2v2Afy1/HrU5VQ40Cuf2OO8xCNXXl6BBOSkpS5FnTSyUnJys5Odmminy3Zo00cqT5hUf1\n2rc3reJXX5V+/WtzzfjnP5ceekhq187u6oKntFTKzjYTvOzb5/nIzjat1TING0odOpQ/Lr7YPMfF\nmVBt0cLzER1t3+flrZKSysO5qKjmR3Gxd/uVPSq+gSj7d8U3EJX9X2VvPir7I13bbdXt43JV/29v\n9wv2sfx13MpUFdCB3t60aeXbfRX0BRwiIiK0YsUKJVYzCidcFnDIyzN/CBcsMC08eO/UKfN1e/RR\n88fwrrukGTNMwIQLyzIt2vR0adeu8ufdu8tv3WrYUIqPl7p393x06WK+Fq1aVf3HCYDzBaUlXFhY\nqM8+++zMyOj9+/dr586datmypTp37hyMEmyRmmreJXM92HeNGkmzZklTpkiPPWZua/rTn6SbbjLX\njhMSQid8LMtM21kWtGVhu2uXVFBg9mnSxLRi+/Y1XVx9+piu+I4ducYJhLOgtIQ//PBDjRgx4px7\nhCdPnqyFCxees3+4tIRnzJDefVfau9fuSkJffr70wgvS3/5m5p/u398slTh+vHTBBXZXVy4317NV\nW/ZcNuo7Olrq3duE7MUXlz936WKuLwKoX1hPOIAuvli64grp2WftriR8lJZKq1dLzz8vrVplum37\n9jXLJQ4bJl1+udS8eeBrOHRI2rPH87F7d/lKUVFRUs+e5UFbFrbdutGyBVCOEA4Qt9sMJlq0yD8j\n6HCuwkIzG9lrr5mu/6+/Nl3UffpIl11mAi8+3jx362a+HzV1YZeUmFbrkSPlj8OHpc8/N4/sbNOz\nceKE2b9RIzMfdq9e5lEWuBdeyGA8ADUjhANk2TIpKcncLtK+vd3VhD/LMt3UGzeaR3q6CcyKaxi7\nXObaa9kjOtqMgq04AvfYMdPSrSg62ixA0bWrdP75poVbFrpdutCyBVB7hHCATJsmrVtnuilhn+PH\nTQt2/35zf21hoWnFFhaaW32iojzvHW3RQmrbVmrTpvzRokXoDAIDEFocfZ9wKFu7VrrqKrurQNOm\n5V3EAOA0jMcMgC+/lLKyzLq5AABUhRAOgPXrzfMPf2hvHQAAZyOEA2DjRnPvan2aahEA4DtCOAA+\n+sjcswoAQHUcHcJJSUlKTEzUkiVL7C7FawUF0n//aybpAACgOo4eHb106dKQu0Vp82ZznyktYQBA\nTRzdEg5FH31kVrbp2dPuSgAATkcI+9nGjaYrmskdAAA1IYT96PRp0x3N9WAAgDcIYT/audNMh8j1\nYACANwhhP/roI7OqzoABdlcCAAgFhLAfbd5sArhRI7srAQCEAkLYj7ZskYYMsbsKAECoIIT95MgR\ns37t4MF2VwIACBWEsJ9s22aeCWEAgLcIYT/ZskVq3VqKj7e7EgBAqCCE/WTrVtMKZpIOAIC3HB3C\nobKAg2WZEGZQFgDAFyzg4Af79klHj3I9GADgG0e3hEPFli3medAge+sAAIQWQtgPtm6VLrjArJ4E\nAIC3CGE/KBuUBQCALwjhOioqktLSGJQFAPAdIVxH6enSqVPSwIF2VwIACDWEcB3t2GHuDb70Ursr\nAQCEGkK4jtLSpB49pKZN7a4EABBqCOE6SkuT+ve3uwoAQCgihOugtFTauZMQBgDUDiFcB599Jh0/\nTggDAGqHEK6DtDTzTAgDAGrD0SHs9AUc0tKkTp3MEoYAAPjKZVmWZXcRZ8vPz1dsbKzy8vIcvYDD\n6NFSo0bSypV2VwIACEWObgk7mWUxMhoAUDeEcC19+aV05AghDACoPUK4lhiUBQCoK0K4ltLSpBYt\npC5d7K4EABCqCOFaSkuT+vUz80YDAFAbtQrhBQsWKD4+XjExMUpISNC2bduq3X/RokXq16+fmjRp\nog4dOuinP/2pjh49WquCnYJBWQCAuvI5hJctW6Z77rlHc+fOVVpamvr27avRo0crNze30v03btyo\nyZMn6+c//7l2796t5cuXa+vWrbrjjjvqXLxdjh2TPv+cEAYA1I3PITx//nxNnTpVkyZNUq9evfTM\nM8+ocePGWrhwYaX7b968WfHx8ZoxY4bOP/98DR06VFOnTtXWrVvrXLxdduwwz4QwAKAufArh4uJi\nbd++XSNHjjyzzeVyadSoUdq0aVOlr7n88st18OBBrV69WpLkdru1fPlyXXfddXUo215paVJ0tNSz\np92VAABCmU8hnJubq5KSEsXFxXlsj4uLU05OTqWvGTp0qF555RXdcsstatiwodq3b6/mzZvrqaee\nqn3VNktLky69VIqMtLsSAEAoC/jo6N27d+uuu+7Sww8/rE8//VRr1qxRdna2pk6dGuhTBwyDsgAA\n/uBTW65169Zq0KCB3G63x3a326127dpV+pp58+Zp2LBhmjVrliSpT58+evrpp/WDH/xAv//9789p\nVVeUlJSkyLOam8nJyUpOTvalbL86eVLKyJDuvNO2EgAAYcKnEI6KitKAAQOUmpqqxMRESZJlWUpN\nTdXMmTMrfc2JEycUFRXlsS0iIkIul0s1rR2xdOlSxy3gsGuXVFJi7hEGAKAufO6OnjVrlp577jm9\n9NJL2rNnj6ZNm6YTJ04oJSVFkjRnzhxNnjz5zP7jxo3T66+/rmeeeUbZ2dnauHGj7rrrLg0ZMqTK\n1rOTpaeb5z597K0DABD6fB5aNGHCBOXm5urBBx+U2+1Wv379tGbNGrVp00aSlJOTo4MHD57Zf/Lk\nyTp+/LgWLFige++9V82bN9fIkSM1b948/30WQZSeLnXrJjVpYnclAIBQx3rCPrr2WrOG8Btv2F0J\nACDUMXe0j9LT6YoGAPgHIeyDb7+VDh+WLrnE7koAAOGAEPYBg7IAAP5ECPsgPd3MktWjh92VAADC\nASHsg/R0M190w4Z2VwIACAeEsA8YlAUA8CdC2EuWRQgDAPyLEPZSTo509CghDADwH0eHcFJSkhIT\nE7VkyRK7S2FkNADA7xy9Iq6TFnBIT5diYqT4eLsrAQCEC0e3hJ0kPV266CKpQQO7KwEAhAtC2Evp\n6cyUBQDwL0LYC6WlZh1hrgcDAPyJEPbCF19IhYWEMADAvwhhLzAyGgAQCISwF9LTpebNpQ4d7K4E\nABBOCGEvlM2U5XLZXQkAIJwQwl743//oigYA+B8hXIPiYmnPHkIYAOB/hHAN9u41QUwIAwD8jRCu\nwa5d5vnii+2tAwAQfhwdwk5YwCEjQ2rTRmrd2rYSAABhigUcarB7t9S7t60lAADClKNbwk6QkWEW\nbgAAwN8I4WqUlEiZmbSEAQCBQQhXIztbOnWKEAYABAYhXI2MDPNMCAMAAoEQrkZGhnTeeVLHjnZX\nAgAIR4RwNTIyTCuYOaMBAIFACFeD25MAAIFECFfBsrg9CQAQWIRwFb78UioooCUMAAgcQrgKu3eb\nZ0IYABAohHAVMjKkRo2k+Hi7KwEAhCtHh7CdCzhkZEg9e0oNGgT91ACAeoIFHKpQdnsSAACB4uiW\nsJ0IYQBAoBHClTh2TPr6a9MdDQBAoBDClcjMNM+EMAAgkAjhSpSFcI8e9tYBAAhvhHAlMjOlTp2k\nJk3srgQAEM5qFcILFixQfHy8YmJilJCQoG3btlW7f1FRkX7zm9+oa9euio6OVrdu3fTCCy/U5tRB\nkZVFKxgAEHg+36K0bNky3XPPPXr22Wc1ePBgzZ8/X6NHj1ZWVpZat25d6WtuvvlmHTlyRM8//7y6\nd++ur776SqWlpXUuPlAyM6Vhw+yuAgAQ7lyWZVm+vCAhIUFDhgzRE088IUmyLEudO3fWzJkzNXv2\n7HP2f+edd3Trrbdq//79at68uVfnyM/PV2xsrPLy8oJ+n3BpqemGnjdPuuuuoJ4aAFDP+NQdXVxc\nrO3bt2vkyJFntrlcLo0aNUqbNm2q9DVvvvmmBg4cqD/+8Y/q1KmTevbsqfvuu0/fffdd3SoPkAMH\npO++Y2Q0ACDwfOqOzs3NVUlJieLi4jy2x8XFKbNsSPFZ9u/frw0bNig6OlorVqxQbm6ufvGLX+jo\n0aP617/+VfvKA4TbkwAAwRLwaStLS0sVERGhxYsXq2nTppKkv/71r7r55pv19NNPq1GjRoEuwSeZ\nmWbhhi5d7K4EABDufArh1q1bq0GDBnK73R7b3W632rVrV+lr2rdvr44dO54JYEnq3bu3LMvSoUOH\n1L179yrPl5SUpMhIzxKTk5OVnJzsS9k+ycqSLriAhRsAAIHnUwhHRUVpwIABSk1NVWJioiQzMCs1\nNVUzZ86s9DVXXHGFli9frhMnTqhx48aSpMzMTEVERKhTp07Vns+OBRwyM+mKBgAEh8/3Cc+aNUvP\nPfecXnrpJe3Zs0fTpk3TiRMnlJKSIkmaM2eOJk+efGb/W2+9Va1atdKUKVOUkZGh9evXa/bs2frp\nT3/quK5oiRAGAASPz9eEJ0yYoNzcXD344INyu93q16+f1qxZozZt2kiScnJydPDgwTP7N2nSRO+9\n955++ctfatCgQWrVqpVuueUWPfLII/77LPyksFA6eJAQBgAEh8/3CQeDXfcJ79gh9e8vffyxdPnl\nQTstAKCeYu7oCrKyzDMtYQBAMBDCFWRmSq1bSy1b2l0JAKA+IIQrYFAWACCYCOEKCGEAQDARwt+z\nLBPCLGEIAAgWQvh7brdUUEBLGAAQPITw91i4AQAQbITw9zIzzXzR1UxlDQCAXzk6hJOSkpSYmKgl\nS5b49bjZ2dLhw57bMjOl+HipYUO/ngoAgCrVuxmzCgulrl2lO++UHnqofPvYsVJEhPTWW349HQAA\nVQr4esJO06SJtGqV1Lt3+ba8PCk1VRoxwr66AAD1T70LYUkaNMjz4+hoqbhYSkiwpx4AQP3k6GvC\n/lJYKB09WvX/Z2eb+4SHDw9aSQAA1I8Q/tvfTOv39OnK/5/bkwAAdqgX3dFJSdL550uRVXy2mZnS\needJ7dqZjzdvlk6e5BoxACCw6kUId+1qHlUpm67S5TIfP/KI1KgRIQwACKx6EcI1ycry7Ip+6SWp\neXP76gEA1A9hfU24qMi7/c5ePalVKzN7FgAAgRTWITxpkvSTn1S/z7ffSkeOMCgLABB8YR3Ct94q\njR9f/T7VjYw+ckSaP18qKfF/bQAAhPU14cTEmvcpC+ELLzz3//bvl37zG2nUKOmSS/xbGwAAjg7h\npKQkRUZGKjk5WcnJyQE5R1aW1KmTmc7ybEOGSF99JcXGBuTUAIB6ztEhvHTpUr8v4HC2swdlnY0A\nBgAESlheE964Ubr9dun48Zr3rSmEAQAIlLAM4aNHpc8/r7yLuaKSEmnv3ppD+PRp6e23zSxaAAD4\nS1iG8Lhx0gcflM+AVZUDB6RTp2oO4exs6frrpffe81+NAAA4+ppwoGVlmeceParf78ILpT176LYG\nAPhXWLaEvZWZaeaI7tKl5n0JYACAv4VdCK9bZ7qPvZGZaVq5TFEJALBD2IXw7bdLTz3l3b61GRnt\ndvteEwAAlQm7EE5Lk+6/37t9fQ3h996T2rc3M2kBAFBXYTcwy9vJNQoLpUOHah6UVdHQodLzz0tt\n29auNgAAKgq7EPbW3r3m2ZeWcJMm0uTJgakHAFD/hF13tLeqWz0JAIBgcHQIJyUlKTExUUuWLKlx\n35Mnza1Gb7/t3bEzM6U2baQWLepYJAAAteTo7mhfFnA4dUqaNEm64ALvjl2XOaNfeUVauNDMygUA\nQG05OoR90by59Oij3u+flSVdemntztWhg9Svn1RUJDVsWLtjAAAQNiHsC8syLeGbb67d66+6yjwA\nAKgLR18TDpScHKmggEFZAAB7hUUInzolPfGEdPiwd/szMhoA4ARhEcLZ2WaWrEOHvNs/M9PMF92t\nW93O+/zz0urVdTsGAKD+qlUIL1iwQPHx8YqJiVFCQoK2bdvm1es2btyoqKgoXXbZZbU5bZV69ZKO\nH5cGDvRu/8xME8B1HVS1dCkjpAEAtedzCC9btkz33HOP5s6dq7S0NPXt21ejR49Wbm5uta/Ly8vT\n5MmTNWrUqFoXW53ISO9XQ8rK8m26yqqsWiX9+c91Pw4AoH7yOYTnz5+vqVOnatKkSerVq5eeeeYZ\nNW7cWAsXLqz2ddOmTdNtt92mhISEWhfrL3W5R7gilkAEANSFTyFcXFys7du3a+TIkWe2uVwujRo1\nSps2barydc8//7yys7P10EMP1b7SKliWVFLi/f5FReYaMoOyAAB28ymEc3NzVVJSori4OI/tcXFx\nysnJqfQ1e/fu1QMPPKBFixYpIsL/48D27ZPOO0/avNn7/UtK/BfClmWWT7Qs/xwPAFB/BHR0dGlp\nqW677TbNnTtX3bt3lyRZfk6rZs3MTFm+TFcp+S+EP/xQuuwyaccO/xwPAFB/+DRjVuvWrdWgQQO5\n3W6P7W63W+3atTtn/4KCAn3yySfasWOHZsyYIckEs2VZatiwod59910NHz68yvMlJSUpMtKzxOTk\nZCUnJ5/5uG1badYs7z+HrCzTcj6rMV9rQ4dKqalSnz7+OR4AoP7wKYSjoqI0YMAApaamKjExUZJp\n2aampmrmzJnn7N+sWTOlp6d7bFuwYIHWrl2r1157TV27dq32fL4s4OCtskFZLpd/jtewIVNYAgBq\nx+e5o2fNmqWUlBQNGDBAgwcP1vz583XixAmlpKRIkubMmaMvv/xSL774olwuly666CKP17dt21bR\n0dHq3bu3Xz4BX/lrZDQAAHXlcwhPmDBBubm5evDBB+V2u9WvXz+tWbNGbdq0kSTl5OTo4MGDfi+0\nMpZlrgfffLOZsMMbmZnS6NGBrQsAAG+4LH+PlPKD/Px8xcbGKi8vr9ru6CNHpEsuMWv7jh1b83GP\nHpVatZKWLZMmTPBjwZJmz5Zyc00tAAB4I6SXMmzTxqyI5O3biKws8+yP2bLOdumlUn6+/48LAAhf\nIR3CZbwdZFV2e9KFF/q/hokT/X9MAEB4C4tVlLyVmSl17iw1aWJ3JQAAhHgIFxX5tj8jowEAThLS\nIdyxo/TEE97vH+gQdrul3/3ODAADAKAmIRvCpaXSn/4kjRjh3f4lJdJnnwVmUFaZoiJp/nwpIyNw\n5wAAhI+QHZgVESFNmeL9/gcOSKdOBbYl3LmzuU2JJQ4BAN4I2Zawr/y9cENVCGAAgLccHcJJSUlK\nTEzUkiVL6nyszEwpOlrq0sUPhQEA4AeO7o6ubgGHZcukyEhp/HjvjrVnj2kFB2BJ43OcOGFm8zr/\n/MCfCwAQuhwdwtV5803T9ettCGdkeD+/dF3ddJOZxWv16uCcDwAQmkI2hF95xbf99+yRqlm62K8e\nfVRq3Dg45wIAhK6QDWFffPutuYc3WC3hyy4LznkAAKHN0QOz/KVsZHSwQhgAAG+EZAifPu39ykmS\n6YqWAjtRBwAAvgrJEH7qKal9e++DeM8eM1I5mNdp//c/adAg6fDh4J0TABBaQvKa8PDh0nnneb+E\n4Z49we+K7tBBuuAC6eTJ4J4XABA6XJblS8ducOTn5ys2NlZ5eXlV3ifsi169pDFjzLzOAAA4RUh2\nR/uiuFjat49BWQAA5wn7EN63zwzkIoQBAE4TciH89dfS738v5eR4t3/ZyGi7QnjNGumdd+w5NwDA\n2RwdwpUt4PDFF9Jf/yoVFHh3jIwMqXlzqW3bABVZg3/8Q3rhBXvODQBwtpAdmGVZ3o2OnjxZysqS\nNm3yc5FeOnGCKSwBAJVzdEu4Ok6+PakiAhgAUJWQDWFvWJb9IQwAQFVCLoRLS73fNydHys93Rgjn\n5Pg21SYAIPyFXAj37m2WCvSG3SOjy6xda6bZ3LvX3joAAM4SctNW/vrXJoi9sWePFBkpdesW2Jpq\nMmiQtHixCWIAAMqEXAhPmeL9vnv2mPmbo6ICV483mjaVkpPtrQEA4Dwh1x3tCwZlAQCcjBAGAMAm\nIRXCGzdKzz/v3b7Hj0sHDnh//TgYfv1r6fHH7a4CAOAUIRXC77/vfYhlZZlnJ7WEIyOliJD6igMA\nAinkpq30drrKRYukiROlY8ek2NgAFQoAQB04ul1W2QIO3k5XmZ4udepEAAMAnMvRtygtXbq0ygUc\narJrl9Snj58LAgDAjxzdEq6L9HRnhvCRI9LChb5NvwkACE8hE8L790stW3q3JGFhoZSdLV18ceDr\n8tWuXdIAuM7KAAAS1ElEQVQddzCFJQDA4d3RFTVtKs2eLXXtWvO+u3ebZye2hK+4QsrNlZo3t7sS\nAIDdQiaE27Y199l6Y9cu8+yke4TLREURwAAAI2S6o32Rnm4WbWjSxO5KAACoWq1CeMGCBYqPj1dM\nTIwSEhK0bdu2Kvf9z3/+o2uuuUZt27ZVbGyshg4dqnfffbfWBXvDqYOyKiotlQoK7K4CAGAnn0N4\n2bJluueeezR37lylpaWpb9++Gj16tHJzcyvdf/369brmmmu0evVqffrppxoxYoTGjRunnTt3+nTe\n114z01Z6Y9cuZw7Kqmj4cOnee+2uAgBgJ59nzEpISNCQIUP0xBNPSJIsy1Lnzp01c+ZMzZ4926tj\n9OnTR0lJSfrtb39b6f9XNmPW4MFSQoL05JPVH/vYMalFCzNj1q23ev95BdvKlVKbNtLll9tdCQDA\nLj4NzCouLtb27dv1wAMPnNnmcrk0atQobfLm3iGZ0C4oKFDLli19KnTLFun06Zr3KxsZ7fSWcGKi\n3RUAAOzmU3d0bm6uSkpKFBcX57E9Li5OOTk5Xh3jz3/+swoLCzVhwgRfTi2Xy4wsrkl6utSggdSz\np0+HBwAg6IJ6i9LixYv1yCOPaOXKlWrdunVAzpGeLl14oRQdHZDDAwDgNz6FcOvWrdWgQQO53W6P\n7W63W+3atav2tUuXLtUdd9yh5cuXa8SIEV6dLykpSZGRniUmJycrOTm5yteEwqCsMhkZZgKSF180\ns4EBAOoXn0I4KipKAwYMUGpqqhK/v6hpWZZSU1M1c+bMKl+3ZMkS/exnP9OyZct07bXXen2+sgUc\n/vMfM5L4v/+t+d7f9HTpF7/w+hS2atZMOn5c+vprQhgA6iOfu6NnzZqllJQUDRgwQIMHD9b8+fN1\n4sQJpaSkSJLmzJmjL7/8Ui+++KIk0wWdkpKiJ598UoMGDTrTio6JifF6haSuXaUJE6TGjavf78gR\nE2ih0hLu2FFau9buKgAAdvE5hCdMmKDc3Fw9+OCDcrvd6tevn9asWaM2bdpIknJycnTw4MEz+z/3\n3HMqKSnRjBkzNGPGjDPbJ0+erIULF3p1zv79zaMmZdNVOn2iDgAApFrcJxwMld0n7I0FC6Rf/cqs\nouTNSGoAAOwUVnNHp6dLvXqFXgCnpUkffGB3FQCAYAuJEH79dTOSuCahMGd0Zf7yF2nuXLurAAAE\nW0iE8M9+ZqZ5rI5lmdHToRjCTzwhpabaXQUAINhCYj3hnJyap6z8/HMpP9+7AVxOE6B5SwAADhcS\nIdywoXlUZ8cO89yvX+DrAQDAH0KiO9obO3ZIbdtKNUzc5WiFhVJxsd1VAACCJaxCuF8/s9BDKDpw\nwHRLc20YAOoPx4fwCy9Iw4bVvF9ZCIeqzp2lxx+XLrnE7koAAMHi6GvCSUlJ+vbbSLVpkyyp6kUb\njh41LclQDmGXS5o+3e4qAADB5OgQLlvAoSY7d5rnUA5hAED94/juaG/s2CHFxEg9ethdCQAA3guL\nEP70U+nSS6UGDeyupO7eeku64QYz+QgAILw5PoTffls6fLj6fT75RBo4MDj1BFrjxmbN5BMn7K4E\nABBojg7hoiLp+uuld9+tep+CAikzM3xC+KqrpMWLTRADAMKbowdmRUVJbre53luVtDTTdTtgQPDq\nAgDAHxwdwi6XmQWrOp98YkK6d+/g1AQAgL84ujvaG598Ym5NinT02wnfrV4tvfee3VUAAAIp5EN4\n+/bwuR5c0VNPSS+/bHcVAIBAcnQIP/usNHFi1f9/9KiUlSUNGhS8moJl2TLpxRftrgIAEEiO7sRt\n0cLMqVyVzZvN8+WXB6eeYGra1O4KAACB5ugQvvlmqbpZKzdtMisPde8evJoAAPAXR3dHJyUlKTEx\nUUuWLKn0/zdvNq3gUF2+0BsHD0rffGN3FQCAQHBZlvMmSMzPz1dsbKzy8vKqXMChpMR0V8+ZYx7h\n6ORJqVUr6Xe/k+691+5qAAD+5uiW8MaNZvBVZXbvNrNlheP14DIxMdKqVdIdd9hdCQAgEBwdwmPH\nSu+/X/n/bdpkFmwIx5HRFQ0fXv11cQBA6HL0wKxPPql60NXHH5uVk5hjGQAQqhzdEr7wQql583O3\nW5a0dq30wx8Gvya7FBdLpaV2VwEA8CdHh3BVPv9cOnBAGjHC7kqC49AhqVMnKTXV7koAAP4UkiG8\ndq25LenKK+2uJDg6dpR++UspPt7uSgAA/uToEP7Vryrfvnat1L9/5V3V4cjlkn77W+mCC+yuBADg\nT44O4cqUXQ+uL13RAIDw5egQnj//3G379kmHDxPCAIDQ5+gQrsw775i1g3/wA7srCb4DB6SRI6XM\nTLsrAQD4Q8iF8FtvmVuT6uMEFm3bmvuijx+3uxIAgD84OoRvuslzAYeCAnM9eNw4mwuzSXS0tHKl\nNGCA3ZUAAPzBkTNmffedeR4/fqmmTi1v8r73nlRUVH9DGAAQXhzZEo78/q3ByJGe2998U7roIqlb\nt+DXBACAvzk6hNu2Ld92+rT09tvS9dfbU5OTfPONdNddZuYwAEDocmQIV+bdd6UjR6SkJLsrsV90\ntLR6tbRnj92VAADqwpHXhCvz0ktSnz5Sv352V2K/Jk3MbUoul92VAADqwpEt4VWrPD8+dkxasUKa\nPJngKcPXAQBCnyND+OhRz49feEEqKZFuvdWWclCNstvH4Gx8n0IH36vQ4K/vkyNDeOLE8n+fOiX9\n5S/SbbdJHTrYV5MTWZb0/PPmerld+INRta++kr7+2nNbcbGZdrW4OLi18H0KHXyvQkNYh3BFL78s\nffmldP/9dlfiPC6XtGiR9MEHdldSf6SmStu3e257/32pVSszcLCiiRPNKPaKMjLM2tBpaZ7b77vP\nvNGs6NQp6amnzh0Fb1m1Lh+Awzg6hN1uac4cacIEqXfvuh8vEO8w7T7mqlXSvHn+P66d7P6aSmae\n7v/7P3MZpKL775eefdbzuN26SffeK0VFee775z9Lv/mN57auXc2tdj17em4fNEgaNsyz1m++ke65\n59y5wu+/Xxo40HObZUmPPSalp5+7PVAC9fPkhO+/XccMlFD5/EPpZ8pfHBnCbrd5nj5dioiQnnzS\nP8cNlR8aX47ZsGFgjmunYH9NLcsM/qto3z7pgQekgwc9t69ZI/39757H7dbNvFk8e33ryy4zI/or\natZMGjtWio313D5hgvSLX3jW2qGDmT3u6qs99x03Trr7bs9tJ0+aVcd27/bcPm+eFB9/7uf8l79I\n27Z5bist9S20Q+kPZqgcM1BC5fMPpZ8pfwn6LUqWZamgoKDK/z99Wrr44nxJ0tq1+Vq61NwXm59f\n93OfPn1a+f44UAgeM1DHDYdjTp4snTgh/fvf5dv69jUBHBXl+bMXFeW5gIYdn3/fvuZx9i6ffWae\nK25PSDDB//bb5ce0LOmJJ8wbuIot8n/+U3r4YdMLEFHh7flDD0mDB0vXXVe+ze2WvvrqtL75Jt+j\nB+DoUalRI3MbXW0F82tqWaa3o0EDzzsO8vLM1+C888q3ffeduTTWqZP52pUd8/PPzTG6dy/f9+RJ\n8yanTx+pZcvy7enpUm6uNHx4+bbSUmnpUmnIEHOMsuOmpZnLFrff7lnzggXmVs0rrijflpFhFre5\n6y7PN+aLFplJjyp+/ocPS6++an7uK9b21lvm87jhBs+vw8svS4mJUpcu5ds3bJC++MLza1paanqJ\nhg+XevUq33fHDvN5TJni+Xm88orp4aw4F35e3mk9/ni+UlI8e5beeMNc8inrLZKknBzTq/TjH0st\nWpRvT001tVR88/rdd6f15JP5uuYaz7FFW7easRsVP+fSUlPb0KHSBReUb9+1y3z/brnFfOzNz+l5\n550nVw23srgsK7hXmPLz8xV7djMAAIAwk5eXp2Y1LPkX9BCuqSUsSV98ka9LL+2sgwcP1vgJwFi9\nWpo926wy1bq13dU406lTUo8e0p13moFQqL2TJ8016w4dPFvNH39sWsF9+5ZvO3TItPSmTDGtmTKL\nF5sFWVJSyrcdPWpa3jNmeLamVqwwLf177y3fZlnSb39rWkIVW1Nbtpg6fvUrz5pfeMG0TCteT8/O\nNvtOmODZ8tqwwfQgVPw8jh2Tdu40r6/Y0t+/39RSsSVcXGx6Utq3l2JiyrcfP256+ypeuihrjUdE\neH4tYY+KiVixEWtZ5uHL98iRLWFvlLWWvXkXAcOyzC/32YOD6rPMTDMQqlGj8m379plrpPyxA+AE\n/CkKEy4XAVzRoUPmetPy5Z7bu3cngAE4R8jMHQ34olMn6Z13pB/+0O5KAKBqtAnC0JNPmttV6otD\nh8xIxp07Pbdfc41nVzQAOA0t4TD01VdmmH19ERcnde5srokDQCihJRyG/vAH6Y9/DNzx586dq4iI\nCI/HRRddFLgTnuXQITM6t0xUlLRsmecI2fpow4YNSkxMVMeOHRUREaGVK1ees8+DDz6oDh06qHHj\nxrr66qv1WdnNxQiamr5PU6ZMOef3a+zYsTZVW3899thjGjx4sJo1a6a4uDjdeOONysrKOme/uv5O\nEcJhKBjLHPbp00dut1s5OTnKycnRRx99FPiTSiookC65xH+zqIWTwsJC9evXT08//XSlt0X88Y9/\n1FNPPaVnn31WW7duVZMmTTR69GgVFRXZUG39VdP3SZLGjBnj8fvl5BmfwtWGDRv0y1/+Ulu2bNH7\n77+v4uJiXXPNNTpZoQXgl98py4Hy8vIsSVZeXp7dpYS8khLLKiry7zEffvhhq3///v49qA/eesuy\n+NGonsvlst544w2Pbe3bt7f++te/nvk4Ly/Pio6OtpYtWxbs8vC9yr5PKSkp1o033mhTRajKkSNH\nLJfLZW3YsOHMNn/8TtESDmOlpdLo0edOWuAPe/fuVceOHdW9e3dNnDhRB8+eZNlP0tPNtHcVXXed\nmUgB3svOzlZOTo5Gjhx5ZluzZs00ZMgQbdq0ycbKUJl169YpLi5OvXr10vTp03X07EXWEXTHjh2T\ny+VSy+/n+vTX7xQDs8JYRISZH/b88/173ISEBL3wwgvq2bOnvvrqKz388MO68sorlZ6eriZ1mTT4\nLJYl/exnZs7aV1/122HrpZycHLlcLsXFxXlsj4uLU05Ojk1VoTJjxozR+PHjFR8fr3379mnOnDka\nO3asNm3aVOPsSwgMy7J09913a9iwYWfGv/jrd4oQDnMTJ/r/mKNHjz7z7z59+mjw4ME6//zz9eqr\nr2rK2bO014HLZcK3XTu/HRJwvAkTJpz598UXX6xLLrlE3bt317p16zRixAgbK6u/pk+frt27d2vj\nxo1+P7aju6OTkpKUmJjIoASHi42NVY8ePeo80vajj8wKMRV16eLbco2oXLt27WRZltxl64R+z+12\nqx3vchwtPj5erVu3ZiS7Te68806tWrVK69atU/v27c9s99fvlKNDeOnSpVq5cqWSk5PtLiUsZGdL\ns2adu1B9XR0/flz79u3z+AGtjXXrzET//q4P5g95u3btlJqaemZbfn6+tmzZoqFDh9pYGWpy6NAh\nffPNN3X+/YLv7rzzTr3xxhtau3atulRcy1H++51ydAjDvzIzzZqhubl1O859992n9evX64svvtDH\nH3+sG2+8UZGRkXV+s3T//dKHH5q1XeG7wsJC7dy5Uzu+H8m2f/9+7dy588ygubvvvluPPvqo3nzz\nTf3vf//TpEmT1KlTJ91QcTFVBFx136fCwkLNnj1bW7Zs0RdffKHU1FT96Ec/Uo8ePTwuAyHwpk+f\nrkWLFmnx4sVq0qSJ3G633G63vvvuuzP7+OV3yn8DuP2HW5QC59Spuh8jKSnJ6tixoxUdHW117tzZ\nSk5Otvbv3+/TMd55x7JGjPD/7VP12bp16yyXy2VFRER4PKZMmXJmn4ceeshq3769FRMTY11zzTXW\n3r17bay4fqru+3Ty5Elr9OjRVlxcnNWoUSMrPj7emjZtmvX111/bXXa9U9n3KCIiwnrxxRc99qvr\n7xRLGcIWO3dK8+ZJTz8ttWhhdzUAYA+6o+uxHTukMWPq3j1dE8uSPv3Uc1vfvtKSJQQwgPqNEK7H\niopMQDZtGtjzvPaaNHCgtHdvYM8DAKGG7mh4sKy6zz19/LhnsBcXSxs2SCNGBGdeawAIFbSE4eGR\nR6SUlNq//p//NDN0VRhAqKgo6aqrCGAAOBshDA9du0p9+nhuO3TIDKDKy/Pc/uyz5rpuRcOHS3/6\nUyArBIDwQXc0avTuu9LYsdKBA1KHDuXbJ06U4uKkxx+3rzYACGWEMLxS9lNClzIA+A8LOMArhC8A\n+J+jrwmzgAMAIJzRHQ0AgE0c3RIGACCcEcIAANiEEAYAwCaOvCZsWZYKCgp03nnnycWwXABAmHJk\nCAMAUB/QHQ0AgE0IYQAAbEIIAwBgE0IYAACbEMIAANiEEAYAwCaEMAAANvn/Zsi71CRQ97EAAAAA\nSUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt_sol1 = plot(sol, [t, 0, 20], linestyle ='-')\n", "plt_sol2 = plot(diff(sol,t), [t, 0, 20], linestyle =':')\n", "(plt_sol1 + plt_sol2).show(figsize=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

伝達関数H(s)を使った解法

\n", "\t

\n", "\t\t次に、ラプラス変換と伝達関数を使って同じ微分方程式を解いてみましょう。\n", "$$\n", "\t\t\\ddot z + \\dot z + z = u(t)\n", "$$\n", "\t\tの両辺をラプラス変換します。zのラブラス変換をZ, uのラブラス変換をUとすると以下の様になります。\n", "$$\n", "\t\ts^2 Z + s Z + Z = U\n", "$$\n", "\t\tここで伝達関数$H(s) = \\frac{Z}{U}$を計算します。\n", "$$\n", "\t\tH(s) = \\frac{1}{s^2 + s + 1}\n", "$$\n", "\t\t単位ステップ応答のラブラス変換は、$\\frac{1}{s}$なので、$Z(s)$は以下の様に求めまります。\n", "$$\n", "\t\tZ(s) = H(s) U(s) = H(s) \\frac{1}{s} = \\frac{1}{s^3 + s^2 + s}\n", "$$\n", "\t\t最後に$Z(s)$をラプラス逆変換すると$z(t)$の解析解が求まります。\n", "\t

\n", "\t

\n", "\t\t上記の処理をSageを使って行ってみます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "-1/3*(sqrt(3)*sin(1/2*sqrt(3)*t) + 3*cos(1/2*sqrt(3)*t))*e^(-1/2*t) + 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(s, t) = var('s t')\n", "P = 1/(s*s + s + 1)\n", "# 単位ステップ応答は伝達関数/sの逆ラプラス変換で得られる\n", "ip = inverse_laplace(P/s, s, t)\n", "show(ip)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

結果のプロット

\n", "\t

\n", "\t\t求まった解をプロットするとdesolveの同じ結果と同じ結果を得ることができました。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFnCAYAAACRo/HLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl8VNX9//H3hAQSAgTZkoggEUUQFCgKwa1SkGCtqbQV\nE6wCasUV+wO1D7RflFqr1u9X6gJVUVzZWuqCVRaNIKgRMCIGQgICsgiJRCSBJECW+/vjGJKBLDPJ\nzNw7M6/n4zGPgZsz936yvufce+45LsuyLAEAgICLsLsAAADCFSEMAIBNCGEAAGxCCAMAYBNCGAAA\nmxDCAADYhBAGAMAmhDAAADZxZAhblqXi4mIxjwgAIJQ5MoQPHTqkuLg4HTp0yO5SAADwG0eGMAAA\n4YAQBgDAJoQwAAA2IYQBALAJIRwCKiulvXulb76Rjh2zuxoAgKcI4SBUWSktWyalpUmnny5FR0td\nu0pnnSXFxJjniROlFSukqiq7qwUA1CfS7gLgndWrpVtukXJzpb59pfR0E8SnnSa1bi1t3y59/bX0\n/vvSCy9IffpIDz8s/eY3kstld/UAgNpclgNnxCguLlZcXJyKiorUrl07u8txBMuSnnlGmjJFSk6W\nnnhCGjKk/mC1LBPYjzwiLV8uDRokPfuseS0AwBk4HR0ESkul66+X7r7bPFasMGHaUM/W5ZIuvdSc\ntl650vz/ooukP/+Z68YA4BSO7glfccUVioyMVHp6utLT0+0uyxYlJdKIEdKGDdKcOeY6cFNUVEiP\nPy499JDUr5+0cKHUq5dPSwUAeMnRIRzup6MrKqRf/1patcr0fs8/v/n7XL/eXEfet0+aN0+68srm\n7xMA0DScjnawe+4x13P/8x/fBLAkDRworV0rXXaZlJoqzZ7tm/0CALxHCDvUu+9KTz0lPfmkNHKk\nb/fdrp305pvSbbeZkdYPPWQGcgEAAotblBxo715pwgTTU73zTv8co0ULM9r61FOlBx4wx5w1S4rk\nJwIAAoY/uQ50991SVJQZiOXPe3tdLun++00Q33yzdPiw9PrrJqABAP5HCDvM8uXSokXS3LlSx46B\nOeb48VJsrBmwFRkpvfwyQQwAgUAIO8jRo+b082WXmUAMpGuuMdNhXnedCeCXXpIiGDEAAH5FCDvI\n009LO3ZIb71lzxSTaWkmiK+/3vSIn3+eIAYAfyKEHaKkxExFeeONZk5ou1x3nbk/ecIEE8SzZjHn\nNAD4CyHsEM8/L/34ozR1qt2VSOPGmR7xTTeZFZqefJIgBgB/IIQdoKzM9IKvv17q0cPuaowbbzTX\nqG+/XTrlFGnaNLsrAoDQQwg7wEsvSd9/b24XcpLbbjO98wceMEF81112VwQAocXRIZyWlhbyCzgc\nPSo99pgZDX3mmXZXc7KpU00QT5pkgvj3v7e7IgAIHY4O4QULFoT8Ag5z55rZqh54wO5K6uZySX//\nuwni8ePNlJepqXZXBQChgVWUbDZ4sNSpk/T++3ZX0rDKSunaa6X//ldautTcywwAaB7uArXRV19J\n69aZRRScrkUL02u/9FLTE/7iC7srAoDgRwjbaPZsKTExeNb0bdXKrL7Ut680apS0ebPdFQFAcCOE\nbVJSIr3xhrkVKCrK7mo816aN9N575s3D5ZdLO3faXREABC9C2Cb/+pd06JCZECPYdOhgFppo1Uoa\nMUIqKLC7IgAIToSwTV54QRo5UkpKsruSpklMlD74wPToU1KkgwftrggAgg8hbIPsbOnzz4NjQFZD\nzjjD9Ih37ZJ+9SuptNTuigAguBDCNpg715zSveoquytpvn79pCVLzEjv3/5WOnbM7ooAIHgQwgFm\nWdKCBdLvfhdcA7IaMmSI9Pbb0kcfmRm1KirsrggAggMhHGBr1pgRxWlpdlfiWyNGSAsXmrWQq5dD\nBAA0jBAOsAULzKCmSy+1uxLfu/pqM+r7zTelsWOl8nK7KwIAZ3N0CKelpSk1NVXz58+3uxSfqKw0\nITVmjJmBKhSNHi39+9+mR0wQA0DDmDs6gFatkn7+c+mzz6ShQ+2uxr/eeUe65hozG9j8+VJ0tN0V\nAYDzOLonHGreesucih4yxO5K/O/XvzanpZcuNfdD//ij3RUBgPMQwgFiWWYE8dVXSxFh8lX/1a/M\niOlNm6SLL5Z277a7IgBwljCJA/tt2CB9+60J4XAydKg5/V5SYv6dnW13RQDgHIRwgLz7rtSuXXiu\nw3v22VJmptS5s3TJJdLKlXZXBADOQAgHyNKlZtWhli3trsQeiYnSxx9LF1xg5pqeN8/uigDAfgEJ\n4dWrVys1NVVdu3ZVRESEFi9eHIjDOsaPP5q5okeNsrsSe7VrZ5ZBTEszE3r88Y/cwgQgvAUkhEtK\nSjRgwADNmjVLLpcrEId0lA8/lKqqTA8w3LVsKb3yivTMM9LMmdLw4dK+fXZXBQD2iAzEQUaNGqVR\nP3UDHXhbst8tXSr17St162Z3Jc7gckl33ikNHGjuJT7vPOnll81oagAIJ1wT9jPLMiEc7qei63LR\nRWbUeHKyWVFq4kSpqMjuqgAgcALSEw5nGzdKe/cSwvXp3FlavFh6/nnp3nul9983//7lL+2uLLhV\nVUkHD0o//FDzOHDAPBcXm7Wfy8rM84n/PnrUvL6xh2VJkZGeP1q2NI+oqJp/1/Vo6OPevjYqypx5\nAZyKEPazpUul1q3NZBWom8sl3XqrCd5bbjFTXY4eLf3f/0lJSXZX5yyWJeXnm3vO9+yRvvvOPKr/\nnZ9vgvbHH03bE7VubQbIxcZKMTHm/9WPNm2kLl2kVq3M3OYREQ0/JDMfekVF44/ychP0x47VPMrL\n3f9/4sNXg/aqg7kp4V8d4r56SDXP1d+f2t+nE7fxMed+7A9/MPPjN5ejQzgtLU2Rke4lpqenKz09\n3aaKvLd0qTRsGHMne6J7d2nJErMk4j33SH36SLfdJv3pT1JCgt3VBU5VlbRjh5SXJ23fLm3b5v5c\nVlbTNiZG6trVPLp3lwYPljp1kjp0kDp2dH906BBcP4eWZQK8voBuLMCb8rHqjx89Kh06ZP5vWb55\nVH9OlnVyINfurXvzMX+3D4ZjV78hDPSx27SRTwR8AYeIiAi9/fbbSk1NrbdNqCzgcPiw+cP35JNm\nIBI8d/iw6QnPmGH+EN52m3TffVJ8vN2V+Y5lmak8N240U3tWP+fk1ARtVJQ5G9Czp3TGGTXPSUlm\noF/79pxuBYJZQHrCJSUl+uabb46PjN6+fbs2bNigDh06qFsIDxlescK8q+Z6sPfatJEefFC6+24T\nxP/4h/TPf5r7i++8U+rf3+4KPVd9CvnEsN20yfS2JPP5nnOO+bzGjpX69ZN69zY93FBd9hJAgHrC\nH3/8sYYNG3bSPcLjxo3TnDlzTmofKj3hO+4wp6O3bbO7kuD344/SrFkmiL/7rias0tLMaVgnsCxp\n/37Tk60dths31qwiFR1tTrP362duW6t+7t49fBb2AFCD9YT9qE8fs37wc8/ZXUnoqKgwI6jnzTOj\nqsvKzP3Gv/iFeVxyidS2rX9rsCzzRiAnxzw2b67594EDpk1UlJkz+8SwPeMMerYAahDCfpKfb+ZL\nnj/f9Nbge4cOmSBevlzKyDDBGBlpwq5PH3N6t08fE4annSbFxTXe26yqMkG6f3/NY/duM1Dq22/N\n844d5pq1ZHq2vXu7H69PH+mss0wQA0BDCGE/WbBASk83UzKG08heu1iW9M03JozXr6/pnf7wg3u7\nNm1MT7ltWzOgqby85nHsmLm3tqrK/TUxMVKPHmYwVPVzdfD26EHPFkDTOfoWpWC2cqX5Q00AB4bL\nZXqfZ53lvn3/fmnLFvNm6OBB03suLq4ZEBUV5f445RQzgUiXLua5c2ezjRHIAPyBEPaTFSvM4gSw\nV3WQAoATMR7TD/buNb2vyy6zuxIAgJMRwn7wySfm+dJL7a0DAOBshLAfZGaawTtcDwYANIQQ9oPP\nPpOGDrW7CgCA0zk6hNPS0pSamqr58+fbXYrHjhwxt8hceKHdlQAAnM7Ro6MXLFgQdPcJZ2WZe07p\nCQMAGuPonnAwysw0a7Oed57dlQAAnI4Q9rHMTOmCC8z0iQAANIQQ9iHLMoOyuB4MAPAEIexDO3ea\nhRu4HgwA8AQh7EOZmeY5OdneOgAAwYEQ9qHMTOnMM5mrGADgGULYh5ikAwDgDULYR0pLpQ0bGJQF\nAPAcIewjX3whVVTQEwYAeI4Q9pHMTKlNG6lfP7srAQAEC0LYRzIzpcGDpRYt7K4EABAsHB3CwbKA\nA5N0AACawtGTKwbLAg7bt0v793M9GADgHUf3hIMFk3QAAJqCEPaBzEzp7LOlDh3srgQAEEwIYR/4\n4gszKAsAAG8Qws1UXm4m6Rg0yO5KAADBhhBuppwc6ehR6Wc/s7sSAECwIYSb6csvJZdLGjjQ7koA\nAMGGEG6mrCwzKKtNG7srAQAEG0K4mbKyuB4MAGgaQrgZKirMoCyuBwMAmoIQbobcXKmsjJ4wAKBp\nCOFmyMoyzwzKAgA0haND2OkLOGRlSWedJQXB9NYAAAdiAYdm+PJLTkUDAJrO0T1hJ6uslNavJ4QB\nAE1HCDfRli1SaSkjowEATUcIN9H69eaZQVkAgKYihJvoq6+k00+XTjnF7koAAMGKEG6ir76SBgyw\nuwoAQDBrUgjPnDlTSUlJiomJUXJystatW9dg+7lz52rAgAGKjY3VqaeeqptuukkHDhxoUsFOYFmE\nMACg+bwO4YULF2rKlCmaPn261q9fr/79+yslJUWFhYV1tv/00081btw4/eEPf1BOTo4WLVqktWvX\n6pZbbml28XbZt0/av58QBgA0j9chPGPGDE2cOFE33HCDevfureeee06tW7fWnDlz6mz/+eefKykp\nSXfccYdOP/10XXjhhZo4caLWrl3b7OLtsmGDee7f3946AADBzasQLi8vV1ZWloYPH358m8vl0ogR\nI5SZmVnna4YOHardu3dryZIlkqSCggItWrRIV155ZTPKtld2tlm68PTT7a4EABDMvArhwsJCVVZW\nKj4+3m17fHy88vPz63zNhRdeqDfeeEPXXnutWrZsqcTERLVv317PPvts06u2WXa21K+fFMGwNgBA\nM/g9RnJycnT33XfroYce0pdffqlly5Zpx44dmjhxor8P7TfZ2dK559pdBQAg2Hk1d3SnTp3UokUL\nFRQUuG0vKChQQkJCna957LHHdPHFF2vy5MmSpH79+mnWrFm65JJL9Mgjj5zUq64tLS1NkZHuJaan\npys9Pd2bsn2qvFzavFm66SbbSgAAhAivQjgqKkqDBg1SRkaGUlNTJUmWZSkjI0OTJk2q8zWlpaWK\niopy2xYRESGXyyXLsho8nhMXcNi6VTp2jJ4wAKD5vD4dPXnyZM2ePVuvvfaacnNzdeutt6q0tFTj\nx4+XJE2dOlXjxo073v6qq67Sm2++qeeee047duzQp59+qrvvvltDhgypt/fsZNnZ5pkQBgA0l9dL\nGY4ZM0aFhYWaNm2aCgoKNGDAAC1btkydO3eWJOXn52v37t3H248bN06HDx/WzJkzdc8996h9+/Ya\nPny4HnvsMd99FgGUnS0lJkodO9pdCQAg2Lmsxs4J26C4uFhxcXEqKipy3Ono1FTp6FFp2TK7KwEA\nBDtusvESI6MBAL5CCHvh0CHp22+l886zuxIAQCgghL2wcaN5picMAPAFQtgL2dlSixZSnz52VwIA\nCAWEsBeys6WzzpKio+2uBAAQCghhLzAoCwDgS4SwhyyLEAYA+BYh7KF9+6QDBwhhAIDvODqE09LS\nlJqaqvnz59tdir7+2jwTwgAAX/F62spActICDtnZUmyslJRkdyUAgFDh6J6wk2RnS337ShF8xQAA\nPkKkeCg7m5myAAC+RQh7oKJC2ryZ68EAAN8ihD2wdatZOYkQBgD4EiHsgexs89yvn711AABCCyHs\ngZwcqUsXqXNnuysBAIQSQtgDmzezaAMAwPcIYQ/k5EjnnGN3FQCAUEMIN6KiQsrLI4QBAL5HCDdi\n+3apvJzT0QAA3yOEG5GTY57pCQMAfM3RIeyEBRxycqT27aWEBNtKAACEKBZwaMTmzaYX7HLZWgYA\nIAQ5uifsBDk5XA8GAPgHIdyAqqqanjAAAL5GCDdg1y6prIwQBgD4ByHcgOqR0ZyOBgD4AyHcgJwc\nKTZW6tbN7koAAKGIEG5A9ZzREXyVAAB+QLw0gDmjAQD+RAjXw7K4PQkA4F+EcD327pWKi+kJAwD8\nhxCux+bN5pkQBgD4CyFcj5wcqVUrKSnJ7koAAKHK0SFs5wIOOTnS2WdLLVoE/NAAgDDBAg71yM2V\neve25dAAgDDh6J6wnfLyTE8YAAB/IYTrUFws5ecTwgAA/yKE67Bli3kmhAEA/kQI1yEvzzz36mVv\nHQCA0EYI1yEvT0pMlGwaEwYACBNNCuGZM2cqKSlJMTExSk5O1rp16xpsf+zYMT3wwAPq0aOHoqOj\ndcYZZ+iVV15pyqEDIi+PXjAAwP+8vkVp4cKFmjJlil544QUNHjxYM2bMUEpKirZs2aJOnTrV+Zpr\nrrlG+/fv18svv6yePXtq3759qqqqanbx/pKXJw0ZYncVAIBQ57Isy/LmBcnJyRoyZIieeuopSZJl\nWerWrZsmTZqk++6776T2S5cu1dixY7V9+3a1b9/eo2MUFxcrLi5ORUVFAb9PuKpKatNG+utfpcmT\nA3poAECY8ep0dHl5ubKysjR8+PDj21wul0aMGKHMzMw6X/Puu+/q/PPP1+OPP67TTjtNZ599tu69\n914dOXKkeZX7yXffSWVljIwGAPifV6ejCwsLVVlZqfj4eLft8fHxyqseUnyC7du3a/Xq1YqOjtbb\nb7+twsJC3XbbbTpw4IBeeumlplfuJ9WfBiEMAPA3v09bWVVVpYiICM2bN09t2rSRJD355JO65ppr\nNGvWLLVq1crfJXglL0+KipJ69LC7EgBAqPMqhDt16qQWLVqooKDAbXtBQYESEhLqfE1iYqK6du16\nPIAlqU+fPrIsS3v27FHPnj3rPV5aWpoiI91LTE9PV3p6ujdleyUvT+rZU4p09KzaAIBQ4FXUREVF\nadCgQcrIyFBqaqokMzArIyNDkyZNqvM1F110kRYtWqTS0lK1bt1akpSXl6eIiAiddtppDR7PjgUc\nmDMaABAoXt8nPHnyZM2ePVuvvfaacnNzdeutt6q0tFTjx4+XJE2dOlXjxo073n7s2LHq2LGjJkyY\noM2bN2vVqlW67777dNNNNznuVLRECAMAAsfrk65jxoxRYWGhpk2bpoKCAg0YMEDLli1T586dJUn5\n+fnavXv38faxsbH64IMPdNddd+mCCy5Qx44dde211+rhhx/23WfhI2Vl0q5dhDAAIDC8vk84EOy6\nTzg7WzrvPOmTT6SLLgrYYQEAYYq5o2vh9iQAQCARwrXk5UkdOkj1zL4JAIBPEcK1sHADACCQCOFa\nGBkNAAgkQvgnlkUIAwACixD+yf79UlERIQwACBxC+CeMjAYABBoh/JO8PMnlks480+5KAADhwtEh\nnJaWptTUVM2fP9+n+/32W6nWpF6STAj36CE5cCZNAECICrsZs0pLpaQkafJk6U9/qtl+1VVSRYW0\nZIlPDwcAQL0c3RP2h9atpTfflG67rWZbcbH0wQfmdDQAAIESlqvmnjgvdFSUdOyYNGiQPfUAAMJT\nWPSEKyqkw4fr//ju3eY+4WHDAlcTAABhEcIvvij17m1OO9eF25MAAHYIi9PRV1xhrvfWN8YrL0+K\njZVOPdX8/4svzNrCl1wSuBoBAOEnLEL49NOliRPr/3j1wg3VA7Puv19q04YQBgD4V1iEcGNOnDN6\n7lzplFPsqwcAEB5C+prw0aOetTsxhDt3liJ5ewIA8LOQDuGxY6Wbb264zcGD0vffMygLABB4IR3C\n6enS1Vc33GbLFvNcVwj/8IM0a5a5fQkAAF8L6ZOuv/td422qb08666yTP5adLU2ZYu4f7tPHt7UB\nAODonrC/FnCoLS/P3JrUtu3JH/v5z6U9ewhgAIB/OLonvGDBAp8v4HCiEwdl1eZySR07+vXwAIAw\n5uiecFN9+KF0yy2ejY5uKIQBAPCnkAzhgwelnTsbXxu4qkraurXxEK6okJYvl8rLfVcjAAAhGcK/\n+520bFnj7Xbvlo4caTyEc3KklBRpxQrf1AcAgOTwa8L+5unCDeedJ23YIJ17rv9rAgCEj7AP4ZYt\nzdzSjTnvPP/XAwAILyF3Onr1amnXLs/a5uVJZ54ptWjh35oAAKhLyIXwDTdITz3lWdumjIwuLPS+\nJgAA6hJyIZyVJd17r2dtvQ3hDz+UEhKk7dubVhsAALWF3DXhDh08a1dSYkZHexPCQ4dKzz9vVlkC\nAKC5Qi6EPfXNN+bZmxCOjZVuusk/9QAAwk/InY72VPXtSb162VsHACB8OTqEvVnA4cgR6YwzpCVL\nPNt3Xp6ZF5q5oQEAdnH06WhvFnA4ckQaM0bq0cOzfTdnzujXX5cWLpT++9+mvR4AAMnhIeyN9u2l\nxx7zvH1eXtNnwOrcWUpKMnNJR0U1bR8AADj6dLS/WFbzesKjRknPPEMAAwCaJyxDuKBAOnSIJQwB\nAPYKiRA+ckR69llp3z7P2jMyGgDgBCERwtu2SZMnm8k3PJGXJ0VESD17Nu+4r71mZtECAKApmhTC\nM2fOVFJSkmJiYpScnKx169Z59LpPP/1UUVFR+tnPftaUw9arb1/p8GFp0CDP2uflmYFVrVo177hz\n5kjLlzdvHwCA8OV1CC9cuFBTpkzR9OnTtX79evXv318pKSkqbGRlg6KiIo0bN04jRoxocrENadnS\n89WQmjMoq7YPPpD+/vfm7wcAEJ68DuEZM2Zo4sSJuuGGG9S7d28999xzat26tebMmdPg62699VZd\nd911Sk5ObnKxvuKrEGZ0NACgObwK4fLycmVlZWn48OHHt7lcLo0YMUKZmZn1vu7ll1/Wjh079OCD\nDza90npYllRV5Xn7Y8ekHTsYGQ0AsJ9XIVxYWKjKykrFx8e7bY+Pj1d+fn6dr9m6davuv/9+zZ07\nVxERvh8HtnOnFBcnffqpZ+23b5cqK303MtqypJwc8wwAgDf8Ojq6qqpK1113naZPn66ePw1Ftnyc\nVq1bS9OmSWed5Vn76tuTfNUTXr7cDAyr3i8AAJ7yatrKTp06qUWLFiooKHDbXlBQoISEhJPaHzp0\nSF988YW++uor3XHHHZJMMFuWpZYtW2r58uW67LLL6j1eWlqaIiPdS0xPT1d6evrx/3fpIt17r+ef\nQ16e1KaNlJjo+Wsacuml0rJlZrQ1AADe8CqEo6KiNGjQIGVkZCg1NVWS6dlmZGRo0qRJJ7Vv166d\nNm7c6LZt5syZWrFihf7zn/+oRyOrLXizgIOnqgdluVy+2V9MjDRypG/2BQAIL14v4DB58mSNHz9e\ngwYN0uDBgzVjxgyVlpZq/PjxkqSpU6dq7969evXVV+VyuXTOOee4vb5Lly6Kjo5Wnz59fPIJeMtX\nI6MBAGgur0N4zJgxKiws1LRp01RQUKABAwZo2bJl6ty5syQpPz9fuz2duqqZKiqkRx+Vxo71fPar\nvDx6rgAAZ3BZvh4p5QPFxcWKi4tTUVFRg6ej9+yRBgyQFiyQPJkD5MABqWNH0/7aa31YsKQpU8zt\nT88849v9AgBCV1CvJ3zaaVJhoee3B23ZYp79sXBDr17m1icAADwV1CFczdNBVv5cPWniRN/vEwAQ\n2kJiFSVP5eWZ3nNsrN2VAAAQ5CFcUeFde0ZGAwCcJKhDuHt36R//8Ly9v0O4sFD629+k77/33zEA\nAKEjaEPYsqS//EX6+c89a19ZKX3zjX9DuKJCeuwxadMm/x0DABA6gnZglssl3Xyz5+137ZKOHvXP\noKxqCQnmNqjIoP2qAgACKWh7wt7y9cIN9SGAAQCecnQIp6WlKTU1VfPnz2/2vvLypFatzHVkAACc\nwNH9toYWcJg3T2rbVrrqKs/2lZtrTkW3aOHDAutx5IgZnEXgAwAa4uiecEMWLZIWL/a8fW6u1Lu3\n/+qpbexYady4wBwLABC8HN0Tbsibb3o+XaVkQviSS/xXT23/8z9cGwYANC6oo8LT6SoPHpTy8wPX\nEx44MDDHAQAEt6A9He2N6pHRNi1hDABAnYIyhL1drSg31zz78x5hAAC8FZQh/PjjUs+enrfPzTUj\nlQO5cENenjR0qLRtW+COCQAILkF5Tfjyy6XERM/bB3JkdLXERBP8x44F9rgAgODhsixvxhgHRnFx\nseLi4lRUVFTvfcLe6NNHGjlSeuopHxQHAICPBOXpaG+Ul5uFGwLdEwYAoDEhH8LbtpnVjQhhAIDT\nBF0I799vBmbl53vWvnpktF0hvGKF9N//2nNsAICzOTqE61rAYedO6W9/MxNweCI3V2rXziwzaIcX\nX5See86eYwMAnC1oB2ZZlmczZo0fb4L48899W6OnSkqk1q09n90LABA+gvIWJcnzULPj9qTaAnlv\nMgAguDj6dHRzWZb9IQwAQH2CLoS9OXleUCAVFTkjhL//XqqqsrsKAICTBF0I9+1rBmZ5wu6R0dXW\nrpXi46UNG+ytAwDgLEF3TXjyZOmcczxrm5tr1vX1Zp5pfxgwQHrtNSkpyd46AADOEnQhfPPNnrfN\nzTUBHBXlv3o80bKldP319tYAAHCeoDsd7Q0GZQEAnIwQBgDAJkEVwqtWSXPneta2pMTMruWkEP7L\nX6S//tXuKgAAThFUIfzee9KMGZ613bLFPPfp4796vNWihRQRVF9xAIA/Bd20lVVVngXZvHnSdddJ\nP/4otW/vp0IBAGgGR/fL6lrAwdOe5MaNUteuBDAAwLkcfYvSggUL6l3AoTGbNpmJPQAAcCpH94Sb\nw6khfOCA9PLL0rFjdlcCALBb0ITwli1Sly7SunWNty0tlbZvl/r1839d3tqxQ7rpJunrr+2uBABg\nN0efjq6tXTtp0iSpW7fG2+bmmoUenNgTHjhQ2r9f6tjR7koAAHYLmhBOSJD+/GfP2m7aZJ49nWM6\nkCIiCGDMcUgqAAASrUlEQVQAgBE0p6O9sWmT1L271Lat3ZUAAFC/JoXwzJkzlZSUpJiYGCUnJ2td\nAxdq33rrLY0cOVJdunRRXFycLrzwQi1fvrzJBXti40ZnnoquzbKkw4ftrgIAYCevQ3jhwoWaMmWK\npk+frvXr16t///5KSUlRYWFhne1XrVqlkSNHasmSJfryyy81bNgwXXXVVdrg5eK6774rrV7tWVun\njoyu7YorpNtus7sKAICdvJ4xKzk5WUOGDNFTTz0lSbIsS926ddOkSZN03333ebSPfv36KS0tTX+u\n5yJvXTNmXXyxdPbZ0ksvNbzvw4fNaeiXX5bGj/f40wq4994zdV56qd2VAADs4tXArPLycmVlZen+\n++8/vs3lcmnEiBHKzMz0aB+WZenQoUPq0KGDV4WuXu3ZvbWbN5tnJ96eVNuVV9pdAQDAbl6dji4s\nLFRlZaXi4+PdtsfHxys/P9+jfTzxxBMqKSnRmDFjvDm0XC6pVavG21WPjHbSwg0AANQloLcozZs3\nTw8//LAWL16sTp06+eUYmzZJSUlSbKxfdg8AgM94FcKdOnVSixYtVFBQ4La9oKBACQkJDb52wYIF\nuuWWW7Ro0SINGzbMo+OlpaUpMtK9xPT0dKWnp9f7mmAYlFVt40bpwQelV17hdioACEdehXBUVJQG\nDRqkjIwMpaamSjLXeDMyMjRp0qR6Xzd//nzdfPPNWrhwoUaNGuXx8aoXcHjvPWnyZDNlZWPrOWzc\nKI0d6/EhbBUdLX3/vVRQQAgDQDjy+nT05MmTNX78eA0aNEiDBw/WjBkzVFpaqvE/DUWeOnWq9u7d\nq1dffVWSOQU9fvx4Pf3007rggguO96JjYmI8XiHp1FOlq65qPKiKi6Xdu4OnJ3zmmZ7fdgUACD1e\nh/CYMWNUWFioadOmqaCgQAMGDNCyZcvUuXNnSVJ+fr527959vP3s2bNVWVmpO+64Q3fcccfx7ePG\njdOcOXM8OubAgebRmOxs83zuuZ5/PgAA2MXr+4QDoa77hD0xa5Z0993mXmFPRlIDAGCnkJo7+uuv\nza1JwRbAX38tffKJ3VUAAAItKEL4zTfNesKN2bBB6t/f//X42vTpZpQ0ACC8OD6ELUu64QYzzWND\nqqrMNeFgDOFZs6QlS+yuAgAQaI5fT9jlkvbvNyHbkG3bpJKS4AzhEyYgAwCECceHsCTFxDTepnpR\npmAMYQBAeHL86WhPbdggJSRIXbrYXUnTlZVJ5eV2VwEACJSQCuHzzrO7iqbbs0fq1En64AO7KwEA\nBIrjQ/j556WUlMbbBevI6Gpdu0qPP85EIwAQThx9TTgtLU0//BCpU09Nl1T/og0HD0q7dgV3CLtc\n0p132l0FACCQHB3C1Qs4NIZBWQCAYOT409GeyMoyI6h797a7EgAAPBcSIfzFF2aBh0hH9+s989//\nSmPG2F0FACAQHB/C778v7dvXcJusLGnQoMDU42/VbyTKyuytAwDgf44O4dJS6corpYyM+tsUF5t5\npUMlhEeNkv71L88mKAEABDdHn8CNjpb27pXatKm/zZdfmudQCWEAQPhwdAhHREiJiQ23ycqSWrdm\nUBYAIPg4+nS0J7KypAEDQmNQVm0ZGcyeBQChLuij64svzHXUUPP00+Z0/OWX210JAMBfHN0TfvFF\ns5ZwfX78Udq6VTr//MDVFChvvCEtWGB3FQAAf3J0T7hdu4bX2l2zxjwPHRqYegKpbVu7KwAA+Juj\nQ3jMGBPE9fn8c6ljR+nMMwNXEwAAvuLo09FpaWlKTU3V/Pnz6/x4ZqaUnGwWPwhVe/dKBw7YXQUA\nwB8c3RNuaAGHqipzOvreewNcVACVlZle/sMPS1Om2F0NAMDXHN0T/vxzM/iqLrm5UlFRaF4PrhYT\nIy1eLN14o92VAAD8wdEhnJIiffhh3R/LzDSTeVxwQWBrCrQRI6RTTrG7CgCAPzj6dPSaNVKvXnV/\nbNUqs34wo4gBAMHK0T3h3r2l9u3r/tjHH0s//3lg67FTZaVkWXZXAQDwJUeHcH2+/VbauVO67DK7\nKwmMb7+VuneXPvnE7koAAL4UlCH88cfmtqRLLrG7ksDo3l0aP15KSLC7EgCALzk6hOu7/ejjj6Vz\nz5U6dAhsPXaJiJAeeUQ66yy7KwEA+JKjQ/jo0bq3r1wZPqeiAQChy9Eh/PTTJ2/bulXasUMaPjzw\n9QAA4EuODuG6LF0qRUVJv/iF3ZUE3rZt0siR0nff2V0JAMAXgi6ElywxA7LatLG7ksDr3NlM11lY\naHclAABfcHQI/+537gs4lJWZ68FXXGFvXXZp187MINa/v92VAAB8wZEzZh05Yp5/+9sFmjixZgGH\nlStNEI8aZU9dAAD4kiN7wpE/vTU48brvokVmVaG+fQNfEwAAvuboEI6Pr9lWXi699ZY0Zkxorx/s\niQMHzNKGu3bZXQkAoDkcGcJ1ycgwyxqOGWN3JfZr2VJ6+20pJ8fuSgAAzeHIa8J1mTfPzBh13nl2\nV2K/Nm3M/dIRQfMWCgBQF0f+GV+yxP3/+/ZJCxZIf/gDp6KrEcAAEPwc+ad8/373/z/9tBQdLd1y\niz31oH7Vt4/B2fg+BQ++V8HBV98nR4bwDTfU/PvQIem556SJE6W4OPtqciLLkl5/XfroI/tq4A9G\n/QoKpO+/d99WXm62l5cHtha+T8GD71VwCOkQrm3OHOnwYWnSJLsrcR6Xy3x9li2zu5LwsWyZtG6d\n+7aMDKlLl5MDd+xY6e673bdt3myWpFy/3n37//yPdNNN7tuOHZOeecasJw0gNDk6hPfulR59VEpL\nk7p1a/7+/PEO0+59LlkiPf647/drJ7u/ppK0e7c0c6ZUWem+/YEHpJdect9vjx4mbFu2dG/797+b\n9rWdfrr0zjtSr17u2888U+rXz73WwkLpnnuk3NyTa7joIvdtliX94x8nt/Unf/08OeH7b9c+/SVY\nPv9g+pnyFUeGcPU14bFjzWINTzzhm/0Gyw+NN/uMjvbPfu0U6K9pVZW5/a22rVul//f/TBjXlpEh\nzZrlvt+ePU0wtm/v3nbQoJpgrRYXJ6Wmntx23DhzvNq1nnqqmT1u5Ej3tsOHSzfe6L6ttFT685+l\nr75y3/6//2vW3j7R7NlSdvbJ270RTH8wg2Wf/hIsn38w/Uz5SsBvUbIsS4cOHar340ePSuecUyxJ\n2ry5WMuXS61bS8XFzT92RUWFin2xoyDcp7/2Gwr7/P3vzc/dv/9ds23AAHMmpmVL9589l8tcHrGr\nVkk6/3zzOLHJ3r3mDUXt7eeeK11/vbRyZc0+LUuaPFmaPt30zKu99JL08MNmqdDadyE8+qh5Q1H7\nzcC+fdLevRUqLCx2OwNQUmLeOJ94VsAbgfyaVlSYNzqxse6f8/79UosWUocONdtKSqTt282tktHR\nNfvMyTFnTGq/2SkpkT74QLrwQnOpotqaNeZrd/XVNduqqqQXXpCGDZPOPrtmv19+KWVlmbtCanvx\nRXOr5uDBNdtyc6X335fuvNP9a//GG+b4tT//PXukhQulCRPcP7933zVfj9Gja7YVFUmvvGK2de9e\ns/3jj6Vvv3X/mlZVSf/8p5npsE+fmrbr10tffnny5ZZXXzWzH55/fu3jVeiJJ4p1443m56ja229L\nHTuaxXuq5edL770n/eY30imn1GzPyDC1XH55zbYjRyr0zDPFuvxy8+a22tq1ZozGVVfVbKuqMrfE\nDh0q9exZsz0nR9q0SbrmGvN/T35O27ZtK1cjt/S4LMuyGmzhY8XFxYpjhBUAIMQVFRWpXbt2DbYJ\neAg31hOWpF27inXuud20e/fuRj8BGK++ak49fvaZ1Lat3dU407FjphczcaJ0//12VxPcjh4106cm\nJLj3ID/6yKz2Vbt3s2ePNHeu6dHV7nm99poZJV67h3TggPTQQ9Ltt0u9e9ds/89/pC1bpKlTa7ZZ\nlhnQNnq06alXW7dO+vxz6a673Gt+4w3TQ6vddtcu0zu9+mr3nte6dWZSnNo9ukOHTA3nnCPFxNRs\nLygwz7Wn2a2sNO3btKmZhre6ZuY6cA7LMo/a8y5YlukNu1zu2ysrzfbaPyeNcWRP2BPVvWVP3kXA\nqKgwf9Bq/3EId7m50hlnuJ+e27jRnPLz5hcJAPzFkQOz4L3ISAK4tt27zfWm2td5JTNQigAG4BRB\nM3c04I1u3cxAlcsus7sSAKgfPeEQ9PTT5vpXuNi3z9y2c+ItNykpUqtW9tQEAJ6gJxyC1q93H4Yf\n6jp1MgNgSkrsrgQAvENPOATNmSM98oj/9j99+nRFRES4Pc455xz/HfAEu3aZySmqRUWZWaiSkwNW\ngiOtXr1aqamp6tq1qyIiIrR48eKT2kybNk2nnnqqWrdurcsvv1zffPONDZWGt8a+TxMmTDjp9+uX\nv/ylTdWGr0cffVSDBw9Wu3btFB8fr9GjR2vLli0ntWvu7xQhHIICcQtEv379VFBQoPz8fOXn5+uT\nTz7x/0ElHTxoBlf9858BOVxQKSkp0YABAzRr1qw6b4t4/PHH9eyzz+qFF17Q2rVrFRsbq5SUFB07\ndsyGasNXY98nSbriiivcfr+cPONTqFq9erXuuusurVmzRh9++KHKy8s1cuRIlZWVHW/jk98py4GK\nioosSVZRUZHdpQS9ykrLOnLEt/t86KGHrIEDB/p2p15YutSyiottO3xQcLlc1jvvvOO2LTEx0Xry\nySeP/7+oqMiKjo62Fi5cGOjy8JO6vk/jx4+3Ro8ebVNFqM/+/fstl8tlrV69+vg2X/xO0RMOYZZl\npmObMsX3+966dau6du2qnj176ve//712nzjJso9kZZ284lBKChOSeGvHjh3Kz8/X8OHDj29r166d\nhgwZoszMTBsrQ11Wrlyp+Ph49e7dW7fffrsOHDhgd0lh7+DBg3K5XOrw04wzvvqdYmBWCHO5zApU\nXbv6dr/Jycl65ZVXdPbZZ2vfvn166KGHdOmll2rjxo2KjY312XEsS7r1VjNz0uuv+2y3YSk/P18u\nl0vxtad1khQfH6/8/HybqkJdrrjiCv32t79VUlKStm3bpqlTp+qXv/ylMjMzG519Cf5hWZb++Mc/\n6uKLLz4+/sVXv1OEcIi7/nrf7zMlJeX4v/v166fBgwfr9NNP17/+9S9NmDDBZ8dxuczE7Sf8jAMh\nbcyYMcf/3bdvX5177rnq2bOnVq5cqWHDhtlYWfi6/fbblZOTo08//dTn+3b06ei0tDSlpqYyKMHh\n4uLi1KtXr2aPtP3oI7OIfW1du7rPvYumSUhIkGVZKqie6PgnBQUFSkhIsKkqeCIpKUmdOnViJLtN\n7rzzTr3//vtauXKlEhMTj2/31e+Uo0N4wYIFWrx4sdLT0+0uJSTs2iXde+/JC9U31+HDh7Vt2za3\nH9Cm+OQT0/OtqvJRYTguKSlJCQkJysjIOL6tuLhYa9as0YUXXmhjZWjMnj179MMPPzT79wveu/PO\nO/XOO+9oxYoV6l57LUf57nfK0SEM39q0SXrzTamwsHn7uffee7Vq1Srt3LlTn332mUaPHq3IyMhm\nv1maOlX68EP3lUvguZKSEm3YsEFfffWVJGn79u3asGHD8UFzf/zjH/XXv/5V7777rrKzs3XDDTfo\ntNNO069//Ws7yw47DX2fSkpKdN9992nNmjXauXOnMjIydPXVV6tXr15ul4Hgf7fffrvmzp2refPm\nKTY2VgUFBSooKNCRI0eOt/HJ75TvBnD7Drco+U9ZWfP3kZaWZnXt2tWKjo62unXrZqWnp1vbt2/3\nah8ffWRZo0ZZ1rFjza8HxsqVKy2Xy2VFRES4PSZMmHC8zYMPPmglJiZaMTEx1siRI62tW7faWHF4\nauj7VFZWZqWkpFjx8fFWq1atrKSkJOvWW2+1vv/+e7vLDjt1fY8iIiKsV1991a1dc3+nWMoQtli/\nXvrLX6TZs820kwAQjjjxF8ays6Vf/Ur68Uf/H+vE2d4GDpTeeosABhDeCOEwdviwVFYmtW7t3+Ms\nXmzu9c3N9e9xACDYcDo6zFmW+1zTVVXNHxh19Kj7EoLHjklLl5peN4OuAKAGfxLD3IkT8PzlLyYs\nm/rWbN486YwzTA+7WsuWUmoqAQwAJ+LPItycf740cqR7OB84YK7fnrhe78yZ0ty57tt+9jNzLzL3\n+gJA45iLCG5+9auTt61eLf3mN9J330m1p4Zeu1ZKSJCuu65mW+/e5gEAaBzXhNEoy5L27zcjmTml\nDAC+Q08YjXK5pC5d7K4CAEKPo/s1LOAAAAhlnI4GAMAmju4JAwAQyghhAABsQggDAGATR14TtixL\nhw4dUtu2beU6cUonAABChCNDGACAcMDpaAAAbEIIAwBgE0IYAACbEMIAANiEEAYAwCaEMAAANiGE\nAQCwyf8H9QegKx4IRygAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt_1 = plot(ip, [t, 0, 20], linestyle ='-')\n", "plt_2 = plot(diff(ip,t), [t, 0, 20], linestyle =':')\n", "(plt_1 + plt_2).show(figsize=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

微分方程式の数値解法

\n", "\t

desolve_system_rk4を使った解法

\n", "\t

\n", "\t\tdesolve_system_rk4は、ルンゲクッタ法を使って効率よく微分方程式の数値解をもとめる関数です。\n", "\t

\n", "\t

\n", "\t\t例題の微分方程式を$\\ddot z, \\dot z$について見てみます。\n", "$$\n", "\t\t\\ddot z + \\dot z + z = u(t)\n", "$$\n", "\t\t以下のような等価な連立微分方式を得ることができます。\n", "$$\n", "\t\t\\left\\{\\begin{array}{l}\n", "\t\t\t\\dot z = \\dot z\\\\\n", "\t\t\t\\ddot z = - z - \\dot z + u(t)\n", "\t\t\t\\end{array} \\right.\n", "$$\n", "\t

\n", "\t

\n", "\t\tこれをdesolve_system_rk4関数を使って定義すると以下のようになります。\n", "\t\t初期条件icsは$[ t_0, z(t_0), \\dot z(t_0)]$として解きます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFnCAYAAACRo/HLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl8VNX9//H3hAQCAcISCIsgEUXWEkUBcakIsmiJ8lUx\nUb8srtSF9geFPtAWpbbfb61WqhWkaN0qm9IK+FVAjaCA7AICgQGEAgoZSIEEApGQ3N8fx5AMZJlJ\nZubemXk9H495THJz595PmIR3zrnnnuOyLMsSAAAIuRi7CwAAIFoRwgAA2IQQBgDAJoQwAAA2IYQB\nALAJIQwAgE0IYQAAbEIIAwBgE0eGsGVZysvLE/OIAAAimSND+MSJE0pMTNSJEyfsLgUAgKBxZAgD\nABANCGEAAGxCCAMAYBNCGAAAm8TaXQBqLidHOnxYSkiQmjSRGjSwuyIAgC9oCYehoiJpyRLp7rul\n5s2lZs2kLl2kdu2khg2lzp2l0aOl5csl7vICAOdyWQ68GTcvL0+JiYnKzc1Vw4YN7S7HUVaskB55\nRMrKMsF7xx1S165Sy5bS6dPSoUPSqlUmpPfulTp2lH7zGykjQ4rhTy4AcBRCOEwUFkrjxkl//avU\nq5f04ovSNddILlf5+xcXS8uWSS+9JC1cKKWmSq+8Il17bUjLBgBUgrZRGMjPl26/XXr1Vekvf5FW\nrpT69Kk4gCXT6r3pJmnBArN/7drS9debID99OnS1AwAq5uiW8ODBgxUbG6uMjAxlZGTYXZYtTp6U\n+veXtm2T/vlPacCA6h2nqMi0nn/7W+nyy82xLr00sLUCAPzj6BCO9u7ooiJp6FBp6VLzuOqqmh9z\nyxZzHdnjkWbOlH72s5ofEwBQPXRHO9iECdJHH0lz5wYmgCWpWzdp3TrTVZ2WZq4xAwDsQQg71Acf\nmO7jKVOkW24J7LETE0139Lhx0pgx5tl5/SEAEPmYrMOBjhwxtyHddpv0xBPBOUdMjPT889LFF5tz\n5OdL06ZxGxMAhBIh7DCWZSbaKC6W/va3ykdAB8Ljj5uZth54QPrhB+n116VatYJ7TgCAQQg7zIIF\n0r/+Jb33npScHJpzjhol1akjDR9ugvidd6RYfjIAIOj4r9ZBzpyRxo83tyHddVdoz33PPSaI09NN\nHbNnS3Fxoa0BAKINVwAdZPp0ac8e6YUX7Dn/HXeYVviCBdL995sucQBA8BDCDnH8uDR5sgm/bt3s\nq2PIEHP/8MyZZuQ0o6YBIHjojnaIKVOkggLpd7+zuxJp2DApL0966CGpUSPp97+3uyIAiEyEsAOc\nOGEmzXj4YbMakhM8+KCUmyv96lfmvuLx4+2uCAAiDyHsADNmmDmix461uxJv48aZbvIJE0wQP/yw\n3RUBQGRxdAinp6dH/AIOP/xgZsa67z6pTRu7q7nQ735nWsSjR0sNG5rR0wCAwHB0CM+ZMyfiF3B4\n913p0CHT2nQil8ssn5ibK/33f0sNGki33mp3VQAQGVhFyUaWJfXoIbVuLX34od3VVO7sWXPv8uLF\n5vHTn9pdEQCEP25RstG6ddLGjaar1+liY80EHn36mNuYNmywuyIACH+EsI2mT5fatpUGDbK7Et/E\nx0vz50udOpmad+ywuyIACG+EsE2OHZPmzDEjjsNpwYQGDaRFi8y81jffLO3bZ3dFABC+CGGb/OMf\nUmGhmSEr3DRpIn3yiVS7ttS/v+Tx2F0RAIQnQtgmb71lrq06ZXIOf7VqJX36qVmHeOBAcz8xAMA/\nhLANtm83A7Luu8/uSmrmkktMi3j/fnPbUn6+3RUBQHghhG0wa5aZgeqWW+yupOa6djXXiDdvNqsw\nnTljd0UAED4I4RCzLBPCd95pRhtHgl69zPKHS5dK995r7ikGAFSNEA6x1avNmsH33mt3JYHVr5/0\n3nvmFqZ77jGDzgAAlSOEQ2zmTDNDViTOOHXbbdK8eSaIMzIIYgCoiqNDOD09XWlpaZo9e7bdpQRE\ncbH0z3+a9XpjHP0vX3233Wa+x4ULpbvv5hoxAFSGuaNDaNUqM+3j8uXSddfZXU1wffSR9F//JQ0e\nLM2dK9WpY3dFAOA8Edoec6YPPpCaNZOuucbuSoLv1ltNt/SSJSaIc3PtrggAnIcQDhHLMiF8223h\nNU1lTQwebO4j3rhRuvFGKTvb7ooAwFkI4RDJypJ275Zuv93uSkLr+utN9/uRI6YrftcuuysCAOcg\nhEPkgw+k+vXNrTzRpmtX6auvzHXhPn3MtXEAACEcMvPnmxmyImWCDn+1bSutWCFdfrm5Peu11+yu\nCADsF5IQXr58udLS0tS6dWvFxMRo4cKFoTitY2RnSxs2SD/7md2V2KtpU+nzz6UHHzRLOD7yiPTD\nD3ZXBQD2CUkI5+fnKzU1VdOmTZPL5QrFKR3lk0/M88CB9tbhBLVrS9OmSa+/blaSuvFG6fvv7a4K\nAOwRG4qTDBo0SIMGDZIkOfC25KBbtEjq0UNq3tzuSpzjgQekbt3MvcQ/+Yn06qtmEhMAiCZcEw6y\noiLTEh482O5KnKdnT2nTJjNY7e67zVSX//mP3VUBQOiEpCUczdavl44elX7sCMB5kpLMjFpDh0qP\nPWZGUr/2GtfPa+rMGenwYSknx0yUcvy49yMvTzp9WiooKH0ueZR8Xlho7m8//1FcbJ4lc897bKxv\nj9q1zSMuLvQfR+o0sQh/hHCQLVokNWpklvtD+Vwu0wr+6U/NoK0hQ8wfLS+8IHXpYnd1zlJcbMJ1\n/37z2LfPPGdnSx5P6ePYsfJfX6+eWcu6YUPzcd26ZsR+fLz5vEmT0s/j4sx743KZEDv/Y8n09Jw9\nW/pc0ePMGenkSRPsZ86Yhy8fFxUF5t+tVq2ah3mtWqX/BtKFH/u67fyvS6V/1JS9WufLtup+LZDH\nirTz+Lr/ww+bFeNqytEhnJ6erthY7xIzMjKUkZFhU0X+W7xYuvlm0xJA5Vq1MnNOz58vjR9vrhU/\n/LA0eXL0XU8/eVLavt1M8rJtm3l2u6UDB7xHlCckmNu/WrUyq3NdeaWUnFz6SEqSGjc2wZuYaAIl\nnBQXl4ayvwEeiI9PnTI9CT/84N0DUNIrUPKxr9vK+3p5oV3Cl23V/VogjxXq81T0x0woa27QQAER\n8gUcYmJiNH/+fKWlpVW4T6Qs4HDsmLkt57XXzEAk+O7MGemVV6Tf/c785/fQQ9IvfmECJ5Lk5Zmw\nLQnakse+faX7tGsnde4sdexoPr74YvPvcPHFppclCm84ACJGSNpn+fn52r1797mR0Xv27NHmzZvV\npEkTtWnTJhQl2OKLL8xfutE4S1ZN1a4tjR0rDR9uuqX/9jfppZekO+8023v2tLtC/xw/7h2yJS3c\n774zX3e5pJQU0/2enm5Ct0sXE7wJCfbWDiB4QtIS/uKLL9S3b98L7hEeMWKE3njjjQv2j5SW8Jgx\n0ocfSnv32l1J+Dt50txX/Je/SN9+awZwDRsm3XWXCSqnOHr0wqDNypIOHjRfj4mR2rc3IVsStJ07\nm5nE6tWzt3YAocd6wkHUrZtpsf3973ZXEjmKiqSPP5bmzJEWLjTh3K2bGU197bVmmcgmTYJfw/79\n0s6d5uF2l16/LVkpqlYt6dJLvYO2c2epQwczGAoAJEI4aA4fNgNj3n1Xuvdeu6uJTKdPm/WK339f\nysw0o4IlE3Y9e5oQvOQS0/JMSTHXT+PiKj9mUZHpOj561LyH331nBkMdOFAavLt3m2vWklmUon17\nqVOnC8O2Tp3gfv8Awh8hHCTvvWcmoPj+ezNyFcFlWdKePWa1ppUrzRrG33574eQfdeua23MaNDBd\nwyW30BQVmVZ1bu6Fx65fX2rTxgyGuuwy03XcoYN5tGkTPetDAwg8bpwJks8/N/9ZE8Ch4XKZFmn7\n9tJ//3fp9txcE8b79pmP8/JKnyVz61itWuaRkGC6skseTZuakE1MZAQygOAghINk6VLpppvsrgKJ\niebe2SuvtLsSALgQk7kFwfffm2uHffvaXQkAwMkI4SBYscI833CDvXUAAJyNEA6CFSvMyNwWLeyu\nBADgZIRwEKxYIV13nd1VAACcztEhnJ6errS0NM2ePdvuUnyWmyt98w0hDAComqNHR8+ZMyfs7hNe\nvdosOEAIAwCq4uiWcDhascIsH9ehg92VAACcjhAOsJLrwUzuAACoCiEcQGfOSGvW0BUNAPANIRxA\nGzeaRQWuvdbuSgAA4YAQDqDVq83KOUyRCADwBSEcQGvWmACuXdvuSgAA4YAQDqDVq6VeveyuAgAQ\nLgjhADlyRNq7lxAGAPiOEA6QNWvMc+/e9tYBAAgfhHCArFkjNW8uXXyx3ZUAAMIFIRwgJdeDmaQD\nAOArR4dwuCzgUFwsrV1LVzQAwD8s4BAAbreUl8egLACAfxzdEg4Xa9aYbuirr7a7EgBAOCGEA2D9\neunyy6UwaLQDAByEEA6ADRukHj3srgIAEG4I4RoqLJQ2bZKuusruSgAA4YYQrqGsLKmggJYwAMB/\nhHANbdhgBmVdcYXdlQAAwg0hXEPr10udOkn169tdCQAg3BDCNcSgLABAdRHCNVBYKG3ezKAsAED1\nEMI1sG2b9MMPtIQBANVDCNfA+vVSTIyUmmp3JQCAcOToEHb6Ag5ffy117CglJNhdCQAgHLGAQw1s\n2sStSQCA6nN0S9jJioqkb76hKxoAUH2EcDXt3i3l59MSBgBUHyFcTZs2mefu3e2tAwAQvgjhatq0\nSbroIikpye5KAADhihCupo0b6YoGANRMtUJ46tSpSklJUd26ddW7d2+tW7eu0v1nzpyp1NRUJSQk\nqFWrVnrggQd09OjRahXsFJs2MSgLAFAzfofw3LlzNW7cOE2ePFkbN25U9+7dNXDgQOXk5JS7/8qV\nKzVixAg99NBDysrK0rx587R27Vo9/PDDNS7eLtnZksdDCAMAasbvEJ4yZYoeeeQRDR8+XB07dtT0\n6dNVr149vfHGG+Xuv3r1aqWkpOixxx7TxRdfrD59+uiRRx7R2rVra1y8XUoGZdEdDQCoCb9CuLCw\nUBs2bFC/fv3ObXO5XOrfv79WrVpV7muuueYaHThwQIsWLZIkeTwezZs3T7feemsNyrbXxo1Sw4ZS\nu3Z2VwIACGd+hXBOTo6KioqUnJzstT05OVnZ2dnlvqZPnz569913dffdd6t27dpq2bKlGjVqpFde\neaX6Vdts82bTFe1y2V0JACCcBX10dFZWln7xi1/omWee0ddff60lS5Zo7969euSRR4J96qD55hvp\nJz+xuwoAQLjza+7opKQk1apVSx6Px2u7x+NRixYtyn3NH//4R1133XUaO3asJKlr166aNm2arr/+\nev3hD3+4oFVdVnp6umJjvUvMyMhQRkaGP2UHVEGB5HZL/+//2VYCACBC+BXCcXFx6tGjhzIzM5WW\nliZJsixLmZmZGjNmTLmvOXXqlOLi4ry2xcTEyOVyybKsSs/nxAUcsrKk4mKpWze7KwEAhDu/u6PH\njh2r1157Te+884527Nih0aNH69SpUxo5cqQkaeLEiRoxYsS5/YcMGaJ//etfmj59uvbu3auVK1fq\nF7/4hXr16lVh69nJtmwxz1272lsHACD8+b2U4bBhw5STk6NJkybJ4/EoNTVVS5YsUbNmzSRJ2dnZ\nOnDgwLn9R4wYoZMnT2rq1Kn61a9+pUaNGqlfv3764x//GLjvIoS++UZq316qX9/uSgAA4c5lVdUn\nbIO8vDwlJiYqNzfXcd3RN99sAviDD+yuBAAQ7pg72k9btjAyGgAQGISwHzwe82BQFgAgEAhhP5QM\nyqIlDAAIBELYD998I9WtawZmAQBQU4SwH7Zskbp0kWrVsrsSAEAkIIT9sGUL14MBAIFDCPuouFja\nvt20hAEACARC2Ef79kmnThHCAIDAcXQIp6enKy0tTbNnz7a7FG3bZp47d7a3DgBA5PB72spQctIC\nDtu2SQ0aSG3a2F0JACBSOLol7CRZWaYV7HLZXQkAIFIQwj7ato2uaABAYBHCPmBkNAAgGAhhH5SM\njKYlDAAIJELYB1lZ5pmWMAAgkAhhH2zbZtYQZmQ0ACCQCGEflAzKYmQ0ACCQCGEfZGXRFQ0ACDxC\nuArFxaX3CAMAEEiEcBX272fOaABAcBDCVSiZM5oQBgAEmqND2AkLOGRlMTIaABAcLOBQBUZGAwCC\nxdEtYSdgUBYAIFgI4UqUjIzmejAAIBgI4Urs3y/l5xPCAIDgIIQrUTJnNN3RAIBgIIQrUTJndNu2\ndlcCAIhEhHAlsrKkTp0YGQ0ACA5CuBLbtnE9GAAQPIRwBSxL2r6d68EAgOAhhCvw3XfSyZOmOxoA\ngGAghCvgdpvnyy+3tw4AQOQihCuwc6cUFyelpNhdCQAgUjk6hO1cwMHtltq3l2IdPbs2ACCcOTpi\n7FzAwe2mKxoAEFyObgnbiRAGAAQbIVyO06elffsIYQBAcBHC5di929wnTAgDAIKJEC4HtycBAEKB\nEC6H2y01aSIlJdldCQAgklUrhKdOnaqUlBTVrVtXvXv31rp16yrd/8yZM3rqqafUrl07xcfH65JL\nLtFbb71VnVOHBIOyAACh4PctSnPnztW4ceM0Y8YM9ezZU1OmTNHAgQO1c+dOJVXQdLzrrrt05MgR\nvfnmm2rfvr0OHTqk4uLiGhcfLG43c0YDAILPZVmW5c8LevfurV69eumll16SJFmWpTZt2mjMmDGa\nMGHCBfsvXrxY99xzj/bs2aNGjRr5dI68vDwlJiYqNzc35PcJW5bUuLH0619LEyeG9NQAgCjjV3d0\nYWGhNmzYoH79+p3b5nK51L9/f61atarc13z44Ye66qqr9Nxzz+miiy7S5ZdfrvHjx6ugoKBmlQfJ\n4cNSbi7d0QCA4POrOzonJ0dFRUVKTk722p6cnCx3yZDi8+zZs0fLly9XfHy85s+fr5ycHP385z/X\n0aNH9fe//736lQcJI6MBAKES9Gkri4uLFRMTo1mzZql+/fqSpBdffFF33XWXpk2bpjp16gS7BL+4\n3VJMjHTppXZXAgCIdH6FcFJSkmrVqiWPx+O13ePxqEWLFuW+pmXLlmrduvW5AJakTp06ybIsfffd\nd2rfvn2F50tPT1fseSsoZGRkKCMjw5+y/eJ2S+3aSQ772wAAEIH8CuG4uDj16NFDmZmZSktLk2QG\nZmVmZmrMmDHlvubaa6/VvHnzdOrUKdWrV0+S5Ha7FRMTo4suuqjS89mxgAO3JwEAQsXv+4THjh2r\n1157Te+884527Nih0aNH69SpUxo5cqQkaeLEiRoxYsS5/e+55x41bdpUo0aN0vbt2/Xll19qwoQJ\neuCBBxzXFS0RwgCA0PH7mvCwYcOUk5OjSZMmyePxKDU1VUuWLFGzZs0kSdnZ2Tpw4MC5/RMSEvTp\np5/qiSee0NVXX62mTZvq7rvv1rPPPhu47yJAzpyR9uwhhAEAoeH3fcKhYNd9wm631LGj9PnnUt++\nITstACBKMXd0GdyeBAAIJUK4DLdbql9fatnS7koAANGAEC6jZFCWy2V3JQCAaEAIl8HIaABAKBHC\nZRDCAIBQIoR/dOyYdOQIIQwACB1C+EeMjAYAhBoh/KOSEL7sMnvrAABED0eHcHp6utLS0jR79uyA\nHjc/33Q/l+V2S23aSAkJAT0VAAAVCvpShjURjAUcioulbt2kO++U/vSn0u07dtAVDQAILUeHcDDE\nxEivvCJ17ly6LT9fWrxYuu46++oCAESfqAthSbrlFu/Pa9UyizekptpTDwAgOjn6mnCgnDljWrsV\nOXhQKiqS+vULXU0AAERFCL/5phn1fOJE+V/n9iQAgB2iojt64EDT5dygQflfd7ul+HipbVvz+ZYt\nZvT0DTeErkYAQPSJihBu10568MGKv+52m5ZyzI/9Ak89ZbqwCWEAQDBFRQhX5fw5o19/XWrc2L56\nAADRIaKvCRcW+rbfzp3eIdy8uRQXF5yaAAAoEdEh/POfS8OGVb7PyZPS998zKAsAEHoRHcK33SYN\nHVr5Pjt3mufyQjg3V/rb38wsWwAABFpEXxMeMqTqfSq7PWn7dumJJ6SePaUrrghsbQAAODqE09PT\nFRsbq4yMDGVkZATlHG63lJwsJSZe+LVevaTvvjPXiAEACDRHh3AwFnA43/kjo8tyuQhgAEDwROQ1\n4eXLzX3BlU1VWaKyEAYAIJgiMoSPHpX27pXq1at8P8u68Pak8hQXS198If3wQ+BqBAAgIkP4ttuk\nzEzTnVyZ7783reWqQnj3bunGG81yhwAABIqjrwkHm68LN3ToIK1fL115ZfBrAgBEj6gP4bg4KSWl\n6n179Ah+PQCA6BJx3dHLl0v79/u2r9sttW8vxUb1nyIAALtEXAiPGiX9+c++7et2m65mfxw/7n9N\nAACUJ+JCeP16aeJE3/b19/akFSukZs1Kp7oEAKAmIq4jtlEj3/Y7fVrat8+/EO7RQ3r5ZRPEAADU\nVMSFsK927zb3CfsTwnXrmpWZAAAIhIjrjvaVr7cnAQAQLI4O4fT0dKWlpWn27NlV7nv2rHTZZdL8\n+b4d2+2WGjeWkpJqWCQAANXk6O5ofxZwOH1auusu6ZJLfDt2yXSVVc2qVZ5586Q335Q++sj/1wIA\nUMLRIeyPBg2k//kf3/d3u6WOHat3roYNpdatzVzSdepU7xgAADi6OzpYLKtmqycNGCDNmEEAAwBq\nJipD+MgRM+kGg7IAAHaKiBA+dcrcv5ud7dv+jIwGADhBRISw2y2NHy8dOuT7/jEx0qWX1uy8778v\nffZZzY4BAIhe1QrhqVOnKiUlRXXr1lXv3r21bt06n163cuVKxcXF6coArwl4xRVmXeDu3X3b3+2W\n2rWr+TXdV1+VFiyo2TEAANHL7xCeO3euxo0bp8mTJ2vjxo3q3r27Bg4cqJycnEpfl5ubqxEjRqh/\n//7VLrYysbGmdeuLmgzKKuvjj6W//rXmxwEARCe/Q3jKlCl65JFHNHz4cHXs2FHTp09XvXr19MYb\nb1T6utGjR+vee+9V7969q11soAQqhOPja34MAED08iuECwsLtWHDBvXr1+/cNpfLpf79+2vVqlUV\nvu7NN9/U3r179fTTT1e/0gpYllRc7Pv+hYXSnj0MygIA2M+vEM7JyVFRUZGSk5O9ticnJyu7gqHJ\nu3bt0pNPPqmZM2cqxtf+Yj9kZ5vJM5Yu9W3/PXvMFJeBCmHLknbtMs8AAPgjqKOji4uLde+992ry\n5Mlq3769JMkKcFrVri1Nnux7qAb69qQvvpA6dJC2bg3M8QAA0cOvaSuTkpJUq1YteTwer+0ej0ct\nWrS4YP8TJ05o/fr12rRpkx577DFJJpgty1Lt2rX1ySef6MYbb6zwfOnp6YqN9S4xIyNDGRkZ5z5v\n2lQaN87378HtlurXl1q29P01lbnmGjNC+se/MQAA8JlfIRwXF6cePXooMzNTaWlpkkzLNjMzU2PG\njLlg/4YNG2rreU3EqVOnaunSpfrnP/+pdu3aVXo+fxZw8JXbbVqu1Vm4oTx16kg//lMAAOAXvxdw\nGDt2rEaOHKkePXqoZ8+emjJlik6dOqWRI0dKkiZOnKiDBw/q7bfflsvlUufOnb1e37x5c8XHx6tT\np04B+Qb8FaiR0QAA1JTfITxs2DDl5ORo0qRJ8ng8Sk1N1ZIlS9SsWTNJUnZ2tg4cOBDwQstTXGxW\nTho2zLRufeF2S0G6VRkAAL+4rECPlAqAvLw8JSYmKjc3t9Lu6OxsqVs36R//kAYNqvq4x45JTZpI\ns2dL6ekBLFjSU09JJ09KL70U2OMCACJXWK8n3KKFWRHJ1z8jgrlww0UXSQUFgT8uACByhXUIl/B1\nkFVJCPvade2Pn/888McEAES2iFhFyVc7d5oWa0KC3ZUAABDmIXz2rH/7MzIaAOAkYR3Cl1wiPf+8\n7/sHO4SPHZP+/Gfp8OHgnQMAEDnCNoQtS3rmGemmm3zbv6jIzPEczBAuLJQmTZI2bw7eOQAAkSNs\nB2a5XNL99/u+//790g8/BDeEmzc3reHatYN3DgBA5AjblrC/gnl7UlkEMADAV44O4fT0dKWlpWn2\n7Nk1PpbbLcXHS23bBqAwAAACwNHd0ZUt4DBnjlk8YehQ347ldkuXXSYFYUnjCxQWmsFZrVsH/1wA\ngPDl6JZwZRYsMA9fhfL2pPvuk+65JzTnAgCEL0e3hCsze7bv01VK0o4d/g3kqonx4wO3VCIAIHKF\nbQhLvgddbq508KDUsWNw6ylx1VWhOQ8AILyFbXe0P0pGRtu0hDEAAOUKyxAuLPSvK3r7dvPMlJUA\nACcJyxB+7jkz0tlXO3aYW5NCuXDDt99KN95oZukCAKA8YXlNeOBA/+733bEjdNeDSyQnS40ascYw\nAKBiYRnCV19tHr7avl0aNCh49ZSnfn1p/vzQnhMAEF7CsjvaH4WFpms41C1hAACqEvEhvHu3WXeY\nkdEAAKcJuxDOyZH++Efp0CHf9t+xwzzb1RJeuVL6+GN7zg0AcDZHh3B5Czjs3y/96U9myUBfbN8u\nNW5slhm0w4wZ0osv2nNuAICzuSzLnztuQyMvL0+JiYnKzc2tcAEHy/Jtxqzhw02X9FdfBbhIH504\nYQZpMY0lAOB8jm4JV8bXUNu+3d5BWQ0aEMAAgPKFbQj7wrLMNWEGZQEAnCjsQtifzvPvv5dOnnTG\n7UnHjklFRXZXAQBwkrAL4S5dpD/8wbd9S0ZG290S3rhRSkqSvv7a3joAAM4SdjNmjRsnde7s277b\nt0u1a0vt2gW1pCp17Sr97W9SSoq9dQAAnCXsQviBB3zfd8cOqUMHKdbm7zIuTnrwQXtrAAA4T9h1\nR/vD7pHRAABUJqJDmJHRAAAnC6sQXrZMmjnTt31zc83Ulk5qCT//vPTss3ZXAQBwirAK4Y8+kl5+\n2bd9nTIyuqwzZ8wDAAApDKet9HW6yrfekkaNMvcJJyQEp04AAGrC0S3h8hZw8HUKyB07pIsvJoAB\nAM7l6FuU5syZU+ECDlVhZDQAwOkc3RKuiR07nBnCeXnSrFlSQYHdlQAA7BY2Ibxzp5n6cd26qvc9\nc0b69ls/xzxWAAASyklEQVRnDcoqsW+fdO+9TGEJAHB4d3RZDRtKY8dKbdtWva/bbRZL6NIl+HX5\nq2tX6eBBqWVLuysBANgtbEK4RQvpySd923frVvPsxBB2uQhgAIARNt3R/ti2TWrVSmrc2O5KAACo\nWLVCeOrUqUpJSVHdunXVu3dvravkQu0HH3ygAQMGqHnz5kpMTFSfPn30ySefVLtgX2zdarp9ne70\nabsrAADYye8Qnjt3rsaNG6fJkydr48aN6t69uwYOHKicnJxy9//yyy81YMAALVq0SF9//bX69u2r\nIUOGaPPmzX6dd8ECafly3/bdts2ZXdFl3XGHmUwEABC9/J4xq3fv3urVq5deeuklSZJlWWrTpo3G\njBmjCRMm+HSMrl27Kj09Xb/5zW/K/Xp5M2b16SN162bW5a3MqVNS/frS669L99/v+/cVah99ZNY6\nvvlmuysBANjFr4FZhYWF2rBhg54sM0LK5XKpf//+WrVqlU/HsCxLJ06cUJMmTfwqdOVKqbCw6v22\nbzdTWzq9O/rWW+2uAABgN7+6o3NyclRUVKTk5GSv7cnJycrOzvbpGM8//7zy8/M1bNgwf04tl8u0\nHKuybZt57tzZr8MDABByIb1FadasWXr22We1cOFCJSUlBeUcW7dK7dqZLmkAAJzMrxBOSkpSrVq1\n5PF4vLZ7PB61aNGi0tfOmTNHDz/8sObNm6e+ffv6dL709HTFxnqXmJGRoYyMjApfEw6Dskrs3Ck9\n/bQ0bRq3UwFANPIrhOPi4tSjRw9lZmYqLS1NkrnGm5mZqTFjxlT4utmzZ+vBBx/U3LlzNWjQIJ/P\nV7KAw7/+Jf3619I330h161b+mq1bpUoy2lHi46Xdu6VDhwhhAIhGfndHjx07ViNHjlSPHj3Us2dP\nTZkyRadOndLIkSMlSRMnTtTBgwf19ttvSzJd0CNHjtTLL7+sq6+++lwrum7duj6vkNSunXTnnVUH\ncF6etH+/8wdllWjb1re5sAEAkcnvEB42bJhycnI0adIkeTwepaamasmSJWrWrJkkKTs7WwcOHDi3\n/2uvvaaioiI99thjeuyxx85tHzFihN544w2fznnlleZRlZLpKsMlhAEA0c3v+4RDobz7hH3x6qvS\nmDHSyZNSnTpBLBAAgACIqLmjN282yxeGWwC73eY+aABAdAmLEH7/fTOSuCqbN0vduwe/nkB7+mnJ\nx8nGAAARxPEhbFnS8OHS4sWV71dcLG3ZEp4h/OKL0mef2V0FACDUHL+esMslHT9uQrYy334r5eeH\nZwi3amV3BQAAOzg+hCXfrvGWLMr0k58EtxYAAALF8d3Rvtq8WUpONo9wdeaMb4tUAAAiQ0SFcDh2\nRZc4ckRq1swscQgAiA6OD+HXX5duuqnq/cI9hJs1k559lu50AIgmjr4mnJ6erqNHY5WcnCGp4gmh\njx4101WGcwhLZqIRAED0cHQIlyzgUJVvvjHP4R7CAIDo4vjuaF9s2GAWd+jY0e5KAADwXUSE8Pr1\n0hVXSLGObtf7ZtEiacQIu6sAAISC40P444+lgwcr32f9eumqq0JTT7CdPSsdOyYVFNhdCQAg2Bwd\nwgUF0q23Sp9+WvE+x45Ju3dHTggPGSItXCjFx9tdCQAg2BzdgVunjnTokJSQUPE+X39tniMlhAEA\n0cPRIexySS1aVL7P+vVS/fpShw6hqQkAgEBxdHe0L9avl668UqpVy+5KAmvFCmnJErurAAAEU0SE\ncCR2Rf/1r9LLL9tdBQAgmBzdHf33v5t7gN99t/yv5+RI//53ZIbwjBmSD/OUAADCmKNDuGFDqWXL\nir++YYN5jsQQTky0uwIAQLA5OoTvuqvy1uDq1VLjxlL79qGrCQCAQHH0NeH09HSlpaVp9uzZ5X59\nxQrp2mulGEd/FzWTk2MWqAAARB5Ht4QrW8Dh7FnTEn7qqRAXFUKFhdKll0oTJkhPPml3NQCAQHN0\nG3LVqopbgd98I508KV13XWhrCqW4OOn996UHH7S7EgBAMDg6hAcNkj77rPyvrVgh1a4dmYOyyrr5\nZql5c7urAAAEg6O7o9eulS67rPyvrVwpXX01cywDAMKXo1vCl18uNWp04XbLKh2UFS2Ki833DQCI\nHI4O4Yrs22eWN4zk68Fl7dtnbsNavdruSgAAgRSWIbxihXnu08feOkKlTRvpv/5LatLE7koAAIHk\n6BD+5S/L3750qdSli9S0aWjrsUtMjPTnP5vueQBA5HB0CJ89e+E2y5IWLzYjpwEACGeODuFXXrlw\n25Yt5nowIQwACHeODuHyLFok1asnXX+93ZWE3r//LQ0ZIh04YHclAIBACLsQXrxYuukmqU4duysJ\nvaZNpbw86fBhuysBAASCo0P4zju9F3DIyzMjo6O1K7pBA+mLL6QePeyuBAAQCI6cMatkQNaQIXP0\nxBOlCzh8/rn52uDBNhUGAEAAObIl7HKZ54EDvbd/9JFZVeiSS0JfEwAAgebIEK5Vyzy3aFG6raDA\nrCh011321OQkx4+bJRz37bO7EgBATTgyhMuzYIGUmyuNGGF3JfaLjZX+8Q+znCMAIHw58ppwed5+\nW7rmGmaNkqT69aW9e0t7DAAA4cmRLeHFi70/P3hQWrJEGjnSlnIciQAGgPDnyBD+z3+8P3/3XSku\nTho2zJ56ULGS28fgbLxP4YP3KjwE6n1yZAjfe2/px6dPS9OmmVWEyltbOJpZlhms9tVX9tXAfxgV\n+/576bvvvLedPWsmWyksDG0tvE/hg/cqPER0CJf13HPSoUPS5Ml2V+JML7wgzZ9vdxXR4/PPpXXr\nvLctWyY1a2Yum5R1//0XrgS2c6eUnCytXeu9fdIk6b77vLedPSu9/LK0Z4/3dsuqdvkAHMbRIbxv\nnwnhsWOlyy6r+fGC8Remncd0uaTMTOlPfwrsce3mhPfpwAHppZcubLE++aQ0fbr3cdu1M2EbH++9\n7wsvSL//vfe2Nm2khQulTp28t3fqVDoTWkmtR49Kv/61lJXlve9vfytdccWFNb/wglngpKxgBnaw\nfp6c8P7bdcxgCZfvP5x+pgLGcqBdu3ItSVbfvrlW69aWdeJEYI47ZMiQwBwoDI8ZrONGwjGLiiwr\nJ8d725dfWladOpa1a5f39qNHzf6+HLe6yh6zuNj7fJZlWcuXW9Ybb3hv++EHy2rSxLLeecd7+/PP\nW1a7dhfW+fLLlrVunfe+xcXVrzOQIuFnymnHjeZjBuu4gTpmyG9RsixLJ06cqPDrRUVS5855kqS1\na/P01ltScbGZN7qmzp49q7xAHCgMjxms40bCMUeNMhOgfPBB6bauXc013bg475+9WrWkkyftq1WS\nfvIT8zh/l717zXPZ7T16SGPGSJ9+6n3M554z2zt0KN337bdNt/jevVJMmT6y3/1OSk2V0tJKtx05\nImVnn1VOTp5q1y7dfuyYVLu2lJDgz3fsLZT/ppZluv1jY0tn6pPMz0NMjNSwdNZcFRSYn4mLLjIL\nyJQcc+9e8//WpZd677tqlXmfmjYt3b51q5SdLfXv713HO+9IvXub96PkuBs3SmvWSKNHe+87Y4bU\npYt07bWl23bulP7v/6RHH/XukZk7V2rc2Pv7P3RIeu89M/YmKal030WLTN1Dh5ZuO3lSeust6Wc/\nk9q1K92+cqX0739f+G/66qvSDTeY+kps3mwu4Tz4oPf3MXOmueX0qqtKt+XmntULL+RpxAjvRXo+\n/NC8Fz/9aem2w4fN9qFDpSZNSrcvXSqdOeM942JBwVm9/HKe+vc371+J9evNuI3bb/eu7R//MO9H\n2R7YHTvM93L33eZzX35OGzRoIFfZH6xyuCwrtFeY8vLylJiYGMpTAgAQcrm5uWpY9i+5coQ8hKtq\nCUvSv/+dp+7d2+jAgQNVfgMwZs+Wnn7a/OXcuLHd1ThTYaH5y/v++6Xf/MbuasJbQYHk8Zhr3GVb\nzV99JdWt633N+vvvTUvvgQek5s1Lt8+aJeXnSw89VLrt+HHz3owebXojSnz4obk2/utfe9fx29+a\nFvrVV5duW7fODJYbP95733fekTp2lHr2LN22b59Zme3OO71bXitXmtZ8amrptrw80xK64gozYU7Z\nYxQVec9pf/as+b6bNzf/HiVOnTJfO/+/teJi739HhI5lefeCSOb9cLm8txcXm0esH/3HjmwJ+6Kk\ntezLXxEwzp41XTD16tldiXNs2ya1b+/dPbdzp9nGZCcAnIC/vSJEbCwBXNbBg+Z63Ny53ts7dCCA\nAThH2MwdDfijVSvps8+8B68AgNPQEo5AM2aY623R4uBBqU8facMG7+19+8pr5C4AOA0hHIFWrJC+\n/truKkKneXNz+0Rxsd2VAIB/COEI9OabZvakYJk8ebJiYmK8Hp07dw7eCc+zd6/3vbqxsablX3aE\nbDRavny50tLS1Lp1a8XExGjhwoUX7DNp0iS1atVK9erV080336zdu3fbUGl0q+p9GjVq1AW/X7fc\ncotN1Uav//3f/1XPnj3VsGFDJScna+jQodq5c+cF+9X0d4oQjkChGHjUtWtXeTweZWdnKzs7WytW\nrAj+SWVuE0lNlV55JSSnCyv5+flKTU3VtGnTyr0t4rnnntMrr7yiGTNmaO3atUpISNDAgQN15swZ\nG6qNXlW9T5I0ePBgr98vR0+7GKGWL1+uJ554QmvWrNFnn32mwsJCDRgwQKdPnz63T0B+pwIy71aA\n5eaaaStzc3PtLiXsFRdbVmFhYI/5zDPPWFdccUVgD+qHzz6zrJMnbTt9WHC5XNaCBQu8trVs2dJ6\n8cUXz32em5trxcfHW3Pnzg11efhRee/TyJEjraFDh9pUESpy5MgRy+VyWcuXLz+3LRC/U7SEI5hl\nmenYzl/JJxB27dql1q1bq3379rrvvvt04MCBwJ9E0saNZmq5svr1q9m0iNFo7969ys7OVr9+/c5t\na9iwoXr16qVVq1bZWBnKs2zZMiUnJ6tjx4569NFHdfToUbtLinrHjx+Xy+VSkx/nyAzU7xS3KEUw\nl8usw9yiRWCP27t3b7311lu6/PLLdejQIT3zzDO64YYbtHXrViUEMB0tS3r8cTPX6/n3+8I/2dnZ\ncrlcSk5O9tqenJys7Oxsm6pCeQYPHqw77rhDKSkp+vbbbzVx4kTdcsstWrVqVZWzLyE4LMvSL3/5\nS1133XXnxr8E6neKEI5wI0YE/pgDy8yM3rVrV/Xs2VMXX3yx3nvvPY0aNSpg53G5pPffN2v1AtFi\n2LBh5z7u0qWLunXrpvbt22vZsmXq27evjZVFr0cffVRZWVlauXJlwI/t6O7o9PR0paWlMSjB4RIT\nE9WhQ4caj7T94gtpyhTvba1amZWMUDMtWrSQZVnyeDxe2z0ej1oEuqsEAZWSkqKkpCRGstvk8ccf\n18cff6xly5apZcuW57YH6nfK0SE8Z84cLVy4UBkZGXaXEhEOHpQmTDCTzQfSyZMn9e2333r9gFbH\n6tXSggWBrw/mP/IWLVooMzPz3La8vDytWbNGffr0sbEyVOW7777Tf/7znxr/fsF/jz/+uBYsWKCl\nS5eqbdu2Xl8L1O+Uo0MYgZWVZa6t1vQS4Pjx4/Xll19q3759+uqrrzR06FDFxsbW+I+lcePMWqDM\n7Vw9+fn52rx5szZt2iRJ2rNnjzZv3nxu0Nwvf/lL/f73v9eHH36oLVu2aPjw4brooot022232Vl2\n1KnsfcrPz9eECRO0Zs0a7du3T5mZmbr99tvVoUMHr8tACL5HH31UM2fO1KxZs5SQkCCPxyOPx6OC\ngoJz+wTkdypwA7gDh1uUgic/v+bHSE9Pt1q3bm3Fx8dbbdq0sTIyMqw9e/b4dYylSy3rxhstq6Cg\n5vXAWLZsmeVyuayYmBivx6hRo87t8/TTT1stW7a06tataw0YMMDatWuXjRVHp8rep9OnT1sDBw60\nkpOTrTp16lgpKSnW6NGjrcOHD9tddtQp7z2KiYmx3n77ba/9avo7xVKGsMXWrdIzz0jTp0tJSXZX\nAwD2oDs6irnd0uDBZnH2YMvK8v68a1dp3jwCGEB0I4SjWH6+GQRVt25wz/PRR1KXLtKWLcE9DwCE\nG7qj4cWyzP25NVFQIMXHl35eWCh9+qk0aJAUw599AHAO/yXCywsvSDWZb2PWLDPDVZk5zhUXJ91y\nCwEMAOfjv0V4ad3ahGhZR46YW5vy8ry3z5olvfee97ZrrpGefpq1fQHAF3RHo0qffioNGGDW8W3X\nrnT7PfdI9etLM2bYVhoAhDVCGFWyLNMKrl+fiTQAIJBYwAFVcrmkxES7qwCAyOPoa8Is4AAAiGR0\nRwMAYBNHt4QBAIhkhDAAADYhhAEAsIkjrwlblqUTJ06oQYMGctV0DkUAABzKkSEMAEA0oDsaAACb\nEMIAANiEEAYAwCaEMAAANiGEAQCwCSEMAIBNCGEAAGzy/wHhrZK+XQjR9wAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 連立一次常微分方程式を解く関数desolve_system_rk4を使った例\n", "(z,z_dot,t) = var('z z_dot t')\n", "P=desolve_system_rk4([z_dot,-z - z_dot + 1],[z,z_dot],ics=[0,0,0],ivar=t,end_points=20)\n", "plt_f = list_plot([(i,j) for i, j, k in P], plotjoined=True, linestyle ='-')\n", "plt_g = list_plot([(i,k) for i, j, k in P], plotjoined=True, linestyle =':')\n", "(plt_f+plt_g).show(figsize=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

連立微分方程式の行列とベクトルを使った表現

\n", "\t

\n", "\t\t先の連立微分方程式を行列とベクトルで表現すると以下のようになります。\n", "$$\n", "\t\t\\left[\\begin{array}{c}\n", "\t\t\t\\dot z\\\\\n", "\t\t\t\\ddot z\n", "\t\t\t\\end{array}\\right]\n", "\t\t = \\left[\\begin{array}{cc}\n", "\t\t\t\t0 & 1\\\\\n", "\t\t\t\t-1 & -1\n", "\t\t\t\\end{array}\\right]\n", "\t\t\t \\left[\\begin{array}{c}\n", "\t\t\t\tz\\\\\n", "\t\t\t\t\\dot z\n", "\t\t\t\\end{array}\\right]\n", "\t\t + \\left[\\begin{array}{c}\n", "\t\t\t\t0\\\\\n", "\t\t\t\t1\n", "\t\t\t\t\\end{array}\\right]\n", "\t\t u(t), \\left[\\begin{array}{c}\n", "\t\t\t\tz\\\\\n", "\t\t\t\t\\dot z\n", "\t\t\t\t\\end{array}\\right]\n", "\t\t\t = \\left[\\begin{array}{c}\n", "\t\t\t\t0\\\\\n", "\t\t\t\t0\n", "\t\t\t\t\\end{array}\\right]\n", "\t\t , u(t) = 1, t \\ge 0\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tここでベクトル$x(t)$を以下のように定義します。\n", "$$\n", "\t\tx(t) = \\left[\\begin{array}{c}\n", "\t\t\t\tz(t)\\\\\n", "\t\t\t\t\\dot z (t)\n", "\t\t\t\t\\end{array}\\right]\t\n", "$$\t\t\n", "\t\t$x(t)$で連立微分方程式を表すと以下のように簡単になります。\n", "$$\n", "\t\t\\dot x(t) = \\left[\\begin{array}{cc}\n", "\t\t\t\t\t0 & 1\\\\\n", "\t\t\t\t\t-1 & -1\n", "\t\t\t\t\t\\end{array}\\right]\n", "\t\t\tx(t) + \\left[\\begin{array}{c}\n", "\t\t\t\t\t0\\\\\n", "\t\t\t\t\t1\n", "\t\t\t\t\t\\end{array}\\right]\n", "\t\t\tu(t), x(0) = \\left[\\begin{array}{c}\n", "\t\t\t\t\t0\\\\\n", "\t\t\t\t\t0\n", "\t\t\t\t\t\\end{array}\\right]\n", "\t\t\t, u(t) = 1, t \\ge 0\n", "$$\t\t\n", "\t

\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

オイラー法を使った解法

\n", "\t

\n", "\t\tオイラー法は、以下のように行列とベクトルで表現された連立一次常微分方程式を差分近似式を使って解く方法です。\n", "$$\n", "\t\t\\dot x(t) = A x(t) + B u(t), x(0) = x_0\n", "$$\n", "\t\tこれを差分近似で表現すると以下のようになります。\n", "$$\n", "\t\tx(t + \\Delta t) \\simeq x(t) + (A x(t) + B u(t)) \\Delta t\n", "$$\t\t\t\t\n", "\t

\n", "\t

\n", "\t\tそれでは、オイラー法を使って先の微分方程式を解いてみましょう。\n", "\t\tA, Bを以下のように定義し、刻み幅$\\Delta t$ = 0.05, 終了時刻$t_f$ = 20で計算します。\n", "$$\n", "\t\tA = \\left[\\begin{array}{cc}\n", "\t\t\t0 & 1\\\\\n", "\t\t\t-1 & -1\n", "\t\t\\end{array}\\right]\n", "$$\n", "$$\n", "\t\tB = \\left[\\begin{array}{c}\n", "\t\t\t0\\\\\n", "\t\t\t1\n", "\t\t\\end{array}\\right]\n", "$$\t\t\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFnCAYAAACRo/HLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl8VPW9//HXZIEAQjAsAQJCBJVNCbKK2jYFQbTEa1VM\nSmuwWEBQoeBSFBdcqlYr117h+pOKOwGlVqFXREVZRCooYQkhYS2bJCUiSQhgQnJ+f3wbkoEsM8nM\nnDOT9/PxmMfIyZkznxCHd77f811clmVZiIiISMCF2V2AiIhIQ6UQFhERsYlCWERExCYKYREREZso\nhEVERGyiEBYREbGJQlhERMQmCmERERGbKIRFRERsohAWERGxidchvGbNGpKSkoiLiyMsLIwlS5bU\neP7f//53hg8fTtu2bYmOjmbIkCF88skndS5YREQkVHgdwkVFRSQkJDB37lxcLlet569evZrhw4ez\nbNkyNm7cSGJiIqNGjWLz5s11KlhERCRUuOqzgUNYWBgffPABSUlJXr2ud+/eJCcnM3PmzLq+tYiI\nSNAL+D1hy7IoLCwkJiamxnMKCgrQBk8iIhLKAh7Czz33HEVFRYwePbracwoLC4mOjqawsDCAlYmI\niARWRCDfbMGCBTzxxBMsWbKE1q1b13p+cnIyERHuJaakpJCSkuKvEkVERAImYCG8cOFCxo8fz+LF\ni0lMTPT4NS1atPBzZSIiIvYISHd0Wloa48aNY+HChVx77bWBeEsRERHH87olXFRUxK5du84Mmtqz\nZw+bN28mJiaGTp06MWPGDL777jveeOMNwHRBjx07lr/85S8MGDCA3NxcAJo0aaJWroiINGheT1Fa\ntWoViYmJ58wRTk1NZf78+dx+++3s27ePzz//HIDExERWr159znXKz69KQUEB0dHR5OfnK6hFRCRk\n1WuesL8ohGu2cSMsWACrVsHOnXDqFLRqBR07ws9/DiNHwtVXgwdrqYiIiI0UwkFk61aYNg0++wxi\nY2HECOjVC6Ki4PvvTSB/9hkcOWKOT58Ov/kNRAR0DLyIiHhK/zwHAcuC2bPhD3+ACy+E996DG2+E\n8PBzzy0rg9Wr4YUX4Le/hf/+b5g7F668MvB1i4hIzbSLksOVlcGUKaZVO2UKbNkCN99cdQADhIXB\nz34GS5bAhg3QuDFcdRX8/vdQXBzQ0kVEpBYKYQcrK4OxY+Gll+Dll+G556BRI89f378/rFtnWtFz\n5pjW8N69fitXRES8pBB2sIcfhrffhrQ0mDChbtcID4epU+Grr+DoURg0CP75T9/WKSIidaMQdqiF\nC+GPf4Rnn4Vbb63/9fr3h6+/hksugcREeP/9+l9TRETqRyHsQNu3m0FVv/413Huv767burUZPX3D\nDTB6tJnmJCIi9nH06OjyDRwa0qYNpaUwbhx06gSvvOL7ub6NG8M770CTJibkS0ogNdW37yEiIp5x\ndAg3xA0c5swxg6nWrDFB6Q/h4fDqqxAZCbffboL4jjv8814iIlI9R4dwQ7NvH8yYAZMnm2lF/hQW\nZkZcR0bC+PFw3nmQnOzf9xQREXcKYQeZORNatICnnw7M+4WFwf/8DxQWmpW1WrYEbXIlIhI4CmGH\n2LLF3KudOxeaNw/c+4aFma7pY8fgl7+ETz/V6loiIoGitaMd4vrrzdrP27aZLuJAO3nStIK3bDFz\ninv0CHwNIiINjaYoOcCaNfDRR/Dkk/YEMJhBYEuWQFwcXHcd/GfbZxER8SO1hB3guuvg0CFITzfd\nw3batw8GD4YLLoAvvoCmTe2tR0QklKklbLNt22DZMrMoh90BDNC5MyxdChkZcNttZv1qERHxDwf8\ns9+wvfACdOjgm6UpfaV/f7Oa1vvvwwMP2F2NiEjoUgjbKCfHbNAwZYp3uyMFwg03mN2Xnn/ezCcW\nERHf0xQlG82da8J3/Hi7K6nalCmwe7dZPKRzZxg50u6KRERCi6MHZo0cOTJk144uLTXBNmoU/O//\n2l1N9UpL4b/+C1atgrVr4dJL7a5IRCR0ODqEQ3l09EcfmbnBGzaYe7BOdvw4XH01fP+92Q6xfXu7\nKxIRCQ26J2yTv/4V+vSBfv3srqR2551nRkyXlkJSEhQV2V2RiEhoUAjbIDfXhNq4cb7fqtBfOnaE\nf/wDMjPNOtOauiQiUn8KYRu89ZbZTnDMGLsr8U7fvpCWBh98AH/4g93ViIgEP4WwDd56ywx2iomx\nuxLvJSWZuc3PPQfz5tldjYhIcFMIB9j27WaThGDeu3fKFJg0Ce68Ez77zO5qRESCl0I4wBYtMnsG\nB/O+vS4XvPgiXHMN3HyzuU8sIiLeUwgHkGWZEL7hBoiKsrua+omIMN/LBReYqVbadUlExHsK4QDa\nuhWysoK7K7qyFi3MiOlTp8y94uPH7a5IRCS4KIQDaNEiOP98GDbM7kp854ILKqYu3XQTFBfbXZGI\nSPBQCAfQ4sVw443O26yhvvr1gw8/hJUrITVVc4hFRDylEA6Q7GzYscNMTQpFP/+5mUP87rtwzz3m\n/reIiNTM0SGcnJxMUlISaWlpdpdSb0uXmsFYQ4faXYn//PKXZtvDOXNg1iy7qxERcT5Hb2W4cOHC\nkNnAYckSM6WnaVO7K/Gv3/0O8vLgwQfN93r//XZXJCLiXI4O4VDx/fdmG8CXX7a7ksD4wx/g5El4\n4AEzp/i+++yuSETEmRTCAfDRR2aw0i9+YXclgeFyme5oyzIt4bAwmD7d7qpERJxHIRwAS5fCgAEN\nax9elwsef9z88nHvvWYu8YMPBs+uUSIigaAQ9rOSEli+vGG2BF0uePJJMyBt5kzTLf/886ZlLCIi\nCmG/W78eCgqCe63o+nC54OGHoVUruOsuOHoU/vpXs+yliEhDF5A2yZo1a0hKSiIuLo6wsDCWLFkS\niLd1hE8/Natk9etndyX2mjQJ3nnHPK6/Ho4ds7siERH7BSSEi4qKSEhIYO7cubga2E3BTz4xc4PD\nw+2uxH4pKaZrfsMGGDwYdu60uyIREXsFJISvvfZaHn/8cW644QasBrSU0rFjpjv6mmvsrsQ5fv5z\n+Ppr89+DBpmeAhGRhkp35vzoiy+gtFQhfLaLLoJ//tO0jEeMgBkz4LHHIDLS7spCh2XBiRPmF8Fj\nx+CHH8xzfr4ZqX7qFPz447n/7cm6340aVTwaN3b/89lfi4oyj8r/ffafGzfWYD1puBTCfvTpp9Ct\nG8TH212J87RsCf/3f/CnP5mR0198Ydae7tzZ7sqcr7AQ9u41j0OH4PBh90dOjlm1rKSk+mtERlYd\nhrXdNrEsc90ffzQ7Zp39KC2t2/dUXo8noe3pOZ7+ufz7DgureNh918yyPHuUldXta/X9uievrU5V\nf7fBdqxVK+jQ4dzz6kIh7EeffgrDh9tdhXOFhZnVtX76U9MqvvRSePppuPNOtYxOnoTt22HbNrMH\n9Z49FY+8vIrzIiKgXTszB719exg40Dy3aWMGBLZs6f5o0QKaNPHf329paUVIlz/KW9tnt77r8ufj\nx833X9trfHHXy+VyD+bycK7qmKeh6clDnO/uu+Evf/HNtRwdwsnJyUScNZclJSWFlJQUmyry3L59\nsGuXaelJza64AjZtMoF8111mBPUrr0Dv3nZX5n+lpSZst241gZuRYZ737KloTcTFmR6VXr1g1Ci4\n8ELziI83YeukX1jCw80jKsq+Gspb694EfVnZuY/yFp0nx10u3zzKQ90X59b09UC9tqqfTSgca9v2\n3GN15egQDuYNHFatMs8/+Ym9dQSLli3N2tpjxsD48ZCQAOPGwaOP+q7bx25lZWY7y2++qXikp5t7\nt2DCtjxoe/Uyj549TetVPOdyVdyX1t+dOF1AQrioqIhdu3adGRm9Z88eNm/eTExMDJ06dQpECQG3\napXpXm3Vyu5KgsvVV5tW8Zw58NRT8NZbZn/i3/8eYmPtrs5zlgW7d7sH7saN5n4umJZt//5w441m\nDnlCgvlFREQaFpcVgDlDq1atIjEx8Zw5wqmpqcyfP/+c8wsKCoiOjiY/Pz9oW8LdusHIkfA//2N3\nJcErP98sczl7tuleTEmBKVOgb1+7K3NnWbB/vwnaDRvM87ffVixI0qWLCdzyx+WXm/u1IiIBCWFv\nBXsIHzoEHTvCe+/BzTfbXU3w++EHePVV8wvN/v0mhJOT4dZbAz+a2rIgN9c9cDdsgCNHzNfj4irC\ndsAA08pt3TqwNYpI8FAI+8GCBebeZm6ub2/gN3SnT5sdqRYsgH/8wwyq6dMHEhPN48orfdf9b1km\n/HfsMIOmMjLM89atFaOTW7c2QVseuP37N6ydskSk/hTCfjBhAqxZA5mZdlcSugoLTRAvX27mGO/f\nb463b29GVffoYQZ0tW9v7iU3bVoxJ/T0aTOn9ccfTdAeOVLxOHjQjEzevdt0h4MZ7dmtm7nHX/7o\n1w8uuMD++aQiEtwUwn7Qvbtpmf3v/9pdScNgWWbhig0bKlqsO3aYhSs82SjC5TIt6DZtTGh37Wqm\nAHXtah49epi5tSIivqYQ9rHcXLN4QlqauW8p9jp50rRwT56smBcaEWFaxI0bQ3S0CWBtsCEidnD0\nPOFgtHateb76anvrEKNJE9NtLCLiRA5aayc0fPWVGbEbF2d3JSIi4nQKYR/76iuzDKOIiEhtHB3C\nycnJJCUlkZaWZncpHjl1yizSMGSI3ZWIiEgwcPQ94WBbO3rjRjP1RSEsIiKecHRLONh89ZWZj3rZ\nZXZXIiIiwUAh7EPr1pn9XCMj7a5ERESCgULYRyzLtITVFS0iIp5SCPvIv/4FOTkKYRER8ZxC2EfW\nrTPPgwfbW4eIiAQPhbCPbNhg1hn21S4+IiIS+hTCPrJ+vRmUJSIi4imFsA+UlJg5wgMG2F2JiIgE\nE4WwD2zbZlbLUktYRES8oRD2gfXrzVZ4ffvaXYmIiAQThbAPbNgAvXub1bJEREQ85egQDpYNHDQo\nS0RE6kIbONRTUZG5J3zXXXZXIiIiwcbRLeFgkJ4OpaVqCYuIiPcUwvX0zTcQFQW9etldiYiIBBuF\ncD1t3Ah9+kCEozv2RUTEiRTC9ZSerqlJIiJSNwrhejh5ErZvVwiLiEjdKITrYetWMyjr8svtrkRE\nRIKRQrge0tPNSlm9e9tdiYiIBCOFcD2kp0PPnmZ0tIiIiLcUwvWgQVkiIlIfCuE6On0atmxRCIuI\nSN05OoSdvHZ0VpbZvlCDskREpK4cvcSEk9eOTk83zwkJ9tYhIiLBy9EtYSdLT4euXcGhvyOIiEgQ\nUAjXkQZliYhIfSmE68CyFMIiIlJ/CuE62LsX8vM1KEtEROpHIVwH5YOy1BIWEZH6qFMIz5kzh/j4\neJo0acLgwYPZsGFDjee/8847JCQk0KxZMzp06MC4ceM4evRonQp2gvR0aN8eYmPtrkRERIKZ1yG8\naNEipk+fzqxZs0hPT6dPnz6MGDGCvLy8Ks9fu3Ytqamp/O53vyMzM5PFixezfv16xo8fX+/i7bJ5\ns6YmiYhI/XkdwrNnz2bChAncdtttdO/enZdffpmmTZsyf/78Ks//5z//SXx8PJMnT6Zz584MGTKE\nCRMmsH79+noXb5etW+HSS+2uQkREgp1XIVxSUsK3337L0KFDzxxzuVwMGzaMdevWVfmaK664ggMH\nDrBs2TIAcnNzWbx4Mddff309yrZPfj7s26cQFhGR+vMqhPPy8igtLSX2rJuhsbGx5OTkVPmaIUOG\n8Pbbb3PrrbfSqFEj2rdvT8uWLXnppZfqXrWNMjLMs0JYRETqy++jozMzM5kyZQqPPfYYGzduZPny\n5ezdu5cJEyb4+639YutWiIiA7t3trkRERIKdV2tHt27dmvDwcHJzc92O5+bm0q5duypf88wzz3DV\nVVcxbdo0AHr37s3cuXO5+uqreeqpp85pVVeWnJxMRIR7iSkpKaSkpHhTtk9t3QqXXAKNG9tWgoiI\nhAivQjgyMpJ+/fqxYsUKkpKSALAsixUrVnDPPfdU+ZoTJ04QGRnpdiwsLAyXy4VlWTW+nxM3cNCg\nLBER8RWvu6OnTZvGvHnzePPNN8nKymLixImcOHGCsWPHAjBjxgxSU1PPnD9q1Cjef/99Xn75Zfbu\n3cvatWuZMmUKgwYNqrb17FSWpRAWERHf8Xorw9GjR5OXl8cjjzxCbm4uCQkJLF++nDZt2gCQk5PD\ngQMHzpyfmprK8ePHmTNnDvfeey8tW7Zk6NChPPPMM777LgLk0CE4dkwhLCIivuGyausTtkFBQQHR\n0dHk5+c7qjv6o4/g+uvN2tFduthdjYiIBDutHe2FrVuheXPo3NnuSkREJBQohL1Qfj/Y5bK7EhER\nCQUKYS9oUJaIiPiSQthDJSWwfbtCWEREfEch7KHsbBPECmEREfEVhbCHtm41zwphERHxFYWwh7Zu\nhbg4OP98uysREZFQ4egQTk5OJikpibS0NLtLYetWuOwyu6sQEZFQ4vWKWYHkpLWjt22Dm2+2uwoR\nEQkljm4JO0VREfzrX9Czp92ViIhIKFEIeyA722zeoBAWERFfUgh7IDPTPPfoYW8dIiISWhTCHsjM\nhE6dzLrRIiIivqIQ9kBmprqiRUTE9xTCHti+XV3RIiLiewrhWvz4I+zapZawiIj4nkK4Fjt2QFmZ\nQlhERHxPIVwLjYwWERF/UQjXIjMT2rWDmBi7KxERkVCjEK5FZqZawSIi4h+ODmEnbOCwfbvuB4uI\niH9oA4calJSYgVmTJtlWgoiIhDBHt4Tttnu3CWK1hEVExB8UwjUoHxmtEBYREX9QCNcgMxNatYI2\nbeyuREREQpFCuAblI6NdLrsrERGRUKQQroFGRouIiD8phKtRWgpZWQphERHxH4VwNf71Lzh1SiEs\nIiL+oxCuhkZGi4iIvymEq5GZCc2bQ4cOdlciIiKhSiFcjcxM0wrWyGgREfEXR4ewnWtHa2S0iIj4\nm9aOroJlmZbw6NEBf2sREWlAHN0StsvBg1BUBN27212JiIiEMoVwFbKyzLNCWERE/EkhXIXsbGjU\nCLp0sbsSEREJZQrhKmRnQ7duEOHoO+YiIhLsFMJVyMqCSy6xuwoREQl1dQrhOXPmEB8fT5MmTRg8\neDAbNmyo8fzi4mIeeughunTpQlRUFBdeeCGvv/56Xd46ILKzdT9YRET8z+sO10WLFjF9+nReeeUV\nBg4cyOzZsxkxYgQ7duygdevWVb7mlltu4ciRI7z22mt07dqVw4cPU1ZWVu/i/aGoCA4cUEtYRET8\nz2VZluXNCwYPHsygQYN48cUXAbAsi06dOnHPPfdw//33n3P+xx9/zK9+9Sv27NlDy5YtPXqPgoIC\noqOjyc/PD/g84fR0uPxyWLcOBg8O6FuLiEgD41V3dElJCd9++y1Dhw49c8zlcjFs2DDWrVtX5WuW\nLl1K//79efbZZ+nYsSOXXHIJ9913H6dOnapf5X6SnW2e1RIWERF/86o7Oi8vj9LSUmJjY92Ox8bG\nkl2eXmfZs2cPa9asISoqig8++IC8vDzuvPNOjh49yquvvlr3yv0kOxvatoXzz7e7EhERCXV+n4RT\nVlZGWFgYCxYs4LzzzgPghRde4JZbbmHu3Lk0btzY3yV4JTtbrWAREQkMr0K4devWhIeHk5ub63Y8\nNzeXdu3aVfma9u3bExcXdyaAAXr06IFlWRw8eJCuXbtW+37JyclEnDVZNyUlhZSUFG/K9kpWFvTv\n77fLi4iInOFVCEdGRtKvXz9WrFhBUlISYAZmrVixgnvuuafK11x55ZUsXryYEydO0LRpUwCys7MJ\nCwujY8eONb5foDdwsCzYsQPGjAnYW4qISAPm9TzhadOmMW/ePN58802ysrKYOHEiJ06cYOzYsQDM\nmDGD1NTUM+f/6le/olWrVtx+++1s376d1atXc//99zNu3DjHdUUfOmSmKKk7WkREAsHre8KjR48m\nLy+PRx55hNzcXBISEli+fDlt2rQBICcnhwMHDpw5v1mzZnz66afcfffdDBgwgFatWnHrrbfyxBNP\n+O678BFt3CAiIoHk9TzhQLBrnvCcOfD738OJE1o3WkRE/E9rR1eijRtERCSQFMKVZGWpK1pERAJH\nIVyJ5giLiEggNcgQ3rsXHngAvv224tiJE7B/v0JYREQCp0GGcIsWsHQpNGlScWzHDvOs7mgREQmU\nBjkEqVUryMx0P7Ztm3lWS1hERAKlwbSET56s+etvv22ei4v9X4uIiAg4PISTk5NJSkoiLS2tXtc5\netS0cFeurP6c6Gizj/BZG0SJiIj4jaO7o321dnRYGNxyC/ToUf05O3eaEBYREQkUR7eEfaVlS/jz\nn6tv5VqWpieJiEjgNYgQrk35xg3lI6O/+QYmTzbhLCIi4i8hHcKehmh2tnkubwkfOwbr18O//+2f\nukRERCDEQ3jhQujf3yzEUZPsbIiMhPh48+dhw0wIa5CWiIj4U0iHcOfOMHw4NG1a83lZWedu3OBy\n+bc2ERERR4+Orq8hQ8yjNhqUJSIidgjplrCnqgvhggK47z7Yvj3wNYmISOhr8CF84gTs21f1mtGN\nG8OyZRVLWoqIiPhSSIbw0aMwfjz861+1n7tzp3muqiXcuDFs3Qo33+zT8kRERIAQDeF9+2D1ajPi\nuTZnT086mwZoiYiIv4TkwKy+fc2IZ09kZ0Pr1hAT49+aREREzubolrCvNnCoiScjo4uLYcEC2LvX\nb2WIiEgD5LIs5y3OWFBQQHR0NPn5+T7ZwKEmAwbAZZfBq69Wf87JkxAXB08+CZMm+bUcERFpQEKu\nOzorC5o0MQt11KZ844Zbbqn5vCZNTCs4Oto3NYqIiEAIhvDMmWbN59Wraz83NxcKCz1bqEMBLCIi\nvhZyIfzKK55vvFA+Mvrii/1Xj4iISHUcPTCrLmJiql54oyrZ2RAeDl27en79zEwoLa1bbSIiIpWF\nXAh7Izvb7JzUqJFn52/aBL16wcqVfi1LREQaiAYfwt5s3NCnD3z0EVx1lf9qEhGRhiOkQviGG2DO\nHM/P37HDuxB2uWDkSLOcpYiISH2FVAj36gUdO3p2bnEx7NmjQVkiImKfkBod/cc/en7unj1mgFVd\n9xG2LK0rLSIi9RNSLWFv1LZxQ03mz4eEBBPEIiIideXoEPbn2tHZ2dC8ObRr5/1re/aEm24yXdoi\nIiJ1FTJrR7/1Flx9NXTp4tl73HEHbN4MGzbUvU4REZH6cHRL2FOFhTBuHHz5peevyc7WoCwREbFX\nSAzMat4cCgq8u0ebnQ3XXOO/mkRERGoTEi1hgKgos9uRJ374AY4cqfvI6HKzZ8PSpfW7hoiINFwh\nE8LeqM/I6Mo++8zcVxYREamLkOiO9taOHeb5oovqd51//ENzhUVEpO6CviVcVgYXXADvvOP5a7Kz\nzcpazZrV770VwCIiUh91CuE5c+YQHx9PkyZNGDx4MBs8nOezdu1aIiMjufzyy+vytlUqLoaJE82S\nlZ7yduMGERERf/A6hBctWsT06dOZNWsW6enp9OnThxEjRpCXl1fj6/Lz80lNTWXYsGF1LrYqUVHw\n4INmBStP+TKEi4th9WqtniUiIt7zOoRnz57NhAkTuO222+jevTsvv/wyTZs2Zf78+TW+buLEiYwZ\nM4bBgwfXuVhfKC2FnTt9F8IrVsBPfwpZWb65noiINBxehXBJSQnffvstQ4cOPXPM5XIxbNgw1q1b\nV+3rXnvtNfbu3cujjz5a90p9ZP9++PFH34VwYiJs3Ajdu/vmeiIi0nB4FcJ5eXmUlpYSGxvrdjw2\nNpacnJwqX7Nz504efPBB3nnnHcLCfD8ObN482LLF8/PLR0b7arWsqCjo21eDtERExHt+naJUVlbG\nmDFjmDVrFl27dgXAm6Wqk5OTiYhwLzElJYWUlBTAdC0/9BA8+ihcdpln18zOhsaNzYhqERERO3kV\nwq1btyY8PJzc3Fy347m5ubSrYjuiwsJCvvnmGzZt2sTkyZMBE8yWZdGoUSM++eQTfvazn1X7fgsX\nLqxxA4fwcMjNhdOnPf8esrPN/ODwcM9f4yntMSwiIt7wqn84MjKSfv36sWLFijPHLMtixYoVDBky\n5JzzW7RoQUZGBps2bWLz5s1s3ryZiRMn0r17dzZv3sygQYPq/Q24XBAZ6fn5/pqeNGYMTJ3q++uK\niEjo8ro7etq0aYwdO5Z+/foxcOBAZs+ezYkTJxg7diwAM2bM4LvvvuONN97A5XLRs2dPt9e3bduW\nqKgoevTo4ZNvwFvZ2XDbbb6/bmIiREf7/roiIhK6vA7h0aNHk5eXxyOPPEJubi4JCQksX76cNm3a\nAJCTk8OBAwd8XmhVvO3+LSqCgwf9s4XhHXf4/poiIhLaXJY3I6UCpKCggOjoaPLz82u8J/yTn8AV\nV8Czz3p23U2bzEjmdevA5unKIiIiwb2Bw7hx0KGD5+f7avckERERXwjqDRxSU+Gaazw/Pzsb2rSB\n88/3Tz07d8L993s3WltERBquoA5hb/l744ajRyEtDfbt8997iIhI6Ajq7mhv7djh+aIedTFwoFkW\nU3OFRUTEE0HbEv7sM3jvPc/Ptyz/t4RdLgWwiIh4LmhD+IMPYO5cz8/PyYHCQg3KEhER5wja7uiX\nXjJrR3sqkCOj8/Lg5Eno1Mn/7yUiIsHL0S3h5ORkkpKSSEtLq/Lr3qz/nJ1tzr/wQh8VV4Of/hSe\neML/7yMiIsHN0S3h2jZw8MaOHRAfD40a+eRyNXrtNejY0f/vIyIiwc3RIVyduuxW5O9BWZUNHBiY\n9xERkeDm6O7o6rz+OrRvDyUlnr8mKwu6d/dbSSIiIl4LyhDu2xfuvdfzLQxPnYK9e8GmjZtERESq\nFJQhnJAA06d7fv6OHVBWFtgQXr4cLr/cu9a6iIg0LEEZwt7KyjLPgeyObt/e7NR0/Hjg3lNERIJL\nUA7M8tb27dC2LcTEBO49L7vMu8VERESk4Qm6lvCPP8J//7d3myRs3677wSIi4jxBF8KHD8ODD8Ke\nPZ6/RiOjRUTEiYIuhLt0gaIisyqVJ0pLzRxhO1rClgVvvgkbNgT+vUVExPmCLoTBLNQR5mHl+/aZ\nKUp2hLDsbwEDAAAUUElEQVTLBX/8I3zySeDfW0REnM/RA7OSk5OJiIggJSWFlJSUOl1j+3bzbFd3\n9KZNEBVlz3uLiIizOTqEfbF2dFYWNGtm345GCmAREalO0HVHX3MNPPCA5+dv325awd6uNS0iIuJv\njm4JVyUlxSyE4SknTE+yLDOqu0MHe+sQERFnCbqW8G9/CyNHenauZVW0hO30/PPQq5cZqS0iIlIu\n6FrC3jhyBH74wf6W8C9/Cb17m18KREREyoV0CJePjLY7hLt2NQ8REZHKgqo7euNGeOcdz8/PyoLw\ncAWgiIg4U1CF8EcfmSUrPbV9O3TrBo0a+a8mERGRugqqEJ45E3bu9Px8J4yMLnfiBNx8M3z6qd2V\niIiIUwRVCIN3rVonbdzQpImZq3z6tN2ViIiIU4TswKyCAti/H3r2tLsSw+WC996zuwoREXGSoGsJ\neyoz0zxfeqm9dYiIiFTH0SGcnJxMUlISaWlpHDkCHTvCF1949tqMDLPTklO6o0VERM7m6O7oyhs4\nfP+9WS2rc2fPXpuRYUZGO20Dhaws2LEDkpLsrkREROzm6JZwZa1aweOPw4UXenb+tm1mlSqneftt\nuPturZ4lIiJBFMLeyshwZgjfd59pCWtXJxERcXR3dF3l5UFOjjNDODra7gpERMQpgqYl/PnnsHat\nZ+du22aenRjCIiIi5YKmJfzMM9CyJVx5Ze3nZmSYRT26dfN/XXVVXGy6pCMj7a5ERETsUqeW8Jw5\nc4iPj6dJkyYMHjyYDRs2VHvu3//+d4YPH07btm2Jjo5myJAhfPLJJ16/57JlMH++Z+dmZMAllzg3\n4A4ehJgYWLHC7kpERMROXofwokWLmD59OrNmzSI9PZ0+ffowYsQI8vLyqjx/9erVDB8+nGXLlrFx\n40YSExMZNWoUmzdv9up9w8PhvPM8O9epg7LKxcXBk09qDrOISEPnsizvJssMHjyYQYMG8eKLLwJg\nWRadOnXinnvu4f777/foGr179yY5OZmZM2dW+fWCggKio6PJz88/M0/YU5ZlWpn33efdjksiIiKB\n5lVLuKSkhG+//ZahQ4eeOeZyuRg2bBjr1q3z6BqWZVFYWEhMTIx3lXrou+/g2DFnt4RFRETAyxDO\ny8ujtLSU2NhYt+OxsbHk5OR4dI3nnnuOoqIiRo8e7fH7fv65WQO6mh5vNxkZ5lkhLCIiThfQ0dEL\nFizgiSeeYMmSJbRu3brW85OTk4mIiKCgAI4fh7FjYcyYFFJSUqp9zbZt0LQpdOniu7r95bXX4PBh\ndZuLiDRUXoVw69atCQ8PJzc31+14bm4u7dq1q/G1CxcuZPz48SxevJjExESP3q/y2tGeysiAXr3M\n5g1Od+gQHDhgdxUiImIXr6IqMjKSfv36saLS3BrLslixYgVDhgyp9nVpaWmMGzeOhQsXcu2119a9\nWg+Uh3AwmDkT/t//s7sKERGxi9ftxWnTpjFv3jzefPNNsrKymDhxIidOnGDs2LEAzJgxg9TU1DPn\nL1iwgNTUVP785z8zYMAAcnNzyc3NpaCgwGffRLmyMudu3CAiInI2r0N49OjRPP/88zzyyCP07duX\nLVu2sHz5ctq0aQNATk4OByr1sc6bN4/S0lImT55Mhw4dzjymTp3q8Xt+/DHs3Fn7eXv3wokTCmER\nEQkOXs8TDoSz5wm3aQNTp8JDD9X8ur/9DW6+2WzecNYAbscqLoYvvoCBA+H88+2uRkREAiko1o7e\ntcuzrf82bYJ27YIngAGOHoVrr4WFC+HWW+2uRkREAikoQtjT7f/S0yEhwb+1+Fq7dmZ/YSdvNiEi\nIv4RBBN5PLdpE/Tta3cV3rvoIs9a+iIiElpCJoSPHDHzboOtJSwiIg2X40N40SJITDQbM9Rk0ybz\nHIwt4XLFxXZXICIigeT4EI6JgZ49a++u3bQJmjWDrl0DU5evPfYY1LDeiYiIhCDHD8y65hrzqE16\nOvTpExzLVVZl6FCIjzctft0fFhFpGBwdWcnJySQlJZGWllbrucE6KKvc1VdDaqoCWESkIQmKxTpq\nc+IENG9u1mG+444AFCgiIuIDjm4JA6xaZbb7q8mWLWbdaI2MFhGRYOLoEC4rg2HD4O9/r/m8b76B\nyEi49NLA1OUvp0/D3XfDRx/ZXYmIiASCowdmuVxm44baeqTXrzet4MaNA1OXv0REwP79kJdndyUi\nIhIIjg/hLl1qP2/9es9GUAeDDz+0uwIREQkUR3dHe+LYMcjONrsQiYiIBJOgD+FvvjHPCmEREQk2\njg7hhQvhtttqPmf9erPL0kUXBaamQDh6FN55p/alOkVEJLg5OoQjIsxSlDVZvx4GDAjelbKqsnGj\n+eVj5067KxEREX8K6sU6LAs6dIDf/haeeiqABfpZSQl8/73Za1hEREJXULcfDx2CnJzQux8cGakA\nFhFpCBwdwrWtHf311+Y51EJYREQaBkfPE77//oX07duC6Oiqv/7ll2Yecfv2AS0rYCwLcnPVKhYR\nCVWObgknJsLSpdV/fdUq+OlPA1dPoD34oNlj2Hl37UVExBcc3RJevRp69Kj6a8eOme0L7747sDUF\n0pgx8POfa49hEZFQ5egQ7tOn+nWjv/zShFMot4R79zYPEREJTY7ujq7JqlXQsSPEx9tdiYiISN0E\ndQj/5CfqphURkeDl6BC+//6qjxcWmlWlQrkrurLJk2HGDLurEBERX3P0PeEffqj6+FdfQWlpwwnh\niy6C886zuwoREfG1oFy28g9/gNdfh8OH1R0tIiLBy9Hd0dVZtgyGD1cAi4hIcAu6ED50CLZsgZEj\n7a5ERESkfhwdwidPnnvs449NC3j48MDXY7dnn4V337W7ChER8RVHhnBxsXkeNuzcDRyWLoVBg6BV\nK5uKs1F6uvYYFhEJJY4cmHX0aAGtWkWzfXs+3btXDMw6fhzatIHHH4f77rOxQJto+UoRkdDiyJZw\nxH8mTnXo4H582TI4dQpuuinwNTmBAlhEJLQ4MoSr8+67kJAAF15odyUiIiL1FzQhfPQoLFkCv/mN\n3ZXY7//+D9assbsKERGpL0eG8Kefnnts4UKzStaYMYGvx2meeQYWLLC7ChERqS9HhvChQ+5/tix4\n+WX4xS8gNtaempzkH/+AuXPtrsKoPHJdnEs/p+Chn1Vw8NXPyZEhPHas+5+XL4etW+H3v7elHMeJ\njnbOIK2G8A+GZUFenhkUWNkrr5glVCs7fhy6dTO/KFW2YAH86lfnXnvxYtixw7f1VqUh/JxChX5W\nwSGkQ7gyy4Inn4QBA8zWhSL+VFoKubnux/bsMVPjvvzS/fipU1BU5H4sKsqM3u/Y0f14ZCQ0a+Z+\nzLLMGIePP3Y//u67MHAglJW5H//663N7iUQkuDk+hNPSYO1aeOqp+rf+/PEbpp3X3L8fVq/2/XXt\n5q86PbnunXfCdde5H+vcGd57D/r0cT9+zz0wZIj7NSMizMpmCQnu595yC8yb537M5TI7hd1xh/vx\njIw0rrwSws76dP7Xf8Ff/+p+bO1a03N09upyx46dG+L+EGqfKbuv6S/B8v3b+dm3iyNDOD/fPO/e\nbf6hu+kmuOaa+l83WP6n8fSaDz0EU6f6/rp2C8QHsbgYpk+Hzz93P+fOO2H2bPdjERFw882mNVzT\nNesiKgqaNnU/tmlT2jk1gNnCc+JE92MFBaalHhXlfnzIEJg2zf3Yrl3wwgumy9xXQu0zZfc1/SVY\nvv+GGMIB30/YsiwKCwur/XpZGXTrVgDAyJEFnH8+PP+8+cemvk6fPk2BLy7kkGs+8ojp4vT07UPt\n+6/vdTdsgEsugf79K87p2tU8O/HvtHyp1spfuvJKc//57I/UE09A27YV554+fZpvving4YchOdm9\nlXzrrdC7Nzz8cMWxnBz47DNISoLKu4mePAmNG1e00oPl53/2NUtLzd9BZGTFOSUl5vtu29Z8j+X2\n7YMff4SLL644VlwMhw+fZufOArfBops3w969pteisldfNbcYLr204lhGBqxaBZMnu5+7Z89pli4t\ncNsvffdu0xszeTI0b15xfNEiaNLE/JzKHTkCb75pxiC0b1/x/f/tbwX88AOMHl1x7qlTMH8+jBhR\n8f8+mM/Gzp3njmN44w3zeenVq+LvdMcOc6vk7OmjS5dCXBxcfnnFsUOHTO/dDTe4//L5+efm77zy\nz+mHH+CTT2DYMPdlitevN/+/Dx1acay4GD74AK64Ajp1cv873rfv3J//+++bn8VFF7n/HW/adO5i\nUJ98Au3awWWXuR+v7f/T5s2b4/Kg+zbgy1aW7xUsIiISqvLz82lR+TfYagQ8hGtrCVsWHDlSwEUX\ndeLAgQMefRMiNSksNL/9//GP0Lev3dUEj5IS0xVf+Zf5VatMK/jqqyuOHT4Mf/4z3H23uX9e7sUX\nTcvnT3+qOFZQYLYhffJJSEysOD5vnlmW9v333WsYNQpSU83tgHKffWZGpi9a5F7b009Dz56mlVVu\nzx74299gwgT3Fv2KFablddVVFccKC00rq18/aNmy4vj+/ebvonJLsazMtJrPP9+0RMudPm2eIwLe\nx2ivqta1P33aHAsPrzhWVmaOR0a6n19SYp4r90yUlZnjkZHu4yNKSszXKvdWWJZp1Tdq5P5+xcWm\n16PyzwjMgMpGjc7tCfnxRzjvvHPPDQ8/95ZPbRzbEvZEeWvZ098kGrqSEjPC9he/cM7UJbscP266\n4saPb3j/EIpI8HHkwCzxzooV5p5QRobdldhvxw4zn/ybb+yuRESkdmoJhwDLgu3bTVdcQ1NVN9j3\n3zfM/aZFJPioJRwCXK6GGcAHDpipONu2uR9XAItIsFAIS9Bq29ZMgSgf1CEiEmzUHR1itm0zo1WH\nDbO7Et8rKzMPDbgSkVChlnCImTXLLPHp3/eYRVhYmNujp5/7w0+fNmuHP/+8X98mqK1Zs4akpCTi\n4uIICwtjyZIl55zzyCOP0KFDB5o2bco111zDrl27bKi0Yavt53T77bef8/m67uy1VMXvnn76aQYO\nHEiLFi2IjY3lxhtvZEcVu63U9zOlEA4xc+eaXaf8rXfv3uTm5pKTk0NOTg5fnr27gY9FRJi5vlde\n6de3CWpFRUUkJCQwd+7cKucnPvvss7z00ku88sorrF+/nmbNmjFixAiKi4ttqLbhqu3nBDBy5Ei3\nz5eTl10MVWvWrOHuu+/m66+/5rPPPqOkpIThw4dzstJC7T75TFkOlJ+fbwFWfn6+3aVIFR577DGr\nb9++fn+fkhK/v0XIcrlc1ocffuh2rH379tYLL7xw5s/5+flWVFSUtWjRokCXJ/9R1c9p7Nix1o03\n3mhTRVKdI0eOWC6Xy1qzZs2ZY774TKklLHWyc+dO4uLi6Nq1K7/+9a85cOCAT68/Zw787GdmtRup\nv71795KTk8PQSgvutmjRgkGDBrFu3TobK5OqrFy5ktjYWLp3786kSZM4evSo3SU1eMeOHcPlchET\nEwP47jOlIS4h7LnnzAYPkyb59rqDBw/m9ddf55JLLuHw4cM89thj/OQnPyEjI4NmZ2+aW0d9+5oF\n3MvK3Jehk7rJycnB5XIRW3m3ASA2NpacnBybqpKqjBw5kptuuon4+Hh2797NjBkzuO6661i3bp1H\nyyCK71mWxdSpU7nqqqvOjH/x1WfK0SGcnJxMREQEKSkppKSk2F1O0PnuO/cdV3xlxIgRZ/67d+/e\nDBw4kM6dO/Puu+9y++231+maZy+6MWSIeYg0NKMrbXPUq1cvLr30Urp27crKlStJrLzgtgTMpEmT\nyMzMZO3atT6/tqO7oxcuXMiSJUsUwHX0wgvw+OP+f5/o6GguvvjiOo+0zc4224RlZ/u4MDmjXbt2\nWJZFbm6u2/Hc3FzatWtnU1Xiifj4eFq3bq2R7Da56667+Oijj1i5ciXty/eGxHefKUeHsNRPoHqu\njh8/zu7du93+B/VGp06QkKD5v/4UHx9Pu3btWLFixZljBQUFfP311wxRl4OjHTx4kO+//77Ony+p\nu7vuuosPP/yQL774ggsuuMDta776TOmfvQbi9Gn48ksz2Km+7rvvPkaNGkXnzp05dOgQjz766Jnb\nBp6q3P3ctCm89Vb962roioqK2LVrF9Z/1t/Zs2cPmzdvJiYmhk6dOjF16lSefPJJunXrRpcuXXj4\n4Yfp2LEjN1Te+0/8rqafU0xMDLNmzeKmm26iXbt27Nq1iwceeICLL77Y7TaQ+N+kSZNIS0tjyZIl\nNGvW7EyLNzo6mqj/7Gvok8+Ur4Zv+5KmKPnevHmW1bixZR08WP9rJScnW3FxcVZUVJTVqVMnKyUl\nxdqzZ49Hry0rs6zx4y1r9uz61yHuVq5cablcLissLMztcfvtt58559FHH7Xat29vNWnSxBo+fLi1\nc+dOGytumGr6OZ08edIaMWKEFRsbazVu3NiKj4+3Jk6caP373/+2u+wGp6qfUVhYmPXGG2+4nVff\nz5SWrWwgSkth82a4/HK7K4GZM80G6XUcwyUiEjIUwuJ3p0/rfq+ISFU0MKuBys+HN94w92b96amn\nYOhQM99XRETcKYQbqL/9DaZOhbNG1/vclVeaNZ+d198iImI/dUc3YIcPgy9nPWRkQHo6/OY3vrum\niEgoU0u4ATs7gJcuhQ8+qPv13n/fdD9rvWcREc8ohOWMDz44d75uQQFs2QJn78z15JPw6qvux+69\nF7Zu1VrPIiKecvSYVa0dHVivvgqVtsoEYO1auO46OHgQ4uIqjh86BGFn/QrXtKn/axQRCSW6Jyw1\nys83azr37QuRkXZXIyISWhzdEhb7RUfDwIF2VyEiEpp0T1hERMQmCmERERGbKIRFRERsohAWERGx\niUJYRETEJo6comRZFoWFhTRv3hxX+c7vIiIiIcaRISwiItIQqDtaRETEJgphERERmyiERUREbKIQ\nFhERsYlCWERExCYKYREREZsohEVERGzy/wHRAXqlGiJzsAAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# オイラー法での解法、インタラクティブにdtを変更できるようにした例\n", "from pylab import arange\n", "A = matrix([[0, 1], [-1, -1]])\n", "B = vector([0, 1])\n", "u = 1\n", "tf = 20\n", "dt = 0.05\n", "\n", "# 微分方程式をEuler法で解く\n", "x = vector([0, 0])\n", "xx = []\n", "for t in arange(0, tf, dt):\n", " xx.append(x.list())\n", " dx = A*x + B*u\n", " x = x + dx*dt\n", "# 結果をプロット\n", "t = arange(0, tf, dt)\n", "x1_lst = [v[0] for v in xx]\n", "x2_lst = [v[1] for v in xx]\n", "lst_plt1 = list_plot(zip(t, x1_lst), plotjoined=True, linestyle ='-')\n", "lst_plt2 = list_plot(zip(t, x2_lst), plotjoined=True, linestyle =':')\n", "(lst_plt1 + lst_plt2).show(figsize=5)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 7.2", "language": "", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.10" } }, "nbformat": 4, "nbformat_minor": 0 }