15.1 微分方程模型.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 微分方程模型\n",
    "\n",
    "**用sympy和scipy**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU9b3/8dcHwhaWIISENQmr7GsIi0tV3NHiLggIiOLSant7f63W5dar3lbb3qpdvMoVAQkguCBc61KK4lYTCPu+k4QQSFgSlpD9+/tj5rZcGzVkJjmZmffz8chjtjM5n29m5s3hnM98jznnEBGR8NLA6wJERCT4FO4iImFI4S4iEoYU7iIiYUjhLiIShqK8LgAgNjbWJSUleV2GiEhIWbNmzRHnXLuqHqsX4Z6UlERGRobXZYiIhBQzy/ymx7RbRkQkDCncRUTCkMJdRCQMKdxFRMKQwl1EJAwp3EVEwpDCXUQkDCncRUQ8kFt4ht8t38nuvFO18vvrxZeYREQiQWWl48s9R0hNy+Sv2/KodI52LZvQI65F0NelcBcRqWUFRaW8teYA89Oz2HfkNG2aN+aei7oxcUQCXdpE18o6Fe4iIrVkQ3YB89Iy+Z8NBykpr2RY4nn8aExPrhnQniZRDWt13Qp3EZEgOlNawf9sOMi8tEw25RQS3bghNw/rzKQRifTt2KrO6lC4i4gEwd78U6SmZfHWmmxOFJfTK74FT43rx41DOtGyaaM6r0fhLiJSQ+UVlfx122FS07L4YvcRohoYV/dvz+SRiaR0bYOZeVabwl1E5BwdPlHMwlVZvLEqm0MniukY05T/d2UvbhvehbiWTb0uD1C4i4hUi3OOr/YcZV5aJn/ZepiKSsfFvdrx9A39ufT8dkQ1rF9fG1K4i4h8i8IzZby95gCp6ZnszT9N6+hGTL+wK3ekJJAU29zr8r6Rwl1EpAqbcwqZ91UmSzfkUFxWyZCE1vznrYMYO7ADTRvVbhtjMHxnuJvZa8B1QJ5zrr//vjbAIiAJ2A/c5pw7br6jBy8C1wJFwFTn3NraKV1EJLiKyyp4b2Mu89Iy2ZBdQLNGDblxSCcmjkikf6cYr8s7J9XZcp8D/BF4/az7HgFWOOeeNbNH/LcfBq4Bevp/RgD/5b8UEam39h85zfz0TN5cc4CCojK6t2vOk9f35cahnYlpVvdtjMHwneHunPvMzJK+dvc44BL/9bnASnzhPg543TnngDQza21mHZxzucEqWEQkGMorKvl4ex7z0jL5fJevjfGqfu2ZODKBUd3aetrGGAw13ece/7+B7ZzLNbM4//2dgOyzljvgv++fwt3MZgAzABISEmpYhojIuck7WcyiVdksWJVFbmEx7Vs15SdX9GL88C7EtaofbYzBEOwDqlX9U+eqWtA5NxOYCZCcnFzlMiIiweCcI33fMealZfLR5kOUVzou6hnLL67vx+V94updG2Mw1DTcD//v7hYz6wDk+e8/AHQ5a7nOwMFAChQRqakTxWUsWZtDalomu/JOEdOsEVNHJzFxZCJd63EbYzDUNNyXAVOAZ/2XS8+6/4dm9ga+A6mF2t8uInVty8FCUtOyWLo+h6LSCgZ1juE3twzk+kEdQ6KNMRiq0wq5EN/B01gzOwD8Al+oLzaz6UAWcKt/8ffxtUHuxtcKOa0WahYR+SfFZRW8vymX1LRM1mYV0CSqAeMGd2TSyEQGdm7tdXl1rjrdMhO+4aExVSzrgB8EWpSISHVlHS1ifnomizOyOV5URrfY5jxxXV9uGdqZmOjQbGMMBn1DVURCTkWl45PteaSmZ/LpznwamHFFn3gmj0pkdPfQb2MMBoW7iISM/JMlLM7IZkF6FjkFZ4hv1YSHLuvJhJQE2seETxtjMCjcRaRec86xev9xUtMy+WBzLmUVjtHd2/L42D5c3jeeRmHYxhgMCncRqZdOFpfx7rocUtOy2HH4JC2bRjFpZCITRyTSI66F1+XVewp3EalXtuWeIDUtk3fX5XC6tIL+nVrx3M0DuH5QR6IbK7KqS38pEakXPtpyiFc/38vq/cdpEtWA6wZ2ZPKoRAZ1jtEB0hpQuIuIp5xz/OHj3fxu+U6S2kbz2LV9uGVYZ85r3tjr0kKawl1EPFNZ6Xj6z1uZ/eV+bh7ameduHhCW87x4QeEuIp4or6jkZ29v5J21Odx1QVceH9uHBg20+yVYFO4iUueKyyr44YJ1/HXbYf71il788LIe2q8eZAp3EalTJ4vLuOf1DNL2HuOpcf24c1SS1yWFJYW7iNSZo6dKmDp7NdtyT/Di+MGMG9zJ65LClsJdROrEwYIzTJqVTs7xM8y8cxiX9Y73uqSwpnAXkVq3J/8Uk19N52RxOfOmjyClaxuvSwp7CncRqVWbcwq587VVNDBYOGMk/TvFeF1SRFC4i0itSdt7lLvnZhDTrBGpd48I+1Pb1ScKdxGpFcu3HuYHC9aS0CaaedNT6BDTzOuSIorCXUSC7p21B/jpWxvp37EVs6el0EZTCdQ5hbuIBNVrX+zjqfe2Mrp7W2bemUyLJooZL+ivLiJB4Zzjhb/u4sUVu7iqXzwvjh9C00YNvS4rYincRSRglZWOp97bypy/7efWYZ351U2aAMxrCncRCUhZRSU/e2sjS9blcM9FXXn02j6aJ6YeCOifVjP7kZltNrMtZvZj/31tzGy5me3yX54XnFJFpL4pLqvgvnlrWLIuh59edb6CvR6pcbibWX/gHiAFGARcZ2Y9gUeAFc65nsAK/20RCTMnisu487VVfLwjj2du6M8PLtXMjvVJIFvufYA051yRc64c+BS4ERgHzPUvMxe4IbASRaS+OXKqhAkz01ibeZzfjx/CpJGJXpckXxNIuG8GLjaztmYWDVwLdAHinXO5AP7LuKqebGYzzCzDzDLy8/MDKENE6tKB40Xc9vJX7Mk/xatTkrl+UEevS5Iq1PiAqnNum5k9BywHTgEbgPJzeP5MYCZAcnKyq2kdIlJ3duedZPKsVZwqKSd1+giSkzQBWH0V0AFV59ws59xQ59zFwDFgF3DYzDoA+C/zAi9TRLy28UABt778FWUVjkUzRinY67lAu2Xi/JcJwE3AQmAZMMW/yBRgaSDrEBHv/W3PESbMTKN5kyjeum8UfTu28rok+Q6B9rm/bWZtgTLgB86542b2LLDYzKYDWcCtgRYpIt75aMshHly4jqS20bx+1wjaxzT1uiSphoDC3Tl3URX3HQXGBPJ7RaR+eDMjm4ff3sjAzq2ZM204raM1AVio0DdURaRKr36+l2f+vI0Le8TyyuRhNNcEYCFFr5aI/B/OOX63fCd/+Hg31/RvzwvjB9MkShOAhRqFu4j8XWWl4xfLtjAvLZPbk7vwy5sG0LCBvnUaihTuIgL4JgD718UbWLbhIPd+rxuPXN1b0wmEMIW7iHCmtIIH5q/hkx35PHx1b+6/pLvXJUmAFO4iEa7wTBl3z11NRuZxfnXTACakJHhdkgSBwl0kguWfLOHO11axO+8kf5wwlLEDO3hdkgSJwl0kQmUfK2LyrHQOnyhh1pThXNyrndclSRAp3EUi0K7DJ5k0K53iskpS7x7BsESdUyfcKNxFIsz67AKmzl5Fo4YNWHTvSHq31zwx4UjhLhJBvth1hBnzMoht0YTU6SNIaBvtdUlSSxTuIhHiw825PLRwPV1jmzNvegpxrTQBWDhTuItEgMWrs3nknY0M7tKa2VNTiIlu5HVJUssU7iJhbuZne/jl+9u5qKdvArDoxvrYRwK9yiJhyjnHbz7awUsr9zB2YAeev20wjaMCOj+PhBCFu0gYqqh0PLF0MwvSs5iQksAzN/TXBGARRuEuEmZKyyv5yeL1vLcxlwcu6c5PrzpfE4BFIIW7SBgpKi3n/tS1fLozn0ev7c2MizUBWKRSuIuEicKiMu6au5p1Wcd57uYB3D5cE4BFMoW7SBjIO1HMna+tYm/+aV6aOJSr+2sCsEincBcJcVlHi5g0K50jp0p4bepwLuwZ63VJUg8o3EVC2I5DJ5k8K53Sikrm3z2CIQmaAEx8Amp6NbN/MbMtZrbZzBaaWVMz62pm6Wa2y8wWmVnjYBUrIv+wNus4t73yFWaw+N5RCnb5P2oc7mbWCXgISHbO9QcaAuOB54DnnXM9gePA9GAUKiL/8PmufCb+dzqtoxvx1n2j6RXf0uuSpJ4J9OtqUUAzM4sCooFc4DLgLf/jc4EbAlyHiJzl/U253DVnNYlto3nzvlF0aaOZHeWf1TjcnXM5wG+BLHyhXgisAQqcc+X+xQ4AnQItUkR8Fq7K4ocL1jKoc2sW3TuKuJaa2VGqFshumfOAcUBXoCPQHLimikXdNzx/hpllmFlGfn5+TcsQiRj/tXIPP39nExf3ase86SOIaaaZHeWbBbJb5nJgn3Mu3zlXBrwDjAZa+3fTAHQGDlb1ZOfcTOdcsnMuuV07nbtR5Js45/jVB9t47sPtXD+oIzMnJ9OscUOvy5J6LpBwzwJGmlm0+SauGANsBT4BbvEvMwVYGliJIpGrotLx6JJNvPLpXiaNTOCF2zWzo1RPIPvc0/EdOF0LbPL/rpnAw8BPzGw30BaYFYQ6RSJOSXkFDy1cx8JV2Tx4WQ+eHqeZHaX6AvoSk3PuF8Avvnb3XiAlkN8rEulOl5RzX+oaPt91hMfH9uHui7p5XZKEGH1DVaSeOXqqhOlzM9h4oIDf3DKQW5O7eF2ShCCFu0g9su/IaabOXsWhwmJenjSMK/u197okCVEKd5F6Yk3mce6euxozY+GMkQzVdAISAIW7SD3wwaZcfrxoPR1imjJnWgpJsc29LklCnMJdxGOzvtjHM3/eypAurXl1ynDaNNdcexI4hbuIRyoqHc/8eSuzv9zP1f3a88L4wTRtpC8nSXAo3EU8UFxWwY/fWM+HWw5x1wVdeWxsH/WwS1Ap3EXq2LHTpdw9dzXrsgt44rq+TL+wq9clSRhSuIvUof3+VsfcwmJeumMo1wzQuU6ldijcRerIuqzjTJ+bgXOOBfeMYFhiG69LkjCmcBepAx9tOcSP3lhHfCtfq2NXtTpKLVO4i9SyOV/u49/f28qgzq2ZNSWZti2aeF2SRACFu0gtqax0/PL9bbz6xT6u7BvPi+OHaB52qTMKd5FaUFxWwU8Wr+f9TYeYOjqJJ67rq1ZHqVMKd5EgO366lHtezyAj8ziPj+3D9Au74jufjUjdUbiLBFHW0SKmzl7FgYIz/OmOoYwdqFZH8YbCXSRI1mcXMH3OaiqcY/7dIxiepFZH8Y7CXSQIlm89zIML19KuZRPmTEuhe7sWXpckEU7hLhKg17/az5PLtjCgUwyzpg4nVq2OUg8o3EVqqLLS8dyH23nls71c3iee308YTHRjfaSkftA7UaQGissq+H9vbuC9jbncOSqRX1zfT62OUq8o3EXOUUFRKTNeX8Oq/cd49Nre3HNRN7U6Sr2jcBc5B9nHipgyexUHjp3hDxOGcP2gjl6XJFKlBjV9opmdb2brz/o5YWY/NrM2ZrbczHb5L3WWXwkLGw8UcONLX3L0VCmpd49QsEu9VuNwd87tcM4Nds4NBoYBRcAS4BFghXOuJ7DCf1skpK3YdpjbX0mjaaOGvH3/aFK6qodd6rcah/vXjAH2OOcygXHAXP/9c4EbgrQOEU+kpmVyz+sZ9IxvwZIHLqBHnHrYpf4L1j738cBC//V451wugHMu18ziqnqCmc0AZgAkJCQEqQyR4KmsdPz6ox28/OkexvSO4w93DFGro4SMgLfczawx8H3gzXN5nnNupnMu2TmX3K5du0DLEAmqkvIKfrxoPS9/uoeJIxJ4ZfIwBbuElGC8W68B1jrnDvtvHzazDv6t9g5AXhDWIVJnCovKmDEvg/R9x3j46t7c9z21OkroCcY+9wn8Y5cMwDJgiv/6FGBpENYhUieyjxVx88t/Y11WAS+OH8z9l3RXsEtICmjL3cyigSuAe8+6+1lgsZlNB7KAWwNZh0hd2XSgkLvmrqakrILXp6cwsltbr0sSqbGAwt05VwS0/dp9R/F1z4iEjE+25/GDBWs5L7oxC+8ZQY+4ll6XJBIQHSGSiLcgPYsnlm6mT4eWvDZlOHGtmnpdkkjAFO4SsZxz/PYvO/jTJ3u49Px2/PGOoTRvoo+EhAe9kyUilZZX8rO3NvDu+oNMSOnC0+P6E9UwWN/pE/Gewl0iTuGZMu6bt4av9h7lp1edzwPqiJEwpHCXiJJTcIapr61i/9HTvHD7YG4Y0snrkkRqhcJdIsbmnELumrOaM2UVzL0rhdHdY70uSaTWKNwlIqzckccP5q8lplkj3r5/NL3i1eoo4U3hLmFv0eosHl2ymfPjWzJ72nDi1eooEUDhLmHLOcfzy3fy+493c3Gvdrw0cSgt1OooEULvdAlLpeWVPPL2Rt5Zl8PtyV145sb+NFKro0QQhbuEnRPFvlbHv+05yr9e0YsfXtZDrY4ScRTuElYOFpxh2uzV7Mk/xX/eOoibh3X2uiQRTyjcJWxsPXiCaXNWUVTia3W8oIdaHSVyKdwlLHy2M58H5q+lZdMo3rx/FL3bt/K6JBFPKdwl5C3OyObRdzbRI64Fc6al0D5GrY4iCncJWc45XvjrLl5csYuLesby0sShtGzayOuyROoFhbuEpNLySn7+zibeXnuAW4d15pc3DVCro8hZFO4Sck4Wl3F/6lq+2H2Ef7m8Fw+NUaujyNcp3CWk5Bb6Wh13553it7cO4ha1OopUSeEuIWNb7gmmzV7NqZJyZk8bzkU923ldkki9pXCXkPDFriPcn7qG5k2iePO+UfTpoFZHkW+jcJd6rbLSMX9VFv++bAs94lowe9pwOsQ087oskXovoHA3s9bAq0B/wAF3ATuARUASsB+4zTl3PKAqJSLtOHSSx5ZsIiPzOBf1jOVPE4fSSq2OItUSaO/Yi8CHzrnewCBgG/AIsMI51xNY4b8tUm1FpeX86oNtjP395+zJP8VvbhnI63elKNhFzkGNt9zNrBVwMTAVwDlXCpSa2TjgEv9ic4GVwMOBFCmRY8W2w/zb0i3kFJzhtuTO/PyaPpzXvLHXZYmEnEB2y3QD8oHZZjYIWAP8CIh3zuUCOOdyzSyuqieb2QxgBkBCQkIAZUg4yC08w5PLtvDRlsP0jGvB4ntHkdK1jddliYSsQMI9ChgKPOicSzezFzmHXTDOuZnATIDk5GQXQB0SwsorKpnzt/08v3wnFc7xs6vP5+4Lu9E4St82FQlEIOF+ADjgnEv3334LX7gfNrMO/q32DkBeoEVKeFqfXcCj72xia+4JLj2/HU+N60+XNtFelyUSFmoc7s65Q2aWbWbnO+d2AGOArf6fKcCz/sulQalUwkbhmTJ+89F25qdnEdeyCf81cShX92+vKQREgijQPvcHgflm1hjYC0zD14Gz2MymA1nArQGuQ8KEc45lGw7y9HvbOHa6hKmjk/jJFb00k6NILQgo3J1z64HkKh4aE8jvlfCz78hp/m3pZj7fdYSBnWOYM204/TvFeF2WSNjSN1SlVpWUV/Dyyr38aeVumjRswFPj+jFxRCING2gXjEhtUrhLrfnb7iM8/u5m9h45zfWDOvLE2D7EtdJZkkTqgsJdgi7/ZAm/fH8bS9blkNg2mrl3pfC9XprBUaQuKdwlaCorHQtXZ/HcB9s5U1bBQ5f14IFLe9C0UUOvSxOJOAp3CYptuSd4dMkm1mUVMLJbG565YQA94lp4XZZIxFK4S0BOl5Tz4opdzPpiHzHNGvG72wZx45BO6lkX8ZjCXWrsL1sO8eSyLRwsLGZCShcevro3raM1yZdIfaBwl3OWU+Cb5Gv51sOcH9+StyYMITlJk3yJ1CcKd6m2sopKZn+5j+eX7wLg59f05q4Lu9KooSb5EqlvFO5SLWsyj/PYkk1sP3SSy/vE8eT3+9H5PE3yJVJfKdzlWxUWlfHsh9tZuCqLDjFNeWXyMK7sG68DpiL1nMJdquSc4931OTzz3jYKzpRx94Vd+ZcretG8id4yIqFAn1T5J3vyT/H4ks18tfcog7u05vUb+9Ovoyb5EgklCnf5u+KyCl5auYeXV+6hSaMGPHNDf+5ISaCBJvkSCTkKdwHg8135PPHuZvYfLWLc4I48PrYv7Vo28bosEakhhXuEyztZzDPvbWPZhoN0jW1O6vQRXNgz1uuyRCRACvcIVVHpWLAqi19/uJ2Sskp+NKYn91/SXZN8iYQJhXsE2pxTyGPvbmZDdgEX9GjL0+P6062dJvkSCScK9whyqqSc55fvZPaX+2jTvDEv3D6YcYM7qmddJAwp3COAc46PthziyWVbOXyymDtSEvjZVb2JidaJqUXClcI9zGUfK+LJZVtYsT2PPh1a8dKkoQxNOM/rskSklincw1RZRSWvfr6PF1fspIEZj4/tw9TRSURpki+RiKBwD0MZ+4/x2JLN7Dh8kiv7xvPk9/vRsXUzr8sSkToUULib2X7gJFABlDvnks2sDbAISAL2A7c5544HVqZUx/HTpTz34XbeWJ1Np9bN+O87k7mib7zXZYmIB4Kx5X6pc+7IWbcfAVY45541s0f8tx8OwnrkGzjneHttDr98fxuFZ8q49+JuPDSmpyb5EolgtfHpHwdc4r8+F1iJwr3W7M47yWNLNpO+7xhDE1rzHzcOoE+HVl6XJSIeCzTcHfAXM3PAK865mUC8cy4XwDmXa2ZxVT3RzGYAMwASEhICLCPyFJdV8MePd/PKZ3uIbhzFr24awO3JXTTJl4gAgYf7Bc65g/4AX25m26v7RP8/BDMBkpOTXYB1RJRPd/om+co6VsRNQzrx6Ng+xLbQJF8i8g8Bhbtz7qD/Ms/MlgApwGEz6+Dfau8A5AWhzohXXlHJiu15zPsqky92H6Fbu+YsuGcEo7trki8R+Wc1Dnczaw40cM6d9F+/EngKWAZMAZ71Xy4NRqGRKu9EMW+szmbhqixyC4vpENOUh6/uzV0XJtEkSpN8iUjVAtlyjweW+OcliQIWOOc+NLPVwGIzmw5kAbcGXmZkcc6RtvcYqWmZfLTlEOWVjot6xvLk9/sxpnecvogkIt+pxuHunNsLDKri/qPAmECKilQnist4Z80BUtOz2J13iphmjZh2QRJ3jEika2xzr8sTkRCiRuh6YHNOIfPTM3l33UHOlFUwqEtrfnvrIK4b2EHzq4tIjSjcPVJcVsH7m3KZl5bJuqwCmjZqwLhBnZg0MpEBnXUyahEJjMK9jmUePc2C9CwWZ2RzvKiMbrHN+bfr+nLz0M6agldEgkbhXgcqKh0fb88jNS2TT3fm07CBcWXfeCaNTGR097Y6WYaIBJ3CvRblnyxhcUY2C9KzyCk4Q3yrJvz48p6MH55A+5imXpcnImFM4R5kzjlW7TtGanoWH27OpazCcUGPtjxxXR/G9ImnkdoYRaQOKNyD5GRxGUvW5ZCalsnOw6do1TSKySOTmDgyge46+bSI1DGFe4C2HjxBanom767Loai0ggGdYvj1zQO5flBHmjVWG6OIeEPhXgMl5RV8sOkQ89IyWZN5nCZRDbh+UEcmj0xkUJfWXpcnIqJwPxfZx4qY729jPHa6lKS20Tw+tg+3DOtM6+jGXpcnIvJ3CvfvUFHp+HSnbzbGlTvzMeDyPvFMHpXIBd1jNX+6iNRLCvdvcOSUr41xfpqvjbFdyyY8eGkPxqck6GTTIlLvKdzP4pxjTeZx5qVl8sGmQ5RWVDKqW1sevbYPV/ZTG6OIhA6FO3CqpJx3/W2M2w+dpGWTKO4YkcCkkQn0iGvpdXkiIucsosN9x6GTpKZlsmRdDqdKyunboRW/umkA4wZ3JLpxRP9pRCTERVyClZZX8sHmXOanZbFq/zEaRzXguoEdmDQykSFdWmueFxEJCxET7geOF7FwVRaLVmdz5FQpCW2iefTa3twyrAttmquNUUTCS1iHe2Wl49Nd+cxPy+Tj7b7zdF/W29fGeFEPtTGKSPgKy3A/drr077MxZh0rIrZFYx64pAcTRiTQSW2MIhIBwibcnXOszSogNS2TP2/KpbS8kpSubfjpVedzVb/2NI5SG6OIRI6QD/fTJeUsXX+Q1LRMtuaeoEWTKMYP78LEEYmc315tjCISmUI63BetzuKZ97ZxsqSc3u1b8h839mfc4E60aBLSwxIRCVjAKWhmDYEMIMc5d52ZdQXeANoAa4HJzrnSQNdTlU6toxnTJ47JoxIZmnCe2hhFRPyCsSP6R8C2s24/BzzvnOsJHAemB2EdVbqwZywvjB/CsMQ2CnYRkbMEFO5m1hkYC7zqv23AZcBb/kXmAjcEsg4RETl3gW65vwD8DKj0324LFDjnyv23DwCdqnqimc0wswwzy8jPzw+wDBEROVuNw93MrgPynHNrzr67ikVdVc93zs10ziU755LbtWtX0zJERKQKgRxQvQD4vpldCzQFWuHbkm9tZlH+rffOwMHAyxQRkXNR4y1359zPnXOdnXNJwHjgY+fcROAT4Bb/YlOApQFXKSIi56Q2vrb5MPATM9uNbx/8rFpYh4iIfIugfNvHObcSWOm/vhdICcbvFRGRmtGEKyIiYcicq7KZpW6LMMsHMmv49FjgSBDLCQUac2TQmCNDIGNOdM5V2W5YL8I9EGaW4ZxL9rqOuqQxRwaNOTLU1pi1W0ZEJAwp3EVEwlA4hPtMrwvwgMYcGTTmyFArYw75fe4iIvLPwmHLXUREvkbhLiIShkIm3M3sajPbYWa7zeyRKh5vYmaL/I+nm1lS3VcZXNUY80/MbKuZbTSzFWaW6EWdwfRdYz5ruVvMzJlZyLfNVWfMZnab/7XeYmYL6rrGYKvGezvBzD4xs3X+9/e1XtQZLGb2mpnlmdnmb3jczOz3/r/HRjMbGvBKnXP1/gdoCOwBugGNgQ1A368t8wDwsv/6eGCR13XXwZgvBaL91++PhDH7l2sJfAakAcle110Hr3NPYB1wnv92nNd118GYZwL3+6/3BfZ7XXeAY74YGAps/obHrwU+wDdt+kggPdB1hsqWewqw2zm31/nOx/oGMO5ry4zDd+Yn8J0JaoyF9rn3vnPMzrlPnHNF/ptp+KZYDmXVeZ0BngZ+DU3+z1oAAAJxSURBVBTXZXG1pDpjvgf4k3PuOIBzLq+Oawy26ozZ4ZtGHCCGEJ863Dn3GXDsWxYZB7zufNLwTZ3eIZB1hkq4dwKyz7pd1Rme/r6M880lX4hvVspQVZ0xn206vn/5Q9l3jtnMhgBdnHPv1WVhtag6r3MvoJeZfWlmaWZ2dZ1VVzuqM+YngUlmdgB4H3iwbkrzzLl+3r9TUGaFrAPVOcNTtc8CFSKqPR4zmwQkA9+r1Ypq37eO2cwaAM8DU+uqoDpQndc5Ct+umUvw/e/sczPr75wrqOXaakt1xjwBmOOc+08zGwXM84+5sornhoOg51eobLkfALqcdbuqMzz9fRkzi8L3X7lv+29QfVedMWNmlwOPAd93zpXUUW215bvG3BLoD6w0s/349k0uC/GDqtV9by91zpU55/YBO/CFfaiqzpinA4sBnHNf4TvbW2ydVOeNan3ez0WohPtqoKeZdTWzxvgOmC772jLL8J35CXxngvrY+Y9UhKjvHLN/F8Ur+II91PfDwneM2TlX6JyLdc4lOd8ZwNLwjT3Dm3KDojrv7XfxHTzHzGLx7abZW6dVBld1xpwFjAEwsz74wj2/TqusW8uAO/1dMyOBQudcbkC/0eujyOdwtPlaYCe+o+yP+e97Ct+HG3wv/pvAbmAV0M3rmutgzH8FDgPr/T/LvK65tsf8tWVXEuLdMtV8nQ34HbAV2ASM97rmOhhzX+BLfJ0064Erva45wPEuBHKBMnxb6dOB+4D7znqN/+T/e2wKxvta0w+IiIShUNktIyIi50DhLiIShhTuIiJhSOEuIhKGFO4iImFI4S4iEoYU7iIiYej/A00ekQ5JJWrrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 一些参数. 来自15.2\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "x0 = np.array([41, 49, 61, 78, 96, 104], dtype=np.float)\n",
    "plt.plot(np.linspace(0, 1, x0.shape[0]), x0)\n",
    "plt.show()\n",
    "lamb = np.array([x0[i] / x0[i+1] for i in range(x0.shape[0]-1)])\n",
    "Theta = (np.exp(-2/(x0.shape[0]+1)), np.exp(2/(x0.shape[0]+2)))\n",
    "OK = [i <= Theta[1] and i >= -Theta[0] for i in lamb]\n",
    "x1 = np.array([x0[:i].sum() for i in range(1, x0.shape[0]+1)])\n",
    "alpha1_x0 = np.array([x0[i] - x0[i-1] for i in range(1, x0.shape[0])])\n",
    "z1 = np.array([(x1[i]+x1[i+1])/2 for i in range(x1.shape[0]-1)])\n",
    "B = np.array(np.c_[-x0[1:], -z1, np.ones(z1.shape[0])])\n",
    "Y = alpha1_x0\n",
    "u = np.matmul(np.linalg.inv(np.matmul(B.T, B)), np.matmul(B.T, Y))\n",
    "a1, a2, b = u[0], u[1], u[2]\n",
    "x1_t = lambda k: (x0[0] - b/a) * np.exp(-a*k) + b/a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "得到带参符号解: Eq(x1_t(t), C1*exp(t*(-a1 - sqrt(a1**2 - 4*a2))/2) + C2*exp(t*(-a1 + sqrt(a1**2 - 4*a2))/2) + b/a2)\n",
      "待定常数:  {C2, C1}\n",
      "关于待定常数的多元方程组:\n",
      "[Eq(41, C1 + C2 + b/a2), Eq(x1_t(6), C1*exp(-3*a1 - 3*sqrt(a1**2 - 4*a2)) + C2*exp(-3*a1 + 3*sqrt(a1**2 - 4*a2)) + b/a2)]\n",
      "解得待定常数: {C1: ((41*a2 - b)*exp(3*a1 + 6*sqrt(a1**2 - 4*a2)) - (a2*x1_t(6) - b)*exp(6*a1 + 3*sqrt(a1**2 - 4*a2)))*exp(-3*a1)/(a2*(exp(6*sqrt(a1**2 - 4*a2)) - 1)), C2: -((41*a2 - b)*exp(3*a1) - (a2*x1_t(6) - b)*exp(6*a1 + 3*sqrt(a1**2 - 4*a2)))*exp(-3*a1)/(a2*(exp(6*sqrt(a1**2 - 4*a2)) - 1))}\n",
      "最终解: x1_t(t) = b/a2 - ((41*a2 - b)*exp(3*a1) - (a2*x1_t(6) - b)*exp(6*a1 + 3*sqrt(a1**2 - 4*a2)))*exp(-3*a1)*exp(t*(-a1 + sqrt(a1**2 - 4*a2))/2)/(a2*(exp(6*sqrt(a1**2 - 4*a2)) - 1)) + ((41*a2 - b)*exp(3*a1 + 6*sqrt(a1**2 - 4*a2)) - (a2*x1_t(6) - b)*exp(6*a1 + 3*sqrt(a1**2 - 4*a2)))*exp(-3*a1)*exp(t*(-a1 - sqrt(a1**2 - 4*a2))/2)/(a2*(exp(6*sqrt(a1**2 - 4*a2)) - 1))\n"
     ]
    }
   ],
   "source": [
    "# sympy的一个例子\n",
    "# 参考: https://blog.csdn.net/and_w/article/details/80160033\n",
    "import sympy\n",
    "\n",
    "sympy.init_printing\n",
    "# 定义符号常量x 与 f(x)\n",
    "t = sympy.Symbol('t')\n",
    "a1, a2, b = sympy.symbols('a1, a2, b')\n",
    "x1_t = sympy.Function('x1_t')\n",
    "# 微分方程ode\n",
    "ode = sympy.Eq(x1_t(t).diff(t, t) + a1*x1_t(t).diff(t) + a2*x1_t(t), b)\n",
    "init1 = sympy.Eq(x1_t(0), x1[0])\n",
    "init2 = sympy.Eq(x1_t(1), x1[1])\n",
    "# 调用dsolve函数,返回一个Eq对象，hint控制精度\n",
    "ode_sol = sympy.dsolve((ode), x1_t(t))\n",
    "print('得到带参符号解:', ode_sol)\n",
    "\n",
    "# 初始条件字典\n",
    "x1_to_value = {x1_t(0): 41, x1_t(x1.shape[0]-1): x1[-1]}\n",
    "# 确定待定常数集\n",
    "known_params = {a1, a2, b, t}\n",
    "free_params = ode_sol.free_symbols - set(known_params)\n",
    "# 将初始条件带入微分方程的后得到的多元方程\n",
    "C_eqs = [sympy.Eq(ode_sol.lhs.subs(t, 0).subs(x1_to_value), ode_sol.rhs.subs(t, 0)),\n",
    "         sympy.Eq(ode_sol.lhs.subs(t, x1.shape[0]).subs(x1_to_value), ode_sol.rhs.subs(t, x1.shape[0]))\n",
    "        ]\n",
    "print('待定常数: ', free_params)\n",
    "print('关于待定常数的多元方程组:', C_eqs, sep='\\n')\n",
    "sol_params = sympy.solve(C_eqs, free_params)\n",
    "print('解得待定常数:', sol_params)\n",
    "final_sol_eq = sympy.Eq(ode_sol.lhs, ode_sol.rhs.subs(sol_params))\n",
    "print('最终解:', final_sol_eq.lhs, '=', final_sol_eq.rhs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df5DV9X3v8eeLBXSNPwBFIgsGo8RGY4ORone8c6/VhKXJnSvJmISkNzItM2SsuW2nudxCO3e05nrVkNRbm8aOqVS0iT+uscokGkqxSZrUoGswAhqEBKOwVNBlEcOC7PK+f3w/R84u5+yes3t+7dnXY+bM+e7nfD/f/Xx3z3ff+/l83+fzUURgZmbWaMbVuwFmZmaFOECZmVlDcoAyM7OG5ABlZmYNyQHKzMwa0vh6N6BRnHHGGTFr1qx6N8Ma3LPPPvt6REytdztGC19XVopi15UDVDJr1iw6Ojrq3QxrcJJ+Ve82jCa+rqwUxa4rD/GZmVlDcoAyM7OG5ABlNgpJapG0UdJ30tdTJK2TtC09T87bd4Wk7ZK2SmrPK79E0qb02h2SlMpPkPRgKt8gaVZencXpe2yTtLh2Z2xjUdUDlC8ks6r4I+DFvK+XA+sjYjawPn2NpAuARcCFwALg65JaUp07gaXA7PRYkMqXAPsi4jzgduC2dKwpwA3ApcA84Ib869es0mrRg/KFZFZBkmYAHwP+Lq/4amB12l4NLMwrfyAiDkfEDmA7ME/SWcCpEfFUZBNy3jugTu5YDwNXpX8K24F1EdEVEfuAdRy7Fs0qrqoByheSWVX8X+B/AkfzyqZFxG6A9HxmKm8DXs3bb2cqa0vbA8v71YmIXmA/cPogx+pH0lJJHZI69u7dO5zzMwOq34PyhWRWQZL+C7AnIp4ttUqBshikfLh1jhVE3BURcyNi7tSp/shYo3t04y4uv/VJzln+XS6/9Uke3bir3k16R9U+B5V/IUm6opQqBcqqfiEBdwHMnTvX647YaHA58F8lfRQ4EThV0j8Ar0k6KyJ2p1GHPWn/ncDMvPozgM5UPqNAeX6dnZLGA6cBXan8igF1vl+5U7N8j27cxcq1W+ns7mH6pFaWtZ/PwouP+z97xN9jxSOb6DnSB8Cu7h5WPLIJoKzvVaytj27cxR8/+Nxx+79868dKOm41P6g7pi6k7m748Y/hX/8VduyAmTPhPe+BWbNg2jSYPBkmTcoeEyZUsyXWzCJiBbACIP3j9z8i4r9JWgksBm5Nz4+lKmuAb0n6S2A62T3cpyOiT9IBSZcBG4Brgb/Oq7MYeAq4BngyIkLSWuD/5N3PnZ9ri1VWpQLHUFau3frO98jpOdLHyrVbS/4+xdra8asu/uEnrxSsM2v5d0sKUlULUGPpQtq8GS67DH796yz4zJwJa9bAoUOF9584EcaPP/ZoaYFx40DKHrnt8ePhxBOzR2tr/+2JE+HwYejpgYMHs8eRIxBx7AHZsSZMOP4hwdGjxx59fYXbKh1rX/5zS8uxtueeI+Dtt7N25R69vceO39eXbeefZ+5cW1rghBOy88o9T5yYHfPQof6P3HFyz0ePHjtWfjsnTDj2M8t/9PVlv6u33sqeDx481t5c+3/wg+wfjFHkVuAhSUuAV4BPAkTEFkkPAS8AvcD1EZH7bV8H3AO0Ak+kB8DdwH2StpP9w7coHatL0peAZ9J+N0VEV7VPbCyqROAoRWd3T1nlhRRr6/0bXi1So3T1mOqo6S6kv/qr7I/k+vVZoDrppOwP65498Ktfwd69sG/fscfBg9kfyd7eY4+I7Bi54HL0aFY+8I/z3r3Z8+HD2R/yk07KAtYZZ2R/0HN//HOPvr4scOUehw9nf5ih/x/1XKAYKFc/PyjkP3p7jz1D1qbcY+LELEhMnNg/cMCxc80FmN7erF25IJELFOPG9Q8uucCVHyylY8fKD1pvv50dM//n19OT1XnXu+Dkk7Pnk07KerYDA2Sji4jvk0YGIuIN4Koi+90M3FygvAP4QIHyQ6TrssBrq4BVw22zlaYSgaMU0ye1sqvAMadPai35GMXa1FeB1dprEqCa+ULq7oZvfQt+93fhyiuPlUvZ0N60adVugZk1m0oEjlIsaz+/3/AcQOuEFpa1n1/yMYq1tUUacZDyTBIjdO+9WY/ouuvq3RIzaxbL2s+ndUJLv7JyA0cpFl7cxi2fuIi2Sa0IaJvUyi2fuKisYcRibf3MpTOL1CidZzMfgQj427+FSy+FD32o3q0xs2aRCxDVzuLLfa+RHHewts59z5SGzeJrej/4Abz4ItxzT71bYmbVVIuU74FGGjhqqVhbR3oODlAj8Pd/n6WPf+pT9W6JmVVLrVK+7Xi+BzUCL70El1ySZdGZWXMaLOXbqssBagQ6O2H69Hq3wsyqqVYp33Y8B6hhOnoUdu92gDJrdsVSuyud8m3Hc4AapjfeyD7A6gBl1txqlfJtx3OSxDB1ptkAHaDMmlstU76tPweoYXKAMhs7RlPKdzPxEN8wOUCZmVWXA9Qw5QLUWWfVtx1mZs3KAWqYOjth6tTRMeu1mdlo5AA1TLt2eXjPzKyaHKCGyR/SNTOrLgeoYXKAMjOrLgeoYejthddec4AyM6smB6hh2LMnm+rIAcrMrHqqFqAknSjpaUk/k7RF0l+k8hsl7ZL0XHp8NK/OCknbJW2V1J5XfomkTem1OyQplZ8g6cFUvkHSrLw6iyVtS4/FlTw3fwbKzKz6qjmTxGHgyoh4S9IE4EeSnkiv3R4RX8nfWdIFwCLgQmA68M+S3hcRfcCdwFLgJ8DjwALgCWAJsC8izpO0CLgN+LSkKcANwFwggGclrYmIfZU4MQcoM7Pqq1oPKjJvpS8npEcMUuVq4IGIOBwRO4DtwDxJZwGnRsRTERHAvcDCvDqr0/bDwFWpd9UOrIuIrhSU1pEFtYpwgDIzq76q3oOS1CLpOWAPWcDYkF76gqTnJa2SNDmVtQGv5lXfmcra0vbA8n51IqIX2A+cPsixBrZvqaQOSR179+4t+bw6O2HcODjzzJKrmJlZmaoaoCKiLyLmADPIekMfIBuuOxeYA+wGvpp2V6FDDFI+3Dr57bsrIuZGxNypU6cOei75Ojth2jQY76l2zcyqpiZZfBHRDXwfWBARr6XAdRT4BjAv7bYTmJlXbQbQmcpnFCjvV0fSeOA0oGuQY1WEPwNl9dLMyUdmA1Uzi2+qpElpuxX4MPDzdE8p5+PA5rS9BliULo5zgNnA0xGxGzgg6bJ0AV0LPJZXJ3eRXAM8me5TrQXmS5qchhDnp7KKeP31bB4+szrIJR99kGwUYoGky9Jrt0fEnPR4HI5LPloAfF1SbvW9XPLR7PTI3ad9J/kIuJ0s+Yi85KNLyf6xvCFviN6s4qo5SHUWsDpdDOOAhyLiO5LukzSHbMjtZeDzABGxRdJDwAtAL3B9yuADuA64B2gly97LZQPeDdwnaTtZz2lROlaXpC8Bz6T9boqIrkqd2MGDMHPm0PuZVVr6B2xYyUfAjnStzJP0Min5CEBSLvnoiVTnxlT/YeBrA5OPUp1c8tH9FTtBszxVC1AR8TxwcYHyzw1S52bg5gLlHcAHCpQfAj5Z5FirgFVlNLlkPT3Q2lqNI5sNLf3T9yxwHvA3EbFB0u+QJR9dC3QAX0wZrG1kH8/IySUMHaHE5CNJZScfkfXMOPvss0d2sjameSaJYejpgZNOqncrbKxq1uSjSnl04y4uv/VJzln+XS6/9Uke3bir5m2wynCAGoaDB92DsvprtuSjSnh04y5WPLKJXd09BLCru4cVj2xykBqlHKCGwUN8Vi/NnHxUCSvXbqXnSF+/sp4jfaxcu7VOLbKR8Cd5ytTXB2+/7SE+q5umTT6qhM7unrLKrbE5QJWpJ73P3YOyemjm5KNKmD6plV0FgtH0Sb5gRyMP8ZUpF6DcgzJrPMvaz6d1Qku/stYJLSxrP79OLbKRcA+qTAcPZs/uQZk1noUXZ1nvK9dupbO7h+mTWlnWfv475Ta6OECVyUN8Zo1t4cVtDkhNwkN8ZfIQn5lZbThAlclDfGZmteEAVSYP8ZmZ1YYDVJk8xGdmVhsOUGXyEJ+ZWW04QJXJQ3xmZrXhAFUmD/GZmdWGA1SZPMRnZlYbDlBl8hCfmVltOECV6eBBGD8eJkyod0vMzJqbA1SZvBaUmVltVC1ASTpR0tOSfiZpi6S/SOVTJK2TtC09T86rs0LSdklbJbXnlV8iaVN67Y60wBppEbYHU/kGSbPy6ixO32ObpMVUiJd7NzOrjWr2oA4DV0bEB4E5wAJJlwHLgfURMRtYn75G0gVkC6NdCCwAvp4WZQO4E1hKthro7PQ6wBJgX0ScB9wO3JaONQW4AbiUbOnrG/ID4Uh4uXczs9qoWoCKzFvpywnpEcDVwOpUvhpYmLavBh6IiMMRsQPYDsxLS1mfGhFPpWWn7x1QJ3esh4GrUu+qHVgXEV0RsQ9Yx7GgNiIe4jMzq42q3oOS1CLpOWAPWcDYAEyLiN0A6fnMtHsb8Gpe9Z2prC1tDyzvVycieoH9wOmDHGtg+5ZK6pDUsXfv3pLOyUN8Zma1UdUAFRF9ETEHmEHWGzpueek8KnSIQcqHWye/fXdFxNyImDt16tRBmnaMh/jMzGqjJll8EdENfJ9smO21NGxHet6TdtsJzMyrNgPoTOUzCpT3qyNpPHAa0DXIsUbMQ3xmZrVRzSy+qZImpe1W4MPAz4E1QC6rbjHwWNpeAyxKmXnnkCVDPJ2GAQ9IuizdX7p2QJ3csa4Bnkz3qdYC8yVNTskR81PZiHmIz8ysNqq55PtZwOqUiTcOeCgiviPpKeAhSUuAV4BPAkTEFkkPAS8AvcD1EdGXjnUdcA/QCjyRHgB3A/dJ2k7Wc1qUjtUl6UvAM2m/myKiqxIn5SE+M7PaqFqAiojngYsLlL8BXFWkzs3AzQXKO4Dj7l9FxCFSgCvw2ipgVXmtHpqH+MzMasMzSZTp4EEP8ZmZ1YIDVJncg7J6atYZWswKcYAqQ4QDlNVdU87QYlaIA1QZDh3Knj3EZ/XSrDO0mBXiAFUGrwVljaAZZ2gxK8QBqgxe7t0aQTPO0GJWiANUGbzcuzWSZpqhxawQB6gyeIjP6q1ZZ2gxK6SaM0k0HQ/xWQNoyhlazApxgCqDh/is3pp1hhazQjzEVwYP8ZmZ1Y4DVBlyPSgP8ZmZVZ8DVBncgzIzqx0HqDI4QJmZ1Y4DVBk8xGdmVjsOUGVwD8rMrHYcoMrQ0wMSnHBCvVtiZtb8HKDKkFvuXYVmJDMzs4qqWoCSNFPSv0h6MS2s9kep/EZJuyQ9lx4fzavT0AureS0oM7PaqeZMEr3AFyPip5JOAZ6VtC69dntEfCV/5wELq00H/lnS+9K0LLmF1X4CPE42OeYT5C2sJmkR2cJqn85bWG0u2WzLz0pak9awGbaeHidImJnVStV6UBGxOyJ+mrYPAC9SYO2YPA2/sFpuiM/MzKqvJveg0tDbxcCGVPQFSc9LWpW3ZHTDL6zmIT4zs9qpeoCSdDLwbeCPI+JNsuG6c4E5wG7gq7ldC1RvqIXVDh70EJ+ZWa1UNUBJmkAWnL4ZEY8ARMRraUXQo8A3gHlp94ZfWM09KDOz2qlmFp/I1pV5MSL+Mq/8rLzdPg5sTtsNv7CaA5SZWe1UM4vvcuBzwCZJz6WyPwM+I2kO2ZDby8DnYXQsrOYhPjOz2qlagIqIH1H4XtDjg9Rp6IXV3IMyM6sdzyRRBgcoM7PacYAqg4f4zMxqxwGqDO5BmZnVjgNUiY4cgb4+Bygzs1pxgCqRFys0M6stB6gSebFCM7PacoAqUa4H5QBlZlYbDlAlyvWgPMRn9dSM66yZFVPNmSSaiof4rEE03TprZsW4B1UiJ0lYI2jGddbMinGAKpF7UNZommWdNbNiHKBK5ABljaSZ1lkzK8YBqkQe4rNG0WzrrJkV4wBVIvegrBE04zprZsU4i69EDlDWIJpunTWzYhygSuQhPmsEzbjOmlkxHuIrkXtQVimSTpV0boHy36xHe8walQNUiQ4ehIkTYZx/YjYCkj4F/Bz4dpoJ4rfyXr6nPq0ya0xV+3M7yJQsUyStS1OlrMv7vEZDT8nS0+PhPauIPwMuiYg5wO+R3ev5RHqt0NCd2ZhVzf5AbkqW9wOXAdenaVeWA+sjYjawPn09cEqWBcDXJbWkY+WmZJmdHrlPr78zJQtwO9mULORNyXIpWbrtDfmBcDi8WKFVSEvKoCMingZ+G/hzSX9Igc8UmY1lQwaoFDgGll0xVL1BpmTJn0ZlNf2nV2nYKVkcoKxCDuTff0rB6gqy9/KF9WqUWSMqpQf1kKQ/VaZV0l8Dt5TzTQZMyTIt7z/I3cCZabeaT8lSDgcoq5DrGDCUl/6BWwD8fl1aZNagSglQl5J9evzfyD7/0En2WYySFJiSpeiuBcqqOiVLOXOGHToEJ5446C5mQ4qIn0XE9gLlRyLim7mvJT1V25aZNZ5SAtQRoIfsw3wnAjvSdCpDKjQlC/Ba7lPv6XlPKq/5lCzlzBnmHpTVmP8dsjGvlAD1DFmA+i3gP5J9Yv3hoSoVm5KF/tOoLKb/9CoNOyWLe1BWY06YsDGvlJkklqRPnAP8O3C1pM+VUK/YlCy3kt3XWgK8QvrEeqNPydLTA2ecMZIjmJlZOYYMUHnBKb/svhLqFZuSBeCqInUadkoW96CskiRdEBEvDCi7IiK+n/uy9q0yayyeF6FEhw75HpRV1FDZsaWMUpg1NQeoEvX0uAdlFTVodmxEbC5Sz2zMcIAqkYf4rMKGnR1rNlY4QJXIaeZWYcPKjjUbS7weVAmOHoW333YPyipquNmxZmOGe1AlOHQoe3YPyipluNmxZmOJA1QJcgHKPSgzs9pxgCqBA5SZWe05QJXAy72bmdWeA1QJ3IMyM6s9B6gSuAdlZlZ7DlAlcA/KzKz2/DmoEjjN3Kz6Ht24i5Vrt9LZ3cP0Sa0saz+fhRePaCFsG+XcgypBbojPPSirN0kzJf2LpBclbZH0R6l8iqR1kral58l5dVZI2i5pq6T2vPJLJG1Kr92R1lsjrcn2YCrfIGlWXp3F6Xtsk7SYCnl04y5WPLKJXd09BLCru4cVj2zi0Y27KvUtbBRygCqBh/isgfQCX4yI9wOXAddLugBYDqyPiNnA+vQ16bVFwIXAAuDrklrSse4ElpItDjo7vQ6wBNgXEecBtwO3pWNNAW4gm+h2HnBDfiAciZVrt9JzpK9fWc+RPlau3VqJw9so5QBVAidJWKOIiN0R8dO0fQB4EWgDrgZWp91WAwvT9tXAAxFxOCJ2ANuBeZLOAk6NiKfSKtT3DqiTO9bDwFWpd9UOrIuIrojYB6zjWFAbkc7unrLKbWxwgCqBe1DWiNLQ28XABmBaROyGLIgBZ6bd2oBX86rtTGVtaXtgeb86EdEL7AdOH+RYA9u1VFKHpI69e/eWdC7TJxX+769YuY0NDlAlcA/KGo2kk4FvA38cEW8OtmuBshikfLh1jhVE3BURcyNi7tSpUwdp2jHL2s+ndUJLv7LWCS0saz+/pPrWnBygSuAelDUSSRPIgtM3I+KRVPxaGrYjPe9J5TvJFkbMmUG2OOLOtD2wvF8dSeOB04CuQY41YgsvbuOWT1xE26RWBLRNauWWT1zkLL4xrmoBStIqSXskbc4ru1HSLknPpcdH815r2EwjByhrFOn9fzfwYkT8Zd5La4Dce30x8Fhe+aJ0vZxDlgzxdBoGPCDpsnTMawfUyR3rGuDJdJ9qLTBf0uSUHDE/lVXEwovb+PHyK9lx68f48fIrHZysqp+Dugf4GtnN13y3R8RX8gsGZBpNB/5Z0vsioo9jmUY/AR4nuyn7BHmZRpIWkWUafTov02gu2fDDs5LWpJu6w9LTAxMnwjj3N63+Lgc+B2yS9Fwq+zPgVuAhSUuAV4BPAkTEFkkPAS+QZQBen64rgOvIrtNWsmvqiVR+N3CfpO1kPadF6Vhdkr5EttgiwE0R0VWtEzWrWoCKiB/m92qG8E6mEbAjXRjzJL1MyjQCkJTLNHoi1bkx1X8Y+NrATKNUJ5dpdP9wz8XLvVujiIgfUfheEMBVRercDNxcoLwD+ECB8kOkAFfgtVXAqlLbazYS9egTfEHS82kIMPcZippnGkHp2UZe7t3MrPZqHaDuBM4F5gC7ga+m8ppnGkHp2UbuQZmZ1V5NA1REvBYRfRFxFPgG2afRocEzjdyDMjOrvZoGqFwabPJxIJfh19CZRu5BmZnVXtWSJCTdD1wBnCFpJ1lm3RWS5pANub0MfB4aP9Po0CH3oMzMaq2aWXyfKVB89yD7N2ymUU9PFXtQR96CA9ug7yD0HYa+Q3A093wk7ZR/Cy1tR6TtyNuupWKJZFCZtgw8/kiPOdityQLHPvvTMPG0EX5PMxsJrwdVgkOH4JRTKnSwro2w8x+hexN0Pw9v/bJCB7aKOvO3HaDM6swBqgQVS5LY/nfQ8QcQfXDKbJhyCbz39+DU98OEU6DlRBh3wrHncRNAuf/883sAaVtK28rbLkWxZMdSFenNRBRp7zCP3+94IzlmXnuLHnPAsU88EzOrLweoEow4SeLoEfjpF+Glv4Z3fwQuvx9OOL1i7TMza0YOUCUYUQ/q8Bvwo0/Ba0/Cb/wJzLkNxvnHbmY2FP+lLMGwe1AR8OSHYf+LcNlqeO+1FW+bmVmzcoAqwbDTzF//Cex7DuZ9w8HJzKxMnp97CBEjSDPfcS+0tMJ7Pl3xdpmZNTsHqCEcOZIFqbIDVN9heOVBmPHxLEPPzMzK4gA1hGEv9975XXh7H5zzuYq3ycxsLHCAGsKwV9PdcR+c+G5494cr3iYzs7HAAWoIw+pBHX4j60HN+qxTys3MhskBagjD6kH96sHsw7ke3jMzGzYHqCEMK0DtuA8mXQSTPliVNpmZjQUOUEMoe4jvzZfgjZ/ArM8NmPPNzMzK4QA1hLJ7UK9+O3ue9dmqtMfMbKxwgBpC2T2ofRvh5HPhpLaqtcnMbCxwgBpC2T2o7udh0m9WrT1mZmOFA9QQyupB9R7MVsd1gDIzG7GqBShJqyTtkbQ5r2yKpHWStqXnyXmvrZC0XdJWSe155ZdI2pReu0PKMg8knSDpwVS+QdKsvDqL0/fYJmnxSM6jrB7U/i0QR2GyA5RVR5Hr6kZJuyQ9lx4fzXutIa8rs1JUswd1D7BgQNlyYH1EzAbWp6+RdAGwCLgw1fm6pJZU505gKTA7PXLHXALsi4jzgNuB29KxpgA3AJcC84Ab8gNhucoKUN3PZ89OL7fquYfjryuA2yNiTno8Do19XZmVomoBKiJ+CHQNKL4aWJ22VwML88ofiIjDEbED2A7Mk3QWcGpEPBURAdw7oE7uWA8DV6X/AtuBdRHRFRH7gHUUvqBLUtYQ377nYfy74ORzhvvtzAZV5LoqpmGvK7NS1Poe1LSI2A2Qns9M5W3Aq3n77UxlbWl7YHm/OhHRC+wHTh/kWMeRtFRSh6SOvXv3FmxweT2on8FpF4F8a89q7guSnk9DgLmeTcNeV2alaJS/pIU+0RqDlA+3Tv/CiLsiYm5EzJ06dWrBhvX0QEsLTJhQ8OX8g2VDfL7/ZLV3J3AuMAfYDXw1lTfsdWVWiloHqNfS8ALpeU8q3wnMzNtvBtCZymcUKO9XR9J44DSyoY9ixxqWkpd779mVLa/h+09WYxHxWkT0RcRR4Btk94igga8rs1LUOkCtAXLZP4uBx/LKF6UMonPIbto+nYYBD0i6LI2DXzugTu5Y1wBPpvH0tcB8SZPTUMf8VDYsJS/3vi+XIOEelNVW7p++5ONALsOvYa8rs1JUbS0ISfcDVwBnSNpJlgF0K/CQpCXAK8AnASJii6SHgBeAXuD6iOhLh7qOLHOpFXgiPQDuBu6TtJ3sP7xF6Vhdkr4EPJP2uykiSr2pfJySl3t/J4PvouF+K7MhFbmurpA0h2zI7WXg89DY15VZKZT9c2Rz586Njo6O48o/+1l45hnYtm2IA/z4M/D6U3D1y1VpnzUGSc9GxNx6t2O0KHZdmeUrdl01SpJEw+rpKXGIz1McmZlVlAPUEEpKkug7BG9udYKEmVkFOUANoaQe1P4XIfp8/8nMrIIcoIZQUg/qza3Z82nvr3p7zMzGCgeoIZSUZn4gZVCcfG7V22NmNlY4QA2hpDTzA9vgpBkw/qSatMnMbCxwgBrCW2/ByScPsdOBbXDK7Jq0x8xsrHCAGkJJAeqtbXDK+2rSHjOzscIBahAR8OtfDxGgDnfB4TfcgzIzqzAHqEEcPgx9fUMEqFyChAOUmVlFOUAN4te/hokTHaDMzOqhapPFNoPTT896UYNOV3hgW7ZA4cnvrVm7zMzGAvegSqBCS7XlHNgGJ70HWk6oWXvMzMYCB6iROvCSh/fMzKrAAWokIvwZKDOzKnGAGonDe+HImw5QZmZV4AA1ErkMvlP9IV0zs0pzgBoJp5ibmVVNXQKUpJclbZL0nKSOVDZF0jpJ29Lz5Lz9V0jaLmmrpPa88kvScbZLukPK8u0knSDpwVS+QdKsqpzImy+BxsO7qnN4M7OxrJ49qN+OiDl569AvB9ZHxGxgffoaSRcAi4ALgQXA1yW1pDp3AkuB2emxIJUvAfZFxHnA7cBtVTmDA9vg5HNgnD9OZmZWaY00xHc1sDptrwYW5pU/EBGHI2IHsB2YJ+ks4NSIeCoiArh3QJ3csR4Grsr1ripq/xY47YKKH9bMzOoXoAL4J0nPSlqayqZFxG6A9HxmKm8DXs2ruzOVtaXtgeX96kREL7AfOL2iZ9B3OOtBnXZhRQ9rZmaZeo1NXR4RnZLOBNZJ+vkg+xbq+cQg5YPV6X/gLDguBTj77LMHb/FAB16C6IXTPlBePTMzK0ldelAR0Zme9wD/CMwDXkvDdqTnPWn3ncDMvOozgM5UPqNAeb86ksYDpwFdBdpxV0TMjYi5U6dOLe8kurdkz5McoMzMqqHmAUrSuySdktsG5gObgTXA4rTbYuCxtL0GWJQy884hS4Z4Og0DHpB0Wbq/dO2AOrljXQM8me5TVc7+zaAWL1RoNSVplaQ9kjbnldUkA1bS4vQ9tknKXV9mVXhCbT0AAAjmSURBVFOPHtQ04EeSfgY8DXw3Ir4H3Ap8RNI24CPpayJiC/AQ8ALwPeD6iOhLx7oO+DuyxIlfAE+k8ruB0yVtB/6ElBFYUfs3Z8HJk8Rabd3DsWzVnKpnwEqaAtwAXEo24nFDfiA0q4aa34OKiF8CHyxQ/gZwVZE6NwM3FyjvAI4bY4uIQ8AnR9zYwXRvgclzqvotzAaKiB8W+Fzf1cAVaXs18H3gT8nLgAV2pH/Y5kl6mZQBCyAplwH7RKpzYzrWw8DXUu+qHVgXEV2pzjqyoHZ/pc/RLKeR0sxHj96D8NYvfP/JGkUtMmCLHes4kpZK6pDUsXfv3hGclo11DlDD8ebPgXAGnzW6SmbAlpQZCyNMPjLL4wA1HN3p/rQ/A2WNoRYZsMWOZVY1DlDDsX8zjJsIp5xX75aYQW0yYNcC8yVNTskR81OZWdV4Ernh6N4Mp/6G5+CzmpN0P1lCxBmSdpJl1t0KPCRpCfAKKUEoIrZIymXA9nJ8Buw9QCtZckR+Bux9KaGiiywLkIjokvQl4Jm03025hAmzavFf2HId7YXX/w3Orm6SoFkhEfGZIi9VPQM2IlYBq0purNkIeYivXG9sgCP74ayBH0UxM7NKcoAqV+f3shkk3l3wH1YzM6sQB6hy7V4Lp18KEyfVuyVmZk3NAaoch16Hrg4P75mZ1YADVDn+fR0QcFb7kLuamdnIOIuvHLvXwsQpMOWSerfEbNR6dOMuVq7dSmd3D9MntbKs/XwWXlxw1iQb4xygShWRBah3fwTGtQy9v5kd59GNu1jxyCZ6jmQfx9rV3cOKRzYBOEjZcTzEV6ru5+HQv8N0338yG66Va7e+E5xyeo70sXLt1jq1yBqZA1SpdqdZXd49v77tMBvFOrt7yiq3sc0BqhS9B+Glv4HT58FJ0+vdGrNRa/qk1rLKbWxzgCrFC1+Gg6/AnC/XuyVmo9qy9vOZ0NJ/5Y4JLWJZ+/l1apE1Mgeoobz1Mrx4G5z9aZj2n+vdGrPRb+AqUgVXlTJzgBpcHIWO/w6Mgw99pd6tMRv1Vq7dypGj/SPSkaPhJAkrqKkDlKQFkrZK2i5pedkHeGtHNnP5nFvgpBlD729mg3KShJWjaT8HJakF+BvgI2SrgT4jaU1EvFDyQU45Fz72ArROq1IrzcaW6ZNa2VUgGDlJwgpp5h7UPGB7RPwyIt4GHgCuLvsoDk5mFbOs/XxaJ/T/oHvrhBYnSVhBzRyg2oBX877emcreIWmppA5JHXv37q1p48zGooUXt3HLJy6ibVIrAtomtXLLJy7yLBJWUNMO8QEqUNbv7mxE3AXcBTB37lznEpnVwMKL2xyQrCTN3IPaCczM+3oG0FmntpiZWZmaOUA9A8yWdI6kicAiYE2d22RmZiVq2iG+iOiV9AVgLdACrIqILXVulpmZlahpAxRARDwOPF7vdpiZWfmaeYjPzMxGMQcoMzNrSIpwdjWApL3Ar4q8fAbweg2b0yh83sd7T0RMrWVjRrO862o0v5fc9uoreF05QJVAUkdEzK13O2rN522VMpp/pm57/XiIz8zMGpIDlJmZNSQHqNLcVe8G1InP2yplNP9M3fY68T0oMzNrSO5BmZlZQ3KAMjOzhuQANYgRLxnf4CS9LGmTpOckdaSyKZLWSdqWnifn7b8i/Sy2SmqvX8vLJ2mVpD2SNueVlX2uki5JP7Ptku6QVGhZF0sa6Rqq1Pu92HtA0gmSHkzlGyTNGkFbq/p+Haytkhan77FN0uLhnkNFRIQfBR5kE8z+AngvMBH4GXBBvdtV4XN8GThjQNmXgeVpezlwW9q+IP0MTgDOST+blnqfQxnn+p+ADwGbR3KuwNPAfyBbb+wJ4HfqfW6N+mi0a6hS7/di7wHgD4C/TduLgAcb9f1arK3AFOCX6Xly2p5cr9+Ze1DFVWbJ+NHnamB12l4NLMwrfyAiDkfEDmA72c9oVIiIHwJdA4rLOldJZwGnRsRTkV3N9+bVseONhmuoku+B/GM9DFw13B52Dd6vxdraDqyLiK6I2AesAxYM5xwqwQGquCGXjG8CAfyTpGclLU1l0yJiN0B6PjOVN+PPo9xzbUvbA8utsEZ7z1Ti/T7Ye+CdOhHRC+wHTq9g+2vR1ob6nTX1chsjNOSS8U3g8ojolHQmsE7SzwfZdyz8PHKKnetY+hlUQqP9vCrxfh/snOp1vpVsa0P9ztyDKq7pl4yPiM70vAf4R7IhmdfS0ADpeU/avRl/HuWe6860PbDcCmuo90yF3u+DvQfeqSNpPHAaxw/TjUQt2tpQvzMHqOKaesl4Se+SdEpuG5gPbCY7x1zmzmLgsbS9BliUsn/OAWaT3YAdzco61zSsckDSZWm8/tq8Ona8hrmGKvV+H+I9kH+sa4An072fSqlFW9cC8yVNTlmC81NZfdQrO2M0PICPAi+RZcX8eb3bU+Fzey9Z5s/PgC258yMbh14PbEvPU/Lq/Hn6WWxllGWvAfcDu4EjZP8lLhnOuQJzyf6w/QL4Gmk2Fj+K/twb4hqq5Pu92HsAOBH4f2RJCk8D723U9+tgbQV+P5VvB36vnu8fT3VkZmYNyUN8ZmbWkBygzMysITlAmZlZQ3KAMjOzhuQAZWZmDckByo4jaZKkP6h3O8xsbHOAskImkc12bGZVJOl7krolfafebWlEDlBWyK3AuWndnJX1boxZE1sJfK7ejWhUDlBWyHLgFxExJyKW1bsxZqOJpN+S9LykE9MUS1skfaDQvhGxHjhQ4yaOGp7N3MysgiLiGUlrgP8NtAL/EBGbh6hmBThAmZlV3k1kk+UeAv6wzm0ZtTzEZ4UcAE6pdyPMRrEpwMlk19GJdW7LqOUAZceJiDeAH0va7CQJs2G5C/hfwDeB2+rcllHLQ3xWUER8tt5tMBuNJF0L9EbEtyS1AP8m6cqIeLLAvv8K/AZwsqSdwJKIqN/6Sw3Gy22YmVlD8hCfmZk1JA/xmZlVkaSLgPsGFB+OiEvr0Z7RxEN8ZmbWkDzEZ2ZmDckByszMGpIDlJmZNSQHKDMza0j/HyeCL6/O8vwdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# scipy.integrate.RK45\n",
    "# 用龙格库塔法解一阶微分方程组\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import RK45\n",
    "\n",
    "\n",
    "def funcs(t, x):\n",
    "    dxdt = np.zeros(2)\n",
    "    a = 1e-8\n",
    "    dxdt[0] = 0.05*x[0]*(1-x[0]/150000)-a*x[0]*x[1];\n",
    "    dxdt[1] = 0.08*x[1]*(1-x[1]/400000)-a*x[0]*x[1];\n",
    "    return dxdt\n",
    "\n",
    "\n",
    "t0 = 0\n",
    "y0 = [5000, 70000]\n",
    "t_bound = 1000\n",
    "\n",
    "rk45 = RK45(fun=funcs, t0=t0, y0=y0, t_bound=t_bound)\n",
    "\n",
    "t = []\n",
    "x = [[], []]\n",
    "\n",
    "while rk45.status == 'running':\n",
    "    rk45.step()\n",
    "    t.append(rk45.t)\n",
    "    x[0].append(rk45.y[0])\n",
    "    x[1].append(rk45.y[1])\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.plot(t, x[0], color='orange')\n",
    "plt.plot(t, x[1], color='b')\n",
    "plt.xlabel(u't')\n",
    "plt.ylabel(u'x')\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.scatter(x[0], x[1])\n",
    "plt.xlabel(u'x_1')\n",
    "plt.ylabel(u'x_2')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}