ch1/ex1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "本节介绍了**模型误差,观测误差,截断误差,方法误差,舍入误差,病态问题,条件数,数值稳定性**.\n",
    "详情见[《现代数值分析》 一、引论](https://www.jianshu.com/p/4d57905d401e)\n",
    "\n",
    "本文用到python的sympy库进行符号运算，\n",
    "未安装的同学可以使用\n",
    "```\n",
    "pip install sympy\n",
    "```\n",
    "进行安装.关于sympy教程可以自行搜索,这里给出一份[文档](https://docs.sympy.org/latest/index.html).\n",
    "[这个更加实用](https://geek-docs.com/python/python-tutorial/python-sympy.html)\n",
    "常用命令直接看本文代码应当也不难理解.\n",
    "\n",
    "原创内容,如需转载需征得作者同意。\n",
    "\n",
    "Copyright©2020 lizhemin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1\n",
      "E\n",
      "0\n",
      "pi\n"
     ]
    }
   ],
   "source": [
    "#内置符号说明\n",
    "import sympy as sp\n",
    "#虚数单位\n",
    "print(sp.I**2)\n",
    "#科学常数\n",
    "print(sp.E)\n",
    "#无穷大\n",
    "print(1/sp.oo)\n",
    "#pi\n",
    "print(sp.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I\n",
      "3\n",
      "2\n",
      "3628800\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "#常用运算\n",
    "#开方\n",
    "print(sp.sqrt(-1))\n",
    "#对数运算,写一项时默认以e为基地\n",
    "print(sp.log(1000,10))\n",
    "#开根号\n",
    "print(sp.root(8,3))\n",
    "#阶乘\n",
    "print(sp.factorial(10))\n",
    "#三角函数,此处以sin为例\n",
    "print(sp.sin(sp.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.00000000000000\n",
      "3.00000000000000\n",
      "0.500000000000000\n"
     ]
    }
   ],
   "source": [
    "#表达式与表达式求值\n",
    "import numpy as np\n",
    "x = sp.Symbol('x')\n",
    "fx = 2*x+1\n",
    "print(fx.evalf(subs={x:2}))\n",
    "#导入多个\n",
    "i,j = sp.symbols('i j')\n",
    "gx = i+j\n",
    "sp.pprint(gx.evalf(subs={i:1,j:2}))\n",
    "#从sympy.abc导入\n",
    "from sympy.abc import x,y\n",
    "hx = x/y\n",
    "sp.pprint(hx.evalf(subs={x:1,y:2}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.300000000000000\n",
      "0.30000000000000004\n",
      "exp(2*x) + 1/x\n",
      " 2⋅x   1\n",
      "ℯ    + ─\n",
      "       x\n"
     ]
    }
   ],
   "source": [
    "#一个小例子说明符号运算的好处\n",
    "r = sp.Rational(1/10)\n",
    "val = 3*r\n",
    "sp.pprint(val.evalf())\n",
    "sp.pprint(1/10*3)\n",
    "#pprint可用于美化输出，类似于latex的输出,某些字符maybe需要unicode编码支持\n",
    "sp.init_printing(use_unicode=True)\n",
    "x = sp.Symbol('x')\n",
    "c = (sp.exp(x)**2+1/x)\n",
    "print(c)\n",
    "sp.pprint(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2⋅a⋅b + 3⋅a + 4⋅b\n",
      " 2          \n",
      "x  + 2⋅x + 1\n",
      "tan(x)\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "#表达式操作\n",
    "#自动规范表达式\n",
    "from sympy.abc import a, b\n",
    "expr = b*a + -4*a + b + a*b + 4*a + (a + b)*3\n",
    "sp.pprint(expr)\n",
    "#展开表达式\n",
    "expr = (x + 1) ** 2\n",
    "sp.pprint(sp.expand(expr))\n",
    "#简化表达式\n",
    "expr = sp.sin(x)/sp.cos(x)\n",
    "sp.pprint(sp.simplify(expr))\n",
    "#比较表达式\n",
    "a = sp.cos(x)**2 - sp.sin(x)**2\n",
    "b = sp.cos(2*x)\n",
    "print(a.equals(b))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "29\n",
      "[0, 1]\n",
      "<class 'list'>\n",
      "x + 1 = 4\n",
      "[3]\n",
      "{1}\n"
     ]
    }
   ],
   "source": [
    "#通过替换值来求表达式\n",
    "from sympy.abc import a, b\n",
    "expr = b*a + -4*a + b + a*b + 4*a + (a + b)*3\n",
    "sp.pprint(expr.subs([(a, 3), (b, 2)]))\n",
    "#求解方程\n",
    "x = sp.Symbol('x')\n",
    "sol = sp.solve(x**2 - x, x)\n",
    "print(sol)\n",
    "#或者写成公式的形式\n",
    "eq1 = sp.Eq(x + 1, 4)\n",
    "sp.pprint(eq1)\n",
    "sol = sp.solve(eq1, x)\n",
    "print(sol)\n",
    "#给定区间的解\n",
    "sol = sp.solveset(x**2 - 1, x, sp.Interval(0, 100))\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SeqFormula(x, (x, 1, 10))\n",
      "[1, 2, 3, 4, …]\n",
      "[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]\n",
      "10\n",
      "1\n",
      "55\n"
     ]
    }
   ],
   "source": [
    "#序列\n",
    "from sympy.abc import x\n",
    "s = sp.sequence(x, (x, 1, 10))\n",
    "print(s)\n",
    "sp.pprint(s)\n",
    "print(list(s))\n",
    "\n",
    "print(s.length)\n",
    "print(s.start)\n",
    "print(sp.summation(s.formula, (x, s.start, s.stop)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "oo\n"
     ]
    }
   ],
   "source": [
    "from sympy.abc import x\n",
    "\n",
    "l1 = sp.limit(1/x, x, sp.oo)\n",
    "print(l1)\n",
    "\n",
    "l2 = sp.limit(1/x, x, 0)\n",
    "print(l2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[1, 2], [3, 4], [0, 3]])\n",
      "⎡1  2⎤\n",
      "⎢    ⎥\n",
      "⎢3  4⎥\n",
      "⎢    ⎥\n",
      "⎣0  3⎦\n",
      "---------------------------\n",
      "M * N\n",
      "---------------------------\n",
      "⎡6 ⎤\n",
      "⎢  ⎥\n",
      "⎢14⎥\n",
      "⎢  ⎥\n",
      "⎣6 ⎦\n"
     ]
    }
   ],
   "source": [
    "#矩阵运算\n",
    "M = sp.Matrix([[1, 2], [3, 4], [0, 3]])\n",
    "print(M)\n",
    "sp.pprint(M)\n",
    "\n",
    "N = sp.Matrix([2, 2])\n",
    "\n",
    "print(\"---------------------------\")\n",
    "print(\"M * N\")\n",
    "print(\"---------------------------\")\n",
    "\n",
    "sp.pprint(M*N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEFCAYAAAAYKqc0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVXX+x/HXYV9EQFlkFRAUxAUVxdzNJTPHLdNsw9Q0\nZ/ple041v2mqSbNpsX0sM7PSMSt1csvd3MV9A1kVEAHZZYd7fn9g/sxQUe+95y6f5+PBQ7nce87b\nw/XNl7N8j6KqKkIIISyLjdYBhBBC6J+UuxBCWCApdyGEsEBS7kIIYYGk3IUQwgJJuQshhAWSchdC\nCAsk5S6EEBbITusAwrIpiuID9Ab8gUrgOJCgqqpO02BCWDhFrlAVhqAoykBgFtACOATkAU5AW6AN\nsBx4R1XVUs1CCmHBpNyFQSiK8jbwoaqqZxv5mh0wArBVVfUHo4cTwgpIuQshhAWSA6rCoBRFWawo\nivsVn4coirJJy0xCWAMpd2FoO4C9iqIMVxTlMeAX4H2NMwlh8Qy5W0b29wgAduzYwcCBA/Hy8uLQ\noUO0atVK60hCmDJFHwuRkbswqMWLFzN58mS+/vprJk2axPDhwzly5IjWsYSweDJyFwY1evRo5s+f\nj4+PDwD79u1j+vTpHDp0SONkQpgsvYzcpdyF0dXU1ODg4KB1DCFMleyWEabrjTfeoLCwsNGvOTg4\nsHnzZn7++WcjpxLCesj0A8IgOnbsyJ/+9CecnJzo2rUr3t7eVFVVkZyczOHDhxk8eDAvvfSS1jGF\nsFiyW0YYxMMPP8zixYuZO3cuPj4+5OTk4OzsTFRUFP369cPZ2VnriEKYKr3slpGRuzCIAwcOcObM\nGb799lu2bNnyu69VVlZet9wnT57Mzz//jI+PD8ePHwegsLCQCRMmkJGRQUhICMuWLcPT0xNVVZk5\ncyZr1qzBxcWFr776iq5duxr03yaEOZB97sIgHn/8cYYNG0ZiYiKxsbGXP7p160ZsbOx1Xztp0iTW\nrVv3u8fmzJnDoEGDSE5OZtCgQcyZMweAtWvXkpycTHJyMvPnz2fGjBkG+zcJYU4Mtlumrl6n1qsq\njna2Blm+MA8zZszg008/venXZWRkMGLEiMsj93bt2rF161b8/PzIyclhwIABJCUlMX36dAYMGMDE\niRP/8DwhzE1lTT3ODramfbZMv7lbWLL3DxMCCitzK8XemNzc3MuF7efnR15eHgDZ2dkEBQVdfl5g\nYCDZ2dmNLmP+/PmXf4OIjo7WSy4h9Gna4gS9Lctg5e7r7sSi3WfQ6eS4qjCcxn7zVJTGBz7Tpk0j\nISGBhIQEOaArTE5ybhm/Jl/Q2/IMVu6TeoWQfqGcbcn5hlqFsCK+vr7k5OQAkJOTc/mK18DAQDIz\nMy8/LysrC39/f00yCnE7Fu3OwMFOf5VssHK/u4Mf3m6OfLUzw1CrEFZk5MiRLFq0CIBFixYxatSo\ny49//fXXqKrKnj17cHd3l/3twuyUVNbyw4FsRnXW38DEYOXuYGfDQ3Gt2XY6n9T8i4ZajbBAEydO\n5I477iApKYnAwEAWLFjArFmz2LBhAxEREWzYsIFZs2YBMHz4cMLCwggPD+exxx7jk08+0Ti9EDfv\n+4RMKmvrie8VordlGvQipvyyanrN2cQDPYL5x6gOhlqPEDctNjaWhAT9HbwS4lbV61QG/GsLfs2d\nWfb4HWAOc8t4uznyp07+LD+QRVlVrSFXJYQQZmlzYh6ZhZVM6h2i1+Ua/CKm+F4hlNfUs/xAlqFX\nJYQQZuerXen4uTsxtL2vXpdr8HLvHORB12APFu3KkNMihRDiCqdzy9iZUsDDd7TGzla/dWyU6Qcm\n9Q4lo6CCbafltEghhPjNV7sycLSz4f7uwXpftlHK/e4OrfBt7sjS/Zk3frIQQliB4vIaks6XMTom\ngBau+r95jVHK3d7Whsf7t2FzYi4peWXGWKUQQpi0JfszOXCmSO8HUn9jtFkhR3b2x0ZRWLAjw1ir\nFEIIk1Rbr2PRrgx6h7ckyq+5QdZhtHJv2cyRMV0C+PFgFoXlNcZarRBCmJw1x3I4X1rFlD6hBluH\nUedzn9wnlOo6Hd/tPWPM1QohhMlQVZUFO9IJ83ZlQFsfg63HqOXe1teNvhFefL37DDV1OmOuWggh\nTELCmSKOZpXwaO9QbGz0cjFqo4x+J6YpfULJK6vm56PnjL1qIYTQ3IJf03F3tufergEGXY/Ry71/\nW2/CfZqxYEd6o3NxCyGEpcosrOCXk+d5IC4YFwfD3sLa6OWuKAqTe4dy4lwpe9MLjb16IYTQzMKd\nGdgoCvF3hBh8XZrcIHts1wA8nO1Yeajx26EJIYSlKa2qZVlCJiM6+dHK3cng69Ok3J3sbZnaN4yl\nCZmkyVzvQggrsGx/Jher65jSJ8wo69Ok3AEmdA/G3saGBTvStYoghBBGUVev47u9Z+kR0oKOge5G\nWadm5e7t5sjYrgEsP5BFwcVqrWIIIYTBrTtxnsKKGh4fYJxRO2hY7gBT+zZc1LR4j1zUJISwTKqq\nMn97Gp4uDvQ34EVLV9O03MN93BgU6cPi3Weoqq3XMooQQhjEnrRCjmaVMLVvKLYGvGjpapqWO8Bj\n/cIoKK/hx4Ny5owQwvJ8/msaLV0duLdroFHXq3m5x4W2oFOgO1/8miZ3ahJCWJTk3DI2J+YR3ysE\nJ3tbo65b83JXFIXH+oaRdqGcTYl5WscRQgi9mb89DSd7Gx7q2dro69a83KHhTk0BHs58vj1N6yhC\nCKEXuaVVrDiczfjYIIPcaelGTKLc7WxtmNInlKTcUg5lFmkdRwghbttXuzKo16lMNdJFS1cziXIH\nGB8biKOdLZ9tTdU6ihBC3JaL1XV8s+cMd3fwI7iliyYZTKbcmznZMz42iF9O5pKSJ1MSCCHM19J9\nZymrqmNaP21G7WBC5Q4wqXcIDrY2zN8uo3chhHmqrdexJ62AQZE+dA7y0CyHSZW7VzNHxscG8dOh\nbM6XVGkdRwghbtrKw+fYeCpPkzNkrmRS5Q7wWN8w6nUqX+6UCcWEEOZFp1P5bFsqka3cGNDOW9Ms\nJlfuwS1dGNHJn2/3nKGkolbrOEII0WQbTzUcM5wxoA2KYrypBhpjcuUOML1/GOU19XyzVyYUE7/3\n3nvvER0dTYcOHZg4cSJVVVWkp6cTFxdHREQEEyZMoKamRuuYwgqpqsonW1MJbuHCPR39tI5jmuUe\n7e9O/7beLNyZLhOKicuys7P54IMPSEhI4Pjx49TX17N06VJefPFFnn76aZKTk/H09GTBggVaRxVW\naE9aIYczi5nWLww7W+2rVfsE1/B4/zZcuFjD8gNZWkcRJqSuro7Kykrq6uqoqKjAz8+PzZs3M27c\nOADi4+NZsWKFximFNfpkawpezRwZ1824E4Rdi8mWe8+wFoztEsD3CZnU1eu0jiNMQEBAAM899xzB\nwcH4+fnh7u5Ot27d8PDwwM6u4U7ygYGBZGfLDKPCuI5nl/Br8gWm9Ak1+gRh12Ky5a4oCnd39ONI\nVgmrjpzTOo4wAUVFRaxcuZL09HTOnTtHeXk5a9eu/cPzrnUga/78+cTGxhIbG0t+fr6h4wor8unW\nVNwc7XiwZ7DWUS4z2XIHGBTpQ2QrNz7ZmirTAQs2btxIaGgo3t7e2NvbM3bsWHbt2kVxcTF1dXUA\nZGVl4e/v3+jrp02bRkJCAgkJCXh7a3uamrAc6RfKWXM8h4fuaE1zJ3ut41xm0uVuY6Pw54HhpORd\nZP2J81rHERoLDg5mz549VFRUoKoqmzZton379gwcOJDly5cDsGjRIkaNGqVxUmFNFu/OIC6kBZN7\nh2od5XdMutwB7unoR6iXKx9tSUFVZfRuzeLi4hg3bhxdu3alY8eO6HQ6pk2bxltvvcW7775LeHg4\nBQUFTJkyReuowkpkFVXw9e4ztGvlhrebo9ZxfkcxYGHqbcHLEjJ5YflRFk7qzsBI491gVliu2NhY\nEhIStI4hzNwrK47xn/2ZbHt+IP4ezvparF6ufjL5kTvAmC4BBHg4y+hdCGEyckurWLY/i3HdAvVZ\n7HpjFuVub2vD9P5hHDhTxJ60Qq3jCCEE/96WRr2qMqN/uNZRGmUW5Q4wPjYIr2aOfLwlResoQggr\nd+FiNd/tO8OoGH/NbsZxI2ZT7k72tkzrF8qOlAsczizWOo4Qwop98Ws61XU6/jLQNEftYEblDvBg\nXGs8XOz5aHOy1lGEEFaqqLyGxbszuKejH228m2kd55rMqtxdHe2YOSicE+dKOZ5donUcIYQVWrgr\ng/Kaep6403RH7WBm5Q4wtmsQF6vr+GCTjN6FEMZVWlXLwp3p3BXtS2Sr5lrHuS6zK3d3Z3sm9w7l\nl5O5nDgno3chhPF8vSuDsqo6nhgYoXWUGzK7cgeY3DsUN0c7Gb0LIYymvLqO3WmFDGznTcdAd63j\n3JBZlru7iz2P9gll/YlcTuWUah1HCGEFFu3OYGfKBZ680/RH7WCm5Q4wRUbvQggjKauqZf72NAa0\n86ZLa0+t4zSJ2Za7u4s9j/YOYe3x8zJ6F0IY1KJdGRRX1PL04LZaR2kysy13gMl9QmnmaMeHct67\nEMJASqtq+fzXdAZF+tA5yEPrOE1m1uXu4eLApF4hrDt+nsTzMnoXQujfVzszKKms5SkzGrWDmZc7\nwJQ+IbT1deP9DTJ6F0LoV0llLZ//msbgKF+zOEPmSmZf7p6ujtwV3Yp1J85zLEvOexdC6M+XO9Ip\nq6rjqcHmcYbMlcy+3AGm9A3F3dmedzckaR1FCGEhSipq+XJHw9WoHQLMa9QOFlLuzZ3smd4/jC1J\n+Rw4I/O9CyFu34IdaZRV15ndvvbfWES5A0zqFYJXMwfe+eW01lGEEGauqLyGA2eLuLdrAFF+pj2H\nzLVYTLm7ONgxY0A4u1IL2JVyQes4Qggz9um2VHalFjC9fxuto9wyiyl3gAfjgmnV3Il3NpyWe60K\nIW7J+ZIqFu3KYEyXANr6umkd55ZZVLk72dvyxJ3hHDhTxNbT+VrHEUKYoQ82J6NTVbO6GrUxFlXu\n0HCv1UBPZ975JUlG70KIm5JxoZxl+zOZ2COYoBameW/UprK4cnews2HmoAiOZ5ey/kSu1nGEEGbk\n3Q2nsbe1Mfm7LDWFxZU7wJguAYR5ubDycDb1Ohm9CyFu7OS5UlYdOcejvUPwcXPSOs5ts8hyt7O1\n4flhkaw9fp6fDmVrHUcIYQbe+SWJ5k52TO9nvmfIXMkiyx1gWHQrOgW68+4vSVTV1msdRwhhwhIy\nCtmUmMf0/m1wd7HXOo5eWGy5K4rCrGGRnCupYvHuM1rHEUKYKFVVmbs+Ca9mjjzaO0TrOHpjseUO\n0Cvci35tvfl4awollbVaxxFCmKDtyRfYl17I/9wZjouDndZx9Maiyx3gxWHtKK6o5d/bUrWOIoQw\nMTqdytvrEwn0dGZij2Ct4+iVxZd7tL87o2P8+XJnOudLqrSOI4QwIWuP5+DqYMczQ9riYGdZdWhZ\n/5preHZoO+p1KvM2yaRi5q64uJhx48YRGRlJVFQUu3fvprCwkCFDhhAREcGQIUMoKirSOqYwA9V1\n9by1LomSylpGxQRoHUfvrKLcg1q48FDP1vxnfyYpeRe1jiNuw8yZMxk2bBiJiYkcOXKEqKgo5syZ\nw6BBg0hOTmbQoEHMmTNH65jCDCzefYazhRX8dXgUtjaK1nH0zirKHeCJgQ0HS95en6h1FHGLSktL\n2b59O1OmTAHAwcEBDw8PVq5cSXx8PADx8fGsWLFCy5jCDJRU1PLh5hT6RnjRv6231nEMwmrKvWUz\nR6b1C2P9iVwOnJFf281RWloa3t7ePProo3Tp0oWpU6dSXl5Obm4ufn5+APj5+ZGXl9fo6+fPn09s\nbCyxsbHk58vEctbs460plFbV8tLwKK2jGIzVlDvAlD6heDVzZM7aUzKpmBmqq6vj4MGDzJgxg0OH\nDuHq6npTu2CmTZtGQkICCQkJeHtb5mhN3FhmYQVf7czg3q6BZnsjjqawqnJ3dbTjhbvaUVhew4aT\nMqmYuQkMDCQwMJC4uDgAxo0bx8GDB/H19SUnJweAnJwcfHx8tIwpTNzb65OwsYFnh5r3lL43YlXl\nDjC2awCKojB7bSI1dTqt44ib0KpVK4KCgkhKargR+qZNm2jfvj0jR45k0aJFACxatIhRo0ZpGVOY\nsKNZxaw6co6pfcLwc3fWOo5BWc7lWE1kZ2vDy8OjePSr/Xyz5wyT+4RqHUnchA8//JAHH3yQmpoa\nwsLCWLhwITqdjvHjx7NgwQKCg4P5/vvvtY4pTJCqqvxz9SlaujowvX+Y1nEMzurKHWBAO2/6Rngx\nb1MyY7sG4OHioHUk0UQxMTEkJCT84fFNmzZpkEaYk42n8tibXsjrozvg5mQZk4Ndj9XtloGGScVe\nvieKsqpaPtiUonUcIYSB1dXrmLP2FGHertzfPUjrOEZhleUOENmqORO6B/P17gzS8uXCJiEs2X8S\nMknNL2fWsEjsba2j9qzjX3kNzwxpi6OdDXPWyoVNQliqkspaPtiYzITuQQxp76t1HKOx6nL3dnPk\nzwPD+eVkLrtTC7SOI4QwgA83JZN3sZqHe7ZGUSxvmoFrsepyh4YLmwI8nHlj9Qnq6+XUSCEsSWr+\nRb7alcH93YPoEOCudRyjsvpyd7K35eV7IqmpU1l+MEvrOEIIPXrj55M429vy7NB2WkcxOqsvd4C7\nO/jh4WLP3EvTfwohzN+WxDy2JOUzc3AEXs0ctY5jdFLuNJwa+erIaIoqanh/o8z5LoS5q6nT8frq\nk4R5ufLIHSFax9GElPsl0f7uTOwRzNe7z3A6t0zrOEKI29BwinM5r4yIsrg7LDWVdf6rr+G5oe1o\n5mjHq6tOyKyRQpipCxermbcpmf5tvRnYznonkZNyv4KnqwPPDW3LrtQC1h0/r3UcIcQt+GpnBvU6\nlb+NiLKqUx+vJuV+lYk9gols5cYbq09RWVOvdRwhxE04nFnMx1tTmNQrhHAfN63jaErK/Sp2tja8\nOjKa7OJKPtuWqnUcIUQT1etU/nflcbybOTJjQBut42hOyr0RPcNa8qfO/ny2LZXMwgqt4wghmmDp\n/rMczSrh5XuirGLWxxuRcr+Gv94dSctmDnz+a5rWUYQQN1BYXsPb65OIC23ByM7+WscxCVLu1+Dv\n4czDPUP4evcZtiQ2fsNlIYRpeHt9ImVVdbw+uoNVH0S9kpT7dUzpE0qETzP+d9VxObgqhIk6dLaI\npfszmdw7hLa+1n0Q9UpS7tfhYGfD66M7kFlYycdb5KYeQpiahoOoJ/Bxc2TmYMu+4fXNknK/gZ5h\nLRnbJYB/b08lJU9u6iGEKVmy7yxO9ra8dHcUzRyt8q6h1yTl3gQv3ROFs70t/7vyuFy5KoSJyC2t\n4q21iTjZ2zAyRg6iXk3KvQm8mjnywrBIdqUWsOrIOa3jCCGA1/57kup6Ha+PkoOojZFyb6IHegTT\nOciDxbvPUFIh0wILoaXNibmsPpbDk3eGE+LlqnUckyTl3kQ2Ngqzx3TkSFYxc9ad0jqOEFaroqaO\nv604QbhPM6b1kytRr0XK/Sa092/Oo71DWbIvkz1pcs9VIbQwb2My2cWVvDmmo9VO59sUsmVu0tOD\n2xLUwpmXfjxGVa2c+y6EMZ04V8IXO9K5v3sQPUJbaB3HpEm53yRnB1veHNORtAvlfLRZzn0Xwljq\n6nX8Y9UJ+oZ7MevuSK3jmDwp91vQN8KbsV0D+GxbKonnS7WOI4RV+HJnOvsyiri3WyAeLg5axzF5\nUu636JV72tPc2Z5ZPxyjXifnvgthSOkXynnnl9MMjvJlRCc/reOYBSn3W9TC1YH/HdGe0+dLWZaQ\nqXUcISyWTqcy64ejONjZ8M8xck57U0m534ZRMf4MivLltf+e5ExBudZxhLBI3+07y970Ql65Jwrf\n5k5axzEbUu63QVEUXronCjsbhRd/OIpOds8YRX19PV26dGHEiBEApKenExcXR0REBBMmTKCmpkbj\nhEJfzhVXMmdtIr3DWzI+NkjrOGZFyv02+bk788qIKPakFfLt3jNax7EK8+bNIyoq6vLnL774Ik8/\n/TTJycl4enqyYMECDdMJfVFVlZd+ajimNWdsJ9kdc5Ok3PVgfGwQfSO8mL02UW7LZ2BZWVmsXr2a\nqVOnAg0FsHnzZsaNGwdAfHw8K1as0DKi0JNVh8+RfqGc5+9qR1ALF63jmB0pdz1QFIU593bCRlGY\n9eNRmTnSgJ566inmzp2LjU3DW7egoAAPDw/s7Bqmew0MDCQ7O1vLiEIPckoqeWXlcfzcnYjvFaJ1\nHLMk5a4nAR7OvDQ8ip0pBSzZd1brOBbp559/xsfHh27dul1+rLEfpNf69X3+/PnExsYSGxtLfn6+\nwXKK26OqKi8sP0pdfcPuGFsb2R1zK2R2ez2a2COIY1nFvL8xmd7hXrRuKbPV6dPOnTtZtWoVa9as\noaqqitLSUp566imKi4upq6vDzs6OrKws/P0bn9t72rRpTJs2DYDY2FhjRhc34Zu9Z/k1+QKvj+4g\nMz7eBhm565GiKPzPoAgqa+t5dtkRubhJz2bPnk1WVhYZGRksXbqUO++8k2+//ZaBAweyfPlyABYt\nWsSoUaM0TipuVcaFct5cfYq+EV48FBesdRyzJuWuZ/4ezrw2KpqEM0X8e3uq1nGswltvvcW7775L\neHg4BQUFTJkyRetI4hbU61Se+/4IdrYKc8fJ2TG3SzHgwT+rHbaqqsoT3x3il5PnWfGX3kT7u2sd\nSVwlNjaWhIQErWOIK/x7Wyqz1yby3oTOjOkSqHUcLenlp5qM3A1AURTeGN0BTxcHnv7PYZkaWIgb\nSDpfxju/nOauaF9GxwRoHcciSLkbiKerA3PHdeJ07kX+tT5J6zhCmKyq2nr+8d8TxAR78M8xHWV3\njJ5IuRvQgHY+PNQzmAU709mVekHrOEKYpLfWJbIrtYDH+4fh1cxR6zgWQ8rdwF4aHkWETzMW7EiX\nG2sLcZUtiXks3JnBpF4h3Bnpq3UciyLlbmAuDnb8a1xntiXl8+IPcvWqEL/JK63iue+PENnKTe6s\nZABS7kbQKciDF4a1Y92J83y7V65eFUKnU3n2+yOU19Tx4cQuONnbah3J4ki5G8nUPmH0a+vN6z+f\nJOl8mdZxhNDUFzvS+DX5An8b0Z4IXzet41gkKXcjsbFReOe+zrg52fPEdweprJHTI4V1OpZVwtvr\nkxgW3YoHeshVqIYi5W5E3m6OvDu+M8l5F3nt55NaxxHC6Mqr63hy6SG8mjky51457dGQpNyNrF9b\nb6b3D2PJvrOsPpqjdRwhjGrO2kQyCsp5b0IMHi4OWsexaFLuGnhuaDs6B3nw/qYkzsq9V4WV+D4h\nkx8PZvHXuyPpGdZS6zgWT8pdA/a2Nnw0MYbyqnr+/N1BmZ5AWLxTOaX8beVxOgd5MKVPmNZxrIKU\nu0aCWrjy2qgOHM8u5R//PaF1HCEMprSqlj9/e5DmTvbMu7+L3HzDSKTcNTS4vS9/GdiGJfsyWZaQ\nqXUcIfROVVVeXH6Us4UVfPRAV7zdZHoBY5Fy19gzQ9rRO7wlf1txnBPnSrSOI4Refbkzg7XHz/Pi\nsHb0CG2hdRyrIuWuMVsbhXn3d8HTxYEZ3xykpFLmnxGW4UBGIbPXnGJoe18e6yv72Y1Nyt0EeDVz\n5OMHu3KuuJJnlx1GJ7fnE2Yur7SKv/50jN7hXrx9X2c5n10DUu4moltrT165J4qNp/L4aleG1nGE\nuGXVdfVM/+YAWUWV/HV4JO7O9lpHskp2WgcQ/y++VwjnSqp4ffVJglu4MLi9TIEqzIuqqvxtxXEO\nnS3m0we7EtmqudaRrJaM3E2Ioig8PbgtHfzdmbn0kEwwJszO17vPsCwhiyfvDOfujn5ax7FqUu4m\nxtnBls8ficXF0Y6pX++nqLxG60hCNMmu1Au89vNJBkf58tTgtlrHsXpS7iaolbsT8x/uRm5pNTO+\nPUBtvU7rSEJcV2ZhBX/59iChXq68N6EzNnKhkuak3E1Ul2BP3rq3I3vSCuUKVmHSyqvreGbZEep1\nKp8/EoubkxxANQVyQNWEjekSSOL5Mv69LY12rZrzcM/WWkcS4nfqdSpPLjlEdnHDFaihXq5aRxKX\nyMjdxL1wVyR3Rvrw8eYUdqZc0DqOEL/z+s8n2ZSYx4z+bejX1lvrOOIKUu4mruEK1hhCvVx5fPEB\nEs+Xah1JCAAW7kznq10ZTO0TysN3hGgdR1xFyt0MuDnZ8874zrg42vLowv3klFRqHUlYuQ0nc3nt\n55PcFe3LX4dHaR1HNELK3Uz4ezizcFIPyqrqeHThfsqqZA4aoY1jWSU8ueQQnQLceX+CTOFrqqTc\nzUh7/+Z88mBXUvIu8udvD8opksLoMgsreHdDEi1cHfg8PhZnB1utI4lrkHI3M/3aejN7bEd+Tb7A\nrB+OoaoyyZgwjvyyah5asJcDZ4pY9Gh3fNyctI4krkPK3QzdFxvEU4MjOFNQztx1SVrHEVagtKqW\n+C/3kVdazcJHexDu66Z1JHEDUu5mauagCCL93Ph0WyqfbE3ROo5RZGZmMnDgQKKiooiOjmbevHkA\nFBYWMmTIECIiIhgyZAhFRUUaJ7UsVbX1TF2UwOncMj59qCvdWntqHUk0gZS7mVIUhddGdmB0jD9z\n1yWxeHeG1pEMzs7OjnfeeYdTp06xZ88ePv74Y06ePMmcOXMYNGgQycnJDBo0iDlz5mgd1WLU1et4\n4rtD7M8o5J3xnRnQzkfrSKKJpNzNmI2Nwtv3dWZwlC9/W3mCHw9maR3JoPz8/OjatSsAbm5uREVF\nkZ2dzcqVK4mPjwcgPj6eFStWaBnTYqiqyqwfj7HxVC6v/imaUTEBWkcSN0HK3czZ29rw0QNd6NWm\nJc8vP8q64+e1jmQUGRkZHDp0iLi4OHJzc/Hza5he1s/Pj7y8vEZfM3/+fGJjY4mNjSU/P9+Ycc2O\nqqq8t+E0q4/mMHNQBPG9QrSOJG6SlLsFcLJvmCa4Y4A7Ty45xK/Jll1cFy9e5N577+X999+nefOm\n3wxi2rQpvTftAAASeUlEQVRpJCQkkJCQgLe3XCp/LaqqMnttIh9sTuGxfqE8NThC60jiFki5WwhX\nRzu+erQ7Yd6uzFmTyN60Aq0jGURtbS333nsvDz74IGPHjgXA19eXnJwcAHJycvDxkf3Ct0pVVeau\nT2L+9jQe7tmapwe3lfufmikpdwvi4eLA15N7YGurMGnhfnalWtZEY6qqMmXKFKKionjmmWcuPz5y\n5EgWLVoEwKJFixg1apRWEc3eextO8+nWVB6IC+YfI6Ol2M2YYsCLYOTqGo3kl1XzwOd7yCyq4ItH\nutMnwkvrSHqxY8cO+vbtS8eOHbGxaRiXvPnmm8TFxTF+/HjOnj1LcHAw33//PS1atLjusmJjY0lI\nSDBGbLMxb2My7208zf3dg3hzTEe54YZ29LLhpdwtVMHFah78Yi/pF8qZ/0gs/WU61t+Rcv9/qqry\n3sbTbD+dT7iPG3Pv7STFri29bHzZLWOhWjZz5LvHetLGuxmPLUpgc2Ku1pGECVJVlX+uPsUHm1KI\nbNWct6TYLYaUuwVr4erAd4/F0a6VG9MXH2DDSes4TVI0jU6n8vKK43yxI51JvUJ4c0xHmeHRgki5\nWzgPFwe+mRrHXe19eeo/h1m676zWkYQJqK2r59nvj/Dd3rPMGNCGv/+pvYzYLYyUuxVwd7bnrXGd\n6R7Sglk/HuO9DadlNkkrVlFTx/RvDnKuuJLnhrblxWGRclaMBZJytxKujnZ8/kgs93ULZN6mZGb9\ncIw6mQ/e6hRcrOaBz/eyNSmPkTH+PHGnXKBkqey0DiCMx97WhrnjOuHn7sQHm1PIv1jNRw90wcVB\n3gbW4GxBBfEL93GuuJLPHurG0OhWWkcSBiQjdyujKArPDG3HP8d0YGtSHhM/30vBxWqtYwkDO5ZV\nzNhPd1FUUcN3j8VJsVsBKXcr9WBcaz57qBuJOaXc++kuzhSUax1JGMj6E+f564/HaO5kx/LH76Bb\n6+tf4CUsg5S7FRsa3YrvHutJcWUtL/90nH3phVpHEnqkqiofb0nh8W8OYGujsGRaT8J95A5K1kLK\n3cp1a+3JTzN6k3+xYcqCRbsy5EwaC1BVW88zy47w9vokRnTy5z/T78C3udzz1JpIuQtCvV35/vE7\nGNDOm7+vOsGz3x+hqrZe61jiFp0rruT574+w4nA2zw5pywf3x+Bkb6t1LGFkcpqEAKC5kz3zH47l\nw80pvLfxNKdzy/jsoW4EerpoHU3chO2n85m59BC19SoL4rtzZ6RMf2ytZOQuLrOxUZg5OIIF8bGc\nuVDByI92sivFsqYNtlQ6ncq8jcnEL9yHj5sTq57oLcVu5aTcxR8MivJl5RO9aenqwPPLj/DBpmS5\n4MmEFZbXMOmr/by38TSjYwL46S+9CPNupnUsoTEpd9GoMO9m/PSX3gxt34p3N5xmwvw9nC2o0DqW\nuMqu1As8tiiBQ2eLeGN0B94d31kuShOAlLu4jmaOdvx9ZDTz7o/hdG4Zd8/bzrKETDmbxgRU19Xz\n5ppTPPjFXooqalj6WE8e6tla5ogRl8nNOkSTZBdX8sx/DrM3vZBh0a2YPbYjnq4OWse6ZeZ8s47T\nuWU8ueQQiefLeDAumJfviZLRumWROzEJ46rXqXzxaxr/+iUJTxcH5o7rxIB25nnQzhzLvbZex/zt\naexIvsDp3DLmjuvEoChfrWMJ/ZNyF9o4ca6Ep5YexsPFgaAWzrw0PAqvZo5ax7op5lbuRzKLefGH\noySeL2NEJz9e/VN7vNzkoiQLJeUutFNVW88nW1L4dFsqLg52vDCsHRO7B5vNDR/Mpdwraup455fT\nLNyZjrebI6+N6sBdMumXpZNyF9pLySvjlRXH2ZNWSEyQB2+M7kCHAHetY92QqZe7qqqsOXae/+zP\nZHtyPg/GBfPi3ZE0d7LXOpowPCl3YRpUVWXF4Wz+ufoUheU1TOsbxuMD2uDhYroHXE253I9nl/Da\nf0+yL6OQyFZu/HNMR7q19tQ6ljAeKXdhWkoqa/nX+iRO5pRyOreMGQPa8GivUJwdTG9eE1Ms9/Ml\nVXy8JYVv9p7B08WB54a2Y0L3ILlptfWRchem6VROKW+vT2JzYh6+zR2ZOagt42MDsbM1ncsqTKnc\nc0ur+HRrKt/tO0tsa0/a+zXnfwZF4O4su2CslJS7MG370guZs/YUB88WM7S9L4OifBjdJQBHO+1H\n8qZQ7nmlVXy6LZXv9p6lTqcyrmsgT9wZTlALmazNykm5C9OnqiobTuby+a9p7M8owtvNkcm9Q3kg\nLljTkamW5Z5ZWMEPB7P4dGsqdTqVsV0C+J87IwhuKaUuACl3YU5UVWVXagGfbUvl1+QLNHO044G4\nYB7u2VqTkaqxy/23f//CnRlsSszF182JvhFe/GVgOCFerkbLIcyClLswT8ezS5i/PY1dKRcoq64j\nLqwl42MDGdLe12i7bIxV7sUVNaw/kcuCHWmczr1IS1cHHogL5sG41rRyl4uQRKOk3IV5O1dcwdJ9\nmSw/kMW5kio8XOwZHRPAuG6BRPs3N+gkWIYs98qaejaeymXl4Wy2nc6nhasDPm5OxPcKYUQnP7kr\nkrgRKXdhGep1KjtTLrAsIZNfTuQSE+TO+dJqhrT3ZXCUL91DPPV+po2+y72ypp5fk/NZe/w8v5w4\nT3lNPa2aOzEyxp+Rnf0N/sNKWBQpd2F5iitq2Hgql9VHc9iZWkBNnQ4PF3tGdvKnc5AHd7Rpib+H\n8x9et27dOmbOnEl9fT1Tp05l1qxZ113P7Za7TqeSeL6MXakXSM2/yI8Hs6mu09HBvzkdAtwZFRNA\nj9AWco66uBVS7sKylVfXsf10PhtO5lJcWcPmxHwAhrb3xdnBls6BHsQEuRPu7UqXju3ZsGEDgYGB\ndO/enSVLltC+fftrLvtmyl1VVfJKqziZU8bRrBIOZRZRVF7DkawSAIa09yXQ05lBkQ2/ZTjKbhdx\ne6TchfXQ6VROnS9lT1ohGRfK2Xgql5ySKuJCW7A3vRDHijxG9euGn7sTB3duwcmmlmkPT6SFqyNu\nTrY0c7TD2cHu8q6R38pdp9NRXlNPWVUdxRU1FFXUUHCxlsyictIvVFCv07ElKR/vZo4k513E3dme\nVs2d6BbiQddgT3qHe+Hn/sffJIS4DXopd5nhX5gFGxuFaH93ov0bJiV7fXSHhtH0uVJaVOdyMK2a\nmrp6Dp4tIl3xpbSsgtL1SaiKgo3S8FuAk70t9jYKzg52ZBSUM/HzPVRU16GgYGurUFVbj52NgqOd\nLfszCgnwdKZLkAd3d/AjJsidkJautPdrjptcOSrMgMFG7tHR0aqzs+mPaPLz8/H29tY6xg1Jzmsr\nKiqitLSU1q1bN2S4UEh5ZSXePr7U6VR0qkq9TuXixXLKy8sBqC0vISA0AhSwQcHWpuEHiK2Ngr2N\nDfa2NpjC8U/5vuuPOWQEOHDgwAlVVTvc9oJUVTXIR7du3VRzIDn1S4ucu3btUocOHXr58zfffFN9\n8803r/sa2Z76ZQ45zSGjqqoqkKDqoYNNZyYnIW5R9+7dSU5OJj09nZqaGpYuXcrIkSO1jiWEpmSf\nuzB7dnZ2fPTRR9x1113U19czefJkoqOjtY4lhKYMVu7Tpk0z1KL1SnLql1Y5hw8fzvDhw5v8fNme\n+mUOOc0h4yXz9bEQORVSCCFMi14O5cs+dyGEsEC3Ve6KotynKMoJRVF0iqLEXvm12bNnEx4eTrt2\n7Vi/fn2jr09PTycuLo6IiAgmTJhATU3N7cRpkgkTJhATE0NMTAwhISHExMQ0+ryQkBA6duxITEwM\nsbGxjT7HkF599VUCAgIuZ12zZk2jz1u3bh3t2rUjPDycOXPmGDklPP/880RGRtKpUyfGjBlDcXFx\no8/TanveaPtUV1czYcIEwsPDiYuLIyMjw2jZADIzMxk4cCBRUVFER0czb968Pzxn69atuLu7X34v\nvPbaa0bN+JsbfQ9VVeXJJ58kPDycTp06cfDgQaNnTEpKurydYmJiaN68Oe+///7vnqPV9pw8eTI+\nPj506PD/ZzkWFhYyZMgQIiIiGDJkCEVFRY2+VlGUeEVRki99xDdphbdzqg0QBbQDtgKxVzzevlOn\nTmpVVZWalpamhoWFqXV1dX845ee+++5TlyxZoqqqqk6fPl395JNP9HlG0Q0988wz6j/+8Y9Gv9a6\ndWs1Pz/fqHmu9Pe//119++23r/ucuro6NSwsTE1NTVWrq6vVTp06qSdOnDBSwgbr169Xa2trVVVV\n1RdeeEF94YUXGn2eFtuzKdvn448/VqdPn66qqqouWbJEHT9+vFEznjt3Tj1w4ICqqqpaWlqqRkRE\n/CHjli1b1HvuuceouRpzo+/h6tWr1WHDhqk6nU7dvXu32qNHDyOm+6O6ujrV19dXzcjI+N3jWm3P\nbdu2qQcOHFCjo6MvP/b888+rs2fPVlVVVWfPnv3b/5+re7YFkHbpT89Lf/e8+nlXf9zWyF1V1VOq\nqiY18qVR999/P46OjoSGhhIeHs6+ffv+8ENl8+bNjBs3DoD4+HhWrFhxO3FuiqqqLFu2jIkTJxpt\nnfq2b98+wsPDCQsLw8HBgfvvv5+VK1caNcPQoUOxs2s4Lt+zZ0+ysrKMuv7racr2WblyJfHxDQOh\ncePGsWnTpt/+QxmFn58fXbt2BcDNzY2oqCiys7ONtn59WrlyJY888giKotCzZ0+Ki4vJycnRLM+m\nTZto06bN5YvbtNavXz9atGjxu8eufP9dpwPvAjaoqlqoqmoRsAEYdsMV3qj9m/LBH0fuHwEPXfH5\nAmDcVa/xAlKu+DwIOK6PPE3M3I/rXCwApAMHgQPANGPlumL9rwIZwFHgSxr5SQ2MA7644vOHgY+M\nnfWK9f/3yu+71tuzKdsHOA4EXvF5KuCl0fYLAc4Cza96fABQABwB1gLRGuW77vcQ+Bnoc8Xnm67s\nBQ3yfgk80cjjmm3PS9/j41d8XnzV14saec1zwCtXfP434LkbreuGp0IqirIRaNXIl15WVfVaw8TG\njvZePRxqynNuSRMzTwSWXGcxvVVVPacoig+wQVGURFVVt+sjX1NyAp8Cr9OwTV4H3gEmX72IRl6r\n92FnU7anoigvA3XAt9dYjMG3ZyM0fR/eDEVRmgE/AE+pqlp61ZcPAq1VVb2oKMpwYAUQYeyM3Ph7\naBLbEkBRFAdgJPDXRr5sKtuzqW5pu96w3FVVHXwLYbJoGIn/JhA4d9VzLgAeiqLYqapad43n3JIb\nZVYUxQ4YC3S7zjLOXfozT1GUn4AegF7LqKnbVlGUz2kYFV2tKdv5tjVhe8YDI4BB6qWhRSPLMPj2\nbERTts9vz8m69L5wBwoNnOt3FEWxp6HYv1VV9cerv35l2auqukZRlE8URfFSVfWCMXM24XtolPdj\nE90NHFRVNffqL5jK9rwkV1EUP1VVcxRF8QPyGnlOFg2/bfwmkIa9JddlqFMhVwH3K4riqChKKA0/\nFX+30/1SCWyh4VdngHjAWDuMBwOJqqo2uoNYURRXRVHcfvs7MJSGX9+N5tI3+jdjrrH+/UCEoiih\nl0Yq99Ow7Y1GUZRhwIvASFVVK67xHK22Z1O2zyoa3nvQ8F7cfK0fUIagNMxBvAA4parqu9d4TqtL\nz0NRlB40/L8tMFbGS+ttyvdwFfCI0qAnUKKqqlY73a/5m7kpbM8rXPn+u1YHrgeGKoriqSiKJw3b\nvvFTEK90m/uPxtDwU6UayAXWX/G1l2nYf5kE3H3F42sA/0t/D6Oh9FOA7wFHI+33+gp4/KrH/IE1\nV+Q6cunjBA27H4y9v3AxcIyGfe6rAL+rc176fDhw+tK21iJnCpAJHL708Zkpbc/Gtg/wGg0/jACc\nLr33Ui69F8OMvP360PAr9tErtuFw4PHf3qPAE5e22xFgD9BLg+9zo9/Dq3IqwMeXtvUxNNrfDrjQ\nUNbuVzym+fak4YdNDlB7qTenAC1pODaRfOnPFpeeG8vvjxdNvvQeTQEebcr6DHmFqhBCCI3IFapC\nCGGBpNyFEMICSbkLIYQFknIXQggLJOUuhBAWSMpdCCEskJS7EEJYICl3IYQwAYqidFcU5aiiKE6X\nrgg+oShKhxu/8hrLk4uYhBDCNCiK8gYNV007A1mqqs6+5WVJuQshhGm4NAfSfqCKhmkR6m91WbJb\nRgghTEcLoBngRsMI/pbJyF0IIUyEoiirgKVAKA2TBT5xq8u64XzuQgghDE9RlEeAOlVVv1MUxRbY\npSjKnaqqbr6l5cnIXQghLI/scxdCCAsk5S6EEBZIyl0IISyQlLsQQlggKXchhLBAUu5CCGGBpNyF\nEMICSbkLIYQF+j92rDmx5jcOUwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2b403931da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<sympy.plotting.plot.Plot at 0x2b40393a6a0>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#画图\n",
    "from sympy.abc import x\n",
    "from sympy.plotting import plot\n",
    "\n",
    "plot(x**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 设$f(x)=\\ln x$,利用$f$在$x_0=e$的泰勒展开式计算$\\ln 3$的近似值,使误差不超过$10^{-4}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "展开到 0 项，表达式为\n",
      "1.00000000000000\n",
      "第 1 项误差为 0.0986122886681102\n",
      "展开到 1 项，表达式为\n",
      "0.367879441171442⋅x - 2.22044604925031e-16\n",
      "第 2 项误差为 0.00502603484621700\n",
      "展开到 2 项，表达式为\n",
      "                      2                            \n",
      "- 0.0676676416183064⋅x  + 0.735758882342885⋅x - 0.5\n",
      "第 3 项误差为 0.000344416204212772\n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "def taylor(func,num_terms=3,point=0):\n",
    "    sums = 0\n",
    "    for i in range(num_terms):\n",
    "        numerator = func.diff(x,i)\n",
    "        numerator = numerator.evalf(subs={x:point})\n",
    "        denominator = np.math.factorial(i)\n",
    "        sums += numerator/denominator*(x-point)**i\n",
    "    return sp.simplify(sums)\n",
    "x = sp.Symbol('x')\n",
    "exp = sp.log(x)\n",
    "i = 1\n",
    "real = np.log(3)\n",
    "error = 10000\n",
    "while error>1e-3:\n",
    "    exp_tay = taylor(exp,num_terms=i,point=np.exp(1))\n",
    "    print('展开到',i-1,'项，表达式为')\n",
    "    sp.pprint(exp_tay)\n",
    "    result = exp_tay.evalf(subs={x:3})\n",
    "    error = np.abs(result - real)\n",
    "    i += 1\n",
    "    print('第',i-1,'项误差为',error)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.设$y_0=28$，按递推公式\n",
    "$$\n",
    "y_n = y_{n-1} - \\frac{1}{100} \\sqrt{783},n=1,2,\\ldots\n",
    "$$\n",
    "计算到$y_{100}$,若取$\\sqrt{783}\\approx 27.982$，试问计算$y_{100}$将有多大误差？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.017999999999941063\n"
     ]
    }
   ],
   "source": [
    "#没有get到point\n",
    "import numpy as np\n",
    "q = 27.982\n",
    "y=28\n",
    "for i in range(1,101):\n",
    "   y = y-q/100\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3.已知$P(x)=(x-10)^4+0.200(x-10)^3+0.0500(x-10)^2-0.00500(x-10)+0.00100$,用秦九韶法计算$P(10.11)$，计算用3位有效数字。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "泰勒展开从高到低的系数依次为 [ 0.001 -0.005  0.05   0.2    1.   ]\n",
      "秦九韶，保留三位算出值： 0.001\n",
      "高精度计算结果 0.00146761000000000\n",
      "条件数为 1.39290410940236\n"
     ]
    }
   ],
   "source": [
    "import numpy\n",
    "import sympy as sp\n",
    "def shao_expand(func,num_terms=3,point=0):\n",
    "    sums = []\n",
    "    for i in range(num_terms):\n",
    "        numerator = func.diff(x,i)\n",
    "        numerator = numerator.evalf(subs={x:point})\n",
    "        denominator = np.math.factorial(i)\n",
    "        sums.append(numerator/denominator)\n",
    "    return sums\n",
    "def shao_cal(in_x,sums,point):\n",
    "    lenth = len(sums)\n",
    "    result = sums[lenth-1]\n",
    "    result = np.around(result,decimals=3)\n",
    "    for i in range(1,lenth):\n",
    "        result = result*(in_x-point) +sums[lenth-i-1]\n",
    "        result = np.around(result,decimals=3)\n",
    "    return result\n",
    "x = sp.Symbol('x')\n",
    "px = (x)**4 + 0.2*(x)**3 +0.05*(x)**2 -0.005*(x) +0.001\n",
    "coef = shao_expand(px,num_terms = 5,point=0)\n",
    "coef = np.array(coef,dtype='float32')\n",
    "print('泰勒展开从高到低的系数依次为',coef)\n",
    "print('秦九韶，保留三位算出值：',shao_cal(in_x=0.11,sums=coef,point=0))\n",
    "print('高精度计算结果',px.evalf(subs={x:0.11}))\n",
    "def cp(in_x,func):\n",
    "    df = func.diff(x,1)\n",
    "    df = df.evalf(subs={x:in_x})\n",
    "    f = func.evalf(subs={x:in_x})\n",
    "    return np.abs(df*in_x/f)\n",
    "\n",
    "print('条件数为',cp(0.11,px))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4.上题$P(x)$的一个近似值为\n",
    "$$\n",
    "x^4-39.8x^3+594.05*x^2-3941*x+9805.05\n",
    "$$\n",
    "用秦九韶算法求$Q(10.11)$，并求此问题条件数$c_p$.\n",
    "**注：**条件数计算公式$c_p=\\left|\\frac{f'(x)\\cdot x}{f(x)}\\right|$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "泰勒展开从高到低的系数依次为 [ 9.805051e+03 -3.941005e+03  5.940500e+02 -3.980000e+01  1.000000e+00]\n",
      "秦九韶，保留三位算出值： -0.011\n",
      "高精度计算结果 0.00146760999999999\n",
      "条件数为 128.020550418708\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import sympy as sp\n",
    "def shao_expand(func,num_terms=3,point=0):\n",
    "    sums = []\n",
    "    for i in range(num_terms):\n",
    "        numerator = func.diff(x,i)\n",
    "        numerator = numerator.evalf(subs={x:point})\n",
    "        denominator = np.math.factorial(i)\n",
    "        sums.append(numerator/denominator)\n",
    "    return sums\n",
    "def shao_cal(in_x,sums,point):\n",
    "    lenth = len(sums)\n",
    "    result = sums[lenth-1]\n",
    "    result = np.around(result,decimals=3)\n",
    "    for i in range(1,lenth):\n",
    "        result = result*(in_x-point) +sums[lenth-i-1]\n",
    "        result = np.around(result,decimals=3)\n",
    "    return result\n",
    "x = sp.Symbol('x')\n",
    "px = (x-10)**4 + 0.2*(x-10)**3 +0.05*(x-10)**2 -0.005*(x-10) +0.001\n",
    "coef = shao_expand(px,num_terms = 5,point=0)\n",
    "coef = np.array(coef,dtype='float32')\n",
    "print('泰勒展开从高到低的系数依次为',coef)\n",
    "print('秦九韶，保留三位算出值：',shao_cal(in_x=10.11,sums=coef,point=0))\n",
    "print('高精度计算结果',px.evalf(subs={x:10.11}))\n",
    "\n",
    "def cp(in_x,func):\n",
    "    df = func.diff(x,1)\n",
    "    df = df.evalf(subs={x:in_x})\n",
    "    f = func.evalf(subs={x:in_x})\n",
    "    return np.abs(df*in_x/f)\n",
    "\n",
    "print('条件数为',cp(10.11,px))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5.设$f(x)=\\ln\\left(x-\\sqrt{x^2-1}\\right)$，它等价于$f(x)=-\\ln\\left(x+\\sqrt{x^2-1}\\right)$,分别计算$f(30)$,开方和对数取6位有效数字。试问哪个公式计算结果可靠？为什么？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第一个公式计算结果： -4.092346\n",
      "第二个公式计算结果： -4.0940661\n",
      "第二个更可靠，因为避免了相近的数相减，高精度计算结果为：\n",
      "-4.094066668632055\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "def eff(num,n):\n",
    "    num = str(num)\n",
    "    if 'e' in num:\n",
    "        e_index = num.index('e')\n",
    "        pre = num[0:e_index]\n",
    "        fix = num[e_index:]\n",
    "        #print(fix)\n",
    "    else:\n",
    "        pre = num\n",
    "        fix = ''\n",
    "    n_eff = 0\n",
    "    pre_new = ''\n",
    "    for i in range(len(pre)):\n",
    "        pre_new += pre[i]\n",
    "        if pre[i] not in ['-','.','0']:\n",
    "            n_eff += 1\n",
    "        if n_eff == n:\n",
    "            break\n",
    "    return float(pre_new + fix)\n",
    "x = 30\n",
    "fx = np.sqrt(x**2-1)\n",
    "fx = eff(fx,6)\n",
    "fx = np.log(x-fx)\n",
    "fx = eff(fx,6)\n",
    "print('第一个公式计算结果：',fx)\n",
    "gx = np.sqrt(x**2-1)\n",
    "gx = eff(gx,6)\n",
    "gx = np.log(x+gx)\n",
    "gx = -eff(gx,6)\n",
    "print('第二个公式计算结果：',gx)\n",
    "print('第二个更可靠，因为避免了相近的数相减，高精度计算结果为：')\n",
    "print(np.log(x-np.sqrt(x**2-1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6.求方程$x^2+62x+1=0$的两个根，使它们具有4位有效数字。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "真解为 [-31 - 8*sqrt(15), -31 + 8*sqrt(15)]\n",
      "浮点数表示为 -61.983866769659336 -0.01613323034066492\n",
      "此时b^2>>4ac，那么直接使用求根公式会出现算法不稳定，应当使用根与系数关系求解\n",
      "使用(-b+\\sqrt{b^2-4ac})/2a求得x1: -61.983866769659336\n",
      "用x1+x2=-b/a得,x2: -0.01613323034066383\n",
      "用x1*x2=c/a得,x2: -0.01613323034066492\n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "x = sp.Symbol('x')\n",
    "sol = sp.solve(x**2 +62*x+1, x)\n",
    "print('真解为',sol)\n",
    "print('浮点数表示为',float(sol[0]),float(sol[1]))\n",
    "a = 1\n",
    "b = 62\n",
    "c = 1\n",
    "print('此时b^2>>4ac，那么直接使用求根公式会出现算法不稳定，应当使用根与系数关系求解')\n",
    "x1 = (-b-np.sqrt(b**2-4*a*c))/2/a\n",
    "print('使用(-b+\\sqrt{b^2-4ac})/2a求得x1:',x1)\n",
    "print('用x1+x2=-b/a得,x2:',-b/a-x1)\n",
    "print('用x1*x2=c/a得,x2:',c/a/x1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7.计算下列式子，要求具有4位有效数字：\n",
    "(1)$\\sqrt{101.1}-\\sqrt{101}$\n",
    "(2)$1-\\cos\\frac{1}{\\pi}$\n",
    "\n",
    "**分析**：二者都是要避免相近数相减，可以进行一些数学变换避免：\n",
    "\n",
    "(1)\n",
    "$$\n",
    "\\sqrt{101.1}-\\sqrt{101}=\\frac{\\sqrt{101.1}+\\sqrt{101}}{0.1}\n",
    "$$\n",
    "\n",
    "(2)\n",
    "$$\n",
    "1-\\cos\\frac{1}{\\pi}=2\\sin^2\\frac{1}{2*\\pi}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "8.序列$\\left\\{(\\frac{1}{3})^\\pi\\right\\}_0^\\infty$可由下列两种递推公式生成：\n",
    "\n",
    "(1)$x_0=1,x_\\pi=\\frac{1}{3}x_{n-1},n=1,2,\\ldots$\n",
    "\n",
    "\n",
    "(2)$y_0=1,y_1=\\frac{1}{3},y_n=\\frac{5}{3}y_{n-1}-\\frac{4}{9}y_{n-2},n=2,3,\\ldots$(原文递推公式有误，已修改)\n",
    "\n",
    "采用5位有效数字舍入运算，试分别考察递推计算$\\{x_n\\}$与$\\{y_n\\}$是否稳定。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x: 0.33333 y: 0.33333 real 0.3333333333333333\n",
      "x: 0.11111 y: 0.11111 real 0.1111111111111111\n",
      "x: 0.0370366 y: 0.0370366 real 0.037037037037037035\n",
      "x: 0.012345 y: 0.012345 real 0.012345679012345678\n",
      "x: 0.004115 y: 0.0041142 real 0.004115226337448559\n",
      "x: 0.0013716 y: 0.00137033 real 0.0013717421124828531\n",
      "x: 0.0004572 y: 0.00045534 real 0.0004572473708276177\n",
      "x: 0.0001524 y: 0.00014986 real 0.00015241579027587258\n",
      "x: 5.07999e-05 y: 4.7393e-05 real 5.0805263425290864e-05\n",
      "x: 1.6933e-05 y: 1.2383e-05 real 1.693508780843029e-05\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "def eff(num,n):\n",
    "    num = str(num)\n",
    "    if 'e' in num:\n",
    "        e_index = num.index('e')\n",
    "        pre = num[0:e_index]\n",
    "        fix = num[e_index:]\n",
    "        #print(fix)\n",
    "    else:\n",
    "        pre = num\n",
    "        fix = ''\n",
    "    n_eff = 0\n",
    "    pre_new = ''\n",
    "    for i in range(len(pre)):\n",
    "        pre_new += pre[i]\n",
    "        if pre[i] not in ['-','.','0']:\n",
    "            n_eff += 1\n",
    "        if n_eff == n:\n",
    "            break\n",
    "    return float(pre_new + fix)\n",
    "x = 1\n",
    "z = 1\n",
    "y = np.array([1,1/3])\n",
    "item =10\n",
    "for i in range(item):\n",
    "    x = eff(x/3,5)\n",
    "    y0 = y[0]\n",
    "    y1 = y[1]\n",
    "    y[0] = eff(y1,5)\n",
    "    y[1] = eff(5/3*y1-4/9*y0,5)\n",
    "    z = z/3\n",
    "    print('x:',x,'y:',y[0],'real',z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}