计算方法聚合版.ipynb · 徐晨哲/数值分析

加入 Gitee

与超过 1200万开发者一起发现、参与优秀开源项目，私有仓库也完全免费：）

该仓库未声明开源许可证文件（LICENSE），使用请关注具体项目描述及其代码上游依赖。

克隆/下载

计算方法聚合版.ipynb 110.31 KB

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# 最小二乘法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pylab as pl\n",
    "import math"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.4349474726647082\n",
      "2.0559131261889667\n"
     ]
    },
    {
     "data": {
      "text/plain": "[<matplotlib.lines.Line2D at 0x25d718b9c30>]"
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD8CAYAAABuHP8oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAArtUlEQVR4nO3deVzVVf7H8ddhBwWUTQTEXXNXwN1MU3MpbbHFyjZttJmmqabRLC2zckqtqX41NVlqjZnWtKrlbmqZKYL7vguCgCL7dpfz++OLhQrKcuFyL5/n48ED7pfv/d7PEXz79dyzKK01QgghHI+LvQsQQghRORLgQgjhoCTAhRDCQUmACyGEg5IAF0IIByUBLoQQDuqaAa6Umq+USlVK7S1xbI5S6qBSardS6lulVINqrVIIIcQVynMH/gkw7LJja4COWuvOwGHgORvXJYQQ4hquGeBa601A+mXHVmutzcUPfwMiqqE2IYQQV+Fmg2uMA74oz4lBQUG6WbNmNnhJIYSoO+Li4s5prYMvP16lAFdKTQXMwKKrnDMBmAAQGRnJ9u3bq/KSQghR5yilTpV2vNKjUJRSDwG3APfrqyyoorWeq7WO0VrHBAdf8Q+IEEKISqrUHbhSahjwLHCD1jrPtiUJIYQoj/IMI1wMbAHaKqUSlVLjgfcAX2CNUmqnUuo/1VynEEKIy1zzDlxrfW8ph+dVQy1CCCEqQGZiCiGEg5IAF0IIByUBLoQQDkoCXAghqlNOGqyYAgWZNr+0LWZiCiGEuJzFBNs+gg2vgykXWtwAbYfb9CUkwIUQwtaOb4AVz0LaQWh5IwybBcFtbP4yEuBCCGErF07B6mlwYCk0bAZjFht33UpVy8tJgAshRFWZ8mHzO/DLW4CCG6dB7yfA3ataX1YCXAghKktr42571TTIPA0d7oCbXgH/mllhWwJcCCEqI/WA0c99YiOEdICHf4Bm/Wq0BAlwIYSoiPwM2DgLtn4InvVh+ByIGQeuNR+nEuBCCFEeVivs/AzWzoC88xD9MNz4AtQLtFtJEuBCCHEtidvhx0mQFA9NesLYryGsq72rkgAXQogyZafAuhmwcxHUD4Xb50Lnu6ttWGBFSYALIcTlLCajj3vD62AugL5PQf9/gKevvSu7hAS4EEKUdGy9Mbrk3GFoNQSGvQ5BrexdVakkwIUQAiD9hDGL8uByaNgc7v0C2g6zd1VXJQEuhKjbivKMGZSb3wEXVxj0IvT+K7h52ruya5IAF0LUTVrD/u+MWZRZidDpLhg8A/zD7V1ZuUmACyHqnpR9Rj/3yZ+hUScY/RE07WPvqipMAlwIUXfkX4CfXoPYj8HLD25+E6IfMbpOHJAEuBDC+VktsGMhrHvZCPGYcTBwKvgE2LuyKpEAF0I4t9NbYcVkSN4JkX1g+Cxo3NneVdmEBLgQwjlln4U102H3EvANg9HzoOPoWjOL0hYkwIUQzsVcBFs/gI2zwVIE/f4O1z9jrBzoZCTAhRDO48haWPksnD8KbYbD0JkQ2NLeVVUbCXAhhOM7fwxWTYXDKyCgJdz/FbQeYu+qqp0EuBDCcRXlws9vwq/vgquHMRGn11/AzcPeldUICXAhhOPRGvZ+DatfgOwk6DwGBr8Efo3tXVmNumaAK6XmA7cAqVrrjsXHAoAvgGbASeBurfWF6itTCCGKnd0DP06G079CaGe4awFE9rJ3VXbhUo5zPgEuX5JrCrBOa90aWFf8WAghqk9eOvzwDHzYH9IOwi1vw4QNdTa8oRx34FrrTUqpZpcdvhUYUPz1p8AG4FlbFiaEEIAxizLuE1j/ChRkQvc/wcDnwLuhvSuzu8r2gTfSWicDaK2TlVIhZZ2olJoATACIjIys5MsJIeqkU1tgxSSj26RpPxgxGxp1sHdVtUa1v4mptZ4LzAWIiYnR1f16QggnkJUEa16EPf8Dv3C4cwF0uN2pZlHaQmUDPEUp1bj47rsxkGrLooQQdZS5ELb8Gza9AVYz9J8E/Z4Gj3r2rqxWqmyALwUeAl4v/vy9zSoSQtRNh1fByimQfhza3mzMogxobu+qarXyDCNcjPGGZZBSKhGYjhHcXyqlxgOngbuqs0ghhBM7f8wI7iOrIbA1jP0aWg22d1UOoTyjUO4t41uDbFyLEKIuKcyBTXOMLhM3L7jpVegxsc7MorQFmYkphKhZWhtvTq55EbKToev9MGg6+Dayd2UORwJcCFFzkncZsygTfoOwbnD3QmjS3d5VOSwJcCFE9cs9b0zEifsEfAJh1LvQdSy4lGcyuCiLBLgQovpYzBC3ANa/CoXZ0PMxGDAFvBvYuzKnIAEuhKgeJ3+BFc9Cyl5o3h+Gz4aQdvauyqlIgAshbCsz0Vjmdd834N8E7v4vtBslsyirgQS4EMI2TAWw5V34+V+grXDDFOj7JHj42LsypyUBLoSoGq3h0ApY9RxcOAntRsJNM6FhU3tX5vQkwIUQlXfuiNHPfWwdBLWFB76DlgPtXVWdIQEuhKi4gizYNBt++wDcfWDoa9DjT+Dqbu/K6hQJcCFE+VmtsPsLWDsdclKg21gY9BLUD67wpWKXfkiT+DmE6DRSVTAJUZPoPmqi7Wt2YhLgQojySdphzKJM3Abh0TBmMUREV+pSsUs/pGPcNLxVESgIJQ3/uGnEgoR4BUiACyGuLvccrHsZ4v8L9YLg1vehy71VmkXZJH6OEd4leKsimsTPgWoKcItVk1dkpsBkpcBkodB88bPl92MmiwY0Vm28N2vVGqs29qGxao2LUri7uuDmonB3c8HdxQU3V+OYu6vC290VH0836nm44uPhhodb9c40lQAXQpTOYobYj+Gnf4IpF3o/DjdMBi//Kl86RKdBKcPCQ/S5az5Xa01ukYVz2YWcyzE+0nKKOJddSEZeEVkFZrLyTWQXmMkqMJGVbyKrwExOobnKdVeUu6vCx8MI9Dfu7kKflkE2vb4EuBDiSic2Gd0laQegxUAYPguC29rs8qkqmFDSSjkehE+BiaSMfJIy8jmTUfD710kZ+SRnFnAup5ACk7XU6/p5ueHv446fl/ERGeCDn3fxY2836nm44eXugqe7K17urni5FX/t5oKXuyturgoXpVAKXJTCRQEYn5VSWLXGbNGYLNbiD43ZYqWo+OsCk4W8IjN5RRbyiizkFpp//xxYz9Nmf34XSYALIf6QkQCrp8H+76BBJNyzCK672aazKLXW7O34LMd2zSdJB3FSh3JSh3JcN+akSxPyX1p9yfnurorG/t409vcipmlDgn09CazvSVB9T4LqexBU35NgX08C6nng7lq3FseSABdCgCkfNv8f/PKW8XjgVOjzBLh7V/qSWmuSMws4lJLN4bPZHErJ5mhqDifScskuDAOmAeCKhcYqnWD/+tzdrjkRDX0Ia+BNWAMvwht4E1TfExcXmYZfGglwIeoyreHgclj1PGSchva3wU2vGHffFZBXZGZ/Uhb7krI4eDabw8WhnV2i3znUz4vWjepze1Q4zQLr0TyoHs2C6hHR0LvO3TnbigS4EHVV2iFjFuXxnyCkPTy0zFg18BouhvWeM5nGR2Imx9JysBqDNWjg406bRr7c1i2ctqG+tA31pU2IL/4+MsnH1iTAhahrCjJhwyzY9iF41DOWeY0ZD65XxoHWmmNpucSfusD2U+nsOJ1xSVgH+3rSKdyfEZ0a0yncn47h/jTy80TJyoM1QgJciLrCaoVdn8Pal4yx3VEPwqAXjbHdxQpMFnYnZrL9VDrxpy4Qd+oCF/JMgHFn3a1JA4Z3akzncH86RfjTyM/LTo0RIAEuRN2QGAcrJsGZOIjoAff/D8K6UWCyEH/sHFuOnefXY+fZnZhRPJkFWgTXY3C7RsQ0a0h004a0CKovbybWMhLgQjiznFRYNwN2fAb1G2G+9QN2NbyJXw9c4NflvxF3+gJFZiuuLopO4f6M69ecmKYBRDdtSEA9D3tXL65BAlwIZ2Qxwba5sOF1tCmfvU0f4n3rHWz8toC8oq0AtG/sx4O9mtKnVSDdmwXg6yVvMjoaCXAhnEzR4XWYlk+iXtYxtrl2Y0r+/Rw/FEazQM3oqAj6tAykZ4tAucN2AhLgQjiB5Mx8tsTFE7n9n8Tk/UKyNYRn9CQKIofwYNsQBrQNoVlQPXuXKWxMAlwIB6S15mhqDqv3p/DT3pNcn7KIia7L0MqFdWETcev3BG+3DsPL3dXepYpqJAEuhIOwWjU7EzNYte8sa/alcPxcDsNdtvG+1+eEuKWR3fpWfG/5J4P8I+xdqqghEuBC1GJWq2ZHwgWW7Urmxz3JpGYX4uaiuDMym0UeH9M4fRsEd4Thn+LbrK+9yxU1rEoBrpR6GngU0MAe4BGtdYEtChOirtJasy8pi2W7kli+O5kzGfl4uLnQzS+Hp9wXc4fLRjzPmrG4esOINyD6kVJnUQrnV+mfulIqHPgb0F5rna+U+hIYA3xio9qEqFOOpGSzbFcSy3Ync+JcLm4uiv5tgvnH0DYEHF9O793P4+Fq+f18k9nMjrNWukt411lV/cm7Ad5KKRPgAyRVvSQh6o5zOYV8t+MMX8ef4UByFi4KercMZGL/FgzrGEoDHw9IiMX0/VTcleWS53orU7VuQSZqv0oHuNb6jFLqDeA0kA+s1lqvvvw8pdQEYAJAZGTFlqgUwhkVmi2sP5DK1/GJ/HQoDYtV0yXCn5dGtmdE58aE+BavL5KdAt++BLs+x01T6S3IhPOqShdKQ+BWoDmQAfxPKTVWa/1ZyfO01nOBuQAxMTG68qUK4bi01uxOzOSruESW7U4iI89EIz9PHr2+OXdGRdC6ke8fJ5uLYOt/YONsMBdAv6dJ/eW/NOL8FddNVUGE1mA7RO1SlS6UwcAJrXUagFLqG6AP8NlVnyVEHZKZZ+KbHYks3naawyk5eLq5MLRDKKOjI+jXKgjXyxeHOroWVkyB80eg9VAY9hoEtuR0XmP84qZdspN7vvYgIXqSBHgdVpUAPw30Ukr5YHShDAK226QqIRyY1pq4Uxf4fNtpftidTKHZSpcmDXjtjk7c3LkxfqWtOZJ+AlZNhUM/QEALuO9LaDP09293HzWRWKBJ/BxC9DlSVRAJ0ZPoLv3fdZrSuvK9GkqpGcA9gBnYATyqtS4s6/yYmBi9fbtkvHBOmXkmvt2RyOJtCRxKyaa+pxu3dQvj3h6RdAjzL/1JRbnw87/g13fBxQ1umAS9/gJutt/BXDgupVSc1jrm8uNVGoWitZ4OTK/KNYRwdHsSM/l0y0mW7Uoy7rYj/Hn9jk6M7BJGPc8y/oppDfu+gdUvQNYZ6HQ3DJkBfmE1W7xwaDKAVIhKMFmsrNx7lk9+PUncqQv4eLgyOjqC+3pE0jG8jLvti87uNfaiPPULhHaC0fOgae+aKVw4FQlwISrgfE4hi7ed5rPfTnM2q4CmgT68eEt77oyJKL1vu6S8dNjwGsR+DF7+cMtbEPUQuMiCU6JyJMCFKIe9ZzL59NeTfL8riSKzletbBzHz9o4MaBty5UiSy1ktEP8prHsFCjKMDYQHPg8+ATVSu3BeEuBClEFrzYbDaczdeJwtx8/j4+HKPTFNeKhPU1qF+F77AgCnf4MfJ8HZ3dC0LwyfZXSbCGEDEuBCXKbIbGXZriQ++vk4B89mE+rnxfMjruOe7pH4e5dz27GsZFg7HXZ/Ab5hRj93x9GgZFNgYTsS4EIUyy4wsWRbAvM3nyA5s4C2jXx5864ujOwShoebS/kuYi6E396HjXPAaoLr/wHX/x08ZDccYXsS4KLOS80qYP7mkyzaeorsAjO9WgTwzzs6MaBNMKoid8xH1hijS9KPQdsRMHSmMSlHiGoiAS7qrDMZ+fxnwzG+iE3AbLUyvGNjJvRvQZcmDSp2ofPHYNXzcHglBLaC+7+G1oOrpWYhSpIAF3XOqfO5fLDhGF/HJwJwZ3QEj93QkqaBFezmKMyBn9+ELe+BqwcMeRl6/hncZLd3UTMkwEWdcSwth3+vP8r3u5JwdVHc2yOSiTe0JLyBd8UupDXs/dqYRZmdBF3uhcEvga8sKyVqlgS4cHqHzmbz7voj/LAnGU83F4Y1yuKx9Nl0iD9G6o5gYqMqsChU8m5YMRlOb4HGXeCuTyCyZ7XWL0RZJMCF0zqamsNbaw/zw+5k6nm4MrF/S2Jyf6bvnud/X5Y1lDT846YRC1cP8bx0WP8qxC0A74Yw8h3o9oDMohR2JQEunE5Ceh5vrz3CtzsS8XZ35a8DWzG+X3Ma1vPg7Eu3XLKmNoC3Kip7azKrxQjt9a9CQRb0mAADphghLoSdSYALp3E2s4B31x/hi9gEXF0U4/s157EbWhJY/4+lWUN0Wvm3Jju52RgWmLIHml1vzKJs1KEaWyBExUiAC4d3PqeQDzYcY+Fvp7BqzZgeTfjrwNaE+ntdcW6qCiaUtFKOl9iaLPMMrHkR9n4FfhFGP3f722QWpah1JMCFw8ouMDF303Hm/XKCApOFO6IieHJQa5oE+JT5nISoSfiXtTWZudAYErjpTbCa4YZnoe9T4FH29YSwJwlw4XBMFiuLt53mnbVHOJ9bxM2dG/P04Da0Cql/zeeWujVZ1D/o3rYp/LsnXDgB191izKJs2Kza2yJEVVRpS7WKki3VRFVorVm17yyzVh7ixLlcerUI4PkR7egc0aDyFz13FFZOgaNrIKiN0c/d8kab1SyELVTLlmpC1JTtJ9P5548HiD+dQeuQ+sx/OIaBbUMqtlZJSYXZsGkObHkf3LzgppnQcyK4lnO1QSFqAQlwUasdS8th9sqDrNqXQoivJ7NGd2J0VARuruVcHfByWsPuL403KXPOQtexMOhF8G1k28KFqAES4KJWysgr4u21R1j42ym83Fx4Zkgbxl/fHB+PKvzKJu00ZlEmbIWwKBizCCKu+F+pEA5DAlzUKubiNyjfXHOYrHwT9/aI5KnBbQj29bz2k8uSex7Wvwxxn4JPIIx6D7reDy6VvIsXopaQABe1xuaj53h52X4OpWTTu0UgL45sT7vGfpW/oMUM2+fDT68aKwf2+gvcMBm8G9isZiHsSQJc2N3p83nM/HE/q/alENHQm/+MjWJoh9DKv0EJcOJnYxZl6j5ofgMMnw0h19muaCFqAQlwYTc5hWb+/dNR5v18AjdXxaShbRnfrzle7lVYICozEVZPg33fgn8k3L0Q2o2UWZTCKUmAixqntWbZ7mReXb6f1OxC7ogK59lh19HI78qp7+VmKoBf3zU2WEDDgOeh79/AvYJrfQvhQCTARY06mprDi9/v5ddj5+kU7s+HD0TTLbIKK/tpDYd+hJXPQcYpaDfKmEXZINJ2RQtRS0mAixqRX2Th3fVH+Ojn43i7u/LKbR25r0ckri5V6NpIOwwrn4Vj6yH4Onjwe2gxwGY1C1HbVSnAlVINgI+BjoAGxmmtt9igLuEktNas2Z/CjGX7OZORz+ioCJ4bcR1B9aswLLAgCzbOgq3/Afd6MOx16P6ozKIUdU5V78DfAVZqre9USnkAsmybACB26Ye4bp/He+aRrLdG0cTHxJcT+9OjeUDlL2q1wu4lsGY65KZBt7EwaDrUD7Zd4UI4kEoHuFLKD+gPPAygtS4Ciq72HFE3bPnuQ7bG/sYHludww8JUt8+4x/wTh/bMgObl3HvycmfijVmUibEQHgP3LYHwaNsWLoSDqcodeAsgDViglOoCxAFPaq1zbVKZcEg7Tl9g2lYXjuk7Ge6ylRfdF9JYpQOUvW3Z1eSkwboZsOMzqBcMt30AncfILEohqFqAuwFRwBNa661KqXeAKcALJU9SSk0AJgBERsrIAGeVW2jmjdWH+OTXkzTSnnzk/gZDXOMvOafUbcvKYjFB7Mfw02tgyoXejxsbLHhVYWamEE6mKgGeCCRqrbcWP/4KI8AvobWeC8wFYz3wKryeqKV+OpjKtO/2kpSZzwO9mvJQ/GRaqsQrzrtk27KrOb7RmEWZdsBYm3vYLAhuY/O6hXB0lQ5wrfVZpVSCUqqt1voQMAjYb7vSRG13LqeQGcv2s2xXEq1D6vPVY72JbhpArPobYWVtW3a1C2achlVT4cBSaNAUxnwObUfILEohylDVUShPAIuKR6AcBx6pekmittNa81VcIjN/PEBeoYWnB7fhsQEt8HQzpsCXum1Z9CS6l9X/bcqHze/AL28BCgZOgz5PgHsVZmYKUQfIlmqiQpIy8pnyzR42HU6je7OGvHZHJ1qF+FbuYlrDgWXGXXfmaehwO9z0KvhH2LZoIRycbKkmqkRrzZfbE3h1+QEsWvPyrR0Y27MpLpWdSZl60BgWeGIjhHSAh5ZD8+ttW7QQTk4CXFxTcmY+U77ew8bDafRqEcDs0V2IDKzknK2CTNjwOmz9EDzrw/A5EDMOXOVXUYiKkr81okxaa/4Xl8gry/djtmhmjOrAA70qeddttcLORcaY7txzEP0Q3Pgi1Au0feFC1BES4KJUZzMLeO6b3fx0KI0ezQOYc2dnmgbWq9zFErfDj5MgKR6a9IT7v4KwrjatV4i6SAJcXOLiCJOXl+/HZLEyfWR7HurdrHJ33TmpsPYl4867fijcPhc63y3DAoWwEQlw8bvzOYU8980eVu9PoUezAGbf2ZlmQZW467aYjD7ujbOMIYJ9n4T+k8CzkqNVhBClkgAXgDGbctJXu8nKNzF1RDvG92teubvuY+uNWZTnDkOrIcZSr0GtbF+wEEICvK7LKzIz84cDLNp6mutCfVk4vkfldoK/cNIYz31wOTRsDvd+AW2GSneJENVIArwO25WQwdNf7OTE+Vwm9G/B34e0qfiGwkV5xgzKze+AiysMehF6/xXcqrBhgxCiXCTA6yCzxcr7G47xzrojNPL15PNHe9G7ZQWH82kN+7+DVdMgKxE63glDXgb/8GqpWQhxJQlwJxO79MPiNUjSSFXBJERdugbJqfO5PPXFTnaczuC2rmHMuLUj/t4V3IosZb8xi/Lkz9CoE4z+CJr2sXFLhBDXIgHuRGKXfkjHi6sAKgglDf+4acRiLDD1dVwiL3y/FzcXxbv3dmNkl7CKvUD+BWN97tiPjXW5b34Toh8xuk6EEDVOAtyJNImfc8kSrgDeqoiAuP/j6fyefLvjDD2bB/DWPV0Ja+Bd/gtbLbBjIax72Qjx6EfgxmngU4X9LYUQVSYB7kRCdBpcNuhjt7U5TxQ9QcLOM/x9SBseH9gK14oMD0zYZsyiTN4JkX1g+Cxo3NmmdQshKkcC3ImkqmBCSQPAqhXzLMOZbR5DANl8MbE33ZtV4I45+6yx+/vuJeDbGEbPg46jZVigELWIBLgTSYiahH/cNHLx4hnTY2y0dmWwy3bui25c/vA2F8HWD2DjbLAUQb+/w/XPGCsHCiFqFQlwJ9J91ETmZXjy/n4PsvFmktuXdO/ehx63TijfBY6shZXPwvmj0GYYDP0nBLas3qKFEJUmAe4kzBYrb645zH8OBNMypD6L7uvGdaG3l+/J6ceNWZSHfoSAlnDf/6DNTdVbsBCiyiTAnUBKVgFPfL6DbSfTGdO9CdNHdsDboxxD+4py4ec34dd3wdUDBs+AXn8BN4/qL1oIUWUS4A5u89FzPLlkB7mFFt66pwu3dyvHfpJaw96vYfULkJ0Ene8xwtuvcfUXLISwGQlwB2W1at776ShvrT1My+D6LP5TFK0blWO51rN7jNUCT22G0M5w1wKI7FX9BQshbE4C3AGdzynk6S93selwGrd1DWPm7Z2o53mNH2VeOvw0E7bPB68GcMvbEPWgzKIUwoFJgDuYuFPpPL5oB+m5Rcy8vSP39YhEXW1sttUCcZ/A+leMDYW7PwoDnpNZlEI4AQlwB6G1Zt4vJ3h9xUHCGnjzzV/60DHc/+pPOrUFVkwyuk2a9jNmUYZ2rJmChRDVTgLcAeQUmvnHl7tYue8sN7VvxJy7ulx9BcGsJGMW5Z4vwS8c7lwAHW6XWZRCOBkJ8FruaGoOExdu5+T5PJ4fcR1/ur5F2V0m5kLY8m/Y9AZYzcY+lP2eBo9K7iYvhKjVJMBrsVX7zvLMl7vwcHNh4fge9GkZVPbJh1fByinGpJy2N8PQmRDQvOaKFULUOAnwWshi1by15jDv/XSUzhH+fDA2mvCyln89fwxWPgdHVkFgaxj7NbQaXLMFCyHsQgK8lsnMM/G3JTvYeDiNu6IjeOW2jqXvU1mYA5vmGF0mbl5w06vQY6LMohSiDqlygCulXIHtwBmt9S1VL6nuOpCcxcSFcSRn5vPqbR25v2cpQwS1hj1fwZoXIDsZutwHg6eDb6h9ihZC2I0t7sCfBA4Afja4llO61j6VAEt3JfHsV7vx9XJjyYTeRDdteOWFknfBj5Mh4Tdo3BXu/i806VEzjRBC1DpVCnClVARwMzAT+LtNKnIy19qn0myxMmvlQT76+QQxTRvy/v1RhPh5XXqRvHRjIs72BeATCKPeha5jwcXFLm0SQtQOVb0DfxuYDJS5CIdSagIwASAyMrKKL+d4ytqnskn8HDKHjOOJxTvYdDiNB3s3ZdrN7fFwKxHKFjPELYD1r0JhNvR8DAZMAe8GNdsIIUStVOkAV0rdAqRqreOUUgPKOk9rPReYCxATE6Mr+3qOqrR9KgFyLG7c9/5mTp/P47U7OnFvj8v+cTv5i7HoVMpeaN4fhs+GkHY1U7QQwiFU5Q68LzBKKTUC8AL8lFKfaa3H2qY051Byn8qLNlk68bjpSdzzTCx6tCc9WwT+8c3MRFjzorHcq38To5+73SiZRSmEuEKlO1G11s9prSO01s2AMcB6Ce8rJURNIl8bQ/u0hvnmYTxsepaGPm58/3jfP8LbVGAMC3yvOxz8AW6YAo9vg/a3SngLIUol48CrWfdRE4kFQuP+xbumkXxpHUj3gHwWPHkb9T3djFQ/tAJWPQcXTkK7kXDTTGjY1N6lCyFqOZsEuNZ6A7DBFtdyRs1vfJg/J3Ui9uQF/jqwFX8f0gYXFwXnjhjT34+uhaC28MB30HKgvcsVQjgIuQOvZvuTsvjTf7dzLqeQd+/txsguYVCQBZtmw28fgLsPDH0NevwJXK+ywqAQQlxGArward2fwt+W7MDPy52vHutDpzBf2LXEeJMyJwW6jYVB06F+iL1LFUI4IAnwarJg8wleWb6fDmH+zHsohpCcAzB/MiRug/BoGLMYIqLtXaYQwoFJgNuY2WLlleX7+XTLKYZ2aMRbt0Tgs2EyxP8X6gXBrf821i+RWZRCiCqSALehnEIzf/08ng2H0ph4fSTPBm7G5cM7oCgXej8ON0wGr2tsgyaEEOUkAW4jSRn5jPskliOpOXzcP4/BJydA7H5oMcCYRRnc1t4lCiGcjAS4DexJzGT8p7H4FaWwpdUyQratgAaRcM9ncN0tMhFHCFEtJMCraPW+s0xeso3HPX9kvPu3uCQpGDgV+jwB7mXsoiOEEDYgAX6Z8qzdDaC1Zt7Px4ldtZBVnp/TyJwC7W+Dm14x7r6FEKKaSYCXcK21uy+yWDUf/O9HOu99jUfd92ANuA5GfAQtbrBf8UKIOkcCvISrrd1NcYAXZF/gl3mTmHjhG8wePliHzMKl+6PgKn+UQoiaJalTQllrd4foc2C1khe7kKJV07nRksHh8Nu47v43jLHdQghhBxLgJZS2djdAuvLHb+4gfM7Gc8jamuwb59J/wE12qFAIIf4g0wFLKLl290Vm7UIQGWSdPc7zPE7BgyskvIUQtYLcgZdwce3uJvGzjW4TpVBKMV+PYpHn3bw/bgBtQ8vc/lMIIWqUBPhlundoA6cD4dw5UkL68sCZ21FBbfhsXHca+8u4biFE7SEBftGFU7B6KhxYhm7YjJWd3uLPsSH0ahHIhw/E4O8ta3ULIWoXCfCiPNj8Dmx+G5QL1oEv8FrGID7aksTILmG8cVdnPN1c7V2lEEJcoe4GuNZwYCmsmgqZCdDhDkyDZjB5TTrf7jjDuL7NmXZzO2PrMyGEqIXqZoCnHoAVk+HEJgjpAA//QH5Ybx7/PJ71B1OZNLQtfxnQEiWLUAkharG6FeD5GbDhddg2Fzx9YcQbEP0ImUWaR+dvZfupC8y8vSP395Qd4YUQtV/dCHCrFXZ+BmtnQN55iH4YbnwB6gWSmlXAg/O3cSwth/fujeLmzo3tXa0QQpSL8wd4QiysmARJO6BJL3jgG2jcBYBT53N5YN42zuUUsuDhHvRrLdPihRCOw3kDPDsF1r4Euz6H+qFwx0fQ6a7fN1fYn5TFg/O3YbFa+fxPvejapIFdyxVCiIpyvgA3F8G2D2HDLDAXQN+noP8/jD7vYrEn0xn3SSz1Pd1YMqE3rUJkdqUQwvE4V4AfXQcrnoXzR6D1TTDsdQhseckp6w+m8OfP4glv6M3C8T0JbyCzK4UQjsk5Ajz9hDGe+9APENAC7vsS2gy94rTlu5N4aslO2of5seDh7gTW97RDsUIIYRuOHeBFufDLW7D5/8DFDQZNh96Pg9uVwfzl9gSmfL2bmKYBzHs4Bl8vmRovhHBsjhngWsO+b2H1C5CVaLw5OeRl8Asr9fRPfz3J9KX7uL51EHMfiMHbQ6bGCyEcX6UDXCnVBPgvEApYgbla63dsVViZUvYZ/dwnf4bQTjD6Y2ja+/dvX74p8f+FvMznp/0Z0r4R793XTdY1EUI4jarcgZuBZ7TW8UopXyBOKbVGa73fRrVdKi8dNrwGsR+Dlz/c/C9jQo7LH4FcclNiDXxmGsDnp/3pF5TH+/dH4e4q+1cIIZxHpQNca50MJBd/na2UOgCEA7YP8J2LYdXzUJABMeNg4FTwCbjitIubEmsNr5jHMt8ygjGu63kieynurnfZvCwhhLAnm/SBK6WaAd2AraV8bwIwASAyMrJyL5CTAiHtYPgso9ukDCE6DQuKqebxLLHcyCOuK3jRbSG6tJ2KhRDCwVU5wJVS9YGvgae01lmXf19rPReYCxATE6Mr9SJ9noC+T/4+i7IsSTRitulOllr78lfXb3nG7X8oBSkEEVqpFxZCiNqrSgGulHLHCO9FWutvbFNSKVyu/cZjodnC0/X+SWy6N5PdFvMXt2UA5GsPEqInSYALIZxOVUahKGAecEBr/S/blVRxBSYLf/4sjth0b8Y1v8Adyb9h1YpUFURC9CS6j5poz/KEEKJaVOUOvC/wALBHKbWz+NjzWusfq1xVBRSYLExcGMfGw2kl1vIeCxjjG+XOWwjhrKoyCuUXsO+7g/lFFiYs3M4vR88xa3Qn7uleyTdJhRDCATnmTEyM8B7/aSxbjp9n9ujO3BXTxN4lCSFEjXLIAM8rMjPuk1i2nUjnzbu6cEdUhL1LEkKIGudwAZ5TaGbcgli2n0rnrXu6cmvXcHuXJIQQduFQAZ5dYOKRBbHsSMjgnTHdGNml9MWrhBCiLnCYAM8qMPHw/G3sTszkvXu7MbyTbD4shKjbHCLAM/NNPDh/G/vOZPLefVEM6yiDA4UQwiECfPr3e9mflMkHY6MZ0r6RvcsRQohawSEC/PkR7RgdHcH1rYPtXYoQQtQaDhHgIX5ehPh52bsMIYSoVWSHAyGEcFAOcQd+0eXbpSVEyUJVQoi6y2ECvOR2aSgIJQ3/uGnEgoS4EKJOcpgulIvbpZXkrYpoEj/HThUJIYR9OUyAh+i0Mo6fq+FKhBCidnCYAE9VpQ8hTFVBNVyJEELUDg4T4AlRk8jXHpccy9ceJERNslNFQghhXw4T4N1HTWRv9KucJRirVpwlmL3Rr8obmEKIOktpXbmN4isjJiZGb9++vcZeTwghnIFSKk5rHXP5cYe5AxdCCHEpCXAhhHBQEuBCCOGgJMCFEMJBSYALIYSDqtFRKEqpNOBUJZ8eBNS1aZfS5rpB2lw3VKXNTbXWV8xmrNEArwql1PbShtE4M2lz3SBtrhuqo83ShSKEEA5KAlwIIRyUIwX4XHsXYAfS5rpB2lw32LzNDtMHLoQQ4lKOdAcuhBCiBIcIcKXUMKXUIaXUUaXUFHvXY2tKqSZKqZ+UUgeUUvuUUk8WHw9QSq1RSh0p/tzQ3rXamlLKVSm1Qym1vPixU7dZKdVAKfWVUupg8c+7dx1o89PFv9d7lVKLlVJeztZmpdR8pVSqUmpviWNltlEp9Vxxnh1SSg2t7OvW+gBXSrkC/waGA+2Be5VS7e1blc2ZgWe01u2AXsDjxW2cAqzTWrcG1hU/djZPAgdKPHb2Nr8DrNRaXwd0wWi707ZZKRUO/A2I0Vp3BFyBMThfmz8Bhl12rNQ2Fv/dHgN0KH7O+8U5V2G1PsCBHsBRrfVxrXURsAS41c412ZTWOllrHV/8dTbGX+pwjHZ+Wnzap8BtdimwmiilIoCbgY9LHHbaNiul/ID+wDwArXWR1joDJ25zMTfAWynlBvgASThZm7XWm4D0yw6X1cZbgSVa60Kt9QngKEbOVZgjBHg4kFDicWLxMaeklGoGdAO2Ao201slghDwQYsfSqsPbwGTAWuKYM7e5BZAGLCjuNvpYKVUPJ26z1voM8AZwGkgGMrXWq3HiNpdQVhttlmmOEOCqlGNOOXRGKVUf+Bp4SmudZe96qpNS6hYgVWsdZ+9aapAbEAV8oLXuBuTi+F0HV1Xc73sr0BwIA+oppcbatyq7s1mmOUKAJwJNSjyOwPgvmFNRSrljhPcirfU3xYdTlFKNi7/fGEi1V33VoC8wSil1EqNb7Eal1Gc4d5sTgUSt9dbix19hBLozt3kwcEJrnaa1NgHfAH1w7jZfVFYbbZZpjhDgsUBrpVRzpZQHRuf/UjvXZFNKKYXRL3pAa/2vEt9aCjxU/PVDwPc1XVt10Vo/p7WO0Fo3w/iZrtdaj8W523wWSFBKtS0+NAjYjxO3GaPrpJdSyqf493wQxns8ztzmi8pq41JgjFLKUynVHGgNbKvUK2ita/0HMAI4DBwDptq7nmpoXz+M/0LtBnYWf4wAAjHevT5S/DnA3rVWU/sHAMuLv3bqNgNdge3FP+vvgIZ1oM0zgIPAXmAh4OlsbQYWY/TxmzDusMdfrY3A1OI8OwQMr+zrykxMIYRwUI7QhSKEEKIUEuBCCOGgJMCFEMJBSYALIYSDkgAXQggHJQEuhBAOSgJcCCEclAS4EEI4qP8H30xw4gZRpQ0AAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pylab as pl\n",
    "import math\n",
    "\n",
    "\n",
    "def phi_k(x, k):  # x^k\n",
    "    y = x ** k\n",
    "    return y\n",
    "\n",
    "\n",
    "def squar_least(x, y, w, n):  #w是权重，n是阶数\n",
    "    G = np.zeros((n + 1, n + 1))  #创建矩阵存放 n sum: xi xi^2\n",
    "    d = np.zeros((n + 1))  #创建行数sum:yi xiyi\n",
    "    for i in range(n + 1):\n",
    "        for j in range(n + 1):\n",
    "            G[i, j] = sum((w * phi_k(x, i)) * phi_k(x, j))  #算n  x的和 x^2的和\n",
    "\n",
    "    for i in range(n + 1):\n",
    "        d[i] = sum((w * phi_k(x, i)) * y)  #算y的和  xy的和\n",
    "\n",
    "    a = np.dot(np.linalg.inv(G), d)  #解方程\n",
    "    return a\n",
    "\n",
    "\n",
    "x = np.array([1, 4, 9, 16, 25, 36, 49, 64, 81])\n",
    "y = np.array([1, 2, 3.1, 3.95, 5.01, 6, 6.97, 8, 9.02])\n",
    "w = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1])\n",
    "\n",
    "a = squar_least(x, y, w, 2)  #2指的是2阶\n",
    "yuce_y = a[0] * 1 + a[1] * 0.6 + a[2] * 0.6 * 0.6\n",
    "print(yuce_y)\n",
    "\n",
    "xxx = [];\n",
    "yyy = []\n",
    "for i in range(0, 1000):\n",
    "    xxx.append(i / 10)\n",
    "    yyy.append(a[0] * 1 + a[1] * i / 10 + a[2] * i / 10 * i / 10)\n",
    "pl.scatter(x, y)\n",
    "pl.plot(xxx, yyy)\n",
    "\n",
    "a = squar_least(x, y, w, 1)  #1指的是1阶\n",
    "yuce_y = a[0] * 1 + a[1] * 0.6\n",
    "print(yuce_y)\n",
    "\n",
    "xxx = [];\n",
    "yyy = []\n",
    "for i in range(0, 1000):\n",
    "    xxx.append(i / 10)\n",
    "    yyy.append(a[0] * 1 + a[1] * i / 10)\n",
    "pl.scatter(x, y)\n",
    "pl.plot(xxx, yyy)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 拉格朗日"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x=2.5 f(x)= 1.5296665994020477\n",
      "x=40 f(x)= 6.199927388682746\n"
     ]
    },
    {
     "data": {
      "text/plain": "[<matplotlib.lines.Line2D at 0x25d7190f1c0>]"
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbsElEQVR4nO3dfWxdd53n8ffX9vWzHduxnQenaZI2TZuWbQNpp0yZ0iVAWUC0zAht0TKbRay6u2JnGDQLame1QjOziErMjkAaZqQuhYnEk6pS0Q7DQkNKmSnDtDj0OWmaNs+OY18/P17fp+/+cY9dJ7HbxPfh3Hv8eUnVPffch/P1afLxL9/zO+eYuyMiItFSFXYBIiJSeAp3EZEIUriLiESQwl1EJIIU7iIiEVQTdgEAnZ2dvmXLlrDLEBGpKAcPHhxy966lXiuLcN+yZQu9vb1hlyEiUlHM7ORyr6ktIyISQQp3EZEIUriLiESQwl1EJIIU7iIiEaRwFxGJIIW7iEgEKdxFRELyraeP848v9hfluxXuIiIh+fa/HOfnhweK8t0KdxGRkIxOp2hvrC3KdyvcRURCMJfOMDWXpqMpVpTvV7iLiIRgbCYFQHuTRu4iIpExOpMEoENtGRGR6BiZzoV7m8JdRCQ6RqdzbZkOtWVERKJjJGjLtOuAqohIdIwGbRlNhRQRiZCR6SQt9TXEqosTwwp3EZEQjM4ki9ZvB4W7iEgoRqaTRWvJgMJdRCQUozNJ2huLczAVFO4iIqEYmUoW7exUULiLiJScuzM8naSzua5o21C4i4iU2HQyw1w6y1qN3EVEomNkKriujMJdRCQ6hqbnANSWERGJkuFg5L62WSN3EZHIGAlG7mrLiIhEyND8yL1JbRkRkcgYmU7SVFtNQ2110bbxtuFuZt8ys0Eze3nRug4z229mR4PH9kWv3W9mr5vZETO7s1iFi4hUquGpOTqK2G+HSxu5/z3woQvW3QcccPftwIHgOWa2E7gHuD74zN+aWfF+NYmIVKDh6WRRWzJwCeHu7v8EjFyw+i5gX7C8D7h70fofuPucux8HXgduKUypIiLRMDyVpLMMRu5LWefu/QDBY3ewvgc4veh9Z4J1FzGze82s18x64/H4CssQEak8w9NzRZ0pA4U/oGpLrPOl3ujuD7r7bnff3dXVVeAyRETKk7szMp1kbRFPYIKVh/uAmW0ACB4Hg/VngCsWvW8TcHbl5YmIRMtEIk0q40W9rgysPNwfB/YGy3uBxxatv8fM6sxsK7AdeDa/EkVEomN4KncCUzHPTgWoebs3mNn3gTuATjM7A3wJeAB42Mw+A5wCPgHg7q+Y2cPAISANfNbdM0WqXUSk4oxMF/8EJriEcHf3Ty7z0p5l3v9l4Mv5FCUiElVDJbiuDOgMVRGRkhoOrisT+jx3EREpnFJcyx0U7iIiJTU8naS1vobamuLGr8JdRKSEhqbmij7HHRTuIiIlNTKdLPocd1C4i4iU1PBUsugzZUDhLiJSUrnryqgtIyISGelMluHpJN0tCncRkcgYmU7iDl0KdxGR6BiczJ3ApHAXEYmQuMJdRCR65sNdPXcRkQiJB5f77dRJTCIi0RGfnKO1vob6WHXRt6VwFxEpkcHJREn67aBwFxEpmfjknMJdRCRqcuFeX5JtKdxFREokPjlXkpkyoHAXESmJ6bk008mM2jIiIlGycAJTCaZBgsJdRKQk5ue4a+QuIhIhpbz0ACjcRURKopSXHgCFu4hIScQn56iuMtobi38XJlC4i4iUxOBkgs7mWqqqrCTbyyvczezzZvaKmb1sZt83s3oz6zCz/WZ2NHhsL1SxIiKVqpRnp0Ie4W5mPcAfA7vd/QagGrgHuA844O7bgQPBcxGRVS0+NVeyaZCQf1umBmgwsxqgETgL3AXsC17fB9yd5zZERCpexYzc3b0P+CvgFNAPjLv7E8A6d+8P3tMPdBeiUBGRSpXOZIlPzrGutTTXlYH82jLt5EbpW4GNQJOZfeoyPn+vmfWaWW88Hl9pGSIiZW9oKknWqYxwB94PHHf3uLungEeB3wUGzGwDQPA4uNSH3f1Bd9/t7ru7urryKENEpLydm0gAsL5Cwv0UcKuZNZqZAXuAw8DjwN7gPXuBx/IrUUSksp0bD8J9TenCvWalH3T3Z8zsEeC3QBp4DngQaAYeNrPPkPsF8IlCFCoiUqnOjc8CFRLuAO7+JeBLF6yeIzeKFxER4NzEHLFqo6NEZ6eCzlAVESm6gYkE3S31JTs7FRTuIiJFd248UdKWDCjcRUSKbmAiUdKZMqBwFxEpKnenfzxR0jnuoHAXESmqiUSa2VSGDWrLiIhEx0BwAtM6hbuISHQsnMCktoyISHSEcekBULiLiBTVQDBy724t3eV+QeEuIlJU/RMJ2htj1MeqS7pdhbuISBENhDANEhTuIiJFdW4iUfJpkKBwFxEpqoGJ0l96ABTuIiJFk0hlGJpKsr61oeTbVriLiBRJfzBTpqdd4S4iEhlnx3I36ehpU7iLiERG36jCXUQkcvrGZjEr7e315incRUSKpG9slu6WOmprSh+1CncRkSI5OzYbSksGFO4iIkXTNzbLRoW7iEh0ZLNO/1gilGmQoHAXESmKoak5kpksmzRyFxGJjjPBHHe1ZUREImThBCa1ZUREomP+BKaKHLmbWZuZPWJmr5rZYTN7t5l1mNl+MzsaPLYXqlgRkUpxdmyWlvoaWutjoWw/35H714Gfuvu1wI3AYeA+4IC7bwcOBM9FRFaVvhDnuEMe4W5mrcDtwEMA7p509zHgLmBf8LZ9wN35lSgiUnn6xhKVGe7ANiAOfNvMnjOzb5pZE7DO3fsBgsfupT5sZveaWa+Z9cbj8TzKEBEpP32jM6EdTIX8wr0GeCfwd+6+C5jmMlow7v6gu+92991dXV15lCEiUl7GZ1NMJNIVO3I/A5xx92eC54+QC/sBM9sAEDwO5leiiEhlOT0yA8CVaxtDq2HF4e7u54DTZrYjWLUHOAQ8DuwN1u0FHsurQhGRCnNyOBfumzuaQquhJs/P/xHwXTOrBY4Bnyb3C+NhM/sMcAr4RJ7bEBGpKCdHpgHYHOLIPa9wd/fngd1LvLQnn+8VEalkp0dmWNtUS3NdvuPnldMZqiIiBXZyeCbUUTso3EVECu7UyAybOxTuIiKRkUxnOTs2y5UKdxGR6OgbmyXrsHlteDNlQOEuIlJQp0bmp0Fq5C4iEhmnhnPTIMM8gQkU7iIiBXVqZIa6miq6mutCrUPhLiJSQCeHczNlqqos1DoU7iIiBXRqZCb0lgwo3EVECsbdOTUywxUhH0wFhbuISMHEJ+eYSWbYEvI0SFC4i4gUzBvx3EyZbV0KdxGRyDg2NAXAtq7mkCtRuIuIFMyx+DT1sSo2tNaHXYrCXUSkUI7Fp9ja2Rz6NEhQuIuIFMyxoemy6LeDwl1EpCCS6SynR2bY1qlwFxGJjFMj02S9PGbKgMJdRKQgFqZBdoY/UwYU7iIiBXGsjOa4g8JdRKQgjsWn6Gqpo6U+FnYpgMJdRKQg3ohPlc3BVFC4i4jkzd05OjDFjvUtYZeyQOEuIpKn/vEEk3Nptq9TuIuIRMZrA5MAXNNdHjNlQOEuIpK3hXCP0sjdzKrN7Dkz+3HwvMPM9pvZ0eCxPf8yRUTK12sDuZky7U21YZeyoBAj988Bhxc9vw844O7bgQPBcxGRyDo6MMk168qnJQN5hruZbQI+Anxz0eq7gH3B8j7g7ny2ISJSzrJZ5+jgFNu7y6clA/mP3L8GfBHILlq3zt37AYLH7qU+aGb3mlmvmfXG4/E8yxARCUff2CwzyUxZTYOEPMLdzD4KDLr7wZV83t0fdPfd7r67q6trpWWIiITqzYOp5dWWqcnjs7cBHzOzDwP1QKuZfQcYMLMN7t5vZhuAwUIUKiJSjl49lwv3cprjDnmM3N39fnff5O5bgHuAJ939U8DjwN7gbXuBx/KuUkSkTB06O8HmjkZay+SaMvOKMc/9AeADZnYU+EDwXEQkkg71T7BzQ2vYZVwkn7bMAnd/CngqWB4G9hTie0VEytnUXJoTw9N8fFdP2KVcRGeoiois0JFzE7jD9RvLb+SucBcRWaFXzk4AsFPhLiISHYfOTtDeGGN9a33YpVxE4S4iskKH+ie4fuMazCzsUi6icBcRWYFUJsur5ybLsiUDCncRkRU5OjBFMp0ty4OpoHAXEVmRF86MAXDTFW2h1rEchbuIyAo8f2qM9sYYmzsawy5lSQp3EZEVeOHMGDde0VaWB1NB4S4ictmm59K8NjDJjZvawi5lWQp3EZHL9FLfOFkv3347KNxFRC7bC6fHAPg3m9aEW8hbULiLiFym50+PsbmjkbXNdWGXsiyFu4jIZXB3Dp4cLeuWDCjcRUQuy6mRGQYn57h5a0fYpbwlhbuIyGV49vgIALdsUbiLiETGb06MsKYhxvbu8roh9oUU7iIil+E3J0a5eUsHVVXlefLSPIW7iMglGpxMcHxomlu2toddyttSuIuIXKLeE6MA3Fzm/XZQuIuIXLJfvzFMY201N/SU78lL8xTuIiKX6J+Pxrl121pi1eUfneVfoYhIGTg9MsOJ4Rnec3Vn2KVcEoW7iMglePr1IQBuv0bhLiISGU8fHWJ9az1XdZX3/PZ5CncRkbeRyTpPvz7Ee7Z3lu3NOS604nA3syvM7BdmdtjMXjGzzwXrO8xsv5kdDR7Lf0KoiMhbePHMGOOzKX5ve2W0ZCC/kXsa+FN3vw64Ffisme0E7gMOuPt24EDwXESkYv388ADVVcZ7r+kKu5RLtuJwd/d+d/9tsDwJHAZ6gLuAfcHb9gF351mjiEio9h8a4JYtHbQ11oZdyiUrSM/dzLYAu4BngHXu3g+5XwBA9zKfudfMes2sNx6PF6IMEZGCOzk8zWsDU7x/57qwS7kseYe7mTUDPwT+xN0nLvVz7v6gu+92991dXZXzTx0RWV32HxoA4APXraJwN7MYuWD/rrs/GqweMLMNwesbgMH8ShQRCc/+QwPsWNfC5rWNYZdyWfKZLWPAQ8Bhd//rRS89DuwNlvcCj628PBGR8AxMJHj2xAh33rA+7FIuW00en70N+EPgJTN7Plj3Z8ADwMNm9hngFPCJvCoUEQnJP7xwFnf42I0bwy7lsq043N39aWC52fx7Vvq9IiLl4vEXznJDTytXl/ldl5aiM1RFRJZwLD7Fi2fGuevGnrBLWRGFu4jIEh57/ixm8NEbN4Rdyooo3EVELpDOZHm49zTvubqTDWsawi5nRRTuIiIXeOpInP7xBP/hdzaHXcqKKdxFRC7wvWdP0dVSx54KO3FpMYW7iMgiZ0Zn+MWRQe65+YqKuJ3eciq3chGRInjo6eNUm/HJWyq3JQMKdxGRBaPTSX7w7GnuuqmHjW2VeSB1nsJdRCSw79cnmE1l+K/v3RZ2KXlTuIuIAOOzKf7+X07w/uu62b6uJexy8pbPtWVERCrej57r46s/O0Lf2CwAN13RFm5BBaKRu4isWj96ro/7H31pIdgBvvGLN/jRc30hVlUYGrmLyJLmR7Rnx2bZ2NbAF+7cwd27KvM6K8v56s+OMJvKnLduNpXhqz87UvE/q8JdpMJls85UMs1sMsNMMsNMsDydzDCbTJPMOO6ee6877uAOtTVVNMSqqY9V01BbRX2smrbGWtY21fLTl89x/6MvLQRf39gs9z/6EkDFh95iZxeN2C9lfSVRuIsUUL6j3UzWGZ9NMTqTZHQ6yehMKnhMMjKTZGw6xcjCa7nXx2aSZL2wP4cBF37lbCrD//7HQ9x2dSddLXWF3WBI1q+pp388cdH6Sp8GCQp3kYKZ798uHu3e98MXiU/O8a4t7YxOJxmZD+rpXGgPz4f0dC68x2dT+DJBXVtdRXtTjPbGWtoba7l2fSttjTE6mmpprY/RWFdNY201DbEaGmuraarLLdfWGGBUGZjlHgFSmSyzySyzqQyJVG7UPz6bZGgqyVd/dmTJGoamktz85Z/T2VzHdRta2LmhlZ0bW3nXle30tDWQu0FbZXB3Nq5puCjcG2LVfOHOHSFVVTgKdylrpez7ujvJTJa5dJa5VJa5dOa85UQqy9RcmslEionZFJOJNJNz6YXlnx8eYC6dPe87E+ksX/7J4Yu2tTioO5pquW5jK+2NMTqa6oLHWtoaa+lorF0I8Mba6pKF5/eeOXXeQcZ5a5tq+W93XMXh/kkO90/w7V+dIJnJ/czrW+vZvaWd3Ve28zvb1nLt+payDvvv/OtJDp4a5cPvWM8Lp8cjd2xB4S4LyuEAWjbrzKWzJFIZHnu+j6/8v1cXArNvbJYvPvIiL/eNc9PmtiB0Lw7hhXUXvJ5ILfW+LHOpN5cvV0Osmpb6Glrqa97y89/+9M10BEHe3lRLUwmDeiW+cOeO8/4VArmf9X99dOd5fyZSmSxHzk1y8OQovSdH+c3xEX78Yj8A3S113H5NF7df08XvXd1Je1NtyX+O5Tx1ZJC/+PEh3ndtN3/zyXdSVVW+/y9Wyny5fwOW0O7du723tzfsMla1C1sKkPvL/JXff8eyAZ/NOtPJNJOJ9Jsj2kTwPJF7Pjn/OJdbN5vKhW4inWsFJFLZ4DFDIp0luYKAXSxWbdTVVFNXU5X7L7bUcjV1sao3l2uqgufLfG7R+1vqYrQ21NBSH6Olvua8C0vd9sCTS452e9oa+NV978vr5wrDSn7ZuztnRmf59RvD/PJonKePDjE+m8IsN3/8gzvX88Hr13FVV3i3res9McKnHnqGbZ3NfP/eW1nTEAutlnyZ2UF3373kawr38lOKEXQqkw0COM1EIsV/+tazDE0nL3pfY201e65bx2QitfD+hbCeSy/bH55XZdBSH6O5Lje6rY9VUx/Lzcyor1m0HMsFaG5dbv2f/8OhZb93/+dvvyika2uqqA5xBLaSX5BRl8k6L5wZ45dH4jz56iAv9Y0DcFVXEx+8fj0f3LmOGze1lWzkvP/QAH/0/d+yYU0DD/+Xd1f8gWGFewV5u4BwdxKp7MIoeX7EPHneSPnNEM69vui14P2J1KWPkLd2NtFSX7MQ0POj1pa6Rcv1MZqD9kTr/PO6mrz6xJU4Ei6H1lY5Ozs2y/5DAzxx6Bz/emyETNbpbqnjAzvX8cHr13Prtg7qaqoLvt10JsvfPvUGX/v5a9zQs4aH9t5c8cEOCveCW+lf4EzWmZp760D+xi9eZ2oufdFnq6uMlvoaphJp0pcw720+iC8K5PnlupogjHPr/+ejLy05cg8zSDUSjrbxmRRPHhngiVcG+OVrcWaSGZpqq7nt6k7+7bXd3LGjqyC3uDt4cpS/+PEhXjg9xsdu3MgDf/AOGmujcbjxrcI9Gj9hCf3ouT7u++GLJBYd5PvCIy/wq9eH2NLZxMRsivHZFGMzucf5/yZmc62MlcpknY/duDEI68VBfcHouS4X6JfbnphNZpYM0jCnhM0HuEbC0bSmMcbHd23i47s2kUhl+NXrQzz56iBPHYnzxKEBAK5d38J7ru7kXVe2864r2+lurb+k755LZ3jqSJzvPXOKX74Wp7O5jq/fcxN33bR6/uysupH7W4263Z3RmRR9o7P0jc3QN5bg3Pgs8ck5hqaSDE3N8drA5FueMFJbU0VbQ4w1wX9tjTFag+XWC0fPiwJ5vpWx5/88xdklTqooxQhaLQUpB+7O0cEpnjoyyJOvDvLbU2MLB9p72hq4qruZrWsb6WlvCOb315BMZ5lKpDg5MsPRgSl6T46QSGVZ11rHH956JZ++bStNddEby66KtsylBNOFo27ItTu2dzeTymQ5O5a46DoT9bEqOpvr6Gqpo7O5jv3BiGIpr/7lh6iP5dcvVCtC5Hxz6QyHzk5w8OQoL54Z5/jQNMeHppdsXzbWVrO1s4mbt3Tw3h25KZg1FXyrvLcT+bbMcmcGDk3Nsam9gUP9k7w+OMkTrwxc1K/OZJ034lPsuXYdd+zopqetgY1tDWxqb6CnrYG2xth5BwTf6iBfvsEOakWIXKiupppdm9vZtbl9YZ27L0wWmEmmqa2uprGumrVNtWV9/kApFW3kbmYfAr4OVAPfdPcHlntvviP3d3/lwJLXh5hXZXDl2iaOD00vXStw/IGPXNK2NLIWkXJR8pG7mVUD3wA+AJwBfmNmj7v78hOXV+DlvnE+94Pn3jLYf/TZ29ixroWG2uplR92Xc5EgjaxFpBIUqy1zC/C6ux8DMLMfAHcBBQ33zuY6tnU1Mzgxt+RMlJ62hvPuqrLcKdWXOyPk7l09CnMRKWvFOtLQA5xe9PxMsG6Bmd1rZr1m1huPx1e0kfVr6vm//3E3f3n3DTRc0O9eKrTv3tXDV37/Hbmr15ELf7VTRCSKijVyX+qIxnnNfXd/EHgQcj33fDZ2Oa0SjbpFZDUoVrifAa5Y9HwTcLZI2wIU2iIiixWrLfMbYLuZbTWzWuAe4PEibUtERC5QlJG7u6fN7L8DPyM3FfJb7v5KMbYlIiIXK9pJTO7+E+Anxfp+ERFZXnTPyxURWcUU7iIiEaRwFxGJoLK4KqSZxYGTeXxFJzBUoHIqnfbF+bQ/zqf98aYo7Isr3b1rqRfKItzzZWa9y108Z7XRvjif9sf5tD/eFPV9obaMiEgEKdxFRCIoKuH+YNgFlBHti/Npf5xP++NNkd4Xkei5i4jI+aIychcRkUUU7iIiEVTR4W5mHzKzI2b2upndF3Y9pWZmV5jZL8zssJm9YmafC9Z3mNl+MzsaPLa/3XdFhZlVm9lzZvbj4Plq3hdtZvaImb0a/Bl59yrfH58P/p68bGbfN7P6KO+Pig33Rfdp/XfATuCTZrYz3KpKLg38qbtfB9wKfDbYB/cBB9x9O3AgeL5afA44vOj5at4XXwd+6u7XAjeS2y+rcn+YWQ/wx8Bud7+B3NVq7yHC+6Niw51F92l19yQwf5/WVcPd+939t8HyJLm/vD3k9sO+4G37gLtDKbDEzGwT8BHgm4tWr9Z90QrcDjwE4O5Jdx9jle6PQA3QYGY1QCO5GwhFdn9Ucri/7X1aVxMz2wLsAp4B1rl7P+R+AQDdIZZWSl8DvghkF61brftiGxAHvh20qb5pZk2s0v3h7n3AXwGngH5g3N2fIML7o5LD/W3v07pamFkz8EPgT9x9Iux6wmBmHwUG3f1g2LWUiRrgncDfufsuYJoItRwuV9BLvwvYCmwEmszsU+FWVVyVHO4lv09rOTKzGLlg/667PxqsHjCzDcHrG4DBsOoroduAj5nZCXItuveZ2XdYnfsCcn8/zrj7M8HzR8iF/WrdH+8Hjrt73N1TwKPA7xLh/VHJ4b7q79NqZkaup3rY3f960UuPA3uD5b3AY6WurdTc/X533+TuW8j9WXjS3T/FKtwXAO5+DjhtZjuCVXuAQ6zS/UGuHXOrmTUGf2/2kDtGFdn9UdFnqJrZh8n1Wefv0/rlcCsqLTN7D/DPwEu82Wf+M3J994eBzeT+UH/C3UdCKTIEZnYH8D/c/aNmtpZVui/M7CZyB5drgWPAp8kN6Fbr/vhz4N+Tm2X2HPCfgWYiuj8qOtxFRGRpldyWERGRZSjcRUQiSOEuIhJBCncRkQhSuIuIRJDCXUQkghTuIiIR9P8B7p31iEpk/gQAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pylab as pl\n",
    "import math\n",
    "\n",
    "\n",
    "def lagelangri(x, y, ans_x):\n",
    "    ans_y = 0\n",
    "    for i_y in range(len(y)):\n",
    "        temp = 1\n",
    "        for i_x in range(len(x)):\n",
    "            if i_x != i_y:\n",
    "                temp *= (ans_x - x[i_x]) / (x[i_y] - x[i_x])\n",
    "        ans_y += y[i_y] * temp\n",
    "    return ans_y\n",
    "\n",
    "\n",
    "x = np.array([1, 4, 9, 16, 25, 36, 49, 64, 81])\n",
    "y = np.array([1, 2, 3.1, 3.95, 5.01, 6, 6.97, 8, 9.02])\n",
    "print(\"x=2.5 f(x)=\", lagelangri(x, y, 2.5))  #x=2.5的时候  y=\n",
    "print(\"x=40 f(x)=\", lagelangri(x, y, 40))  #x=40的时候  y=\n",
    "xxx = []\n",
    "yyy = []\n",
    "for i in range(0, 900):\n",
    "    xxx.append(i / 10)\n",
    "    yyy.append(lagelangri(x, y, i / 10))\n",
    "\n",
    "pl.scatter(x, y)\n",
    "pl.plot(xxx, yyy)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 分段线性插值"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x=2.5 f(x)= 1.5\n",
      "x=40 f(x)= 6.298461538461538\n",
      "x=2.5 sqrt= 1.5811388300841898\n",
      "x=40 sqrt= 6.324555320336759\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeVklEQVR4nO3deXxU1f3/8dfJvhASIGFJSAhr2BQCARVwAbXgWrTutS61tVaq1Vqs2n7ro6399fcttj/tpqKo7QMXoEaqRYtW0ZKqlECQPQIJkA0ISxIIWWfO748ECxgkhEzunTvv5z+SmSF+PM68c3LuuedjrLWIiIh7hTldgIiIfDkFtYiIyymoRURcTkEtIuJyCmoREZeLCMQ3TU5OtpmZmYH41iIinrRq1aq91tqUtp4LSFBnZmaSn58fiG8tIuJJxpgdJ3pOSx8iIi6noBYRcTkFtYiIyymoRURcTkEtIuJyAdn1ISLihMUFZcxZWkh5VR2pSbHMnp7FzOw0p8s6bQpqEfGExQVlPJy7jromHwBlVXU8nLsOIOjDWksfIuIJc5YWfh7SR9Q1+ZiztNChijqPglpEgt76smrKqurafK78BI8HEy19iEhQqj7cxOI1ZSxYWcLGipoTvi41KbYLqwoMBbWIBA2/3/JJ8T4WrCzh7fW7aGz2Myq1Oz//6igiwgy/+PumY5Y/YiPDmT09y8GKO4eCWkRcb1d1Pa+tLmXByhJ27j9MQkwE1+ekc/2EdEanJX7+urioCH702loamv2kadeHiEhgNfn8vL95DwtXlrCscA9+C2cP6sn9Fw/lktH9iIkM/8LfmZmdRm5BGTV1TSyeNdmBqgNDQS0irrKt8hAL80t4bVUZew810DshmrvOH8x1OelkJsc7XZ4jFNQi4rjDjc28tW4XC1buZOX2A4SHGaYN7831OelckJVCRHj7N6iZANbpFAW1iDjCWsva0moW5JfwxppyDjU0MzA5nh/NGM7XxqXRu3tMx793J9bpBgpqEelSVYcbeb2gZVvd5l0HiYkM49Iz+nF9TjoTB/bEGC/OiU+PglpEAs7vt3y0bR8L8ktYuqFlW90ZaYk8NnM0V45NpXtMZKf9u4wBrLfm1ApqEQmYiuo6FuWXsjC/hNIDdSTGRnLTxAyuy0lnZGp3p8sLGgpqEelUjc1+3tu0mwX5Jfzrs0r8FiYN7sXs6VlMH9W3zW11ncmgNWoRkTZt3XOQBStLyF1dxr7aRvp2j2HW1CFcOz6djF5xTpcX1BTUItJhtQ3NLFlbwYL8ElbtOEBEmOHCEb25YUIG5w1LITys6y8MGmO8tkStoBaRU2OtZU1JFQtWlvDmp+XUNvoYlBLPw5cM5+px/UlJiHa6RM9RUItIu+yvPbKtbief7T5EbGQ4l53Zj+snpJMzoIe21QWQglpETsjvt+Rt3cuC/BLe3bCbRp+fMelJ/J+rzuCKMf1I6MRtdZ2l5WKit9Y+FNQi8gVlVXUsyi9hUX4pZVV1JMVF8vWzM7h+QjrD+2pbXVdTUIsIAA3NPv65cQ8L8ktYvqUSa+Hcock8dMlwvjKqD9ERgd1W11mM8dz9Lu0LamPM/cC3aNmeuA643VpbH8jCRKRrfLb7yLa6Ug4cbqJfYgz3TBvKteP7k95T2+rc4KRBbYxJA+4FRlpr64wxC4EbgBcDXJuIdJLFBWXMWVpIeVUdqUmx3DNtCAAL8kso2FlFZLjh4pF9uC4nnXOHOrOtrvOE7va8CCDWGNMExAHlgStJRDrT4oIyHs5d93mLqrKqOh7KXQfA0N7d+MllI7gqO41e3bStzq1OGtTW2jJjzOPATqAOeMda+07AKxORTvHrf2w+po/gEcndonjn/vM8t63OGO/dQn7S07iNMT2ArwIDgVQg3hhzcxuvu9MYk2+Mya+srOz8SkXklDQ2+8ldXUp5dduXk/YdavRcSHtVe5Y+LgKKrbWVAMaYXGASMP/oF1lr5wJzAXJycrz2A00kaOyvbeTlFTv4y8c72HOwgYgwQ7P/ix/J1KRYB6qTjmhPUO8EzjbGxNGy9HEhkB/QqkTklG3dc5B5edvJXV1KQ7Ofc4cm8+trzuRAbSOPvL7+mOWP2MhwZk/PcrDawGk5jtpbc8X2rFGvMMb8FVgNNAMFtM6cRcRZ1rbcOTgvr5gPCiuJigjj6uw0vjllIMP6JHz+OmPMMbs+Zk/PYmZ2moOVy6lo164Pa+2jwKMBrkVE2qm+ycff1pTxfN52CncfJLlbND+4eBhfPyujzd0bM7PTQiaYvbjsrjsTRYJI5cEG5n+yg/mf7GBfbSPD+yYw55ozuXJsatDcOSinTkEtEgQ276rh+bxiFheU0+jzM214b+6YMpBJg3tp58ZxTAjf8CIiXczvt3z4WSXz8orJ27qXmMgwrpvQn9snD2RwSjeny5MupKAWcZm6Rh+5BaU8n1fMtspa+nSPZvb0LG6amEGP+CinywsKOuZURAJid009f/l4Oy+t2EnV4SbOSEvkievHcukZ/YiKOOm9adLKiytBCmoRh60vq2ZeXjF/X1tOs99y8Yg+fOvcQUzIVNcUaaGgFnGAz295b9Nu5uUVs6J4P/FR4Xz9rAHcPjmTAb3inS4vqIXsedQi0jlqG5pZlF/CCx9tZ8e+w6QlxfLjS0dw3YR0EmPd19ZK3EFBLdIFyqvq+PNH23n5Pzs5WN9MdkYSs6dnMWNUXyLCtf7cmQzGY5cSFdQiAVWw8wDz8op5e/0uAGaM7ssdUwYyLqOHw5VJMFFQi3SyZp+fdzbu5rnlRazeWUVCdATfnJzJrZMy6d9Dra0CzoTgoUwi0j419U0sXFnCC//eTllVHRk943j0ipFcm5NOt2h91KTj9O4ROU079x3mhY+KWZRfyqGGZiZm9uSnV4zkohF9grz3YPDy1nxaQS3SIdZa8nccYN7yYt7ZuIswY7j8zH7cMWUQZ/RPdLq8kObFH40KapEvcXz37h9cPIyIcMO8vGLWllaTGBvJd84fzK3nZNI3McbpcsWjFNQiJ9BW9+4fLvoUCwxKjucXM0fztXFpxEXpY+QmxoPdbfUOEzmBOUsLv9C92wI946P45w/OJ0zrz9JFtNNepA0Hahspq6o74XMKafcyeG5CrRm1yNH21NTz7PIiXlqx84SvUfdu6WoKahGgZP9hnvnXNhbml9Ls83PlmFSy+ibwu/e2hkz3bi/RDS8iHrKt8hB/WraNxWvKCDNwzfj+3HX+4M9PsOuXGKvu3UHGiyfDKqglJG0or+ZPy7bx1voKoiPCuOWcAdx53iD6JR67rBFK3bvFvRTUElJW7TjAH5dt5f3Ne0iIjuC75w/mm1MGktwt2unSpJPoYqJIELLW8tG2ffzh/a18XLSPHnGRPHDxMG6ZlKkzoCUoKKjFs6y1vLdpD39YtpU1JVX0TojmJ5eN4MaJGcTrkCTPMsaow4uI2/n8lrfWVfDHZVvZvOsg/XvE8tjM0Vwzvj8xkeFOlydyyhTU4hlNPj+vF5Tx9AfbKNpby+CUeH5z7RiuHJtKpLqohIyWNWpvTakV1BL06pt8LMwv4ZkPiyirqmNUanf+9PVxzBjVV3cQiicoqCVoHWpo5qVPdvDs8mL2HmogZ0APHrtqNBcMS2k5mEdCltaoRRxWdbiRFz/azgv/3k51XRPnDk1m1tRszhrYUwEtnjyQWkEtQaPyYAPP5RUx/+Md1Db6uHhkH2ZNHcLY9CSnSxMJKAW1uF5ZVR1zP9zGqytLaPL5ufzMVO6eOpjhfbs7XZq4kEHb80S6TPHeWp76YCu5q8swBq7O7s9dFwxmYHK806WJdCkFtbjO5l01/HHZNpasLScyPIybzx7At88bRJqOF5V28OJlCgW1uMaakir+8P5W/rlpN/FR4dx53mDumDKQlASdwyGhrV1BbYxJAp4DRtNy3sk3rbUfB7Au8ajjm8X+8CvD6JsYyx+XbSVv616S4iK5/6Jh3DYpk8Q4ncMhHROq51E/CfzDWnuNMSYKiAtgTeJRbTWLfWDRp/gtpCRE88ilw7nprAF00zkccho8uPJx8qA2xnQHzgNuA7DWNgKNgS1LvKitZrF+C4mxkSx/cKrO4ZBO4635dPua2w4CKoEXjDEFxpjnjDFfuOxujLnTGJNvjMmvrKzs9EIl+J2oWWxNXZNCWjqNFy8mtieoI4BxwFPW2mygFnjo+BdZa+daa3OstTkpKSmdXKYEs7pGH79997MTPq9msSJfrj2LgaVAqbV2RevXf6WNoBY5nrWWN9dW8Ku3NlFRXc+4jCQ2lNfQ0Oz//DVqFiudLSRveLHW7jLGlBhjsqy1hcCFwMbAlybBbH1ZNT97cwMrtx9gdFp3fndjNhMye35h14eaxYqcXHsvr98DvNS646MIuD1wJUkw23uogceXFrIgv4Re8VH879fO4Jrx6YS3HjeqZrESaMaE6HnU1to1QE5gS5Fg1tjs5y8fb+fJf26hrsnHt6YM5J4Lh9I9RnuhRU6XNqzKaVtWuIdf/H0jRZW1XJCVwv9cPpLBKd2cLktCWMitUYucSFHlIX7x940sK6xkUHI8L9w2ganDeztdloQ4L27PU1DLKaupb+L3723hxY+2ExMRzo8vHcGtkzKJilBfQnEHj02oFdTSfn6/ZdGqEuYsLWRfbSPXjU/nh9OzdGiSuIz3ptQKammX/O37+dmbG1lXVk3OgB68cNtEzuif6HRZIiFBQS1fqqK6jl+9tZk3Pi2nb/cYnrxhLFeOSVVvQnEtY3QxUUJEfZOPuf8q4qkPtuGzlnunDeGuCwYTF6W3jEhX06dOjmGt5e31u/jlkk2UVdVx6Rl9efiSEaT31Mm2Eky8NaVWUMvnNlXU8LM3N/BJ0X6G903g5W+fxaTByU6XJXJKvLgop6AW9tc28pt3CnnlPzvpHhvJL2aO5sYJ6USEa7udBCetUYtnNPn8zP9kB//v3c+obfRxyzmZ3HfRUJLiopwuTaTDvHidW0EdopZvqeTnb25ky55DTBmSzE+vGMmwPglOlyXSKTw2oVZQh5od+2p5bMkm3t24m4yeccz9xnguHtlH2+3EM4wHV6kV1B51/LnP90wbwo79h5m3vJiIcMODM7K4Y8pAoiPUAkvE7RTUHtRWt++HctcBcPW4NH40Yzh9usc4WaJIwLTc8OKtxQ8FtQe11e0bIKVbNL+9bmzXFyQip0X7rzyo/ATdvvceaujiSkSc4a35tILac/x+S7eYtn9RUrdvCQXeu5SooPaUfYcauO3FlRysbyb8uF0c6vYtocRjS9Rao/aK/O37+d7LBeyvbeSXV40mLjKcx9/5TN2+JeR4caupgjrIWWt5dnkR//uPQtKSYsm9exKj01rOib5qXH+HqxORzqCgDmLVh5t4YNGn/HPTbmaM6suvrz1TXb9F0PY8cYlPS6qY9fJqdlXX89PLR3L75ExP/sonIgrqoGOt5S8f7+CxJRvpnRDDwrvOYVxGD6fLEnEVb82nFdRB5WB9Ew+9to4l6yqYNrw3v7l2DD3iddKdyNG8+IulgjpIbCyvYdbLq9m5/zAPXTKcO88dRFiYB9+RIp3BY1NqBbXLWWtZsLKER9/YQFJcJK98+2wmDuzpdFkirqXT86RLHW5s5ievrye3oIwpQ5J54oaxJHeLdrosEdfz2IRaQe1WW3Yf5O6XVrO18hD3XzSM700bQriWOkROSmvU0iVeLyjlkdz1xEeHM/+Os5g8RA1mRUKZgtpF6pt8/OzNDbzynxImDuzJ72/M1rnRIqfIoBteJECK99Zy90ur2VRRw90XDOYHFw9TF3ARARTUrrBkbQU/em0tEeGGF26bwNThvZ0uSSSoeWs+raB2VEOzj1+9tZkXP9pOdkYSf7hpHGk6M1rktOhionSasqo67p6/ik9Lq/nWlIE8OGM4URFa6hDpDB5bom5/UBtjwoF8oMxae3ngSvK+5VsqufeVApp9lqdvHs+M0X2dLknEM7x4ONmpzKi/D2wCugeoFk9aXFDGnKWFlFfV0S8xhuyMHry9voKhvRN4+hvjGZgc73SJIp5jPbZK3a7ftY0x/YHLgOcCW463LC4o4+HcdZRV1WGB8up6lqyrYGx6Eq/PmqSQFgkA782n298z8QngQcB/ohcYY+40xuQbY/IrKys7o7agN2dpIXVNvi88vqu6nrgoXR4QkfY5aVAbYy4H9lhrV33Z66y1c621OdbanJSUlE4rMJiVV9W1+XhFdX0XVyISWrx2MbE9M+rJwJXGmO3Aq8A0Y8z8gFblAY3NfuKiwtt8LlVb8EQCx4NrHycNamvtw9ba/tbaTOAG4H1r7c0BryyI7a6p56ZnP6G20feFg5RiI8OZPT3LocpEQoPHJtTaR93ZVhTtY9bLBRxubOb3N2bj89vPd32kJsUye3oWM7PTnC5TxLNC/jxqa+0HwAcBqSTIWWuZl1fMr97ezICecbz87bMY1icBQMEs0tU8NqXWjLoT1DY08+Bra1mytoLpo/rw+LVjSIiJdLoskZDkwftdFNSna+ueQ9w1fxVFlYd46JLhfOe8QZ68M0okmHjthhcF9Wn4x/oKfrhoLdERYcy/4ywm6YB/Ecd5cZqkoO6AZp+fOe8U8syHRYxJT+Kpr4/TljsRCRgF9Snae6iBe14u4OOifdx8dgb/c/lIoiPa3i8tIs7w2g0vCupTsHrnAe6ev5oDhxt5/NoxXDO+v9MlichxvHiJSEHdDtZa5q/Yyc/f3EDfxBhy757EqNREp8sSkRPw2IRaQX0ydY0+frx4Hbmry5ialcIT12eTGKetdyJuFfI3vISaHftquWv+ajbvquG+i4Zy77ShhIV5700g4jXqQh4i3t+8m/teXYMxhudvm8DULDWcFQkGWqMOAT6/5cn3tvC797YwKrU7T988nvSecU6XJSKnwFvzaQX1Ma2y+naPITEuks27DnLN+P48NnM0MZHaeicizgrpoD7SKutIF5aKmnoqauq5dnx/fn3NmboVXCQIefFT295WXJ50olZZH23bp5AWCWIeu5YY2kF9olZZJ3pcRIKABydZIR3U/RJj2nxc53aIiJuEdFCfPajXFx5TqyyR4Oa9+XQIB/WW3QdZsq6Ckf26k5oYgwHSkmL51dVnqCOLiAd46aaXkNz10djs5/uvrqFbdAR//uZEUhKinS5JRDqJB5eoQzOof/vuZ2ysqOHZW3IU0iLieiG39PFJ0T6e+dc2bpyYzsUj+zhdjogEiIdWPkIrqKvrmnhg4acM6BnHTy4b6XQ5IhIAOj0vyD36t/Xsqqnnr3edQ3x0SP2ni4QcD02oQ2dG/can5SxeU86904aSndHD6XJEJEC8eDExJIK6vKqOn7y+juyMJGZNHex0OSLSBby0Pc/zQe33Wx5Y+CnNfssT148lItzz/8kiIc2DE2rvB/W8vGI+LtrHo1eMZECveKfLEZEu4p35tMeDelNFDXOWFvKVkX24Lifd6XJERDrEs0Fd3+TjvlfXkBgXyf/9ms6WFgkVXvyoe3aP2pylhRTuPsiLt0+gZ3yU0+WISBfz0LVEb86o87bsZV5eMbecM4AL1JRWJKR48bdnzwV11eFGHli0hsEp8Tx8yQinyxERh1gPXU701NKHtZYfv76efYcamXfrBGKj1JhWRIKfp2bUrxeUsWRdBT/4yjBGpyU6XY6IOEhr1C5Usv8wP/3bBiZm9uQ75+nuQ5FQ5cEl6pMvfRhj0oG/AH0BPzDXWvtkoAtrj8UFZcxZWkhZVR1R4WEYA7+5bgzhYR78PyUiIas9M+pm4AFr7QjgbGCWMcbxM0IXF5TxcO46ylo7hjf6/PitZdWOAw5XJiLSuU4a1NbaCmvt6tY/HwQ2AY43FZyztJC6Jt8xjzX5LHOWFjpUkYi4gRfPoz6lNWpjTCaQDaxo47k7jTH5xpj8ysrKTirvxMpbZ9LtfVxEQktIXkw0xnQDXgPus9bWHP+8tXautTbHWpuTkpLSmTW2KTUp9pQeF5HQ4MWLie0KamNMJC0h/ZK1NjewJbXP7OlZRBx30TA2MpzZ07McqkhE3CSkbngxLfdjzgM2WWt/G/iS2ueikX2IDA8jItzS0OQnNSmW2dOzmJnt+PK5iDjIgxPqdt2ZOBn4BrDOGLOm9bFHrLVvBayqdnhlxU7qmnz8bdZkxqQnOVmKiLiQl9aoTxrU1to8XPZDqqHZx3N5RUwa3EshLSKeF5R3Ji4uKGN3TQPfvUB3IIrIsY5cTPTQhDr4gtrntzzzYRGj07ozZUiy0+WIiARc0AX1Oxt2UbS3lu+eP8ST586KyOkJ+RtenGat5akPt5HZK44Zo/s6XY6IuJj10NXEoArqj7btY21pNd85f7AOXhKRNnnxF+2gCuqnPthG74Rorh6nvdIi8uW8M58OoqBeW1pF3ta93DFlINER6twiIqEjaIL66Q+3kRATwU1nZThdiogEAQ8tUQdHUBdVHuLt9bu45ZwBJMREOl2OiEiXCoqgnvuvIqLCw7ht0kCnSxERl/Pitl3XB/Wu6npeW13KdTnppCREO12OiAQLLX10nef/XYzfwp3nDXK6FBEJAt6bT7s8qKsPN/HSJzu4/Mx+pPeMc7ocEQkiIXUetROO7i4OkNUnweGKRCRYeHCJ2n0z6uO7iwP8/v2tLC4oc7AqEQk22p4XQG11F69r8qm7uIiELNcFtbqLi8jpOLLy4aEJtfuCWt3FRUSO5bqgnj09i9jIY8/yUHdxEWkvL97w4rpdH0e6iD+Su47DTT7S1F1cRDrAS+dRuy6ooSWsPyjcw5qSKj6YPdXpckQkiHhwQu2+pY8jmv1WzQFEpMO8M592cVD7FNQi0gFeTA3XBnXLjNq15YmIy3loidq9Qe3zWyI0oxYRcW9Qa41aRDqk9Wqilw5lcm1Q+/x+zahFRHBxUDf7NKMWkVPnxdRwbVD7/JaIcC8OuYh0Ce+sfLg3qLXrQ0Q6Qje8dCHt+hCR0+GhCbV7g1q7PkREWrg2qLXrQ0Q6wrReTtQNL11AM2oRkRauDWqtUYtIRxy5mBhyN7wYY2YYYwqNMVuNMQ8Fuig4so/atT9HRES6zEmT0BgTDvwRuAQYCdxojBkZ6MI0oxaRjvBiarSnccBEYKu1tgjAGPMq8FVgY2cXc8Xv86hv7UBeeaiBMAW1iHTQTc+u6PLJXo+4KBbedU6nf9/2BHUaUHLU16XAWce/yBhzJ3AnQEZGRoeKGZwST6PPD8CwPglcpfZbInKKJg9JZubY1M+zpCt1j4kMyPdtT1C39SPpC6v01tq5wFyAnJycDq3iP3FDdkf+mojI59J7xnkuS9pzta4USD/q6/5AeWDKERGR47UnqFcCQ40xA40xUcANwBuBLUtERI446dKHtbbZGPM9YCkQDjxvrd0Q8MpERARo3xo11tq3gLcCXIuIiLRBd5SIiLicglpExOUU1CIiLqegFhFxOWMDcGirMaYS2NHBv54M7O3EcoKdxuO/NBbH0ngcK9jHY4C1NqWtJwIS1KfDGJNvrc1xug630Hj8l8biWBqPY3l5PLT0ISLicgpqERGXc2NQz3W6AJfRePyXxuJYGo9jeXY8XLdGLSIix3LjjFpERI6ioBYRcTnXBLUTDXTdxBiTboxZZozZZIzZYIz5fuvjPY0x7xpjtrT+s4fTtXYVY0y4MabAGPP31q9DeSySjDF/NcZsbn2PnBPi43F/6+dkvTHmFWNMjJfHwxVB7VQDXZdpBh6w1o4AzgZmtY7BQ8B71tqhwHutX4eK7wObjvo6lMfiSeAf1trhwBhaxiUkx8MYkwbcC+RYa0fTcvzyDXh4PFwR1BzVQNda2wgcaaAbMqy1Fdba1a1/PkjLBzGNlnH4c+vL/gzMdKTALmaM6Q9cBjx31MOhOhbdgfOAeQDW2kZrbRUhOh6tIoBYY0wEEEdL1ynPjodbgrqtBroh29nWGJMJZAMrgD7W2gpoCXOgt4OldaUngAeBozuUhupYDAIqgRdal4KeM8bEE6LjYa0tAx4HdgIVQLW19h08PB5uCep2NdANBcaYbsBrwH3W2hqn63GCMeZyYI+1dpXTtbhEBDAOeMpamw3U4qFf609V69rzV4GBQCoQb4y52dmqAsstQa0GuoAxJpKWkH7JWpvb+vBuY0y/1uf7AXucqq8LTQauNMZsp2UZbJoxZj6hORbQ8vkotdauaP36r7QEd6iOx0VAsbW20lrbBOQCk/DweLglqEO+ga4xxtCyBrnJWvvbo556A7i19c+3An/r6tq6mrX2YWttf2ttJi3vhfettTcTgmMBYK3dBZQYY7JaH7oQ2EiIjgctSx5nG2PiWj83F9JyTcez4+GaOxONMZfSsi55pIHuL52tqGsZY6YAy4F1/Hdd9hFa1qkXAhm0vEGvtdbud6RIBxhjLgB+aK293BjTixAdC2PMWFourEYBRcDttEy0QnU8fgZcT8tuqQLgW0A3PDoerglqERFpm1uWPkRE5AQU1CIiLqegFhFxOQW1iIjLKahFRFxOQS0i4nIKahERl/v/3z1IOXOSWn0AAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pylab as pl\n",
    "import math\n",
    "\n",
    "\n",
    "def fenduan(x, y, ans_x):  #  ans_x :f(ans_x)\n",
    "    ans_y = 0\n",
    "    for i in range(len(x) - 1):\n",
    "        if x[i] <= ans_x <= x[i + 1]:\n",
    "            ans_y = y[i] * ((ans_x - x[i + 1]) / (x[i] - x[i + 1])) + y[i + 1] * ((ans_x - x[i]) / (x[i + 1] - x[i]))\n",
    "\n",
    "    return ans_y\n",
    "\n",
    "\n",
    "x = np.array([1, 4, 9, 16, 25, 36, 49, 64, 81])\n",
    "y = np.array([1, 2, 3.1, 3.95, 5.01, 6, 6.97, 8, 9.02])\n",
    "\n",
    "xxx = []\n",
    "yyy = []\n",
    "for i in range(0, 900):\n",
    "    xxx.append(i / 10)\n",
    "    yyy.append(fenduan(x, y, i / 10))\n",
    "\n",
    "pl.scatter(x, y)\n",
    "pl.plot(xxx, yyy)\n",
    "\n",
    "#精度验证\n",
    "print(\"x=2.5 f(x)=\", fenduan(x, y, 2.5))  #x=2.5的时候  y=\n",
    "print(\"x=40 f(x)=\", fenduan(x, y, 40))  #x=40的时候  y=\n",
    "print(\"x=2.5 sqrt=\", math.sqrt(2.5))\n",
    "print(\"x=40 sqrt=\", math.sqrt(40))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 复合梯形公式"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "梯形公式 0.13090909090909092\n"
     ]
    }
   ],
   "source": [
    "#被积函数\n",
    "def fun(x):\n",
    "    return 1 / x\n",
    "\n",
    "\n",
    "#复合梯形\n",
    "def tx(a, b, n):\n",
    "    h = (b - a) / n\n",
    "    x = a\n",
    "    s = fun(x) - fun(b)\n",
    "    for k in range(1, n + 1):\n",
    "        x = x + h\n",
    "        s = s + 2 * fun(x)\n",
    "        result = (h / 2) * s\n",
    "        return result\n",
    "\n",
    "\n",
    "a = 1\n",
    "b = 5\n",
    "n1 = 40\n",
    "n2 = 20\n",
    "t = tx(a, b, n1)  #起始值、结束值、分成几段\n",
    "print(\"梯形公式\", t)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 复合辛普森求积公式"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "辛普森公式 1.6094411643415028\n"
     ]
    }
   ],
   "source": [
    "def fun(x):\n",
    "    return 1 / x\n",
    "\n",
    "\n",
    "#复合辛普森\n",
    "def xps(a, b, n):\n",
    "    h = (b - a) / n\n",
    "    x = a\n",
    "    s = fun(x) - fun(b)\n",
    "    for k in range(1, n + 1):\n",
    "        x = x + h / 2\n",
    "        s = s + 4 * fun(x)\n",
    "        x = x + h / 2\n",
    "        s = s + 2 * fun(x)\n",
    "    result = (h / 6) * s\n",
    "    return result\n",
    "\n",
    "\n",
    "a = 1\n",
    "b = 5\n",
    "n2 = 20\n",
    "p = xps(a, b, n2)  #起始值、结束值、分成几段\n",
    "print(\"辛普森公式\", p)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 高斯求积公式"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.60943789295249\n",
      "1.948161032494511e-08\n"
     ]
    }
   ],
   "source": [
    "from math import exp\n",
    "from math import sqrt\n",
    "from math import log\n",
    "\n",
    "f = lambda x: 1 / x\n",
    "\n",
    "\n",
    "def G_quad(f, a, b, N):\n",
    "    h = (b - a) / N;\n",
    "\n",
    "    t = [-sqrt(3 / 5), 0, sqrt(3 / 5)]\n",
    "    A = [5 / 9, 8 / 9, 5 / 9]\n",
    "    x_fromt = [0, 0, 0]\n",
    "    Sum_G = 0\n",
    "    for k in range(N):\n",
    "        for i in range(3):\n",
    "            x_fromt[i] = h / 2 * t[i] + a + (k + 1 / 2) * h\n",
    "        Sum_G = Sum_G + h / 2 * (A[0] * f(x_fromt[0]) + A[1] * f(x_fromt[1]) + A[2] * f(x_fromt[2]))\n",
    "    return Sum_G\n",
    "\n",
    "\n",
    "quadG = G_quad(f, 1, 5, 15)  #函数体、起始值、结束值、分成几段. 这里是1，5之间分成15段\n",
    "\n",
    "print(quadG)\n",
    "print(log(5) - quadG)  #截断误差\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 高斯列主元"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "高斯列主元：(1, array([0.        , 0.        , 0.        , 0.50819672]))\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def Gauss(A, b):\n",
    "    n = len(b);\n",
    "    index = 1;\n",
    "    x = np.zeros(n)\n",
    "    for k in range(n):\n",
    "        #＇选主元＇\n",
    "        a_max = 0\n",
    "        for i in range(k, n):\n",
    "            if abs(A[i][k]) > a_max:\n",
    "                a_max = abs(A[i][k])\n",
    "                r = i\n",
    "            if a_max < 0.00000001:\n",
    "                index = 0\n",
    "                return\n",
    "                #＇交换两行＇\n",
    "            if r > k:\n",
    "                for j in range(k, n):\n",
    "                    z = A[k][j];\n",
    "                    A[k][j] = A[r][j];\n",
    "                    A[r][j] = z\n",
    "                z = b[k];\n",
    "                b[k] = b[r];\n",
    "                b[r] = z;\n",
    "                #＂消元计算＂\n",
    "            for i in range(k + 1, n):\n",
    "                m = A[i][k] / A[k][k]\n",
    "                for j in range(k + 1, n):\n",
    "                    A[i][j] = A[i][j] - m * A[k][j]\n",
    "                b[i] = b[i] - m * b[k]\n",
    "            #回代过程\n",
    "            if abs(A[n - 1][n - 1]) < 0.0000001:\n",
    "                index = 0\n",
    "                return\n",
    "            for k in range(n - 1, 0 - 1, -1):\n",
    "                for j in range(k + 1, n): b[k] = b[k] - A[k][j] * x[j]\n",
    "                x[k] = b[k] / A[k][k]\n",
    "                return index, x\n",
    "\n",
    "\n",
    "A = [[15, 2, 3, 2], [3, 15, 6, 5], [5, 6, 18, 3], [6, 5, 4, 13]]\n",
    "b = [22, 28, 13, 15]\n",
    "print(\"高斯列主元：{}\".format(Gauss(A, b)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 雅可比迭代"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 1, array([ 6. ,  7.5, -6. ]))\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "#A:左边的3*3的矩阵 b:右边3*1的矩阵 x0:x的初值\n",
    "# it_max:限制的最大迭代次数\n",
    "#ep:设置精度\n",
    "def Jacobi(A, b, x0, it_max, ep):\n",
    "    D = np.diag(np.diag(A))\n",
    "    U = -np.triu(A, 1)\n",
    "    L = -np.tril(A, -1)\n",
    "    B = np.dot(np.linalg.inv(D), (L + U))\n",
    "    f = np.dot(np.linalg.inv(D), b)\n",
    "    x = np.dot(B, x0) + f\n",
    "    k = 1\n",
    "    index = 1\n",
    "    while np.linalg.norm(x - x0) >= ep:\n",
    "        x = x0\n",
    "        x = np.dot(B, x0) + f\n",
    "        k = k + 1\n",
    "        if k > it_max:\n",
    "            index = 0\n",
    "            break\n",
    "        return k, index, x\n",
    "\n",
    "\n",
    "A = [[4, 3, 0], [3, 4, -1], [0, -1, 4]]\n",
    "b = [24, 30, -24]\n",
    "x0 = [0, 0, 0]\n",
    "print(Jacobi(A, b, x0, 100, 0.000001))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 高斯赛德尔迭代"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 1, array([ 6.  ,  3.  , -5.25]))\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "#A:左边的3*3的矩阵 b:右边3*1的矩阵 x0:x的初值\n",
    "# it_max:限制的最大迭代次数\n",
    "#ep:设置精度\n",
    "def G_S(A, b, x0, it_max, ep):\n",
    "    D = np.diag(np.diag(A))\n",
    "    U = -np.triu(A, 1)\n",
    "    L = -np.tril(A, -1)\n",
    "    B = np.dot(np.linalg.inv(D - L), U)\n",
    "    f = np.dot(np.linalg.inv(D - L), b)\n",
    "    x = np.dot(B, x0) + f\n",
    "    k = 1\n",
    "    index = 1\n",
    "    while np.linalg.norm(x - x0) >= ep:\n",
    "        x = x0\n",
    "        x = np.dot(B, x0) + f\n",
    "        k = k + 1\n",
    "        if k > it_max:\n",
    "            index = 0\n",
    "            break\n",
    "        return k, index, x\n",
    "\n",
    "\n",
    "A = [[4, 3, 0], [3, 4, -1], [0, -1, 4]]\n",
    "b = [24, 30, -24]\n",
    "x0 = [0, 0, 0]\n",
    "print(G_S(A, b, x0, 100, 0.000001))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 二分法"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 20, 0.25753116607666016)\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "#原函数\n",
    "f = lambda x: x ** 2 - 3 * x + 2 - math.e ** x\n",
    "\n",
    "\n",
    "#迭代法\n",
    "#f:原函数 起始值 结束值 设置精度\n",
    "def b(f, a, b, ep):\n",
    "    fa = f(a)\n",
    "    fb = f(b)\n",
    "\n",
    "    if fa * fb > 0:\n",
    "        return 0\n",
    "    k = 1\n",
    "    while abs(b - a) / 2 > ep:\n",
    "        x = (a + b) / 2\n",
    "        fx = f(x)\n",
    "        if fx == 0:\n",
    "            return 1, k, x\n",
    "        elif fx * fa < 0:\n",
    "            fb = fx\n",
    "            b = x\n",
    "        else:\n",
    "            fa = fx\n",
    "            a = x\n",
    "        k += 1\n",
    "\n",
    "    x = (a + b) / 2\n",
    "    return 1, k, x\n",
    "\n",
    "\n",
    "ans_b = b(f, 0, 1, 10 ** -6)\n",
    "print(ans_b)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 牛顿法"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 4, 0.257530285439827)\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "#原函数\n",
    "f = lambda x: x ** 2 - 3 * x + 2 - math.e ** x\n",
    "#原函数的导数\n",
    "fd = lambda x: 2 * x - 3 - math.e ** x\n",
    "\n",
    "\n",
    "#f:原函数 导数 起始值 结束值 设置精度\n",
    "def n(f, fd, x, ep, it_max):\n",
    "    index = 0\n",
    "    k = 1\n",
    "    while k < it_max:\n",
    "        x1 = x\n",
    "        fv = f(x)\n",
    "        fdv = fd(x)\n",
    "        if abs(fdv) < ep:\n",
    "            break\n",
    "        x = x1 - fv / fdv\n",
    "        if abs(x - x1) < ep:\n",
    "            index = 1\n",
    "            break\n",
    "        k += 1\n",
    "    return index, k, x\n",
    "\n",
    "\n",
    "ans_n = n(f, fd, 1.5, 10 ** -6, 20)\n",
    "print((ans_n))\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 欧拉法"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[1;31mKeyboardInterrupt\u001B[0m                         Traceback (most recent call last)",
      "Input \u001B[1;32mIn [46]\u001B[0m, in \u001B[0;36m<cell line: 6>\u001B[1;34m()\u001B[0m\n\u001B[0;32m      3\u001B[0m \u001B[38;5;28;01mimport\u001B[39;00m \u001B[38;5;21;01mpylab\u001B[39;00m \u001B[38;5;28;01mas\u001B[39;00m \u001B[38;5;21;01mpl\u001B[39;00m\n\u001B[0;32m      5\u001B[0m h \u001B[38;5;241m=\u001B[39m \u001B[38;5;241m0.01\u001B[39m\n\u001B[1;32m----> 6\u001B[0m n \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;43mint\u001B[39;49m(\u001B[38;5;241m1\u001B[39m \u001B[38;5;241m/\u001B[39m h)\n\u001B[0;32m      7\u001B[0m f \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;01mlambda\u001B[39;00m x, y: \u001B[38;5;241m-\u001B[39m\u001B[38;5;241m50\u001B[39m \u001B[38;5;241m*\u001B[39m y \u001B[38;5;241m+\u001B[39m \u001B[38;5;241m50\u001B[39m \u001B[38;5;241m*\u001B[39m x \u001B[38;5;241m*\u001B[39m\u001B[38;5;241m*\u001B[39m \u001B[38;5;241m2\u001B[39m \u001B[38;5;241m+\u001B[39m \u001B[38;5;241m2\u001B[39m \u001B[38;5;241m*\u001B[39m x\n\u001B[0;32m      9\u001B[0m \u001B[38;5;66;03m#f:函数体\u001B[39;00m\n\u001B[0;32m     10\u001B[0m \u001B[38;5;66;03m# x0初值\u001B[39;00m\n\u001B[0;32m     11\u001B[0m \u001B[38;5;66;03m# y0初值\u001B[39;00m\n\u001B[0;32m     12\u001B[0m \u001B[38;5;66;03m# h 间距\u001B[39;00m\n",
      "Input \u001B[1;32mIn [46]\u001B[0m, in \u001B[0;36m<cell line: 6>\u001B[1;34m()\u001B[0m\n\u001B[0;32m      3\u001B[0m \u001B[38;5;28;01mimport\u001B[39;00m \u001B[38;5;21;01mpylab\u001B[39;00m \u001B[38;5;28;01mas\u001B[39;00m \u001B[38;5;21;01mpl\u001B[39;00m\n\u001B[0;32m      5\u001B[0m h \u001B[38;5;241m=\u001B[39m \u001B[38;5;241m0.01\u001B[39m\n\u001B[1;32m----> 6\u001B[0m n \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;43mint\u001B[39;49m(\u001B[38;5;241m1\u001B[39m \u001B[38;5;241m/\u001B[39m h)\n\u001B[0;32m      7\u001B[0m f \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;01mlambda\u001B[39;00m x, y: \u001B[38;5;241m-\u001B[39m\u001B[38;5;241m50\u001B[39m \u001B[38;5;241m*\u001B[39m y \u001B[38;5;241m+\u001B[39m \u001B[38;5;241m50\u001B[39m \u001B[38;5;241m*\u001B[39m x \u001B[38;5;241m*\u001B[39m\u001B[38;5;241m*\u001B[39m \u001B[38;5;241m2\u001B[39m \u001B[38;5;241m+\u001B[39m \u001B[38;5;241m2\u001B[39m \u001B[38;5;241m*\u001B[39m x\n\u001B[0;32m      9\u001B[0m \u001B[38;5;66;03m#f:函数体\u001B[39;00m\n\u001B[0;32m     10\u001B[0m \u001B[38;5;66;03m# x0初值\u001B[39;00m\n\u001B[0;32m     11\u001B[0m \u001B[38;5;66;03m# y0初值\u001B[39;00m\n\u001B[0;32m     12\u001B[0m \u001B[38;5;66;03m# h 间距\u001B[39;00m\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_310_64.pyx:1180\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_310_64.SafeCallWrapper.__call__\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_310_64.pyx:621\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_310_64.PyDBFrame.trace_dispatch\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_310_64.pyx:930\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_310_64.PyDBFrame.trace_dispatch\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_310_64.pyx:921\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_310_64.PyDBFrame.trace_dispatch\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32m_pydevd_bundle\\pydevd_cython_win32_310_64.pyx:318\u001B[0m, in \u001B[0;36m_pydevd_bundle.pydevd_cython_win32_310_64.PyDBFrame.do_wait_suspend\u001B[1;34m()\u001B[0m\n",
      "File \u001B[1;32mD:\\app\\apps\\PyCharm-P\\ch-0\\213.7172.26\\plugins\\python\\helpers\\pydev\\pydevd.py:1147\u001B[0m, in \u001B[0;36mPyDB.do_wait_suspend\u001B[1;34m(self, thread, frame, event, arg, send_suspend_message, is_unhandled_exception)\u001B[0m\n\u001B[0;32m   1144\u001B[0m         from_this_thread\u001B[38;5;241m.\u001B[39mappend(frame_id)\n\u001B[0;32m   1146\u001B[0m \u001B[38;5;28;01mwith\u001B[39;00m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39m_threads_suspended_single_notification\u001B[38;5;241m.\u001B[39mnotify_thread_suspended(thread_id, stop_reason):\n\u001B[1;32m-> 1147\u001B[0m     \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_do_wait_suspend\u001B[49m\u001B[43m(\u001B[49m\u001B[43mthread\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mframe\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mevent\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43marg\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43msuspend_type\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mfrom_this_thread\u001B[49m\u001B[43m)\u001B[49m\n",
      "File \u001B[1;32mD:\\app\\apps\\PyCharm-P\\ch-0\\213.7172.26\\plugins\\python\\helpers\\pydev\\pydevd.py:1162\u001B[0m, in \u001B[0;36mPyDB._do_wait_suspend\u001B[1;34m(self, thread, frame, event, arg, suspend_type, from_this_thread)\u001B[0m\n\u001B[0;32m   1159\u001B[0m             \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39m_call_mpl_hook()\n\u001B[0;32m   1161\u001B[0m         \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mprocess_internal_commands()\n\u001B[1;32m-> 1162\u001B[0m         \u001B[43mtime\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43msleep\u001B[49m\u001B[43m(\u001B[49m\u001B[38;5;241;43m0.01\u001B[39;49m\u001B[43m)\u001B[49m\n\u001B[0;32m   1164\u001B[0m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mcancel_async_evaluation(get_current_thread_id(thread), \u001B[38;5;28mstr\u001B[39m(\u001B[38;5;28mid\u001B[39m(frame)))\n\u001B[0;32m   1166\u001B[0m \u001B[38;5;66;03m# process any stepping instructions\u001B[39;00m\n",
      "\u001B[1;31mKeyboardInterrupt\u001B[0m: "
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import pylab as pl\n",
    "\n",
    "h = 0.01\n",
    "n = int(1 / h)\n",
    "f = lambda x, y: -50 * y + 50 * x ** 2 + 2 * x\n",
    "\n",
    "# f: 函数体\n",
    "# x0 初值\n",
    "# y0 初值\n",
    "# h 间距\n",
    "# n 不用你管\n",
    "def E(f, x0, y0, h, N):\n",
    "    x = np.zeros(N + 1)\n",
    "    y = np.zeros(N + 1)\n",
    "    x[0] = x0\n",
    "    y[0] = y0\n",
    "    for i in range(N):\n",
    "        x[i + 1] = x[i] + h\n",
    "        y[i + 1] = y[i] + h * f(x[i], y[i])\n",
    "    return x, y\n",
    "\n",
    "\n",
    "[xx, yy] = E(f, 0, 1 / 3, h, n)\n",
    "print(xx)\n",
    "print(yy)\n",
    "\n",
    "xxx = np.arange(0, 1 + h, h)\n",
    "yyy = []\n",
    "for i in xxx:\n",
    "    yyy.append((math.e ** (-50 * i) / 3) + i ** 2)\n",
    "    # yyy.append(xxx)\n",
    "pl.plot(xx, yy, label=\"Euler\", linestyle=\":\")\n",
    "\n",
    "pl.plot(xxx, yyy, color=\"red\", label=\"right\", alpha=0.3)\n",
    "pl.legend()\n",
    "pl.show()\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 改进欧拉法"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1  0.11 0.12 0.13\n",
      " 0.14 0.15 0.16 0.17 0.18 0.19 0.2  0.21 0.22 0.23 0.24 0.25 0.26 0.27\n",
      " 0.28 0.29 0.3  0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4  0.41\n",
      " 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5  0.51 0.52 0.53 0.54 0.55\n",
      " 0.56 0.57 0.58 0.59 0.6  0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69\n",
      " 0.7  0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8  0.81 0.82 0.83\n",
      " 0.84 0.85 0.86 0.87 0.88 0.89 0.9  0.91 0.92 0.93 0.94 0.95 0.96 0.97\n",
      " 0.98 0.99 1.  ] [0.33333333 0.20845833 0.13064896 0.0823306  0.05251912 0.03434945\n",
      " 0.02353091 0.01738182 0.01422614 0.01301633 0.01309771 0.01406107\n",
      " 0.01565067 0.01770667 0.02012917 0.02285573 0.02584733 0.02907958\n",
      " 0.03253724 0.03621077 0.04009423 0.0441839  0.04847744 0.0529734\n",
      " 0.05767087 0.0625693  0.06766831 0.07296769 0.07846731 0.08416707\n",
      " 0.09006692 0.09616682 0.10246676 0.10896673 0.1156667  0.12256669\n",
      " 0.12966668 0.13696668 0.14446667 0.15216667 0.16006667 0.16816667\n",
      " 0.17646667 0.18496667 0.19366667 0.20256667 0.21166667 0.22096667\n",
      " 0.23046667 0.24016667 0.25006667 0.26016667 0.27046667 0.28096667\n",
      " 0.29166667 0.30256667 0.31366667 0.32496667 0.33646667 0.34816667\n",
      " 0.36006667 0.37216667 0.38446667 0.39696667 0.40966667 0.42256667\n",
      " 0.43566667 0.44896667 0.46246667 0.47616667 0.49006667 0.50416667\n",
      " 0.51846667 0.53296667 0.54766667 0.56256667 0.57766667 0.59296667\n",
      " 0.60846667 0.62416667 0.64006667 0.65616667 0.67246667 0.68896667\n",
      " 0.70566667 0.72256667 0.73966667 0.75696667 0.77446667 0.79216667\n",
      " 0.81006667 0.82816667 0.84646667 0.86496667 0.88366667 0.90256667\n",
      " 0.92166667 0.94096667 0.96046667 0.98016667 1.00006667]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoiklEQVR4nO3dd1yVdf/H8deHISAgoOBEBLeW41bUhnlrZo4y07sypXGXaaVmO2daao5M08pxmzmampp7jzS34h7kxIGaIO4t8P39AfUjQznqgeuMz/Px8BHnXNc51/sb+Obye64hxhiUUko5Pw+rAyillLIPLXSllHIRWuhKKeUitNCVUspFaKErpZSL8LJqw6GhoSYyMtKqzSullFPauHHjSWNMWFbLLCv0yMhIYmNjrdq8Uko5JRE5dLNlOuWilFIuQgtdKaVchBa6Ukq5CMvm0LNy/fp1EhISuHLlitVRcpyvry/h4eF4e3tbHUUp5SIcqtATEhIIDAwkMjISEbE6To4xxpCcnExCQgJRUVFWx1FKuYhsp1xEZKyIJIrIjpssFxH5QkT2icg2Eal2p2GuXLlCgQIFXLrMAUSEAgUKuMW/RJRSuceWOfTxQKNbLG8MlMn40w4YeTeBXL3M/+Qu41RK5Z5sC90Y8xtw6harNAO+NenWAsEiUsReAZVSypV8+X0ntu1dmSPvbY859GLAkUyPEzKeO37jiiLSjvS9eCIiIuywaaWUch5TFg5j4dwvSb12lcplatv9/e1x2GJWcwdZ3jXDGDPaGBNtjIkOC8vyzFXLeXp6UrVq1b/+DBgwwOpISikXsOvAOn78oTMloqrS8YUvc2Qb9thDTwCKZ3ocDhyzw/taws/Pjy1bttz261JTU/H09LR/IKWU0zt/8TQDhz5NHp+89Hh7Gl5eeXJkO/Yo9JlARxGZCNQCzhpj/jHdcrvemv8WW/7Ycrdv8zdVC1dlaKOhdnu/yMhIXn75ZRYuXEjHjh159tln7fbeSinXYIzhk+EtOZ2cQOd3p1E4NDLHtpVtoYvIT0BdIFREEoBegHdG0FHAXKAJsA+4BLyUU2Fzw+XLl6latepfj7t27UrLli1vur6vry8rV+bMBxxKKec3YfrH7NyyiMebvsuD1Zrl6LayLXRjTKtslhugg90SZbDnnvTtuN0pl1uVvVLKva3fvoBfpvahXMWHaNtyYI5vT6/lcpf8/f2tjqCUckCJpxMYMjyGoKBC9HhzKh4eOf8Zm0Od+q+UUq4gNTWFPp835+qlc3TrsYTgwNw5qk8L/QY3zqE3atRID11USt2WL77twMF9sTz33KdULvtQrm1XC/0GqampNq978ODBnAuilHJK81dOYOmi0dS4rwUtm7yfq9vWOXSllLKTfYe3Mvqb9hQpVo4PXv0217eve+g2aN68OfHx8X97buDAgTRs2NCiREopR3Px8jk++bw5nh4e9Hh7Or4+uX/AhBa6DaZNm2Z1BKWUAzPG8MlXz5B8Ip533ppERNHyluTQKRellLpL46b1ZPvmBTR5/B3q1nzGshxa6EopdRdWbZnFjKn9qHBvXV59dpClWbTQlVLqDiUk7mPY8OcJzl+UHp2mIB7WVqoWejaaNGnCmTNnbrlO3bp1iY2N/cfzW7ZsYe7cuTmUTCllpSvXLtF78BOkXr9Kt3emkS+ggNWRtNBvxRjD7NmzCQ4OvqPXa6Er5ZqMMfQfGcPxI3G0efkrykVFWx0J0EL/h4MHD1KhQgXat29PtWrV8PT05OTJkwD06dOH8uXL06BBA1q1asVnn3321+smT55MzZo1KVu2LCtWrODatWv07NmTSZMmUbVqVSZNmmTVkJRSdvbD7H5sWjedRx59jSZ12lgd5y+Oe9jizp1w9qx93zMoCO65J9vVdu/ezbhx4xgxYgSRkZEAxMbGMnXqVDZv3kxKSgrVqlWjevXqf70mJSWF9evXM3fuXD7++GMWL15M7969iY2N5auvvrLvOJRSllm7fR6TJ/WiTLn76fh8ztx56E45bqFbqESJEtx3331/e27lypU0a9YMPz8/AJo2bfq35S1atACgevXqekkApVzU0cT9fP5Fa4KCC/Hh29Px9HSsCnWsNJnZsCedU7K6JG76Zd9vzsfHB0i/J2lKSkqO5FJKWefKtUt8PLgpKdev8FHXBYTkK2h1pH/QOXQb1a5dm1mzZnHlyhUuXLjAnDlzsn1NYGAg58+fz4V0SqmcZIyh34hWHD8SxysvD6dCyZpWR8qSFrqNatSowRNPPEGVKlVo0aIF0dHRBAUF3fI19erVY9euXfqhqFJObsL0j9m8fiaPNmxP4zovWx3npiS7qYScEh0dbW48djsuLo4KFSpYkscWFy5cICAggEuXLlGnTh1Gjx5NtWrV7vj9HH28SilYvvEXhgx5mrIVajOw29JcufPQrYjIRmNMlsdJOu4cugNq164du3bt4sqVK7z44ot3VeZKKccXf3QnX454kfxhEfR8e7rlZZ4dLfTb8OOPP1odQSmVSy5cOkPvz5oiBj58dyaB/iFWR8qWw82hWzUFlNvcZZxKOaO0tFQ++rwZyYkH6fDaWEoWr2R1JJs4VKH7+vqSnJzs8mVnjCE5ORlfX1+royilsvDld2+we+dvNG/Rnbo1n7Y6js0casolPDychIQEkpKSrI6S43x9fQkPD7c6hlLqBjOWjmTxgpFUr9GM/7bobXWc2+JQhe7t7U1UVJTVMZRSbmrz7uWMH/8WxUtUomuHnxARqyPdFoeaclFKKascOxnPwCEt8PMPptd7s/HJ42d1pNumha6UcnuXr16k16dNuHblIl3f+YVCBSKsjnRHtNCVUm7NGEPvL1rwR8LvtG0zgkplHrQ60h3TQldKubURP77Njs0Leezxdxz6tH5baKErpdzWzGX/Y/6cYVSp1phXW32W/QscnBa6UsotbYxbytixb1C8RCU+7DTV6Y5oyYpNhS4ijURkt4jsE5EuWSwPEpFZIrJVRHaKyEv2j6qUUvaRkLiPT4c+RUBAfj5+f65THtGSlWwLXUQ8geFAY6Ai0EpEKt6wWgdglzGmClAXGCwieeycVSml7tqFS2foObARKVev0P3dGYTld50T/GzZQ68J7DPGHDDGXAMmAs1uWMcAgZL+b5YA4BSgt+1RSjmU1NQUeg5pSvIf8XR8fRwVStWyOpJd2VLoxYAjmR4nZDyX2VdABeAYsB140xiTduMbiUg7EYkVkVh3OL1fKeVYBo15mb27VvLU072oV6ul1XHszpZCz+qTghuvntUQ2AIUBaoCX4lIvn+8yJjRxphoY0x0WFjYbUZVSqk7993Mvqxa/h0P1onh+Sd7Wh0nR9hS6AlA8UyPw0nfE8/sJeAXk24fEA+Ut09EpZS6O0vWTWTypJ6UKf8A77cdb3WcHGNLoW8AyohIVMYHnc8CM29Y5zBQH0BECgHlgAP2DKqUUndi+77VDB/5EgULl6L3e3Pw9HSoaxLaVbYjM8akiEhHYAHgCYw1xuwUkdcylo8C+gDjRWQ76VM0nY0xJ3Mwt1JKZetY0gH6DW6Gj48/H3eeT0DeYKsj5SibflUZY+YCc294blSmr48Bj9o3mlJK3bk/D0+8dvkCH3VbRLGCpayOlOP0TFGllMtJSb3Oh4MfI+n4ftq/NpZKZWtbHSlXaKErpVyKMYZ+I1uzL241LVv1pf59rayOlGu00JVSLmX0z53ZsHoKD9d/hdaPd7U6Tq7SQldKuYxpi79i9oxBVKr6KG++NCr7F7gYLXSllEtYuXkG48e/RYmoqvR6azoeHp5WR8p1WuhKKacXF7+Bz79oTf7QcPp0XuAyV0+8XVroSimndjRxP70/bYK3tw8ffzCfkHwFrY5kGS10pZTTOnshmR79G3Dt8kW6vTeDiKLufcURLXSllFO6eu0y3Qc+yqmkQ3Tq8C2Vyz5kdSTLaaErpZxOWloqHw17kkP7N/Hii5/z7xpPWR3JIWihK6WcijGGgV//lx2bF9K06Xu0aNDJ6kgOQwtdKeVUvp7chdXLv+fBOjG0ffZTq+M4FC10pZTTmLxwKLOmf0qlqo/yQbsJpN/1Uv1JC10p5RQWrfmB7ye8S1Sp6m574lB2tNCVUg5v3fb5jBj5MmFFSvFJl0Vue+JQdrTQlVIOLS5+A599/jQBQaH067aEQP8QqyM5LC10pZTDOnx8N70HNMbTy5s+XRZRMH/x7F/kxrTQlVIOKfF0Aj36P8L161fo8cFsIotVtDqSw9NCV0o5nHMXT9Htk4e5cCaR99+azL2lH7A6klPQQldKOZTLVy/SpX99kv7YT/vXx1KrcmOrIzkNLXSllMO4nnKN7oMakXBgC/99cSiP3B9jdSSnooWulHIIaWmp9BrajL27VvLMM71p3uANqyM5HS10pZTl0m/sHMP2TfN5/PF3eO7JD62O5JS00JVSlvt83KusWzWJuvVeol2rz6yO47S00JVSlhr103v8uvhrat3/FO+88o1en+UuaKErpSzz7Yw+zJk1mMr/akS3DhO1zO+SFrpSyhKT5g9m8qSelKv4EB/pxbbsQgtdKZXrZiwdyQ/fvU/J0jXo+8F8vL19rI7kErTQlVK5asGqbxk7tiPhJSrRv9sSfPPktTqSy9BCV0rlmiXrJjJyVBsKFS3LgG5LyesbaHUkl6KFrpTKFb9tmsaXw18gtGAkA7svI19AAasjuRybCl1EGonIbhHZJyJdbrJOXRHZIiI7RWS5fWMqpZzZ2m1zGTqsFSH5izGgxzJCggpZHckleWW3goh4AsOBBkACsEFEZhpjdmVaJxgYATQyxhwWkYI5lFcp5WQ27FzEoCFPEZgvjAE9lhEaUszqSC7Llj30msA+Y8wBY8w1YCLQ7IZ1WgO/GGMOAxhjEu0bUynljDbF/Ur/z54kr38w/Xr8SqHQElZHcmm2FHox4EimxwkZz2VWFggRkWUislFEXsjqjUSknYjEikhsUlLSnSVWSjmFLXt+o9+gpvj75aP/h8spVqi01ZFcni2FntWpW+aGx15AdeAxoCHwoYiU/ceLjBltjIk2xkSHhYXddlillHPYtnclfQc+Rh4ff/r2WEp44TJWR3IL2c6hk75HnvlGfuHAsSzWOWmMuQhcFJHfgCrAHrukVEo5je17V9F7QGN8vP3o230JJYpWsDqS27BlD30DUEZEokQkD/AsMPOGdWYAD4mIl4jkBWoBcfaNqpRydH+WeR5vX/p0X0JU+L1WR3Ir2e6hG2NSRKQjsADwBMYaY3aKyGsZy0cZY+JEZD6wDUgDxhhjduRkcKWUY/mzzL29fejbfSkli1eyOpLbEWNunA7PHdHR0SY2NtaSbSul7Gvb3pX0GdAEb28f+nRfQqnila2O5LJEZKMxJjqrZXqmqFLqrmzds+KvaZa+3ZdqmVvIlg9FlVIqS5t/X06/Tx/Hx8ePvt2X6py5xXQPXSl1R2J3LqbvwCb4+PrzSfdftcwdgO6hK6Vu27rt8/l0cIu/zgAtXqSc1ZEUWuhKqdu0asssBg99hoDAAvTvsUzPAHUgWuhKKZst2zCFYV/GEBRSmP49llEkLMrqSCoTLXSllE0WrvmeESNeIn9YBAN6LKNg/uLZv0jlKi10pVS2Zi3/mjFfv06hIqUZ0GMZ+YMKWx1JZUELXSl1S5MXfM53375LeMS99O+2lKDAUKsjqZvQQldK3dS3M/sweWJPokpH07/LYvzzBlkdSd2CFrpSKkv/m/QBs2cMolzFh+j7/jx8ffytjqSyoYWulPobYwyfj3uVXxd/TeWqDfno7Rl4e/tYHUvZQAtdKfWXtLRU+o2MYd2qSdS8/ym6tf8JT0+tCWeh3ymlFADXr1+l17An2b5pPvXqvczbr4xBJKsblilHpYWulOLSlfP0GNSYvXGraNr0Pdq1GmR1JHUHtNCVcnNnzifRdcAjHI3fRqtWn9C6aTerI6k7pIWulBs7fvIg3fvV51TiIdq+MpKmD79mdSR1F7TQlXJT+49so9eAhly6cJq3Ov1I3ZrPWB1J3SUtdKXc0NY9K/hkUFNIS6PHB3Oodk99qyMpO9BCV8rNLI+dyrCvnsPXLx+9Os+lXGR1qyMpO9FCV8qNTF8ygnHjOlEgLIK+3RZTNKyk1ZGUHWmhK+Umxkztzoyp/SgRVZW+XRYSHBhmdSRlZ1roSrm4tLRUPv36JVYt/457Kj3MR+/M1OuyuCgtdKVc2JVrl+j1eTN2bV3Mg3Vi+KDdBDw8PK2OpXKIFrpSLurUuRN0H/goCfHbaN68Ky8/3c/qSCqHaaEr5YLij+7ko4GNOXv6OG3bjOCJ+q9bHUnlAi10pVzMprhfGTCkOSY1lc7vTuP+qo9bHUnlEi10pVzInN++YcyY9gQE5OfD7nMoG1nN6kgqF2mhK+UCjDGMmdKNWdMGULR4RXp3nk/B/MWtjqVymRa6Uk7ueso1+o+KYcPqKdxbqT49356On2+A1bGUBbTQlXJiZ84n0fOzx4jfu4F69dvw1kv/08MS3ZiHLSuJSCMR2S0i+0Skyy3WqyEiqSLylP0iKqWyEn90J516VOfQ/o08/9wg3mkzRsvczWW7hy4insBwoAGQAGwQkZnGmF1ZrDcQWJATQZVS/2/VllkM/TIGjOGDd37hwWrNrI6kHIAte+g1gX3GmAPGmGvARCCrn543gKlAoh3zKaVuMHHuID79rDl58wbRv9cKLXP1F1vm0IsBRzI9TgBqZV5BRIoBzYGHgRo3eyMRaQe0A4iIiLjdrEq5tZTU6wz5pi0rlk2gZOkafPz+HL3AlvobWwo9q9t+mxseDwU6G2NSb3WXcGPMaGA0QHR09I3voZS6iTPnk+g1uCkH9qzjgdqteL/deLy88lgdSzkYWwo9Ach8QGs4cOyGdaKBiRllHgo0EZEUY8x0e4RUyp3tObSZvoOf4Nyp47Ru1Y9WTbtaHUk5KFsKfQNQRkSigKPAs0DrzCsYY6L+/FpExgOztcyVunuL1vzAyNHt8PL0ovN70/U0fnVL2Ra6MSZFRDqSfvSKJzDWGLNTRF7LWD4qhzMq5XaMMYyc9B7zZg6hUJHS9Hx3FhFFy1sdSzk4m04sMsbMBebe8FyWRW6M+e/dx1LKfV24fJY+XzzFrq2LqVT1UXp0mkJe30CrYyknoGeKKuVA9idsp8/gJzideIgnnnifV1oO5FYHGiiVmRa6Ug5i8bqfGDnqFTw8PHj7zYnUrfmM1ZGUk9FCV8piaWmpfPXjWyye+xUFi5SmxzsziCxW0epYyglpoStloVPnTtB7aHP2/76GatUfp0uHn/RKieqOaaErZZHNu5cz6IuWXDp7kpbP9CamWQ+dL1d3RQtdqVxmjOHHuQOZPPFD/P1D6N55DjUqNbQ6lnIBTlfol69c4NDxOEqFV8bb28fqOErdlouXz9F/VAxbN8ymVJma9HjrF0JDilkdS7kIm66H7kjmrRjH+91rsj9hu9VRlLotcfEbaN+tMts2zKZhow4M7rlKy1zZldPtoRctVAqAo4n7KB8VbXEapWwzddEXfP/9B+Tx9uXdtyfz7xp6Dxhlf05X6OGFygDwR1K8xUmUyt7Fy+f49H8vsGn9DIpHVqbH29MoGlbS6ljKRTldoRcpEImHhwcnTh60OopStxQXv4GBw57mVOIhHnn0NTo8N0wveatylNMVuqeXN0H5CpJ08rDVUZTKkjGGH+cMYPLPvfD18ee9t6dQp8Z/rI6l3IDTFTpAUEhhkk8dtTqGUv9w+lwi/Ue0Im7bUkqVqUnXN36mUGgJq2MpN+GUhZ4/fzhx+9ZaHUOpv1mzdQ5f/u9lLp47yRPNPqDN0/3w8PC0OpZyI05Z6AULRLBx02yupVwlj5cei66slZJ6neE/vMWSBSMJCi5Mz64LqH7PI1bHUm7IOQs9tARiIOHkAUoWrmB1HOXGDiTsYOBXz3Ls8E6q1WjKe+0mEOgfYnUs5aacstCLFEw/7CvhxF4tdGUJYww/zx/MxIk98PTwpG2bETxR/3WrYyk355SFXqxgaQCOndhvcRLljpJOH+XTkc/x+45lRJasRueOEwkvXMbqWEo5aaEXKo2AHouuct2C1d/xzbg3uHblon7wqRyOUxZ6Hr8Agvzzk5h8yOooyk2cv3iawd+0YePaaRQsXJK3O8/l3tIPWB1Lqb9xykIHCA4pQnJygtUxlBtYsWk6I8e8ysWzSTxc/xXaP/8FPnn8rI6l1D84baGHhBRl9/EdVsdQLuzi5XMMm/A6a377kZACxejRZZ5et1w5NKct9AL5wzn3+xJS01Lx1DlMZWertsxi5JhXOXvqOA/Wbk2nl0eR1zfQ6lhK3ZLTFnrB0BJ4XU/j+LmjhAdHWB1HuYgLl88ybNxrrF05keCQInR5bzoPVmtmdSylbOK0hV44LAqAIyf2aqEru1i6bhJjJnTiwplEHqwTwxsvDMc/b5DVsZSymdMW+p/Hoiec2Avl6lucRjmzU+dO8Pk3bdmyYRYFQovTvctcalVubHUspW6b8xZ6ofRC1xtdqDtljGHGslH8+GNXrl46T/1H2vFazBB8ffytjqbUHXHaQvfPV4B8PoEkntRj0dXtO3x8N0PGvMz+uNUULVaeju/PplLZ2lbHUuquOG2h4+lJSL5CnEw+YnUS5USup1xj/LRezJv9OQAtWnTnheYf4enpvH8VlPqTU/8Uh4QUZc9p3UNXttmwcxEjx75O0vH9lKtQm05tviaiaHmrYyllN05d6AXyh3PmyFqMMYiI1XGUgzp17gRfTujIxjVTyBuYn9df/ZrGddroz4xyOTYVuog0AoYBnsAYY8yAG5bHAJ0zHl4AXjfGbLVn0KyEhhbH4+o1ki4lUdC/YE5vTjmZtLRUfp4/hKnT+nLt0gUeqvsir8V8rtcrVy4r20IXEU9gONAASAA2iMhMY8yuTKvFA/82xpwWkcbAaKBWTgTOrHBYFF5pEJ+4h4JRWujq/23evZyR49pz/PAuipeoRPs2o/RiWsrl2bKHXhPYZ4w5ACAiE4FmwF+FboxZnWn9tUC4PUPeTLnwqgDsPBRLrSg9QkGlX6t8xPdvEbtmCv7+wbz44hD+0+BNxMPD6mhK5ThbCr0YkPlQkgRuvffdBpiX1QIRaQe0A4iIuPuzO8OLlifA25/dB2Pv+r2Uc0tJvc73Mz9h9qzPuH71Mg/WieH1mKEEBYZaHU2pXGNLoWf1yZHJckWReqQXepa7y8aY0aRPxxAdHZ3le9wOCQykZIHS7Dm4+W7fSjmxpesm8e1PnUlOPERU6WjavzSC8lE1rI6lVK6zpdATgOKZHocDx25cSUQqA2OAxsaYZPvEy4aHB8XCK7Li98lcS71GHs88ubJZ5Rh2H9rEqG87sS9uFcEhRWj/2jc0euglPXpFuS1bCn0DUEZEooCjwLNA68wriEgE8AvwvDFmj91T3kJUiSrk3fQTu5J2UbVw1dzctLJI0umj/O/Hd9mwejJeeXxo9mRnnn+yl950Qrm9bAvdGJMiIh2BBaQftjjWGLNTRF7LWD4K6AkUAEZk7B2lGGOicy72/ytXqibeabD1wBotdBd3+epFJkz/iIXzR5B67Qo173+Kdq0HE5Y/Vz6DV8rh2XQcujFmLjD3hudGZfr6FeAV+0azTWREZfy8fPn9wHp44HUrIqgclpqawtTFXzBten8unD1J2QoP0u65oZSLypV9BqWchlOfKQrgERRMVEhJ4g9tsTqKsjNjDAtXf8dPU3qRfOJg+kW02o7RG04odRNOX+h4elK8WAXWxM/W29G5kNVbZ/PtxK4cPbSDkALFaPfKKB6v21aPJ1fqFpy/0EmfdvHZMZXdybupGFbR6jjqLmyMW8qESV2J37OegMD8tGr1CU83ehdvbx+roynl8Fyi0MuWrIlvCmw5tE4L3Ult3bOCbyf3YM/O38ibNx9PPtmF1k90x883wOpoSjkNlyj0kpFV8fHMQ9z+dRD9ktVx1G3YsX8NEyZ35/dtv+Lr60/jxp14vnkvAgPyWx1NKafjEoXuFVKAqOAoDugZo05j654VfD+lJ7/vWEYeHz8aNurAc817ERwYZnU0pZyWSxQ63t4UL1yO9QlLSDNpeIh+cOaoNuxcxE/TerN310p8ff1p2KgDMc0+JCSokNXRlHJ6rlHoQImISnjumcn+U/spU6CM1XFUJsYYfts0jckz+nFo30Z8/QJo3KQTMc0+1ItnKWVHLlPoFcs8gP88+HXvIi10B5GamsKcleOYMWcwiQm7CQjIz5PNu9LysQ8IyBtsdTylXI7LFHrpyGoUCSjM4s1TaHdfe6vjuLVLVy8wdeEwFiwcwdnkY4TkL0rr1v1p3uANfH38rY6nlMtymUKX0FBqFb+fL3fN4uyVswT5Blkdye38kXyIn+cMYuVv33H50jmKhpfn2bY9aVynDZ6eLvOjppTDcp2/ZV5eRFdpTOiUaczdO5dWlVpZnchtbNu3ismzBrJ90zxSU1MoV6E2/2namfuqPKaXslUqF7lOoQOVqjSgyPRg5m/6WQs9h11Puca8VROYt3A4CfFb8fb2odZ9T/H0450pXaKq1fGUcksuVegeRYpSK7wW47Ys5ErMFXy9fK2O5HKOnzzILwuGsmLF91w8l0xwcGGefLILLRq+pYceKmUxlyp08uSh+j0NmDJvAUsOLOGxso9ZncglpJk0VmyaztzFI4nbvhSTlkapMjVpHPMZjzzwnM6PK+UgXO5vYrVqj1FwzofM3fyzFvpdOnHqCNMXf8WKFd9zNvkYvr4B1H4ohhaN3tFpFaUckMsVuk94CaKLVmfqptmkPq2X071dKanXWbp+Egt/HcOeXSswaWkUL1GJFk0/oMm/X9HDDpVyYC5X6Pj5Ub1CfWYtX8myg8uoX7K+1YmcQlz8BmYvGUnshulcOn+avP5B1KnzAk0bdNA7AynlJFyv0IEHav6HqOXDGLSktxb6LRw/eZC5y75m9drJJB7bi4eHB+UqPES9h17gkftj9BrkSjkZlyx0v4iSPHXP03TZ8jW/xv9Kvah6VkdyGKfPJ7FgxXhWrf2ZQ/s3YoyhaLHy/OepD3msbju94bJSTswlCx1/fxr+uw1T4n6hz4Ju1H11tVuf4HL2QjKLVn/P6vVT2P/7GtLSUgkOKUL9Bq/SuO4rlI2sbnVEpZQduGahAz73VKblPc/QbetIFh1YxKOlHrU6Uq46efY4S1b/wNrY6cTvWUdqagr58oXx4EOteaT2i/yr4sNu/UtOKVfksoWOnx+P1m/HlF1T6DevGw06NHD5Aos/toula35k4+Y5HI3fRppJIyioIA/UbkW9B2Kofs8jeOhRP0q5LNctdMC7fEVaVm5Fj21f8N2273ihygtWR7KrlNTrrN+xgNWx09m+YzGnThwCIKxgJA8/0pa697eicrk6Lv+LTCmVzqULnTx5qP9IO5bvX0LXiW0pV6ActcJrWZ3qrhw+sYdVG6ezedsC9u9ey7Wrl/D08CQiqir/fiCGeve3Iir8XqtjKqUs4NqFDniXKUeXhn35Y+brxIxvyq9vbKR4UHGrY9ks+ewfrN02l807FrH795WcSUoAIDCwAFWqNKR61cY8FN2CfAEFLE6qlLKayxc6Xl4E1W3IRxf68Ob8N3l6wmMsenUVgT6BVifLUtKZY6zfPp/tu5axZ+8ako7tA8Dby4cSpf5FvdovcF+1plQoWUunUpRSf+P6hQ7g50eJRi358PwJuvzWk/uGV2V8y4nUKFbD0lipaansPrSJLb//yu97VnPwwCZOJx0BwMvLm/CIe4lu3IlqlR+lWsX65PHWq0cqpW7OPQodIDCQ6k++zhAPbz5dPYhnh9xPm2Yf80Htznh55Pz/hjSTxoGjO9mxbxV7928g/uBmjifEkXL1CgB5fQOJiKzCg7We5l/3PkLV8nXJk8cvx3MppVyHGGMs2XB0dLSJjY3N/Q1fvMi5DSsZNb8vM5JWcq5EIZrXbssr1dsSERRx129vjOHo6UPsPbiJ+ITtHD6yg6NH40g8vo+0q1cB8PL0omChUkRGVqVs6VpULleH0iX+hXh43PX2lVKuTUQ2GmOyvMCS+xV6BnP4MOsXjWferhmsTdzEH4FCwYgKVIisQaWS9xMZWorQvKGE5g0lj2cejDEYDBevXeTM+ZOcO5/EyVMJJCUf5lTyUU4mH+Zk0mHOnjoKV67+tR0/nwAKFylNsWIViCxRhYql76dcyRq6962UuiN3Xegi0ggYBngCY4wxA25YLhnLmwCXgP8aYzbd6j2tLnQA0tIgMZFju2P5dc1P7E7cxd7kvVy4fpEUD0jxgNSMzx09DXgY8EpL/29mfl5+BAcVJCQ0nLCwSIoUKk1k8UqULvEvCoZF6p63UspublXo2U4ei4gnMBxoACQAG0RkpjFmV6bVGgNlMv7UAkZm/NexeXhA4cIULfw4MQ81gYsXMefPcyRhF4lnj3L2QjJnLiaTkpYCnh4gHvj6BRIQkJ+AgPwUCClGobBIAoMLgqeegamUspYtnwbWBPYZYw4AiMhEoBmQudCbAd+a9N39tSISLCJFjDHH7Z44p3h4QGAgEhhIRNGi3P1sulJK5S5b5gKKAUcyPU7IeO5210FE2olIrIjEJiUl3W5WpZRSt2BLoWd19sqNE++2rIMxZrQxJtoYEx0WFmZLPqWUUjaypdATgMznyocDx+5gHaWUUjnIlkLfAJQRkSgRyQM8C8y8YZ2ZwAuS7j7grFPNnyullAvI9kNRY0yKiHQEFpB+2OJYY8xOEXktY/koYC7phyzuI/2wxZdyLrJSSqms2HTOuzFmLumlnfm5UZm+NkAH+0ZTSil1O/SMF6WUchFa6Eop5SIsu5aLiCQBh+7w5aHASTvGcQY6ZvegY3YPdzPmEsaYLI/7tqzQ74aIxN7sWgauSsfsHnTM7iGnxqxTLkop5SK00JVSykU4a6GPtjqABXTM7kHH7B5yZMxOOYeulFLqn5x1D10ppdQNtNCVUspFOHShi0gjEdktIvtEpEsWy0VEvshYvk1EqlmR055sGHNMxli3ichqEaliRU57ym7MmdarISKpIvJUbubLCbaMWUTqisgWEdkpIstzO6O92fCzHSQis0Rka8aYnfqaUCIyVkQSRWTHTZbbv7+MMQ75h/QLge0HSgJ5gK1AxRvWaQLMI/167PcB66zOnQtjfgAIyfi6sTuMOdN6S0m/ptBTVufOhe9zMOl3BYvIeFzQ6ty5MOZuwMCMr8OAU0Aeq7PfxZjrANWAHTdZbvf+cuQ99L9ufWeMuQb8eeu7zP669Z0xZi0QLCJFcjuoHWU7ZmPMamPM6YyHa0m/9rwzs+X7DPAGMBVIzM1wOcSWMbcGfjHGHAYwxjj7uG0ZswECM246H0B6oafkbkz7Mcb8RvoYbsbu/eXIhW63W985kdsdTxvSf8M7s2zHLCLFgObAKFyDLd/nskCIiCwTkY0i8kKupcsZtoz5K6AC6TfH2Q68aYxJy514lrB7f9l0+VyL2O3Wd07E5vGISD3SC712jibKebaMeSjQ2RiTmr7z5vRsGbMXUB2oD/gBa0RkrTFmT06HyyG2jLkhsAV4GCgFLBKRFcaYczmczSp27y9HLnR3vPWdTeMRkcrAGKCxMSY5l7LlFFvGHA1MzCjzUKCJiKQYY6bnSkL7s/Vn+6Qx5iJwUUR+A6oAzlrotoz5JWCASZ9g3ici8UB5YH3uRMx1du8vR55yccdb32U7ZhGJAH4BnnfivbXMsh2zMSbKGBNpjIkEpgDtnbjMwbaf7RnAQyLiJSJ5gVpAXC7ntCdbxnyY9H+RICKFgHLAgVxNmbvs3l8Ou4du3PDWdzaOuSdQABiRsceaYpz4SnU2jtml2DJmY0yciMwHtgFpwBhjTJaHvzkDG7/PfYDxIrKd9OmIzsYYp72sroj8BNQFQkUkAegFeEPO9Zee+q+UUi7CkadclFJK3QYtdKWUchFa6Eop5SK00JVSykVooSullIvQQldKKRehha6UUi7i/wA0QNcgvZKHNwAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import pylab as pl\n",
    "\n",
    "h = 0.01\n",
    "n = int(1 / h)\n",
    "f = lambda x, y: -50 * y + 50 * x ** 2 + 2 * x\n",
    "\n",
    "#f:函数体\n",
    "# x0初值\n",
    "# y0初值\n",
    "# h 间距\n",
    "# n不用你管\n",
    "def E_r(f, x0, y0, h, N):\n",
    "    x = np.zeros(N + 1)\n",
    "    y = np.zeros(N + 1)\n",
    "    x[0] = x0\n",
    "    y[0] = y0\n",
    "    for i in range(N):\n",
    "        x[i + 1] = x[i] + h\n",
    "        y_b = y[i] + h * f(x[i], y[i])\n",
    "        y[i + 1] = y[i] + h / 2 * (f(x[i], y[i]) + f(x[i + 1], y_b))\n",
    "    return x, y\n",
    "\n",
    "\n",
    "[xx_r, yy_r] = E_r(f, 0, 1 / 3, h, n)\n",
    "print(xx_r, yy_r)\n",
    "\n",
    "xxx = np.arange(0, 1 + h, h)\n",
    "yyy = []\n",
    "for i in xxx:\n",
    "    yyy.append((math.e ** (-50 * i) / 3) + i ** 2)\n",
    "    # yyy.append(xxx)\n",
    "\n",
    "pl.plot(xx_r, yy_r, label=\"E_r\", color=\"green\", linestyle=\"-\")\n",
    "\n",
    "pl.plot(xxx, yyy, color=\"red\", label=\"right\", alpha=0.3)\n",
    "pl.legend()\n",
    "pl.show()\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 龙格库塔法"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1  0.11 0.12 0.13\n",
      " 0.14 0.15 0.16 0.17 0.18 0.19 0.2  0.21 0.22 0.23 0.24 0.25 0.26 0.27\n",
      " 0.28 0.29 0.3  0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4  0.41\n",
      " 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5  0.51 0.52 0.53 0.54 0.55\n",
      " 0.56 0.57 0.58 0.59 0.6  0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69\n",
      " 0.7  0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8  0.81 0.82 0.83\n",
      " 0.84 0.85 0.86 0.87 0.88 0.89 0.9  0.91 0.92 0.93 0.94 0.95 0.96 0.97\n",
      " 0.98 0.99 1.  ] [0.33333333 0.2023572  0.12312403 0.07536562 0.04678383 0.02991649\n",
      " 0.02023579 0.01499437 0.01252523 0.01181687 0.01225555 0.01346886\n",
      " 0.01523085 0.01740439 0.01990631 0.02268612 0.02571319 0.02896894\n",
      " 0.03244209 0.0361258  0.04001592 0.04410992 0.04840628 0.05290407\n",
      " 0.05760273 0.06250192 0.06760142 0.07290112 0.07840094 0.08410083\n",
      " 0.09000077 0.09610072 0.1024007  0.10890069 0.11560068 0.12250067\n",
      " 0.12960067 0.13690067 0.14440066 0.15210066 0.16000066 0.16810066\n",
      " 0.17640066 0.18490066 0.19360066 0.20250066 0.21160066 0.22090066\n",
      " 0.23040066 0.24010066 0.25000066 0.26010066 0.27040066 0.28090066\n",
      " 0.29160066 0.30250066 0.31360066 0.32490066 0.33640066 0.34810066\n",
      " 0.36000066 0.37210066 0.38440066 0.39690066 0.40960066 0.42250066\n",
      " 0.43560066 0.44890066 0.46240066 0.47610066 0.49000066 0.50410066\n",
      " 0.51840066 0.53290066 0.54760066 0.56250066 0.57760066 0.59290066\n",
      " 0.60840066 0.62410066 0.64000066 0.65610066 0.67240066 0.68890066\n",
      " 0.70560066 0.72250066 0.73960066 0.75690066 0.77440066 0.79210066\n",
      " 0.81000066 0.82810066 0.84640066 0.86490066 0.88360066 0.90250066\n",
      " 0.92160066 0.94090066 0.96040066 0.98010066 1.00000066]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoO0lEQVR4nO3deXhU5d3/8fc3kz1kISRAQggECDsESNh3hMoqiCs+dQNKcelirRVra21t+2hFa1UoP7SotSpocQFJFVH2PWEPiAQSIGxJCBCSELLdvz+S+sQYyEAmc2Ym39d1ccnknMz53Mj14eSec+4jxhiUUkq5Py+rAyillHIMLXSllPIQWuhKKeUhtNCVUspDaKErpZSH8LbqwBEREaZt27ZWHV4ppdxSampqrjEmsrZtlhV627ZtSUlJserwSinllkTk6JW26ZSLUkp5CC10pZTyEFroSinlISybQ69NaWkpWVlZFBcXWx2lwfn7+xMTE4OPj4/VUZRSHsKlCj0rK4vg4GDatm2LiFgdp8EYYzh79ixZWVnExcVZHUcp5SHqnHIRkUUiki0i+66wXUTkZRFJF5E9ItLnesMUFxfTrFkzjy5zABGhWbNmjeInEaWU89gzh/4mMPYq28cB8VW/ZgF/r08gTy/z/2os41RKOU+dhW6MWQfkXWWXycA/TaUtQJiIRDkqoFJKeZKSnfsgP79B3tsRV7m0Ao5Xe51V9bXvEZFZIpIiIik5OTkOOLRSSrmP/X9bSfLov3J+y/4GeX9HFHptcwe1PjXDGLPQGJNkjEmKjKz1zlXL2Ww2evXqRffu3Zk0aRLnz5+/4r5r1qxh4sSJzgunlHJbOVuPsPfxd/Bp3ZKQkX0b5BiOKPQsoHW11zHASQe8ryUCAgLYtWsX+/btIzw8nHnz5lkdSSnl5i6fK2TDbS/hFeDLsI9+jpePrUGO44jLFpcBD4vIYqA/cMEYc8oB78ubI9783te63d6Nvg/2pbSolHfGv/O97b3u60Wv+3pRlFvE+7e+/51t962575qOP3DgQPbs2WPXvtu3b2fWrFksXbqUdu3aXdNxlFKebfej/6Q46yyDP/olwXENNztRZ6GLyHvACCBCRLKA3wE+AMaYBUAyMB5IB4qA+xsqrDOVl5fz5ZdfMmPGjDr33bRpEz/5yU/45JNPiI2NdUI6pZTbyMggcUIULYZMJ3Zy7wY9lFj1kOikpCRTc7XFAwcO0KVLF0vy/JfNZqNHjx5kZmaSmJjIypUrsdlq//FozZo1zJgxg4CAAFauXEl0dPQ1HcsVxquUajg5m74hNCsN37gYSEoCB1yuLCKpxpik2rbpWi41/HcO/ejRo5SUlNQ5hx4VFYW/vz87d+50UkKllDsoPHGOtTe9wLrnt0GvXg4p87pooV9BaGgoL7/8MnPnzqW0tPSK+4WFhbFixQp+/etfs2bNGucFVEq5rIqyctZO+Svl+Zfo8dfp4KQ1m7TQr6J3794kJCSwePHiq+7XokULli9fzkMPPcTWrVudlE4p5aq2PvRPLqSk0+1Pd9JiSLzTjutSi3O5goKCgu+8Xr58+RX3HTFiBCNGjAAgNjaWtLS0hoymlHID6f/cQObCL4ieOoDuj4136rG10JVSylEKC4n2OkObyQn0e3u20w+vhW6HvXv3cvfdd3/na35+fjq9opT6VmlBMbZt2wiMCGLQez+FAF+nZ9BCt0OPHj3YtWuX1TGUUi7KGMO6W17Blnua4cmPIwEBluTQD0WVUqqedv32Q7JX7iRoRCLSorllObTQlVKqHo4t28nX//shESO6kzR3mqVZtNCVUuo65adns/We+fi3asawpT+z/ME1Wuh1GD9+/FWX0IXKyxdrLmMAsGvXLpKTkxsomVLKUhUVlK7bhH+ID0M+fAS/8CCrE2mhX40xhk8//ZSwsLDr+n4tdKU82N69NIu0MWHzb4hMco2HvWuh15CZmUmXLl148MEH6dOnDzabjdzcXACeeeYZOnfuzJgxY5g2bRpz58799vs++OAD+vXrR8eOHVm/fj0lJSU89dRTLFmyhF69erFkyRKrhqSUcrA9f1rOzqc/oSKuPV6trm1RvobkupctpqXBhQuOfc/QUOjWrc7dDh48yBtvvMH8+fNp27YtACkpKSxdupSdO3dSVlZGnz59SExM/PZ7ysrK2LZtG8nJyfz+979n1apV/OEPfyAlJYVXX33VseNQSlkm6z97SPvdEpr2i4fOnayO8x2uW+gWatOmDQMGDPjO1zZs2MDkyZMJqLq+dNKkSd/ZPnXqVAASExPJzMx0Sk6llHPlH85m812v4N+iKcM/fgQv74Z58tD1ct1Ct+NMuqEEBX3/w4261o338/MDKtdTLysra5BcSinrlBWVsHbSXCqKSxny+RMENA+xOtL36By6nYYMGcLy5cspLi6moKCAFStW1Pk9wcHBXLx40QnplFINLef9L7l0+BS9500nsp9rPmZSC91Offv25aabbiIhIYGpU6eSlJREaGjoVb9n5MiR7N+/Xz8UVcrdZWQQ1ayMCV8+Qsfpw6xOc0X6CLprUFBQQJMmTSgqKmLYsGEsXLiQPn36XPf7ufp4lVJw9KNUytdtpN1dAx32GLn6uNoj6Fx3Dt0FzZo1i/3791NcXMy9995brzJXSrm+c2kn2HbPfHwigon94yy8LS7zumihX4N3333X6ghKKScpOV/EuolzMRiGLXsU7yB/qyPVyeXm0K2aAnK2xjJOpdyRKa9gzeS/UnQ0m76LZhPeo7XVkeziUoXu7+/P2bNnPb7sjDGcPXsWf3/X/xdfqcbo+JurOLsujc5PTiXutn5Wx7GbS025xMTEkJWVRU5OjtVRGpy/vz8xMTFWx1BK1XTyJLHNL+P71r20uPsHVqe5Ji5V6D4+PsTFucYiN0qpxuf02oN4p24lon97Wk4YaPkVLdfKpQpdKaWscjEjhw23vIi3vw83HbwNLy+XmpG2i/slVkopBystKGbN+OcpLyhm4JKf4RVkzTNB60sLXSnVqP33Ac8FX2fR5+8zaTE43upI100LXSnVqB16aQXZK3fS4RcTib9/qNVx6kXn0JVSjdfJk3SIK0f+PIUOc26zOk29aaErpRqlM+sPEnxkD4Edoom/yf2uaKmNXVMuIjJWRA6KSLqIzKlle6iILBeR3SKSJiL3Oz6qUko5Rn56NuunvMjav2zFJCaCG17RUps6RyEiNmAeMA7oCkwTka41dnsI2G+MSQBGAC+IiK+DsyqlVL2VnC9izbjnKL90maRFDyIedMe2Pf8s9QPSjTFHjDElwGJgco19DBAsIgI0AfIAfWyPUsqlVJSVs3rSCxQePk2/RQ8Q2d81H1Rxvewp9FbA8Wqvs6q+Vt2rQBfgJLAX+JkxpqLmG4nILBFJEZGUxnB7v1LKtez61bvkbThAl6dvI+7O/lbHcTh7Cr22Twpqrp51I7ALiAZ6Aa+KyPceuGeMWWiMSTLGJEVGRl5jVKWUqofMTLr3C6Tn76bQ66kpVqdpEPZc5ZIFVF87MobKM/Hq7geeNZXLJKaLSAbQGdjmkJRKKVUPp1fuIeLiEXzbtabbHX2tjtNg7DlD3w7Ei0hc1QeddwLLauxzDLgBQERaAJ2AI44MqpRS1+PMxkOsnfwiWxbthz59POLyxCup8wzdGFMmIg8DnwM2YJExJk1EZldtXwA8A7wpInupnKJ53BiT24C5lVKqTheP5LBhygvYgvzo+fJM8PbsW2/sGp0xJhlIrvG1BdV+fxJwr4WDlVIereR8EavHPkfZxUuMWPUbQto3tzpSg/OMq+mVUqo6Y9g07RUKD5+i7xsP0GKI+y64dS08++cPpVTjtHcviVNac2JcL9pNG2B1GqfRQldKeZSTi9cRFXie4GG96dyli9VxnEqnXJRSHuPAvFWsnbaAA1suQOfOVsdxOi10pZRHOLZsJ7t//hYhCXF0/M0dHn154pVooSul3F7u9gy23PUK/tHhjFr5ON6BjXNtQC10pZRbK88vZP3kuXj5ejNi5RwCmn9v1ZFGQz8UVUq5r9JSbDtTGPBQIjJsCGGdoqxOZCktdKWUWyq7VEL268uJjvMjatYk0AX/dMpFKeV+THkFa2/+G+t+9iHnQtpomVfRQldKuRVjDBvvf43sz3fS4dFJNB3Ww+pILkMLXSnlVnY88QHH315L6x8OJ/Evd1gdx6VooSul3Mbp5BS+ee5jIsckMPjNHyGN8Frzq9EPRZVS7iE7m5YVp+n3h/G0+eXtiE3PR2vSQldKubwTK/cRcHA34Qmtaf/EWI9f1/x66Z+KUsql5W7PYOMtLxIYFcaEPbcgWuZXpD+zKKVc1vmDp1gz7lm8fLwZ/PGjiL+/1ZFcmha6UsolFZ44x+rRf6b8UgnDPv0VTbu2sjqSy9NCV0q5nrIy0n7xDy5nX2DwB4/QfFAHqxO5BZ2MUkq5lvJy2LaNxDvjiZ3xA1r+oKfVidyGnqErpVxG+eVStt4/n8sZJ7H1S9Iyv0Za6Eopl2DKK1g79WWOvL2ZI9lB0ErnzK+VFrpSynLGGNbfvYAzyal0eGQiXR4ZZ3Ukt6SFrpSy3JbZb3LivQ3E3jeSpBemWR3HbWmhK6UsdXnnfnKXbybqlgEMWjRT12epB73KRSllGZORgV/WYUb/6178hg3QMq8nPUNXSlki7YXP2DL9NSoimhMwYiBe3jarI7k9LXSllNN9/fev2PPY25wv9KEioTd4aRU5gv4pKqWcKv2fG9n18CJCerThhlVz8A70tTqSx9BCV0o5zZElW0mZsYCgjtGMXv0kviEBVkfyKFroSinnOHuWoIw0Qrq0YvSa3+AXHmR1Io9jV6GLyFgROSgi6SIy5wr7jBCRXSKSJiJrHRtTKeXOCr85Dlu30mJAO8Zt/z0BLUKsjuSR6ix0EbEB84BxQFdgmoh0rbFPGDAfuMkY0w24zfFRlVLu6MQX+1jR+ykOrzsBAwcifn5WR/JY9pyh9wPSjTFHjDElwGJgco197gI+NMYcAzDGZDs2plLKHZ1afYANU17AOzSQyJmTQcu8QdlT6K2A49VeZ1V9rbqOQFMRWSMiqSJyT21vJCKzRCRFRFJycnKuL7FSyi2cXneQ9ZOexyc4kBvW/JaQDs2tjuTx7Cn02m7dMjVeewOJwATgRuC3ItLxe99kzEJjTJIxJikyMvKawyql3ENR5mnWT/wLtiA/Rq1+ktCOLa2O1CjYc+t/FtC62usY4GQt++QaYwqBQhFZByQA3zgkpVLKfeTnE3hwNwkz+hA5cwphXaKtTtRo2HOGvh2IF5E4EfEF7gSW1djnE2CoiHiLSCDQHzjg2KhKKVeXvSmdnH8sAy8vOv7xPpp20zXNnanOM3RjTJmIPAx8DtiARcaYNBGZXbV9gTHmgIh8BuwBKoDXjTH7GjK4Usq1ZG9KZ+24Z/EN9WfSvv/FK0ivM3c2MabmdLhzJCUlmZSUFEuOrZRyrDMbD7F23LN4+Xoz8ssnaZYQa3UkjyUiqcaYpNq26fK5Sql6ObPhEOvGPYvN34eRXz5JeM/WdX+TahBa6Eqp63fhAhl/fgevAF9GfvlrwntomVtJC10pdV3MuXPI1q30/elAunZIIKSDXppoNV2cSyl1zbI+38vnA5/mUkEZtuHDtMxdhBa6UuqaHF++i41TXuDSxQpKe/WDAF0C11VooSul7Jb57+1svPVF/CJDGb3+Kb2d38XoHLpSyi5H/72NLdNewb91BGPW/pag1uFWR1I1aKErpep26hTNi48RNbwjff/1UwJbhlqdSNVCC10pdVVHF28kxj+XgLgohv9nMvj4WB1JXYHOoSulrmjPn5azado89n1+CgYM0DJ3cXqGrpSqVeqc9/nmuY8JH9yFbi9MB5vN6kiqDlroSqnvMMawZfabZC78gsgxCYxc9gts/npm7g600JVS/8cYLq7cRNbbq4meOoChSx7Cy1vPzN2FFrpSCgBTVo7s2U1ISR5jls4mdOxARGp7YJlyVVroSilKC4r5auzztOnkT+c5UwiLj7c6kroOepWLUo1cce5FVg7+A3kbD1DSrjNombstPUNXqhEryMzly9F/5lJGNr3nz6TzA6OsjqTqQQtdqUaqNPc8Xwx+mpKzF+n/3k+Ju72f1ZFUPWmhK9UYXbiAz45tdL+tM0ETRxI9upvViZQDaKEr1chkfpiC34E9RPWPJf6P90GTJlZHUg6iH4oq1Yh8Pf9LNt/+EjvePYAZPFjL3MPoGbpSjcTO3yzl6z8tJTQhjlFfzEH0wRQeRwtdKQ9nyivYOP01jv9zLREjezDy00fxDvS1OpZqAFroSnmyigrYuQPvE8eIuWsog9+apbfyezAtdKU81KUz+ZSu20SIfyn95t+HxHfQW/k9nBa6Uh7o3P4TrBn7HL6UMW7DE3jFtrY6knICLXSlPMzJr/az8Za/YsrKSXr3p1rmjYgWulIe5NAb60md/Rq+TYMZ+sUviUyKszqSciItdKU8hEk/zNHn/01Q+5aM/OxxmsQ2szqScjItdKXcXEVJGWU7duObc5IhL9+GrX8SPsH+VsdSFtBCV8qNFedeZM2EF/AtvsDId3+Ef9cuoFeyNFp23fovImNF5KCIpIvInKvs11dEykXkVsdFVErV5tz+E3yW9FvOp6QT+T+jkW5dtcwbuTrP0EXEBswDxgBZwHYRWWaM2V/Lfs8BnzdEUKXU/8lK3sPmaS9jyisYuPQR2kxJtDqScgH2nKH3A9KNMUeMMSXAYmByLfv9BFgKZDswn1KqhvIjmeyYMQ/vkABu2PC0lrn6lj1z6K2A49VeZwH9q+8gIq2Am4FRQN8rvZGIzAJmAcTGxl5rVqUatYqycth/ANvRDIb+dQr+owYT0DzE6ljKhdhT6LVNypkar18CHjfGlF/t1mJjzEJgIUBSUlLN91BKXUFx7kXW3PQikRGGxD/dQtNu3cBLV79W32VPoWcB1W81iwFO1tgnCVhcVeYRwHgRKTPGfOyIkEo1Zmd3H2PdpLlcPpFH9DN3QI8eVkdSLsqeQt8OxItIHHACuBO4q/oOxphvb0cTkTeBT7XMlaq/jMVb2TZjAWLzYtBHjxJ7U2+rIykXVmehG2PKRORhKq9esQGLjDFpIjK7avuCBs6oVKN0afteUu6fh3+rCIYv/yVhXaKtjqRcnBhjzVR2UlKSSUlJseTYSrmy8ksl2PbvhZMnOXVGaHb7DfiG6NOFVCURSTXGJNW2TT9VUcqF5O09zvJuj3P4zfXQpQtRMydqmSu76a3/SrmIjMVb2T7z/4GADBkIHTpYHUm5GS10pSxmyivY/ot3OfxyMoHtWjJs2aM07dbK6ljKDWmhK2Wl0lJOL0rm8MvJtBifyNDFD+lKieq6aaErZZHio2fwT99HVCsbw5c+TNTNA/WZn6pe9ENRpSyw9y8rWNblV+SknYFBg4ieOkjLXNWb+52hl5dDUREEBemtz8rtlF4sZsMP/87pZdsJ6xtP4C3jILyp1bGUh3C7Rry48xDf/Oo1CjJ0UUflXnJTM0nu+QSnl20n7oEbuXHTUwS10jJXjuN2hX4hu4TUBSlkpxy1OopS9jtxglN/e5+ScxcZsOTnDJh/L17eNqtTKQ/jdlMuTdpHAlCQkWtxEqXqVnqxmPzPNtHMv5BuMwcS98x0mrSJsDqW8lBuV+jBbSNAhKJjWujKteVuz2DD7X/DnM9n4leP4NOrO030g0/VgNyu0G1+Pvg2C+ZSVp7VUZSqlTGGvc8ls/93S7AF+JH02mx8euuSt6rhuV2hA/i1bMqlE1royvWUXbzE6qkvk7tqN2F94xmy5GGC4yKtjqUaCbcs9OHPjsWv/JLVMZT6rrw8bDt2EHDpHPGP3kTic7cjNre77kC5Mbcs9OC4SMjIAGNA5ySVxSpKy0n95Xt07m4juF1zBi/7FRIebnUs1Qi55elDztECdr2xg0unL1gdRTVy5/afILnPb0h/OZlDu4tg+HAtc2UZtyz0C9klHFh6gHMHaj7aVCnnMMaQ9uJnrEx6ksIjp+n1ynT6vDoDvN3yh17lIdzyb1+Tqg+ZCo7kwCiLw6jGp6SEr3+/hD1//g8hveIYvPhhwjpFWZ1KKfcs9JD4FgAUZOZYnEQ1NpePHMcv/QAd+gRT+vStdH9yst7xqVyGWxZ6QMtQvPx8KDp21uooqpG4fK6QzT/6BwVb0hj7j9vxGTOKnreEWB1Lqe9wy0IXEfyah3H5zHmro6hG4NgnO0n58etcPnOONvePgmFDIcDX6lhKfY9bFjrAhLduxdtmrI6hPFhZ4WW2zH6T4/9ai390M4YmzyFmXE+rYyl1RW5b6D7hIXDmjNUxlKc6exavHTsp2r6PVtOGMGDB/fiGBFidSqmrcttCP74rl6w319B36Ai89cdf5SAl54tIffQdeo8Mwz8qnFFf/Qbv6BZWx1LKLm55HTpAfm4pmWsyuXhEr3RRjpHxwTY+7fwYmYtWk3lUYMQILXPlVtz2DD0ornJN6Yvp2TTt1sriNMqdXcrOZ8uP/sHpZdvxbx3B0BVziBmvc+XK/bhtoQe3bw5AQYaeoat6OHmS3TNf5/TnB2kz8wb6/e1uvAN1Ck+5J7ct9JD/FvpRfdCFunYXDp1G9u0jxPsSPX/Uj7Zz7qDlsE5Wx1KqXty20H2a+BPQMhSKiqyOotxIRWk5u37/EYfmLieyWySj3n+AwHbtCNRVO5UHcNtCB5iyZBr4+VkdQ7mJk1/uJ+WBRRQeOkn4oM70eX0mtI+2OpZSDuPWhU5AABQWWp1CubrSUjLmJbPlkSX4hAeTuPDHxM8chuhZufIwdl22KCJjReSgiKSLyJxatv+PiOyp+rVJRBIcH/X7Dq46xlePrsAYvWNUfZ8pr+Di9v2wejUxMRW0+/FoJh16kY4/Gq5lrjxSnWfoImID5gFjgCxgu4gsM8bsr7ZbBjDcGHNORMYBC4H+DRG4usKLFZzZcYLLuQX4RwY39OGUGzm97iApD76Byctj3L/uwucHA+h/qy6mpTybPVMu/YB0Y8wRABFZDEwGvi10Y8ymavtvAWIcGfJKgtpUXot+4eBpLXQFQNGp86T8/B1OfLAJ75AAuv72VmzDh4E+21M1AvYUeivgeLXXWVz97HsG8J/aNojILGAWQGxsrJ0Rryw8ofI9zu7IpMWQ+Hq/n3JjxpD31U5WTXmV8qLLxEwbQt+Xfqj/0KtGxZ5Cr22ysdZJaxEZSWWhD6ltuzFmIZXTMSQlJdV74rtZnzZ4+XiTtyOjvm+l3FjB3gyanD1K04J8Ym7sRvxjNxPZv53VsZRyOnsKPQtoXe11DPC9h3mKSE/gdWCcMcYpT57w8vUmcmB7ArxKnHE45WLO7jpGyk/e4uLOdCa+fTv+A/oyaPJNVsdSyjL2FPp2IF5E4oATwJ3AXdV3EJFY4EPgbmPMNw5PeRWjXrkZTurDohuTopPnSX3sPU4s3oj4etP+4XF4jx2jD51QjV6dhW6MKRORh4HPARuwyBiTJiKzq7YvAJ4CmgHzqy4HKzPGJDVc7GpCQ+HoUUxhIRIU5JRDKotUVFC4fR8rbniJ8qLLRN8ygMQX7qJJbDOrkynlEuy6scgYkwwk1/jagmq/nwnMdGw0+5zLLWfd9E/o8mwoHWcMsyKCamCmvILTyalEBVwg6NIlOt7Tn1b3jSGyn86TK1Wd21/L1aRTNEV5l8hL1Q9GPY0xhsP/2sSyjr9k7ZSXyM8phoED6TX/x1rmStXCvW/9p2qRrthILuw+ZnUU5UBZ/9nD7jmLyd+TiV9UOL1emU7w7SP1enKlrsLtCx0gpGtr8rY49bNY1VAuXKBo0w42THod75Aguj5zB90fG4/Nz8fqZEq5PI8o9LBebTiTnErBsbP6AZmbyt6cTtZbX9JnQjSBPj4MfH060bcMxifY3+poSrkNjyj06FGdKd/VAXP+Amihu5WcbUfY9eQH5K7ajS3Aj453PUyTgT1p46Nn5EpdK48o9JbDOtGyqC8EVlgdRdnp4sETbH3wbXK+2oNXgC9xD9xIwtM3E9BcF9BS6np5RKHj40OFnz9Fh07QpEMHq9Ooqyg+ehr/7OP4ZxynIC2jssh/dzMBLbTIlaovzyh0YPXzO8jfn8XNJ4ZbHUXVYIzh+LJdpP3pYypOnWH8/5uMT6/u3JQ5AS9/feKUUo7iMYUe3L012at2czmvEL9wvWPUFZjyCg69tYGDc1dQcOA43mFBtJs1hooRo7AF+rn/TRBKuRiPKfSIvnEcBk6tOUDbqc5ZdUBdQUUFZGVx7I2vSH16JX5RTen2x2l0+dkYfJroVStKNRSPKfSY8Qls87aR9VGqFrpFCo/nkfZ8MqEVeXS6MY6YoW1IXDiLDvcNxcvHZnU8pTyexxS6b1ggTQd05MzKXRhj9JmRTpS9OZ205z7lzIpUTFk5rSb2ptOgQdiaNaPjKKvTKdV4eEyhA/R89AdI2l44fwGahlkdx7MZA6dPk/LEvzn01mbE15voqf3p+quJRCS2tTqdUo2SRxV61Pje4J0DZ05roTeQgmNnOfDS53TuHUBwmDetE8KRxybT9RdjCWgZanU8pRo1jyp0fH3Jy4OsF1bS87XOVqfxGMYYji/fxaH5q8hetRvKK2jy2Gi6/GI8LSa2oIVObynlEjyr0IFjR8o48PoWomffpD/611dJCRWZx/h03DwKj5zBK9CPmDsH0+WRcfpnq5QL8rhCb/c/Aznw+/fJeG+zls51qCgr5+i/U8hZmUq/W9rgVVFB7Kj2+Dx0Ix1/NEIXy1LKhXlcoYfEt6BJlxhOf7oT5k6zOo7byEnJ4NCCrzj50TZK8y7iHRJIj+n9CEjoRK9JwVbHU0rZweMKHSBqQh8OzV3GhW9OE9qxpdVxXNfly3DiBJnvbGTzU/8BLyF8UGfi7h1G+x8OwuavKx4q5U48stDb3zuEzL9/xvk1OwntOM7qOC7lcl4h6W+u5/iSLbRPCiV+bAei+kTR6cmpdJw1UteTV8qNeWShN+0ew+QvHsInLxsuXYKAAKsjWausjK8XruX4+1s5u/EApqwcvxZNKZnYA0aMwC84mD4TrA6plKovjyx0AJ9e3TFffknOxxtoPm2M1XGc7lJ2PmdW7qZtJ3/Izibr1ZVcOFNMqzsG0e6eoUSP6aZ30yrlYTy20AkIYPf6Cxx4PpnRHToQ2TfO6kQN7vzBU2S+u4WTK3aSv/MIYIj6YBp+Hdsy4MNHCYqPQfQhy0p5LM8tdKDjzyfwzSur2PnYe/xgza+tjuNwFWXlmLN52M6d5ZtF60l9fjUAAbGRxE4fSdtpA/Ed2QVEaGJxVqVUw/PoQg+MCiPux6NJf2kFGR9sJ+62vlZHqrf8w9kc/3gHpz7fTd7mg/R/oA9thrWlZd/WdHx8Cm3vHECzXrFWx1RKWcCjCx0g4embOfHRdrbdO5/AqF/TYki81ZGuiblcgpw/R/HhLD67620uHc0GwKdpEyJH9sBn6AC4MYEQHx8Sb7M4rFLKUh5f6L6hgYz8Yg4bJzyH995U6NMKAgOtjnVFl7LzOfHZXs6s3k/elkM0a+HNoEcH4idCWJeWtL5zEDGTetN8UAf9UFMp9R1ijLHkwElJSSYlJcVpxzP5+cimTVTYvClsn0Bwu0inHfuKmYzhYvppQgLKIS+PNT/9kFNrvwFAvG0Ed21NqwkJ9Hp0NDRtCl76gaZSjZ2IpBpjan2Kj8efof+XhIRAv358/dgi9i5+h85P3EzPJyY69aqPgowcTq35mrPbDnN+Zyb5aceQsjJuee8WvPx8ad4vlsCEDrS8oRvRo7vhHejrtGxKKffXaAodgPBwWv78DjL3vMb+3y7m+JLNdH1yCm2mJmLzddwfxeW8QnJTMzi3+zjn9xyj7w874Wcuc+SNFNKWpIGXENCmOc3HJNCsfwcqBg/EK6IpXcfpFIpS6vo1mimX6owxHHhpJWlPf0BZfhGxozow+K+3UhHWlOyMIsJ7x+IbUvvdpcYYLucWUHTqPAWZuRRm5tK6dwRNQrw4/sXXbP3fryjNu/jt/uLrzeh5NxPRtx0X8qHgko3mgzrow5KVUtel3lMuIjIW+BtgA143xjxbY7tUbR8PFAH3GWN21Ct1AxIRuj5yIx1njyTzg200qciHrCwKN+9j9exPAfDy98UW6IfN35feM3vTdkhrzn6dw8pHkqG84jvvF/irwTS5IZ6gyECaDe5Mk/gWNO0ZS3jvNoR1jcbLu/IByaFVv5RSqiHUWegiYgPmAWOALGC7iCwzxuyvtts4IL7qV3/g71X/dWneAb50uGdI5YuKCnxP5ZH4Wisu7Mvicl4B5YWXKSu8jHfLCIiMxNcE0nb6KPwigvGPDCGwdTjB7ZsT2qklBPoRDoy819IhKaUaMXvO0PsB6caYIwAishiYDFQv9MnAP03l/M0WEQkTkShjzCmHJ24oXl74tYqg48zhV9wlGBj4g8HOy6SUUtfAnks8WgHHq73Oqvrate6DiMwSkRQRScnJybnWrEoppa7CnkKv7dKLmp+k2rMPxpiFxpgkY0xSZKT114ErpZQnsafQs4DW1V7HACevYx+llFINyJ5C3w7Ei0iciPgCdwLLauyzDLhHKg0ALrjV/LlSSnmAOj8UNcaUicjDwOdUXra4yBiTJiKzq7YvAJKpvGQxncrLFu9vuMhKKaVqY9d16MaYZCpLu/rXFlT7vQEecmw0pZRS10JXe1JKKQ+hha6UUh7CsrVcRCQHOHqd3x4B5DowjjvQMTcOOubGoT5jbmOMqfW6b8sKvT5EJOVKi9N4Kh1z46Bjbhwaasw65aKUUh5CC10ppTyEuxb6QqsDWEDH3DjomBuHBhmzW86hK6WU+j53PUNXSilVgxa6Ukp5CJcudBEZKyIHRSRdRObUsl1E5OWq7XtEpI8VOR3JjjH/T9VY94jIJhFJsCKnI9U15mr79RWRchG51Zn5GoI9YxaRESKyS0TSRGStszM6mh1/t0NFZLmI7K4as1uvCSUii0QkW0T2XWG74/vLGOOSv6hcCOww0A7wBXYDXWvsMx74D5XrsQ8Atlqd2wljHgQ0rfr9uMYw5mr7fUXlmkK3Wp3bCf+fw6h8Klhs1evmVud2wph/DTxX9ftIIA/wtTp7PcY8DOgD7LvCdof3lyufoX/76DtjTAnw30ffVffto++MMVuAMBGJcnZQB6pzzMaYTcaYc1Uvt1C59rw7s+f/M8BPgKVAtjPDNRB7xnwX8KEx5hiAMcbdx23PmA0QXPXQ+SZUFnqZc2M6jjFmHZVjuBKH95crF7rDHn3nRq51PDOo/BfendU5ZhFpBdwMLMAz2PP/uSPQVETWiEiqiNzjtHQNw54xvwp0ofLhOHuBnxljKpwTzxIO7y+7ls+1iMMefedG7B6PiIykstCHNGiihmfPmF8CHjfGlFeevLk9e8bsDSQCNwABwGYR2WKM+aahwzUQe8Z8I7ALGAW0B74QkfXGmPwGzmYVh/eXKxd6Y3z0nV3jEZGewOvAOGPMWSdlayj2jDkJWFxV5hHAeBEpM8Z87JSEjmfv3+1cY0whUCgi64AEwF0L3Z4x3w88ayonmNNFJAPoDGxzTkSnc3h/ufKUS2N89F2dYxaRWOBD4G43Plurrs4xG2PijDFtjTFtgX8DD7pxmYN9f7c/AYaKiLeIBAL9gQNOzulI9oz5GJU/kSAiLYBOwBGnpnQuh/eXy56hm0b46Ds7x/wU0AyYX3XGWmbceKU6O8fsUewZszHmgIh8BuwBKoDXjTG1Xv7mDuz8//wM8KaI7KVyOuJxY4zbLqsrIu8BI4AIEckCfgf4QMP1l976r5RSHsKVp1yUUkpdAy10pZTyEFroSinlIbTQlVLKQ2ihK6WUh9BCV0opD6GFrpRSHuL/A2DXhj2610UZAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import pylab as pl\n",
    "\n",
    "h = 0.01\n",
    "n = int(1 / h)\n",
    "f = lambda x, y: -50 * y + 50 * x ** 2 + 2 * x\n",
    "\n",
    "# f: 函数体\n",
    "# x0 初值\n",
    "# y0 初值\n",
    "# h 间距\n",
    "# n 不用你管\n",
    "def R_k(f, x0, y0, h, N):\n",
    "    x = np.zeros(N + 1)\n",
    "    y = np.zeros(N + 1)\n",
    "    x[0] = x0\n",
    "    y[0] = y0\n",
    "    for i in range(N):\n",
    "        x[i + 1] = x[i] + h\n",
    "        k1 = f(x[i], y[i])\n",
    "        k2 = f(x[i] + h / 2, y[i] + h / 2 * k1)\n",
    "        k3 = f(x[i] + h / 2, y[i] + h / 2 * k2)\n",
    "        k4 = f(x[i] + h, y[i] + h * k3)\n",
    "        y[i + 1] = y[i] + h / 6 * (k1 + k2 * 2 + k3 * 2 + k4)\n",
    "    return x, y\n",
    "\n",
    "\n",
    "[xx_R_k, yy_R_k] = R_k(f, 0, 1 / 3, h, n)\n",
    "print(xx_R_k, yy_R_k)\n",
    "\n",
    "xxx = np.arange(0, 1 + h, h)\n",
    "yyy = []\n",
    "for i in xxx:\n",
    "    yyy.append((math.e ** (-50 * i) / 3) + i ** 2)\n",
    "    # yyy.append(xxx)\n",
    "\n",
    "pl.plot(xx_R_k, yy_R_k, label=\"R_k\", linestyle=\"--\", color=\"purple\")\n",
    "pl.plot(xxx, yyy, color=\"red\", label=\"right\", alpha=0.3)\n",
    "pl.legend()\n",
    "pl.show()\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtPElEQVR4nO3deVhU5f/G8fczw64Isrgi4b6FkoJppV9bXEtNs1xS0zKytMxWNdfM0tIs00Qzs0Uzy8y1cinNXXEXSMUdV3BXRGB4fn9g/hBRRhk4M8PndV1eceacmbmf0JvDM2dRWmuEEEI4PpPRAYQQQtiGFLoQQjgJKXQhhHASUuhCCOEkpNCFEMJJuBj1xgEBATokJMSotxdCCIe0efPmJK11YE7rDCv0kJAQoqOjjXp7IYRwSEqpQ7daJ1MuQgjhJKTQhRDCSUihCyGEkzBsDj0naWlpJCQkkJKSYnSUAuPh4UFQUBCurq5GRxFCODi7KvSEhAS8vb0JCQlBKWV0nHynteb06dMkJCRQvnx5o+MIIRxcrlMuSqlpSqlTSqldt1ivlFLjlVLxSqkdSqk6dxsmJSUFf3//QlHmAEop/P39C9VvJEKI/GPNHPp0oPlt1rcAKl/7EwlMykugwlLm/yls4xVC5J9cC11r/Q9w5jabtAG+05nWA75KqdK2CiiEEM7ki69fY8fe1fny2rY4yqUscCTLcsK1x26ilIpUSkUrpaITExNt8Na2ZzabCQsLu/5n1KhRt91++vTp9OnTp4DSCSEc2bfvfojnq4dY/uv3+fL6tvhQNKc5gxzvmqG1ngJMAQgPD7fLO2t4enqybdu2fHv99PR0XFzs6rNoIUQBWDh7AYzdzuWyHvR+fXy+vIct9tATgHJZloOAYzZ4XbsSEhJCUlISANHR0TRu3PimbRITE3nqqaeIiIggIiKCNWvWADBs2DAiIyNp2rQp3bp1K8jYQgg7kHj8OAde/JZU1wyeWDgUN3f3fHkfWxT6fKDbtaNd6gPntdbHbfC6dJi8jp+jM2dz0iwZdJi8jrlbEwC4kmqhw+R1LNie+bPjQkoaHSav449dmW995nIqHSavY1nsSQBOXbTuSJIrV67cMOXy008/WZ23b9++9OvXj02bNjFnzhx69ux5fd3mzZuZN28eM2fOtPr1hBCOT2vNmGldSSy7H8uwR6kYWiPf3ivX3/2VUj8CjYEApVQCMBRwvRY0ClgMtATigWSgR36FLQh5mXJZtmwZsbGx15cvXLjAxYsXAWjdujWenp62iCiEcCDjf3iX2JjltB76Di926JWv75VroWutO+WyXgO9bZYoi59eanD9a1ez6YZlTzfzDcvFPFxvWPYr4nbDcglvjzxlcXFxISMjA+CWx41nZGSwbt26HIu7SJEieXp/IYTjmf75NDzf+5fgjo3p+cztD7CwBbmWi5VCQkLYvHkzAHPmzMlxm6ZNmzJhwoTry/n54aoQwr7ti40lbeAfpJnc6fvOtAI550QKPZvsc+j9+/cHYOjQofTt25eGDRtiNptzfO748eOJjo6mVq1a1KhRg6ioqIKMLoSwE8nJV1jQYhjuV03UnfkiVaoUzKU9VOaMScELDw/X2W9wERcXR/Xq1Q3JY6TCOm4hnNWgB7pRc106V96swfNjBtn0tZVSm7XW4Tmtkz10IYSwod+WfEnxg4eJrWOmxyfvFeh7S6ELIYSNxO7fyHcz3iC+62kGrooq8Gs1SaELIYQNHDp4hJ+bDsTzqi+D35qHp1fBH9kmhS6EEHlkSU/nlxZvU2lfIA+GDadMYAVDckihCyFEHk3q9CZl/jVxqXswr/R/ybAcUuhCCJEH44d8it8vpzgeZuLFrz8yNIsUuhBC3KX4w9th0jJOFk+jw6JPMZmMrVS5jms2ZrOZ0NDQ68sdO3a8fnKREEL85+Lls4z89EmuNjvPWy8vpmyZEkZHkkLP7m4vzmWxWG55BqkQwrlYLBaGPtmFs/5H6PfmLOpE1Dc6EmDHhf76H6+z7cQ2m75mWKkwPmv+mc1eLyQkhOeff54lS5bQp08fOnbsaLPXFkLYr/Gd3iBimQ/uT7/I/yLaGx3nOrstdKP8dy2X/wwYMIAOHTrccnsPDw9Wr86f+wMKIezPr59PpsTPpzgWamLEzAm5P6EA2W2h23JP+k7c6ZTL7cpeCOFclixcxvm3/yI5QPPC3xNwcbGvaVY5yiWP5DrnQhQO5y+dZlfPCWQoTdjsvvj6+xsd6SZS6EIIkYuMDAvDP2vLlogVuH3SlAcfftDoSDmy2ykXo2SfQ2/evDmjRuX/nUaEEPZrQN+e7E1aRfsXhtL1yZ65P8EgUujZWCwWq7c9ePBg/gURQtiFqQM+ImzCVVJbtKVbm6FGx7ktKXQhhLiFf35biOsn2zhVJoPBM74p8Mvh3ikpdCu0bduWAwcO3PDY6NGjadasmUGJhBD5bef2XcR2/RaTB7ReOgK/4j5GR8qVFLoV5s6da3QEIUQBSrmazO8th1Iq2QX/756hfI1qRkeyihzlIoQQWWit+WhSZ/ZVXkfSa6E8/uxTRkeymhS6EEJkMfKTN9iycR6hz7fjjXG2vcFzfpNCF0KIa7756AtC+h+n/JUmvNLlc6Pj3DGZQxdCCGDDkuWoIas545fBW59Px2Syr9P6rSF76NmYzWbCwsK49957adWqFefOnbvltitWrOCJJ54ouHBCiHwRF7ObLU9Hke6qafrHQIKDyxgd6a5IoWfz38W5du3ahZ+fHxMnTjQ6khAiH6VcSWZ+0/coeslMUNRTVAu/z+hId82up1ymN55+02M1n6lJxCsRpCWnMaPljJvWh3UPI6x7GMlJycxuP/uGdd1XdL+j92/QoAE7duywattNmzYRGRnJnDlzqFDBmDt+CyHujNaaj6I6c6XEbjJaPU3Xbo599VTZQ78Fi8XC8uXLad26da7brl27ll69ejFv3jwpcyEcyOiJr7Nl0zyC+z3KgKj3jY6TZ3a9h367PWpXL9fbrvcK8LrjPXL4/4tzHTx4kLp169KkSZPbbh8XF0dkZCRLliyhTBnHnHcTojD6YsAoSo89yj1dm9C7y3ij49iEVXvoSqnmSqndSql4pdRNd0xWSvkopRYopbYrpWKUUj1sH7Vg/DeHfujQIVJTU3OdQy9dujQeHh5s3bq1gBIKIfJq5dz5FPtkK2eKw6CRM1Am55isyHUUSikzMBFoAdQAOimlamTbrDcQq7WuDTQGxiql3GyctUD5+Pgwfvx4xowZQ1pa2i238/X1ZdGiRQwcOJAVK1YUXEAhxF3ZuC6aPV1+4KpHBm3+ep9SpQKNjmQz1vxYqgfEa633a61TgVlAm2zbaMBbZV6KrChwBki3aVID3HfffdSuXZtZs2bddruSJUuyYMECevfuzYYNGwoonRDiTiWdOMHqJ0bhmmri3h97UqFm9n1Tx2bNHHpZ4EiW5QTg/mzbTADmA8cAb6CD1joj+wsppSKBSIDg4OC7yZvvLl26dMPyggULbrlt48aNady4MZA5npiYmPyMJoTIA4slnQ8ntcffx8LZXk/zQCvnu1qqNXvoOV0AWGdbbgZsA8oAYcAEpVSxm56k9RStdbjWOjww0Hl+zRFC2L+hn3Zh79413DPucfqNfMPoOPnCmj30BKBcluUgMvfEs+oBjNJaayBeKXUAqAZstElKg+3cuZOuXbve8Ji7u7tMrwjhIN7v8Drlfr9Expud6NLGsS64dSesKfRNQGWlVHngKNAR6Jxtm8PAo8AqpVRJoCqw35ZBjRQaGsq2bduMjiGEuAuzR39BxdknOVixCEPe/cboOPkq10LXWqcrpfoAfwJmYJrWOkYp1eva+ihgBDBdKbWTzCmad7XWSfmYWwghcjX3u9lceW81F0pqem+YhIeHu9GR8pVVJxZprRcDi7M9FpXl62NAU9tGE0KIu7dx1TqSXvqZdA944q8R+Pr7Gx0p3znH0fRCCJHFmQsnmfJtVy4UP0eJqd0c5hZyeSWFnouWLVve9hK6kHn4YnR09E2Pb9u2jcWLF+fwDCFEfjl//gL9RzQj6eohGix+jac6tjI6UoGRQr8NrTULFy7E19f3rp4vhS5EwbKkpzOuwYtUneHL0x3G8EBY4SlzkEK/ycGDB6levTqvvPIKderUwWw2k5SU+fnuiBEjqFatGk2aNKFTp06MGTPm+vN+/vln6tWrR5UqVVi1ahWpqakMGTKEn376ibCwMH766SejhiREofFF695UjTNzvnZ1nn2ir9FxCpz9Xm0xJgbOn7fta/r4QM2auW62e/duvvnmG7788ktCQkIAiI6OZs6cOWzdupX09HTq1KlD3bp1rz8nPT2djRs3snjxYoYPH86yZct4//33iY6OZsKECbYdhxDiJuO6v02p3y9y8iEPhi0qnDemsd9CN9A999xD/fr1b3hs9erVtGnTBk9PTwBatbrxV7l27doBULduXQ4ePFggOYUQmb4YMIoS3yZwuFIGbyyLwuQkV0+8U/Zb6FbsSeeXIkWK3PRY5kmwt+bunnl8q9lsJj3d4a9LJoTDWLd9EZt3jKVayAN0XfkVru4OfaHXPCmcP8buwkMPPcSCBQtISUnh0qVLLFq0KNfneHt7c/HixQJIJ0ThtGzVcsZ+1gFdwYMXt0ZRtkwJoyMZSgrdShEREbRu3ZratWvTrl07wsPD8fHxue1zHn74YWJjY+VDUSHyQdzW7extEUWFrdV4f8AS/H1LGx3JcCq3qYT8Eh4errMfux0XF0f16tUNyWONS5cuUbRoUZKTk2nUqBFTpkyhTp06eX5dex+3EPbmxOEj/FL3DbzPmTFHtaPLC88YHanAKKU2a63Dc1pnv3PodigyMpLY2FhSUlJ47rnnbFLmQog7c/r0WWY++CYlzrjg83kTWhWiMs+NFPodmDlzptERhCjU0tJTmfRAJBUTXLn0dihd+jxvdCS7Yndz6EZNARmlsI1XiLultWbklx05HrCR+PYlePHjAUZHsjt2VegeHh6cPn260JSc1prTp0/j4eFhdBQh7N57Q55n8/q5VOnZgsE/f2Z0HLtkV1MuQUFBJCQkkJiYaHSUAuPh4UFQUJDRMYSwax8/+wY1f7zC5R5teK37JKPj2C27KnRXV1fKly9vdAwhhB35fuAoys08weEq8OH4WSiV022OBdjZlIsQQmQ1YeinqFHbOBVsoc+mryhSRKYnb0cKXQhhl/6YO4tiH2wksbiFZ9d+RpFi3kZHsntS6EIIu7Nj72qmzH+Bf+8/TOOlwwgoK2eBWsOu5tCFEGLJwmV8/01XPAOL8Mqv3xJUqrLRkRyGFLoQwm7Ebd7C/o6Tud8lgvrrBkqZ3yEpdCGEXTi4ew9/P/YRXldNFJnYlfDq9XN/kriBFLoQwnDx8QdY8OBA/C+YKTn1SZo997TRkRySfCgqhDDUxctnmdH+FQLPuHJhUAOa9ehkdCSHJYUuhDDMlauX6T/qUbZXXU76qAj6DH/d6EgOTaZchBCGuHjhIiOadOR48E66RI6lXZPXjI7k8KTQhRAFLvXqVaIejOS+XT641nhVytxGZMpFCFGg0lJTmfBgT4J2Kc50KcPIbz41OpLTkEIXQhSYdEs679/fndKbLZxpH0jv78cYHcmpSKELIQqE1prRH3clZHcyOx905+WfxhkdyelIoQsh8l1GRgYfRnVl/c5ZnB1RnBErp2IySf3YmlX/R5VSzZVSu5VS8Uqp/rfYprFSaptSKkYptdK2MYUQjiojI4Ohj7xA2rg4Iup35M03pmE2m42O5ZRyLXSllBmYCLQAagCdlFI1sm3jC3wJtNZa1wTkNC8hBFprvmjXhxorr5LhVZUBvb6XG1TkI2v20OsB8Vrr/VrrVGAW0CbbNp2BX7XWhwG01qdsG1MI4YjGPtmbEvPOcaKeK0Ojv8PVTY6Uzk/WFHpZ4EiW5YRrj2VVBSiulFqhlNqslOqW0wsppSKVUtFKqejCdN9QIQqj4a16UXb+OY7cZ+LV1VMxu0iZ5zdrCj2n3490tmUXoC7wONAMGKyUqnLTk7SeorUO11qHBwYG3nFYIYRjmDizH8cvLGBPLU2fNVNxcXU1OlKhYE2hJwDlsiwHAcdy2OYPrfVlrXUS8A9Q2zYRhRCO5L3Bvfhj4WcEtgpj0Jbv8fKU+4AWFGsKfRNQWSlVXinlBnQE5mfbZh7QUCnlopTyAu4H4mwbVQhh70a360OND84R4tKCYf3mYTbLNEtByvX/ttY6XSnVB/gTMAPTtNYxSqle19ZHaa3jlFJ/ADuADGCq1npXfgYXQtiXL555jeC5ZzgUamZ01K9S5gZQWmefDi8Y4eHhOjo62pD3FkLYjtaaoc0jqb7kMifquNBn3VRc3dyMjuW0lFKbtdbhOa2TU7WEEHdNa82Yt1+m+pLL7K6JlLnBpNCFEHdFa83YaS/yz9HJ7O1qoX/0N1LmBpNJLiHEHbNYLAx+tDtHvBZQv+UzDHhlJiaTnM5vNNlDF0LcEUuGhc8ffZHQlRbKnH+cAa/8KGVuJ6TQhRBWS716lc8e6EHplSmcfKwoI1d+J1dNtCPynRBCWOVqSgof1nmOMhvSSXy8OK/9GYWLi+yZ2xMpdCFErq6mXmHQxy0pdvokO5t48+qCCbJnbofkOyKEuK3ExEReH/Aw//77N8XHP8SHSybLJXDtlBzlIoS4paQTJ5hWtw81rhYjPGoEPdoPMjqSuA3ZQxdC5OjogYPMqPMqQcddudS+AS9Imds92UMXQtxkc/Q21jQbid85M24jH2DAgFeNjiSsIHvoQogbHDgaw19PvofPeRMZox/haSlzhyF76EKI62L2rWfEx4/j3sCFx5t+xjMvdjI6krgDUuhCCAC+m/Qte8dPgoeg/6j5VK94v9GRxB2SQhdCMPfzKei3llPKoxwPdx0rZe6gpNCFKOQmv/MhRcbu5JKfpvlfw6kYWsPoSOIuSaELUYh9+Gw/Qmae4EQZC8+u/5SS5YKMjiTyQI5yEaIQ0loz7puXOHDsR/ZXSefZLV9KmTsB2UMXopC5knyFAS91Yp9lHvXaPsWA3j/K/T+dhOyhC1GInDt9mi9qP0e9HzypWbY7g179WcrciUihC1FIHNqzl29r9SIo3szxbsGM+uQbuciWk5EfzUIUAn8tXUFs+wn4XXLB5f36vDm4r9GRRD6QPXQhnFx0zDIWDXsNzxQTKZ804xkpc6cle+hCOLHpv0Qxb35fioT50e2zydSOaGB0JJGPZA9dCCc15rm3ofMy/NOrMvb9DVLmhYDsoQvhZCzp6Xzx+CuUXXKJQ+VNDBn1ByX9yxgdSxQAKXQhnEhiYiJRD/Wh0h4zJxt68ubSSbi6uxkdSxQQmXIRwkmcOH2Iz9t1oOIeE/FPBvDaiq+kzAsZKXQhnMDmuDW8MbgecSHruDr6fgbPHY/JJP+8CxuZchHCwU3+4HNSP10CLWDwoN8Jq97Y6EjCIFLoQjgorTWTXxqA99RDJBX3pNeL8wirXt/oWMJAVv1OppRqrpTarZSKV0r1v812EUopi1Kqve0iCiGyS06+wvsNuuHz1WESqyi6xn5Jo8ZS5oVdroWulDIDE4EWQA2gk1LqpivgX9tuNPCnrUMKIf7fmQsnGf1we6pssLDvIU/67PwWv5IljI4l7IA1e+j1gHit9X6tdSowC2iTw3avAnOAUzbMJ4TIIvbAZvq+V4d/g1dy4pWKDFr1NS6urkbHEnbCmjn0ssCRLMsJwA03HFRKlQXaAo8AEbd6IaVUJBAJEBwcfKdZhSjUJr7/ORe/mkfqY8m8OeQX6oU2NzqSsDPW7KHndH1NnW35M+BdrbXldi+ktZ6itQ7XWocHBgZaGVGIwi0jI4OJXd/Ed+h6PJL9eLnHYilzkSNr9tATgHJZloOAY9m2CQdmXbu2cgDQUimVrrX+zRYhhSiskhKTmPzYq1TYoThR00T3FVEUDwgwOpawU9YU+iagslKqPHAU6Ah0zrqB1rr8f18rpaYDC6XMhcibw8d380PDt6iwz5sjTxTnrbnjMbmYjY4l7Fiuha61TldK9SHz6BUzME1rHaOU6nVtfVQ+ZxSi0Fm5ZR4TJnXDu6YXRdu/wTuj3jY6knAAVp1YpLVeDCzO9liORa617p73WEIUThkZGYzo+Dp6+1q8mvrQ//PfqBJSx+hYwkHIxR6EsBPnTp9mXEQ3qvychCmtAsP7r5UyF3dETv0Xwg4s/3MFMV0mUDrJlTPtS/HezE8wu8o/T3Fn5G+MEAZb/Pc3JLRbQNE0F1xGN6L3Oy8bHUk4KCl0IQySkpLC5z/0YfVfX1O1eT2aRY6hSbOGRscSDkwKXQgD7I+LY3aTwZwss42ILu14p9f3eLh5GR1LODgpdCEK2OJvfuBon/mUSTGT8lh7Br/6EddOyhMiTxzuKJd/on+l19u1OHR0r9FRhLgj6WlpDH68F2dfWITFFSr91oNh00dJmQubcbhC37bwMI3G1GT5L6uMjiKE1RLPHmVgr6bUWHyBfZU07WM+54FWcj0WYVsOV+iVQ6uhUKQePmp0FCGsMuOXH3i1fy32pK7m3AcVGRg7g4CypY2OJZyQw82hN2gYzu98R/qx00ZHEeK20lJT+eyZfpRacJqiTwQy4OO51K7ayOhYwok5XKH7lgggxc3ClYSzRkcR4pb2bN/J4ic/IOigmf01TAwYt5SKFcrl/kQh8sDhCh3gXBELqUcvGx1DiBxNHDoW99Eb8Es3kfJaVQaMG4zJ5HCzm8IBOWShn4y4wlWL3OlO2JcrVy/z6bQXufDbNkp71aDCVz1o/VQLo2OJQsQhCz2gXQD/rF2K1loO+RJ2Yd6s+SyeOYhjRXfy0Etd6PnsRIr6FDM6lihkHLLQ/YqVo8hZN44nHaFMoNybVBjHkmFhauR7eH1zkHu9y9P817do+0g3o2OJQsohJ/b2/G6m9YLW7Nmw3egoohCL3riFcZU7U+zrwyRVgAeWjZQyF4ZyyEK/N6IqAEdj4w1OIgqrObMnsr3RKEocVFx+qQJ9d8+kbvi9RscShZxDFnrjxhEAnNt33OAkorA5euoY741pyfS5fThc6wzFZj1Hz6gP5CgWYRcccg49sHIIGUpz6XCS0VFEITLj08mcHf4HB/63hmbt+/BSpzG4urobHUuI6xyy0PeeTuaCl4WrRy8aHUUUAmdOnWJ6+3cpveoqqpiZpk98SfduzxodS4ibOGSh3+PvxQ+NznHBS/bQRf6aNu4rUob8SalLrpx61Jses0dTzK+40bGEyJFDTvx5e7ji38Kfwz5xaK2NjiOc0KUr5xn5ZScOT4pCo0gb14S+y6ZImQu75pB76ADuuhSBh705c/k0/kUDjI4jnMiUjyexfu2nJHnF82DvLjz9zBhKlC5pdCwhcuWwhb5rITRd9hh7Y3bhf39jo+MIJ3AyIYEZHQZRem0qlctWovu8T2lUt5XRsYSwmkNOuQDUa1AdgCM79xicRDiDzwZ+zLyq/Si59ionHitK143fSpkLh+Owhd7of3UASIw/YnAS4ciOJR1gcNcnKPnRNlLcNWVndaDf0ikElSlhdDQh7pjDTrl4hdwDwMWDctVFcefS09IY+8UwNu78jIyMNLw6vsjLk0bh6+ttdDQh7prDFvrCw+fxcMkgNeGC0VGEg1n/x1LWvDCZEkkZ+EdW4523vqbSPWFGxxIizxy20FuElibqiROc9pdCF9Y5dTKRb58bROmlFyjuauJIh0pMHDsbVzeH/WcgxA0c9m/yPf5F8H3Ij517txgdRdg5rTWLFk/jcMc/CLrkyolwV56eOZKgyhWMjiaETVn1oahSqrlSardSKl4p1T+H9c8qpXZc+7NWKVXb9lFvdDXdgj5bgjKxxbl4VS4BIHK2fusG3vigEZNn9CSxwmnOfNCIfpu+kzIXTinXPXSllBmYCDQBEoBNSqn5WuvYLJsdAP6ntT6rlGoBTAHuz4/A/7lwJZ09SxRNNj3IvoQ9hFWsm59vJxzM+bNnmNztPQL/TOJUyzie7Nyfbt8Ok4tpCadmzR56PSBea71fa50KzALaZN1Aa71Wa3322uJ6IMi2MW/mV8SNhg+FArBrc3R+v51wEBkZGUwfPpYf73mJcgvPk1DBhV79FvHCMx9JmQunZ02hlwWyHuydcO2xW3kB+D0voaxhNikeeawBAAc27czvtxMOIGbvekZXeAb3YZuxuIBPVHMG//sjDf+Xr78sCmE3rPlQNKe7MOd4RSyl1MNkFvpDt1gfCUQCBAfn/V6gp0NKodFc3nosz68lHNfe/fv4cfFANq/7hdBiddj15AO8+91IinkXNTqaEAXKmkJPAMplWQ4CbmpQpVQtYCrQQmt9OqcX0lpPIXN+nfDw8DxfJnHcusPU91G470tBa41SOf3sEc7qSvJlpr4yDK8fD3Go0UrqtWlPr4mf4e9b2uhoQhjCmkLfBFRWSpUHjgIdgc5ZN1BKBQO/Al211gV2cZWhrWqw8khpNsf+yv6z+6noV7Gg3loYKCMjg6nDxpL+xQYCzrlxpBw06zGFZzq3yf3JQjixXAtda52ulOoD/AmYgWla6xilVK9r66OAIYA/8OW1veR0rXV4/sXOVKmEN1eaNuK3PZ+y9fBGKfRCIDpmGUtaT6D8fi/OFle4f9KQt9/sJb+dCYGVJxZprRcDi7M9FpXl655AT9tGy13C2WT2ngwkNOZedi/dAGGdCjqCKCDL//yL+RtGsP/fFVQtHcqVuo/xzrQReBctYnQ0IeyGw54pCrDp4BkGLT/O4B2hbPv7CLxtdCJha/ExMczrNYaSa66Q0SCB5r1e5bmn3qeol6/R0YSwOw5d6I9ULcn8Ic1YNeNXXPdekQ9Gncj++APM6DWS4BUXKZVh4mC4O50/+YkGD9QxOpoQdsuhC93HyxUfL1dWV/Cl5N5kjlw4QrBP3g+HFMa5ePks3/02nPTX91ApyZejoa40ndibZxs2MDqaEHbPoQsd4O/dp0gJKU3AlrNsil1HcAMpdEd06mQiX/Qazu6iM7hiOUetVk0JeKgdbz3f3uhoQjgMhy/0GesPgUcgjU0Wdq/bBA06GB1J3IHzZ88w8/WPcJ19iBopLpxvEkb3cYOoU/NRo6MJ4XAcvtA/bBdKkVZV6V/0Iy6b5O5FjuLi5XOM6zqY0r+fwDfFlZMhCt93WvB5ry7yOYgQd8nhC72EtwfgwT2lKvP7oS3ywaidSzhxjN+Wf8LKv77mgbXhJPr7ETCwOa+/0s3oaEI4PIcv9JQ0C1+vPoDnrvI02J7C8fePU8a7jNGxRDaH9uxldt8xFP/7DKuaLSK4zn3Un/cW9eu1kB/AQtiIwxe6m9nEl3/H8/z5QKrtv8DK3X/TKfxZo2OJa/5ZvprVg78iaONVgiwm4iuZaNt5Bh07tDU6mhBOx+EL3WRSRA9qwqFpbmxdfoBV8xdJoduBLXF/88tPo6k+wpdgrUgMd6fRyB50atLY6GhCOC2HL3QATzcz5R4PYwualKVHuDrkKu4ucjODgpaelsbXH3zB/gV/EFttKW7unhTp3JGGPV+my/8ijI4nhNNzikL/98QFvl5/hNqV/Qj99wrLDyynZeWWRscqNI4dO8HcIeNQc/ZR/JwbFbx88Hm+H326v0exov5GxxOi0HCKQr981cLfu0/xcKsGnJjzFYs2zJZCLwB7D2/jl08/peSXV/FLM5NUykTGe2F0fOdlfIp5Gx1PiELHKQq9TrAv0YOaQHIyu1wns2nXAiwZFswms9HRnE5aaipR748nbs8Cjrj8g2eaJ7Uqt6H0i03p0/c5OWJFCAM5RaFfLxEvL8KrPcry37ew9shaGt7T0NhgTmTntp2sGDUV18XHCLjoSnAZN6p+2JsOj79FqYAQo+MJIXCSQgeYvekIq+KT6La/Ct1nt2Fh819o2EMKPS8slnRWRP/CP2/8ROV17gRoE6eCTaS+Vodeb7+Er49MqwhhT5ym0C+kpHHi/BWCWzUgfsRCTszehu4uZ43ejc1r1jN32FT2lVnIpbSTVHKtTmy9CB5+rwudWjUxOp4Q4hacptB7NqxAz4YVALCU8qTi9iJsPbGVOqXl+tnWOH78BHM+nET6gn8pcchEDeDcE9V49s33afpgN9xcPYyOKITIhdMU+n/SLRkEPR6BmraScQs+5vvIWUZHslvpljRWRs9hxe8/UHFkEfzTzZwvlsHZdqWo93onOjesb3REIcQdcKpC/2H9Ib74ay+/dGnM8a//IeX7f9nRZge1StYyOprdsGRY+GvWHNZPXETa6X1srbMGD4+iuP7vCQKbPULkGz0wm+XoICEckVMVeqUSRXmsekmKRFSg1rDH+HXfAoYvHcScLvONjmYorTULfprDxokLCdp+GZ+LroSYMthXsRTPPjeJtg8/h7ubp9ExhRB55FSFXr+CP/UrZJ6Z6NevPU9NXM7wdTOJfiSa8DLhBqcrWGmpqfz41UziU5cSv3cZ5f8pSfXYGhy/x4RXZE1a9OtBYFm5KqUQzsSpCv0/+xMv4VfEg3r+7Xh99EVGVxjOz28uMDpWvjt79iwLp8wgcUE0RbdexjvZlZONV1G6cSVqDXqSe2u2oUtoFaNjCiHyidMV+pEzyTwydiWDHq9O+0YN8L70K2W+OsrSJ5fSpKLzHXJ34Nhe1m6dw46/llJ9fADu6Wb8XSycqeLG+caV6ffaSKpVLW90TCFEAXC6Qi/n58Un7WvRqEogvsU8KNvhAerNXM27Y16m/Mg/qORXyeiIeXL50kVWzJrL3rlr0RuTuOh7iE0R0fh4l6BovSZ4/e8+evWPpIhcS0WIQkdprQ154/DwcB0dHZ3v75N89CwL732bc6nn+eaNvSwe+DfFPYvn+/vaSobOYMfeNWza8TvHP4gjOMYN93QzGUpzsrSFk/WK8+zH3QmtVF9OohKiEFBKbdZa5/ihoNPtof9n78mLjPr9Xz7tEEbjRf3567ERPLDARMfyT7HwuT9xNbsaHTFHaamprP19CXELV3NuzSGKHL/EkuYLQcF9Lg3Ze28wVVrX4fHITgSULW10XCGEHXHaQr+QkkbMsQscOZPMvQ9UovmqIfiv+YINqz/nsfTGTH/6B8oXN35u+cTx42zbv4r4g2s4/uNuKi4phmeaCz6AqYiFhOCiNG85kiebPE3ZUpWNjiuEsGNOW+h17/Fj5TuNcXfJPEnGt25FmpR5F9N5Rdy7B3h+2RM8/15/utTqUmBTFefPn+XPub9zdM0OLLuO4R6fjF+SK0ubLOFMwBkqu93Hgap++DWsQvNurah2fx2ZRhFCWM3p59C11szYcBg3s4lnIspxavW/rOo8gdQjZ9hWPpGdz6TQ9uludArtRFG3ojZ5z9T0VGK2bWbv6i3E/RNDkt8+zhWJwyUmhabLHwMgxc3ChXImkoK9CepwH506P00x7wCbvL8Qwnndbg7d6Qs9OTWd7tM2Uc7Pi7HP1AbAkpJG9Luz2PflH6h0zVH/C8zouoLq1R+kRsl63Fe9PhWKVyDAKwA/T78bbpSRkp7CuZRzJJ0+xbHYfRz99yCJacdJ8T3K6cOHCf4qEL8zLnik//8vP1vrxKPaehIUUJO0f4MIfawhj7Z8EJPJlO/jF0I4lzwXulKqOfA5YAamaq1HZVuvrq1vCSQD3bXWW273mgVV6AAZGZqr6Rl4upnZmXCesUt3M+SJGpRMSWXvtH84EbeZvZU2EJcUR+2ZtXFN8eKSZxoprumkuWRwouwJdoRtw6Q1T/3SGu9kd9ws/1/yeyvGE/tgDMV9y1B+fh0uFveg0v0VCAqvRpmIOlSqWRFXs5S3ECLv8nSUi1LKDEwEmgAJwCal1HytdWyWzVoAla/9uR+YdO2/dsFkUni6ZRbw6ctXOXwmGR9PV4oEFmV7kzDGZHiyqu/b+FquMi11Dvs3xVK1qBnL5RQunEvmTLoHHWvfCy4mzsQo9qcqwu4tQ7FyJbjk50fl8uUZ0rY+mEyZP9aEEMIA1nwoWg+I11rvB1BKzQLaAFkLvQ3wnc7c3V+vlPJVSpXWWh+3eeI8aly1BI2rlri+XLWUN+3Dg/D29QZTMcq91ZZdu+rzdNtQPFwz9+j3J13i8dDSuJhN6JflphlCCPtkTaGXBY5kWU7g5r3vnLYpC9xQ6EqpSCASIDg4+E6z5ousF/QCaFKjJE1qlLy+HBrkQ2iQz/VlKXMhhL2yZmI3pwbLPvFuzTZoradorcO11uGBgYHW5BNCCGElawo9ASiXZTkIOHYX2wghhMhH1hT6JqCyUqq8UsoN6Ahkv2PEfKCbylQfOG+P8+dCCOHMcp1D11qnK6X6AH+SedjiNK11jFKq17X1UcBiMg9ZjCfzsMUe+RdZCCFETqw69V9rvZjM0s76WFSWrzXQ27bRhBBC3Ak520UIIZyEFLoQQjgJKXQhhHAShl2cSymVCBy6y6cHAEk2jOMIZMyFg4y5cMjLmO/RWud4Io9hhZ4XSqnoW12cxlnJmAsHGXPhkF9jlikXIYRwElLoQgjhJBy10KcYHcAAMubCQcZcOOTLmB1yDl0IIcTNHHUPXQghRDZS6EII4STsutCVUs2VUruVUvFKqf45rFdKqfHX1u9QStUxIqctWTHmZ6+NdYdSaq1SqrYROW0ptzFn2S5CKWVRSrUvyHz5wZoxK6UaK6W2KaVilFIrCzqjrVnxd9tHKbVAKbX92pgd+iJ/SqlpSqlTSqldt1hv+/7SWtvlHzKv7LgPqAC4AduBGtm2aQn8TuYNNuoDG4zOXQBjfgAofu3rFoVhzFm2+4vMi8S1Nzp3AXyffcm8zWPwteUSRucugDEPBEZf+zoQOAO4GZ09D2NuBNQBdt1ivc37y5730K/fy1RrnQr8dy/TrK7fy1RrvR7wVUqVLuigNpTrmLXWa7XWZ68trifzZiKOzJrvM8CrwBzgVEGGyyfWjLkz8KvW+jCA1trRx23NmDXgrTLv81iUzEJPL9iYtqO1/ofMMdyKzfvLngv9VvcpvdNtHMmdjucFMn/CO7Jcx6yUKgu0BaJwDtZ8n6sAxZVSK5RSm5VS3QosXf6wZswTgOpk3u1sJ9BXa51RMPEMYfP+sup66Aax2b1MHYjV41FKPUxmoT+Ur4nynzVj/gx4V2ttcZKbdFszZhegLvAo4AmsU0qt11rvye9w+cSaMTcDtgGPABWBpUqpVVrrC/mczSg27y97LvTCeC9Tq8ajlKoFTAVaaK1PF1C2/GLNmMOBWdfKPABoqZRK11r/ViAJbc/av9tJWuvLwGWl1D9AbcBRC92aMfcARunMCeZ4pdQBoBqwsWAiFjib95c9T7kUxnuZ5jpmpVQw8CvQ1YH31rLKdcxa6/Ja6xCtdQjwC/CKA5c5WPd3ex7QUCnlopTyAu4H4go4py1ZM+bDZP5GglKqJFAV2F+gKQuWzfvLbvfQdSG8l6mVYx4C+ANfXttjTdcOfKU6K8fsVKwZs9Y6Tin1B7ADyACmaq1zPPzNEVj5fR4BTFdK7SRzOuJdrbXDXlZXKfUj0BgIUEolAEMBV8i//pJT/4UQwknY85SLEEKIOyCFLoQQTkIKXQghnIQUuhBCOAkpdCGEcBJS6EII4SSk0IUQwkn8H0CN1PcY5ZT+AAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xxx = np.arange(0, 1 + h, h)\n",
    "yyy = []\n",
    "for i in xxx:\n",
    "    yyy.append((math.e ** (-50 * i) / 3) + i ** 2)\n",
    "    # yyy.append(xxx)\n",
    "pl.plot(xx, yy, label=\"Euler\", linestyle=\":\")\n",
    "pl.plot(xx_r, yy_r, label=\"E_r\", color=\"green\", linestyle=\"-\")\n",
    "pl.plot(xx_R_k, yy_R_k, label=\"R_k\", linestyle=\"--\", color=\"purple\")\n",
    "pl.plot(xxx, yyy, color=\"red\", label=\"right\", alpha=0.3)\n",
    "pl.legend()\n",
    "pl.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 一阶常微分"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

一键复制编辑原始数据按行查看历史

提交于 2022-05-11 22:32 . 呜呜呜幸苦死我了

马建仓 AI 助手

尝试更多

代码解读

代码找茬

代码优化

https://gitee.com/changerxu/numerical-analysis.git

git@gitee.com:changerxu/numerical-analysis.git

changerxu

numerical-analysis

数值分析

master

徐晨哲/数值分析

简介

发行版

贡献者

近期动态

徐晨哲/数值分析 .gitee-modal { width: 500px !important; }

简介

发行版

贡献者

近期动态

搜索帮助

徐晨哲/数值分析