{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as pp\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd8jef/x/HXJxtBkFgxYpMGQWpUrapSo0pLUW3pUF26VXere7eq2tKhC9UqVZvaNUMosbfUSGxB9vX74z7681WE5OTcZ3yej0ceMm6539VzPuc6133dn0uMMSillPItfnYHUEop5Xpa/JVSygdp8VdKKR+kxV8ppXyQFn+llPJBWvyVUsoHafFXSikfpMVfKaV8kBZ/pZTyQQF2B7iY8PBwExUVZXcMpZTyKKtWrTpkjInI7Ti3Lf5RUVHEx8fbHUMppTyKiOy+nON02kcppXyQFn+llPJBWvyVUsoHafFXSikfpMVfKaV8kBZ/pZTyQVr8lVLKB7ntOn+llM1OH4ENv4NfAFRqBqWqgYjdqZSTaPFXSv2/nGzYPg8SfoDN0yA74/9/VjgcKjV1fDSDcvXBP9C+rCpftPgrpeDIDkj4EdaMhZP7oFBJiLsHYvuAfxDsXQZ7lsGepbBpivV3SlaF3uMgopa92VWeaPFXytet+xUm3g8mB6pfDze+DTU7QEDw/x9TujY06md9fvIA7FwEM5+DUW3h1m+g5g22RFd5pxd8lfJlq7+HCfdCxSbweCLc/gtEd/3fwn++omWhXg8YMA9KVoExPeGvYWCM63KrfNPir5SvWv4lTH4Eql0Ht/8Kxcpf2d8vXgHunmG9WMx+ESY9AJlpBZNVOZ0Wf6V80aIPYfpgqN0Zeo+FoMJ5+z1BRaDHaGjzPKwdC991tqaFlNvT4q+ULzEG/nwN/nwV6vawCvelpnguhwi0Ggw9f4CDidZ1gJMHnRJXFRwt/kr5CmOsi7SL3oeGd0K3L527VDP6Jug/DU4fht/us5aNKreV7+IvIhVFZJ6IbBSRRBF59ALHtBaR4yKyxvHxUn7Pq5S6QkuHw7IR0OQB6DIM/Pydf47yDaDje7BzASx83/m/XzmNM5Z6ZgFPGmNWi0hRYJWIzDbGbDjvuEXGmM5OOJ9S6kodTIQ/h1pz/B3eKtg7dRv0hV2LYf5b1g1hVVsV3LlUnuV75G+M2W+MWe34/CSwEYjM7+9VSjlJVjpMuA9CwqDLJwXfokEEOn0A4TWsZaQ6/++WnDrnLyJRQANg+QV+3ExE1orIdBG5ypnnVUpdwtzXITkRug6HIuGuOWdwKPT4DtJPwoR7dP7fDTmt+ItIKDABeMwYc+K8H68GKhtj6gOfApMu8jsGiEi8iMSnpKQ4K5pSvmvXYljyKTTqDzXbu/bcZaKh0/uwaxEseNe151a5ckrxF5FArML/kzHmt/N/bow5YYxJdXw+DQgUkf8MQYwxI40xccaYuIiICGdEU8p3pR2HiQOtu3Dbv2FPhgZ9oX4fWPCO1TBOuQ1nrPYR4GtgozHmw4scU9ZxHCLS2HHew/k9t1LqEqY/Ayf2QfdR1s1Ydun0PoTXtJZ/puo7enfhjJF/c+AO4LpzlnJ2FJGBIjLQccytwHoRWQsMA3oZo41AlCowiZOsO25bPgUV4uzNcvYu4DPH4M9X7M2i/pXvpZ7GmMXAJZcPGGOGA8Pzey6l1GU4sR+mPAblG0LLp+1OYykTDU0HwpLhEHc3RDayO5HP0zt8lfI2M5+1Gqx1H+lem620HAxFIqzpqJwcu9P4PC3+SnmTPcshcSI0f9RaZ+9OQopBu1chaSX8/bPdaXyeFn+lvMXZ3j2hZaH5ILvTXFi9XtaUz5yXrXsAlG20+CvlLdZPgH/ioe1L9q7uuRQ/P7jxPUg9CAvfszuNT9Pir5Q3yDwDc16BsvWgfm+701xahUYQ2xeWjoBD2+xO47O0+CvlDZZ9Dsf3Wjdz+XnA07rtSxAQYl2cVrbwgEeJUuqSUlOsnblqdYQqLe1Oc3mKloHWz8DWWbBlpt1pfJIWf6U83fw3IesMtBtqd5Ir0/h+KFUDZjxrdR5VLqXFXylPlrwRVo2GuHvcb2lnbgKCoMPbcGQ7rBhldxqfo8VfKU826wUILgqth9idJG9qXA9V28DiDyE91e40PkWLv1Keatsc66PlYChc0u40eXfdC9a+v8u/sDuJT9Hir5QnysmB2a9AiShofJ/dafKnQhzUvBGWDLOavymX0OKvlCfa9AccXAetn4OAYLvT5F+b56z9B5aNsDuJz/C+4p+TA+t/g9Rku5MoVTBycmD+29ZKmbq32p3GOcrVg+iu1o1fp3SrD1fwvuJ/bJe1Z+iSYXYnUapgbJgEyRusi7x+/nancZ7Wz0FGKiz5xO4kPsH7in/JqlC3J6z4SncNUt4nJ9sa9UfUhqu62Z3GuUrXhro9YPlIOHnQ7jRez/uKP1gbWGSn6whCeZ/EiXBoM7R6xrtG/We1HgLZGbD4I7uTeD3vLP7h1XX0r7zP2VF/6WiIvtnuNAWjVDWI7QPxX8PxJLvTeDVnbOBeUUTmichGEUkUkUcvcIyIyDAR2SYif4tIw/yeN1f/jv517l95iXW/wuGtjrl+7xy3AdBqsLU3wcL37U7i1ZzxCMoCnjTG1AGaAg+JSPR5x9wI1HB8DAA+d8J5Ly28ujV/uFJH/8oLZGfBgrehTAzU7mJ3moIVVgka3QUJP8DRXXan8Vr5Lv7GmP3GmNWOz08CG4HI8w7rCnxvLMuAMBEpl99z56rl05CVpqN/5fnWjYcjO6D1s9496j+rxVPgFwALdMOXguLUR5GIRAENgOXn/SgS2HvO10n89wUCERkgIvEiEp+S4oTRengNHf0rz5edCQvesTZqqd3J7jSuUawcxN0Na8fq6L+AOK34i0goMAF4zBhz4vwfX+CvmP98w5iRxpg4Y0xcRESEc4Lp6F95urXjrALY5jmQCz2VvNQ1g6zR/6IP7E7ilZxS/EUkEKvw/2SM+e0ChyQBFc/5ugKwzxnnzlV4DYi51Rr9nzrkklMq5TTZmbDwXSjfAGp2sDuNaxUrZ839rxkDx/bYncbrOGO1jwBfAxuNMR9e5LDJwJ2OVT9NgePGmP35Pfdl09G/8lR//2wVvlZDfGvUf1bzx0D8dN1/AXDGyL85cAdwnYiscXx0FJGBIjLQccw0YAewDRgFPOiE816+iJrW6H/FKB39K8+RnWUtdyxXH2q2tzuNPYpHQoO+sPoHXffvZM5Y7bPYGCPGmHrGmFjHxzRjzBfGmC8cxxhjzEPGmGrGmLrGmPj8R79COvpXnmb9r3B0p9Wv3xdH/Wdd+zhg4C+9Y9+ZfGDNmENETYi5xbrrV7sGKneXkw0L37PW9dfqaHcae4VVsu76XfUdnHDdbLG3853iD9ba4czTsOwzu5ModWmJE+HwNusdqy+s68/NtU9ATpa+c3ci33pUla4NV91sdQ08fcTuNEpdWE4OLHjX6txZ5ya707iHklWgfi+I/0Y7fjqJbxV/sOZPM07CsoLvMKFUnmz83ercqaP+/9XiSavj59JP7U7iFXzvkVUm2hpNLf9C9wtV7icnx2ppUKqG9/Xrz69S1Rx37H+tq/acwPeKP1gjqvQTsPxLu5Mo9b82T4XkRMeo3wv79edXi6cg8wwsHW53Eo/nm8W/XD2o1cm68Jt23O40SlmMseb6S1a1Vqap/4qoCTHdret2umovX3yz+IPVMzztOKwYaXcSpSxbZsKBv63RrX+A3WncV8vB1qo9nfvPF98t/uVjrV4pSz+D9JN2p1G+zhirc2dYZajX0+407q10bR39O4HvFn+wRhBnjlptH5Sy09bZsG81tHgC/APtTuP+dPSfb75d/Cs0gurXWxeP0lPtTqN8lTEw7w3rTtb6fexO4xlK17aui+joP898u/gDtHoGTh+2bh5Ryg6bp8P+NdZoNiDI7jSeo5Vj9K93/eaJFv+KjaFqG6tplI7+lasZA/PfhBKOO1jV5Yuo5ejXpd1680KLP1g7JJ0+pCt/lOttmgIH1lnvQHWu/8q1esYx+te5/yulxR+s0X+N9tboX9f9K1fJyYF5b0Gp6tadq+rKRdSEurpXR15o8T+rzXOQdgyWjrA7ifIVG3+37uZtNUTX9edHy8GQdUbn/q+QFv+zysdaPX+WfqYdP1XBy8mG+W9DeC1rzbrKO92pL0+ctYH7NyKSLCLrL/Lz1iJy/JxtHl9yxnmdrs1zkJGqOwapgpc4EVI2Qesh2sPHGVoN1p36rpCzRv6jgQ65HLPonG0ehzrpvM5Vuo41f7j8S+0ZrgpOdhbMfwtKR0P0zXan8Q7hNazrJitGQWqy3Wk8glOKvzFmIeAdcyWthlg9wxd/ZHcS5a3W/2rt0tX6We3X70ytnrGeuwvetTuJR3DlI6+ZiKwVkekicpULz3tlwqtDbG+I/xqOJ9mdRnmb7Cyrh0/ZulC7s91pvEupatDwTlj1LRzZYXcat+eq4r8aqGyMqQ98Cky60EEiMkBE4kUkPiUlxUXRLqDlYOvmm4Xv25dBeae1Y6zCpKP+gtHqGfAPgrmv253E7bnk0WeMOWGMSXV8Pg0IFJHwCxw30hgTZ4yJi4iIcEW0CytRGRrdBQk/wJGd9uVQ3iXjFMx7EypcDbU62p3GOxUtC00fhPUTYN8au9O4NZcUfxEpKyLi+Lyx47zu3Y2pxVPgFwAL37M7ifIWy0bAyf3Q7jWwng6qIDQfBIVKwp+v2p3ErTlrqedYYClQS0SSROQeERkoIgMdh9wKrBeRtcAwoJcxxjjj3AWmWDm4+l5YOxYObrA7jfJ0qSmw+BNrnr9yM7vTeLeQ4tZm79vnwo75dqdxW+KuNTguLs7Ex8fbG+L0ERgWa71N7zvB3izKs019yuoc+9Bya1miKliZaTA8DoqEw33zfOqdloisMsbE5XacXnG6lMIlrYu/2+bA1jl2p1Ge6tA2awVKo35a+F0lMMS6aXNfAmy44PoSn6fFPzeN77Pa7c56wVqmp9SV+vMVCAix7uZVrlPvNutGuj+HQnam3Wncjhb/3AQEQ7uhkLIREr63O43yNHuWwcY/oPmjEFra7jS+xc8f2r5sLa1drc/d82nxvxx1ukCla2DuG5B2wu40ylMYA7NehNCy0Owhu9P4pprtoVIz68a6jFN2p3ErWvwvhwi0f8Pa8GXxh3anUZ5i42RIWmHNPQcVsTuNbxKB61+F1IPasuU8WvwvV2RDqN/b6vd/dLfdaZS7y86EOa9ARG2Ivd3uNL6tUhOr6dtfn8Dh7XancRta/K/EdS+C+OnNIyp38d9Yc83thupGLe6g3WvgHwwzhljTcUqL/xUpHmndPbh+AuxdYXca5a5OHrR6y1RpBTVusDuNAuumzdZDYOss2Dzd7jRuQYv/lbpmkHUBb+ZzOoJQFzbreWtjkU4f+tTNRW6vyf0QUQdmPAOZZ+xOYzst/lcqOBTavghJK63WD0qda/s8WPcLXPuE1R5cuQ//QOj4HhzbA4s/tjuN7bT450X9PlCxiTX6T7Wx9bRyL5lpMPVJKFkVrn3c7jTqQqq0gJhbrJU/Pt6xV4t/Xvj5QZdhkJ4KM5+1O41yF399DEe2Q8f3rfYCyj3d8Lr1LmCGbz93tfjnVenaVufAdb/A1tl2p1F2O7wdFn1gjSqrt7U7jbqUYuWtDd+3TIctM+1OYxst/vnR4gkIrwVTnrDeBSjfZAxMfcLq39P+TbvTqMvR5AEIrwnTB1vTdT5Ii39+BARDl0/g+B5rhyblm9ZPsPrGt33J2klKub+AIOvi79Fd1js2H6TFP78qN4O4e2D55/DPKrvTKFc7c8y68F++AcTdbXcadSWqtrbu2l/0gU8+d7X4O8P1L0NoGZg8SFvH+po/h8KpFOj8kdVFUnmWDm9D0XLw2/2QcdruNC6lxd8ZQopbKzwOrocln9qdRrnK1tkQ/zU0GWiN/JXnKRQGN38Gh7f6XNsWZ+3h+42IJIvI+ov8XERkmIhsE5G/RaShM87rVup0hjo3Wa1jtXmU90tNhkkPQOmrrJ7xynNVbW29gC//wqf2/HXWyH800OESP78RqOH4GAB87qTzupeO71kXgSfcC1kZdqdRBcUYmPSgtbfDLV/pmn5v0PZlKFXD+v965pjdaVzCKcXfGLMQOHKJQ7oC3xvLMiBMRMo549xupWhZ6PoZ7FsNc3Q06LWWfwnbZls3C5WJtjtNgcnJMew7dob9x8+QnpVtd5yCFVQYun8JJw/A9GfsTuMSruo1GwnsPefrJMf39rvo/K5Tpws0vh+WjYCoFlC7o92JlDMdTITZL0GN9tb+zl4gKzuH+N1H2Zqcyq5Dp9h9+BS7Dp9mz5HTZGTl/Htc0ZAAwkODKVUkiFKhQVSLCKVddBnqVwjDz88LGthFNoKWT8OCt63nbXRXuxMVKFcV/ws9Mv7TElNEBmBNC1GpUqWCzlRwbngN9i6z5oQHLoIwD/5vUf8v8wz8eo91gb/rZx7fsXPj/hP8tjqJSWv2kXIyHYCQQD+iShWhWkQR2tYuTaVShRGEw6npHD6VwaHUdA6nZrAj5RR/bkxmxPztlC4aTLvoMtxwVVmaVS1FUIAHryNp+RRsmQF/PAYVm0LRMnYnKjBinNSWWESigCnGmJgL/OxLYL4xZqzj681Aa2PMRUf+cXFxJj4+3inZbHF4O3zZymoD0X+61UtEebZpT8OKkdB3AlS/3u40eZJ8Io3f1+xjwuokNh04SaC/0KZWabo1iKRh5RKULhqMXOaL2vHTmczbnMysDQeYvzmF0xnZhAYH0C66DA+1qUb10kUL+L+mgKRshi9aQNS1cPsvHreEV0RWGWPicj3ORcW/E/Aw0BFoAgwzxjS+1O/z+OIPsP43+LU/NH/U2tFJea7NM2DsbdD0Qejwlt1prtiRUxl8MGszY1fsIcdAbMUwbmkYSed65SlRJCjfvz8tM5sl2w8xK/Egf6zdx5nMbG5pWIHH2tUkMqyQE/4LXCz+W5jymEc+d11a/EVkLNAaCAcOAi8DgQDGmC/EGkoMx1oRdBrob4y5ZGX3iuIPMOVxa0u/Pr9ATd3VySMd3QWjroOi5eG+P60VXR4iMzuHH5bu5uM5WziVkc0dTStzZ7PKVI0ILbBzHk5NZ8T87fyw1Nrrum/TyjzUphqlQj3n3w2wenbFfw3dR0G9nnanuWwuH/k7m9cU/8wz8NX1cGIfDFxsbQWpPMeZY/BNezi5H+79E8Jr2J3osi3cksLQKRvYlpxKixrhvNQ5mhplXDcV88+xM3wyZwu/rkqiUKA/97aoygOtqxES6CHTKNmZ8H1Xq/VD/+kQ6Rm3J2nxdyeHtv7//P9df0BQEbsTqcuRnQk/9YBdi+COiVClpd2JLkvS0dO8MjmRORuTqVyqMC90iub6OqUvey7f2bYln+SDWVuYvv4AtcsWZXifBp5zPeDUIRjZBnKyYMA8j2jcd7nF34Mvy3uQ8BpwyyjYlwC/9IfsLLsTqdwYY13g3THP6tzqIYV//uZkOn+6mKXbD/NMh9rMerwl7aLL2Fb4AaqXLsrnfRvxbb+rST6ZTpdP/+LnlXtw14Hn/ygSDr3HQNox+PkOyEq3O5HTaPF3ldqdrP4/W2daF5I84YHvy5Z+Bqu+tbZjbNDX7jS5yskxfDxnC/1Hr6RssRCmDmrBA62rERzgPlMsbWqXZvqjLWhQKYxnJqzjkbEJnEjzgEaIZevCzZ9D0grrOoCXPHe1+LvS1fdYN5Ek/ADz37Y7jbqYjVNg1gvWTT7XvWR3mlwdO51B/9Er+XjOVro1iGTig82JCnfPqcUyxUL44Z4mPN2+FtPXH6DTsEWs2esB7RSuuhlaDoY1P1o9gLyAFn9Xa/M8xPa17iKM/9buNOp8+xLgt/usi3vdvrT2a3Zj6/85/u80z+s3x/BBj/oUCnKf0f6F+PsJD7Wpzvj7m5KTA7d+voSxK/bYHSt3rZ+FWp2svX/XjrM7Tb659yPbG4lAl4+hejtr679N0+xOpM46shPG9ILC4dB7HAS69/r031Yn0f3zJeTkGMYPbEbfppVtndu/Uo0ql2TaoBZcWyOcZ39bx/C5W937OoCfn9XIr0oL6+79v3+xO1G+eGXxX7Al5X96krgd/0Do+R2Ui4Vf74a9K+xOpFK2wLcdITsd+vwMoaXtTnRJo//ayRPj1xJXuQRTBrUgtmKY3ZHypHjhQEbdGUf3BpG8P2sLr/6xgZwcN34BCCoMvX+Gys1h4gBY96vdifLM64r/tuRU+n+7ggd/Wu3eLwBBRaDPeChWDn66FfYsszuR7zqYCKM7Qk4m9Jvq9p06R8zfxit/bKD9VWX4tv/VlHTCHbp2CvT34/0e9bn32iqMXrKLx35e4+bP3cLWAKFSM2uKcP1vdifKE68r/tVLh/Jq1xjmbDzIgz+tcu9WtKERcOfvUCTCuplk8wy7E/mefQkwuhP4BVo38pS5yu5EF2WM4f2Zm3l3xma6xpZneJ+GbrWaJz/8/ITnO9XhmQ61mbx2H/d+H8/pDDdeEn128FaxibV/R+IkuxNdMa8r/gB3NK3M6zfHMGdjMg/8uNq9XwDCKsHdM6F0HRjXB9aMsTuR79izHL67CYKLQv9pbn33rjGG16ZsZPi8bfS6uiIf9owl0N+7nr4iwgOtq/HOLXVZvDWFPqOWc/SUG2+KFBxqNX6rEAcT7oGNf9id6Ip416PnHH2bVuaNbjHM3eQBLwBFwq07f6OutS4k/TXM7kTeb+ci+KGb9a6r/3QoWcXuRBeVnWN4buI6vvlrJ/2uieKt7nXx94b++Rdx29WV+LxvIzbsP0Gfr5Zz/Iwb3wsQXBRu/xXKN4Rf+ll9vDyE1xZ/gNubVObNbnWZuymZgT+sIi3TjV8Agotao4irusHsF2HWi15zM4nb2TDZus4SVskq/MUr2J3oonJyDE//spaxK/byYOtqvNwl2qNW9ORV+6vKMurOOLYln+Se0Ss5k+HGz92QYlab76qtrUaOU56wWoO4Oa8u/gB9mlTizW51mbc5hYE/uvkLQEAw3PI1XH0vLBlm7SfqRbeT2y47E2Y+D+PvsOb2+0116806jDEMnbKB3xL+4Yl2NRncobZPFP6zWtWM4OPbGrBqz1Ee+GmVe18EDilmXQO4ZpDVCfT7rlZfIDfm9cUfrBeAt7rXZf7mFB4ek0C2Oy8l8/O32kC0fhbWjoGv2lqN4VT+nDxgze8vHW69uPafDkVK2Z3qkj5fsJ3RS3Zxd/MqPHJddbvj2KJTvXK82c167j4xfo37P3dveM1qAf3PKqsh3IF1dqe6KJ8o/gC9G1fi1ZuuYs7Gg7wyOdG9byYRgdZDrBuNjv9jdQRN+EmngfJq12JrZ6b9a6D7V9DpA7fvyf9L/F7enbGZm+qX54VOdXxqxH++3o0r8eyNtZny935e/H29ez93wer933+61Qn06xvcdiWQzxR/gLuuiWJAy6r8sGw3IxfusDtO7mrdCA/8ZbUa+P1Ba0lZ2gm7U3kOY2Dxx9aIP6Q43DcX6vWwO1Wu5m46yJDf1nFt9XDe71HfOzZHz6f7W1XjgdbVGLN8D+/N3Gx3nNxFNoQB86FMDPxyl3UtwM2euz5V/AGGdKhNp3rleGv6Jv5Yu8/uOLkrVt66F6DNC5A4Eb5sAUmr7E7l/g5vhx+7w5yXoU4Xqxd76Tp2p8pVwp6jPPjTauqUK8oXdzTy7M3QnWxw+1r0aVKJEfO3M3Lhdrvj5K5oGeg3BZo+BKtGw2dNYPN0u1P9y+ceWX5+wgc96nN1VAmeHL+WFTuP2B0pd37+0Oppay16TjZ8c4PVdTLtuN3J3E/mGZj7BoxoCntXWtdPeoy2VlO5ue0pqdw9eiVlioXwbb/GhAYH2B3JrYgIr3WNsa4DTNvEjPUH7I6Uu4Bg6PAm3DMHCoXB2F7Wnh6pKXYnc07xF5EOIrJZRLaJyJAL/LyfiKSIyBrHx73OOG9ehQT6M+rOOCqULMR938ezLfmknXEuX6WmMHAR1O8FS4bDsIbWumLdHMayeYY1ulr4rtWO+ZF4aHyfdQ3FzaWcTOfOr1fg7yd8f3djIoq69zUJu/g7Bm+xFcN4/Oc1JO7zkAFQhUYwYIH1Dn7TFPjsalgz1tbrePku/iLiD3wG3AhEA71F5ELNUX42xsQ6Pr7K73nzK6xwEN/1b0ygv9Dv25Ukn0yzO9LlKVQCun5mzSdG1LLmEr9sCdvn2Z3MPkd3wdjeMPY2CAixbpi75SuP2HIPICMrhwd+XMXhU+l8268xlUu5Zy9+dxES6M/IOxsRVjiQ+76LJ+WkhyyHDgiy3sEPXAzhtWDSQOuC8LY5trwIOGPk3xjYZozZYYzJAMYBXZ3wewtcxZKF+abf1RxOzeC+7+Ld+x6A85WPtdap9/weMlLhh5thzG1uvbTM6Q5ugIkD4dNGsGMBtBtqPbE8ZMtFsNbyvzx5PfG7j/J+j/rUrVDc7kgeoXTREEbdGceR0xnc/4OHPXcjalmrgbp8Aif2wY+3wFfXw9bZLn0RcEbxjwT2nvN1kuN757tFRP4WkV9FpOKFfpGIDBCReBGJT0lxzZxYvQphfNwrlrVJx3lhkgcsIzuXiDW98dAKuP5V2PUXfHEtfNvJ6jOS40FPiMtljPXf+VNP+LwZbPgdrr4PHl4JzR+1Rlce5Mfle/69e7dzvfJ2x/EoMZHF+bBnLKv3HOO5ies867nr5weN+sGg1dD5Y0hNtu46H3UdbJnpkhcBye8/mIj0ANobY+51fH0H0NgY88g5x5QCUo0x6SIyEOhpjLnuUr83Li7OxMfH5yvblfhw1maGzd3Ga12v4o5mUS47r1OdOQqrf4AVo+D4HiheCRrfCw3vtKaLPFlWOmyZAUs+haSVULgUNBlo3bBVuKTd6fJk+Y7D3P7VclrWjGDUnXFe3a+nIH0yZysfzdnCkBtrM7BVNbvj5E1WBqwdC4veh2N7IKqFNX2Zh+tVIrLKGBOX63FOKP7NgFeMMe0dXz8LYIx56yLH+wNHjDGXfH/r6uKfk2O457uVLNp6iLEDmnJ1lGcWFMAa8W+eBsu/hF2LILAy7x+NAAAbOUlEQVQwRN8MdTpD1TZWP3JPkJ0JO+Zb/dI3TYX04xBWGa55BGJv95z/jgtIOnqarsP/onjhQCY91JxiIYF2R/JYxhgeGZvA1HX7GXVHHNdHu2/LjlxlZ8LfP0N6KjQdmKdf4criHwBsAdoC/wArgT7GmMRzjilnjNnv+Lwb8Iwxpumlfq+riz/A8TOZdB2+mNT0bKY8ci1li4e49PwF4sB6WPGlNT2SdhwCCkH1tlC7M9Rs736j5sw02LPUuqdh42Tr3UxwMStvTHfrxcvfs5dAnsnI5pbPl7D36GkmPdScahGhdkfyeGcysun55VJ2pKTy+8PNqV7a/Zf2FhSXFX/HyToCHwP+wDfGmDdEZCgQb4yZLCJvATcBWcAR4AFjzKZL/U47ij/A5gMn6TbiL2qVLcq4AU29ZrMMsjNh91+wcYo1ij65D8QfKja2+pFHNrLa0oZVcu3SyNNHYO9yq+DvWWZtrpKdAYFFrDucY7pDtbYQ6AUvxPzvKPWbflfTppZ7bxfpSfYfP0PnYYspWSSI3x9uTuEgzx4k5JVLi39BsKv4A0xbt58Hf1r9b0dQr2OMVWQ3TbGmVQ6sswouWJuXRzayVhOFVYawilbL42KRee+HYwykHoQjOxwfO60/kzdAimMM4BcI5RtY9zJUama1x/XgaZ2LGbVwB29M2+jZ89NubPHWQ9zxzXK61i/PR7fF+mRPpMst/r750piLjnXLMbBVNb5YsJ26kcXp3biS3ZGcS8TqPRLZENq+ZF1sSk60OhH+s9r6c+ss4NyBgUBoGSgeCUGh1nr6gGDrz8AQ8A+2XkAyUq35yvSTkHHS+jz1IGSePudX+VvvMMJrQt0eVrGPbAiBhVz9L+FSq3Yf4e0Zm7gxpiz3t6xqdxyvdG2NcB6/viYfzt5CXFRJ+jatbHckt6XF/yKebl+LxH3Hefn3RGLKF/fu9dcBQdaou3wDuNrxvcw0OPEPHE8652Ov9b2M05B2zFqBk5Vm/Zl5BvyDrDYKwaHWPH2xCtbXRcKhZFVrt6ySVaF4RfD3rQucR05l8PCYBCqUKMQ7t9bzyRGpqzzcpjqr9xxl6B8bqFehOPUqhNkdyS3ptM8lHD2VQcdhiwgK8OOPR67VFRkqT3JyDP1Hr2TpjsP89sA1xER68UDCTRw9lUGnYYsQEaYOupawwp51/0d+XO60j881drsSJYoE8WnvBiQdPcOzEzzsJhLlNj5fsJ0FW1J4uUu0Fn4XKVEkiBF9G5F8Mo0nxq8lx503gbGJFv9cxEWV5MkbajJ13X5+XL7H7jjKwyzbcZgPZlmbsvTxtmtHbi62Yhgvdo5m7qZkPl/gAS2gXUyL/2UY2LIarWpG8NqUDZ7TRVDZLuVkOo+MTSAqvAhvdq+r8/w2uKNpZbrUL88HszazdPthu+O4FS3+l8HPT/iwZ31KFA7k4TEJpKZrC2V1adk5hkfHJXAyLZMRtzfU3vw2ERHe7l6XqFJFeOznBI6cyrA7ktvQ4n+ZSoUGM6xXA3YfPsVzv+n8v7q0T+duZcn2wwztGkPtssXsjuPTigQHMKx3A46eymTwr2v1ueugxf8KNKlaiifa1WTy2n2MW7k397+gfNKKnUcY9udWujeIpGfcBRvYKheLiSzOsx1rM2djMqOX7LI7jlvQ4n+FHmxdnRY1wnllciKbDrjXhszKfsdOZ/DYuAQqlSzM0Jtj7I6jztHvmiiur1Oat6ZtYv0/eu1Oi/8V8vMTProtlmKFAhk0NsGzNpFQBcoYw5AJ60hJTefT3jrP725EhHdvrU+JIoE8MjaBUz5+7U6Lfx6EhwbzQY/6bDmYyutTN9gdR7mJMSv2MCPxAIPb1/buO8I9WMkiQXziuHb30u+Juf8FL6bFP49a1ozgvhZV+HHZHmYlHrA7jrLZloMnGfrHBlrWjOCea6vYHUddQtOqpXj4uhpMWJ3EpIR/7I5jGy3++fB0+9rERBZj8IS/OXDcQzaAV06XlpnNI2MSKBoSwAc96uOnO3K5vUHXVadxVEmen7iOnYdO2R3HFlr88yEowI9hvRqQkZXD4z+vIVtvIfdJb07byOaDJ/mgZywRRfPY9lq5VIC/Hx/3iiXA349HxyWQkZVjdySX0+KfT1UjQnnlpqtYuuMwX+gt5D5nVuIBvl+6m/taVKFVzQi746grUD6sEO/cUpe/k47z0ZwtdsdxOS3+TtCjUQU61yvHh7O3kLDnqN1xlIscOJ7G4Al/UzeyOE+3r213HJUHHWLK0btxRb5YsJ0l2w/ZHcelnFL8RaSDiGwWkW0iMuQCPw8WkZ8dP18uIlHOOK+7EBHe6FaXssVCGOS4pV95t5wcwxPj15CemcMnvWIJCtBxlKd6sXM0VcKL8MTPaznqQ+0f8v2IFRF/4DPgRiAa6C0i0ecddg9w1BhTHfgIeCe/53U3xQsFMqx3LPuOpfHipPV2x1EFbOSiHSzZfphXb7qKqroBu0crHBTAsF4NOHwqnWd9qHWLM4YrjYFtxpgdxpgMYBzQ9bxjugLfOT7/FWgrXtjisFHlkgy6rgaT1uzz6SVk3m5d0nHen7mZjnXL0iOugt1xlBPERBbn6fa1mJF4gJ99pHWLM4p/JHDuv1aS43sXPMYYkwUcB0o54dxu56E21YirXIIXJq1n75HTuf8F5VFOpWcxaFwCEUWDeaubbsfoTe69tirXVg/n1T82sD0l1e44Bc4Zxf9Cj/7z3zddzjGIyAARiReR+JSUFCdEc70Afz8+ui0WAR4dl0BWtu8tIfNmQ//YwK7Dp/jotliKF9ZtPb2Jn5/wQc/6hAT6MWhsAulZ3t26xRnFPwk4t3VhBWDfxY4RkQCgOHDk/F9kjBlpjIkzxsRFRHjusrmKJQvzRve6rN5zjGFzt9kdRznJtHX7+Tl+Lw+2rkbTql75xtXnlSkWwru31idx3wk+mOXdyz+dUfxXAjVEpIqIBAG9gMnnHTMZuMvx+a3AXOPlV1Vuql+e7g0jGT53Kyt3/ed1TnmYfcfOMGTC39SvGMZj19e0O44qQO2iy9C3aSVGLtzBwi2eOQNxOfJd/B1z+A8DM4GNwHhjTKKIDBWRmxyHfQ2UEpFtwBPAf5aDeqOhXWOoUKIwj41bw/EzuvzTU2XnmH/v4P7ktlgC/XVZp7d7oVM0NcuE8sT4tRxKTbc7ToFwyqPYGDPNGFPTGFPNGPOG43svGWMmOz5PM8b0MMZUN8Y0NsbscMZ53V1ocACf9IrlwIk0Xpi03meWkHmbz+dvY/nOI7zaNYao8CJ2x1EuEBLoz7DeDTiRlsngX//2yueuDmEKWINKJXiiXU3+WLuPCat1+aenWbX7CB/N2UrX2PLc0vD8RWzKm9UuW4znO9Zh7ibv3P1Li78LDGxVjSZVSvLS7+vZ4QNLyLzF8TOZDBq7hvJhIbx+c4wu6/RBdzar/O/uXxv2edfOfVr8XcDfT/jY0QLgER9YQuYNjDE899s6Dp5IY1ivBhQN0WWdvujs7l9hhQN5ZOxqzmR4z3NXi7+LlCteiPccS8jemb7Z7jgqF+Pj9zJ13X6euKEmDSqVsDuOslHJIkF8dFssOw6dYugU79m5T4u/C7WLLkO/a6L45q+d/LnxoN1x1EVsSz7JK5M30Lx6KQa2rGZ3HOUGmlcP5/6W1Ri7Yg8z1u+3O45TaPF3sSE31qZOuWI89cta3f3LDaVlZvPI2DUUCvLnw56xuiuX+tcT7WpSr0JxnpmwjqSjnt+6RYu/i4UE+jO8TwPSMnN47OcE3f3Lzbw9fRMb95/gvVvrUaZYiN1xlBsJCvDj094NyMkxPDI2gUwPb92ixd8G1SJCGdr1KpbtOMKIedr+wV3M3nCQ0Ut20e+aKNrWKWN3HOWGKpcqwtu31CNhzzHem+nZ1+60+Nvk1kYV6Bpbno/mbNH2D25g75HTPDl+DTGRxRhyo+7KpS6uU71y3NG0MiMX7mDuJs+9dqfF3yYiwus3x1CxZGEGjU3giA/tIORu0rOyeXjMagwwok8jQgL97Y6k3NzzneoQXa4YT4xfy75jZ+yOkyda/G1UNCSQ4b0bcjg1g8ccvWOU6701bRNrk47z3q31qVSqsN1xlAcICfTns9sbkpmV47Hz/1r8bVa3QnFevimahVtSGK7tn11u2rr9jF6yi7ubV6FDTFm74ygPUiW8CG/dUo9Vu496ZPtnLf5uoE/jSnRvEMnHf25h0VbvbSHrbnYdOsXgX/8mtmKYzvOrPLmpfnl6N67EFwu2M29zst1xrogWfzcgIrzeLYYapUN5dNwaj51D9CRpmdk8+NNqAvyFz25vSFCAPhVU3rzcJZraZYvypIfN/+sj3k0UDgrg876NSM+0Lj5mZHneHKInGTplAxv2n+DDnvWJDCtkdxzlwUIC/Rlxe0MysnIY+OMq0jI9o/+PFn83Ui0ilHdurcfqPcd4e/omu+N4rYkJSYxZvoeBrapxXW1dz6/yr2pEKB/2rM/fScd50UP27tDi72Y61yv/b/+fqX97Rw8Rd/J30jGembCOplVL8tQNuh2jcp4brirLoLY1+GVVEj8u32N3nFxp8XdDz3WsQ4NKYQz+dS1bDp60O47XSD6Zxv0/rCIiNJgRtzciQLdjVE72WNsaXFe7NK9OTnT7mzfz9egXkZIiMltEtjr+vGDvWxHJFpE1jo/zN3dX5wkK8GPE7Q0pHBzAvd/F6w1gTpCelc0DP67m2OlMRt7ZiJJFguyOpLyQn5/w0W2xVChRiAd+XO3WzRvzO/QZAvxpjKkB/MnFN2Y/Y4yJdXzcdJFj1DnKFS/EyDsaceBEGg/8uEovAOeDMYaXf09k1e6jvN+jPleVL253JOXFihcKZOSdcZzOyOKBn1a57eZN+S3+XYHvHJ9/B9ycz9+nztGgUgneuaUuy3ce4eXJiR5xEckd/bhsN+NW7uWhNtXoVK+c3XGUD6hZpijv96hPwp5jvPqHe24Ak9/iX8YYsx/A8WfpixwXIiLxIrJMRC76AiEiAxzHxaek6M1OAN0aVOCB1tYmEt954SbSBW3p9sO8+scG2tYuzZPtatkdR/mQjnXLMbBVNcYs38MPS3fZHec/AnI7QETmABe67/35KzhPJWPMPhGpCswVkXXGmO3nH2SMGQmMBIiLi9NhrsPTN9RiW3IqQ6dsoGpEKC1rRtgdySPsPXKah8aspnKpwnzUSzdmUa73dPtabD14kpcnJ1K2eCHaRbvP0uJcR/7GmOuNMTEX+PgdOCgi5QAcf17w/mZjzD7HnzuA+UADp/0X+ICzF5FqlinKQ2NWsy051e5Ibu/46Uzu+W4lmdk5jLozjmK6Abuygb+f8GmfBsREFueRsatZu/eY3ZH+ld9pn8nAXY7P7wJ+P/8AESkhIsGOz8OB5oB7ToK5sdDgAL66K44gfz/u/W4lx07rCqCLScvM5t7vV7Lr0Gm+7NuIqhGhdkdSPqxwUABf33U14aHB3PPdSvYcdo8tIPNb/N8G2onIVqCd42tEJE5EvnIcUweIF5G1wDzgbWOMFv88qFCiMF/c0Yh9x9LoP3olp9Kz7I7kdrKyrRa78buP8uFt9bmmerjdkZQiomgwo/s3JjPb0G/0Co66wfJtcdcVJHFxcSY+Pt7uGG5pxvoDPPjTKppXD+eru+IIDtDNR8Ba0vncxHWMXbGXV7pE0695FbsjKfU/Vu46wu1fLadeZHF+vLdJgWwcJCKrjDFxuR2ntzh6oA4xZXnnlnos2nqIR8euIcsDN5IoCB/N2crYFXt5sHU1LfzKLV0dVZIPe9YnfvdRnhy/lhwbN3DS4u+hesRV5MXO0cxIPMBzE9f5/D0APy7bzbA/t9KjUQWebq9LOpX76lyvPM91rM3Udft59Q/77t/Jdamncl/3XFuF46czGDZ3G8ULBfJcxzqI+N5yxhnrD/DS7+u5rnZp3upe1yf/DZRnua9FVVJOpjNq0U4AXrnpKpc/brX4e7jH29Xk2JlMRi3aSVjhIB5qU93uSC41fd1+Bo1LoH7FMD7r01CbtSmPICI817EOAKMW7cQAr7r4BUCLv4cTEV7pchUnzmTy3szNBAf4cW+LqnbHcomJCUk8OX4tsRXD+LZ/YwoF6YVv5TnOvgCICCMX7sAYGNrVdS8AWvy9gJ+f8F6P+qRl5vD61I0cO53JkzfU9OrpjzHL9/D8pHU0rVKKr+6Ko0iwPpSV5xERnr2xNiLw5YIdGAxDb4pxyd3o+ozxEoH+fgzv04AXJq1n+LxtHD6Vwes3x+DvhS0Nvlq0g9enbqRNrQg+79uoQJbLKeUqIsKQDrURhC8WbMcYeK1rwb8AaPH3IgH+frzVvS6lQoP4bN52jp7K4ONesV5VHIfP3cr7s7ZwY0xZPunVQDdeV15BRHimQy1E4PP52zHA6wX8AqDF38uICE+3r03JIsG8NmUD/b9dycg7G1HUw3vb5OQY3pu1mc/nb6d7g0jevbWeXtxVXkVEGNy+FgKcSMukoGdt9Q5fLzYp4R+e+mUttcoWZXT/xkQUDbY7Up4cP5PJk+PXMmfjQfo0qVTgIyKl7HS2Juf1mp3e4au4uUEko+6KY3tKKl2HLybezfcUvZCN+0/Qdfhi5m9O5uUu0bxxsxZ+5d1ExCWLNbT4e7k2tUoz/v5mBPj7cdvIZXz651aybbyl/EpMTEii24i/OJ2RzbgBTenfvIpXr2BSypW0+PuAehXCmDroWjrVLccHs7fQ96vlHDzhvhtLp2dl8+Kk9Tz+81rqVwhjyqBriYsqaXcspbyKFn8fUTQkkE96xfLurfVYs/cYN36yiLmbDtod6z+2Jady25fL+GHZbga0rMpP9zahdNEQu2Mp5XV0tY8PERF6xlWkYaUSPDxmNXePjqffNVE8dn0NwgoH2Zrt+JlMhv25le+W7KJQkD8jbm9Ix7q62bpSBUVX+/iotMxs3p6+ie+W7iI0OICBrarRv3kUhYNcOx7IzjGMj9/L+zM3c+R0Br2ursiTN9QiPNQzVyYpZbfLXe2jxd/HbTpwgvdnbmbOxmTCQ4MZ1LY6va6u5JKbp1bsPMKrfySSuO8EjaNK8lKXaGIiixf4eZXyZi4p/iLSA3gFa6vGxsaYC1ZrEekAfAL4A18ZY97O7Xdr8XetVbuP8M6MzazYeYSKJQvxaNuadKxb1unvBE6kZTJ93X5+W/0Py3ceoXzxEJ7tWIfO9crpSh6lnMBVxb8OkAN8CTx1oeIvIv7AFqw9fpOAlUDv3Pbx1eLvesYYFmxJ4b2Zm0ncd4LgAD9a1IjghugytK1TmlJ5nIrJzM5h0dYUflv9D7M3HCQ9K4eqEUXoGVeRu5pFaTdOpZzocot/voZ1xpiNjpNd6rDGwDZjzA7HseOAroBu4u5mRITWtUrTskYEy3YeZlbiQWZvOMicjQfxE4irXJLro0tTNTyUUqFBhIcGUyo06N93B8YYjp7OZNfhU+w+fIpdh06z6/Ap/tp2iEOpGZQoHEivqyvSvWEF6lUoriN9pWzkiqt7kcDec75OApq44Lwqj/z8hGuqhXNNtXBe7hJN4r4TzN5wkFkbDvLmtE3/Ob5QoD8liwRxIi2Tk2lZ/35fBCLDCtG4Skm6NahAq5oR2ohNKTeRa/EXkTlA2Qv86HljzO+XcY4LDe8uONckIgOAAQCVKlW6jF+tCpqIEBNZnJjI4jzeribJJ9LYfzyNw6fSOZSaweHUDA6npnP4VAahwQFEhRchqlRhKpcqQsWShQgO0CkdpdxRrsXfGHN9Ps+RBFQ85+sKwL6LnGskMBKsOf98nlcVgNLFQihdTG+6UsrTueI9+EqghohUEZEgoBcw2QXnVUopdRH5Kv4i0k1EkoBmwFQRmen4fnkRmQZgjMkCHgZmAhuB8caYxPzFVkoplR/5Xe0zEZh4ge/vAzqe8/U0YFp+zqWUUsp5dOmFUkr5IC3+Sinlg7T4K6WUD9Lir5RSPkiLv1JK+SC3beksIinA7nz8inDgkJPiOJPmujKa68porivjjbkqG2MicjvIbYt/folI/OV0tnM1zXVlNNeV0VxXxpdz6bSPUkr5IC3+Sinlg7y5+I+0O8BFaK4ro7mujOa6Mj6by2vn/JVSSl2cN4/8lVJKXYTXF38ReUpEjIiE253lLBF5TUT+FpE1IjJLRMq7Qab3RGSTI9dEEQmzO9NZItJDRBJFJEdEbF2ZISIdRGSziGwTkSF2ZjmXiHwjIskist7uLOcSkYoiMk9ENjr+Hz5qdyYAEQkRkRUistaR61W7M50lIv4ikiAiUwryPF5d/EWkItbG8XvsznKe94wx9YwxscAU4CW7AwGzgRhjTD1gC/CszXnOtR7oDiy0M4SI+AOfATcC0UBvEYm2M9M5RgMd7A5xAVnAk8aYOkBT4CE3+TdLB64zxtQHYoEOItLU5kxnPYrV/r5AeXXxBz4CBnORbSPtYow5cc6XRXCDfMaYWY69FwCWYe245haMMRuNMZvtzgE0BrYZY3YYYzKAcUBXmzMBYIxZCByxO8f5jDH7jTGrHZ+fxCpqkfamAmNJdXwZ6Piw/XkoIhWATsBXBX0ury3+InIT8I8xZq3dWS5ERN4Qkb3A7bjHyP9cdwPT7Q7hhiKBved8nYQbFDJPISJRQANgub1JLI7plTVAMjDbGOMOuT7GGrDmFPSJ8rWZi90utbk88Bxwg2sT/b/cNr43xjwPPC8iz2LtdPay3ZkcxzyP9Vb9p4LOc6XZ3IBc4Hu2jxY9gYiEAhOAx85752sbY0w2EOu4vjVRRGKMMbZdMxGRzkCyMWaViLQu6PN5dPG/2ObyIlIXqAKsFRGwpjBWi0hjY8wBO7NdwBhgKi4o/rllEpG7gM5AW+PiNcBX8O9lpySg4jlfVwD22ZTFY4hIIFbh/8kY85vdec5njDkmIvOxrpnYecG8OXCTiHQEQoBiIvKjMaZvQZzMK6d9jDHrjDGljTFRxpgorCdtQ1cV/tyISI1zvrwJ2GRXlrNEpAPwDHCTMea03Xnc1EqghohUEZEgoBcw2eZMbk2s0dfXwEZjzId25zlLRCLOrmgTkULA9dj8PDTGPGuMqeCoWb2AuQVV+MFLi78HeFtE1ovI31hTU+6w/G04UBSY7ViC+oXdgc4SkW4ikgQ0A6aKyEw7cjguiD8MzMS6cDneGJNoR5bzichYYClQS0SSROQeuzM5NAfuAK5zPK7WOEa2disHzHM8B1dizfkX6NJKd6N3+CqllA/Skb9SSvkgLf5KKeWDtPgrpZQP0uKvlFI+SIu/Ukr5IC3+Sinlg7T4K6WUD9Lir5RSPuj/ACOpWmKaJZiAAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def g( x ):\n", " \n", " return np.exp( x )\n", "\n", "def f( x ):\n", " \n", " return np.sin( x )\n", "\n", "x = np.linspace( -4.0, 4.0 )\n", "pp.plot( x, f( x ) )\n", "pp.plot( x, g( f( x ) ) )\n", "\n", "pp.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Seja $\\def\\vect{\\boldsymbol}$$\\vect x^* \\in \\mathbb R^n$ um minimizador local de $f$, então, existe $\\epsilon > 0$ tal que:\n", "$$\n", "\\|\\vect y - \\vect x^* \\| \\leq \\epsilon \\Rightarrow f( \\vect y ) \\ge f( \\vect x^* ).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mas, como $g$ é crescente, temos\n", "$$\n", "f( \\vect y ) \\ge f( \\vect x^* ) \\Rightarrow g\\bigl( f( \\vect y ) \\bigr) \\ge g\\bigl( f( \\vect x^* ) \\bigr)\n", "$$\n", "além disso como $g$ é estritamente crescente,\n", "$$\n", "f( \\vect y ) > f( \\vect x^* ) \\Rightarrow g\\bigl( f( \\vect y ) \\bigr) > g\\bigl( f( \\vect x^* ) \\bigr).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Portanto:\n", "$$\n", "\\|\\vect y - \\vect x^* \\| \\leq \\epsilon \\Rightarrow f( \\vect y ) \\ge f( \\vect x^* ) \\Leftrightarrow g\\bigl( f( \\vect y ) \\bigr) \\ge g\\bigl( f( \\vect x^* ) \\bigr).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f( \\vect x ) = \\| A\\vect x- \\vect b \\| = \\sqrt{ \\sum_{i = 1}^m\\bigl( (A\\vect x)_i - b_i \\bigr)^2 }\n", "$$\n", "$$\n", "F( \\vect x ) = \\frac12f^2( \\vect x ) =\\frac12\\sum_{i = 1}^m\\bigl( (A\\vect x)_i - b_i \\bigr)^2\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\nabla F( \\vect x ) = A^T( A\\vect x - \\vect b ).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "A^T( A\\vect x^* - \\vect b ) = \\vect 0 \\Leftrightarrow A^TA\\vect x^* = A^T\\vect b.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "A \\in \\mathbb R^{m \\times n}\n", "$$\n", "$$\n", "A^T \\in \\mathbb R^{n \\times m}\n", "$$\n", "portanto,\n", "$$\n", "A^TA \\in \\mathbb R^{n \\times n}\n", "$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import sympy as sym" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def f( x_1, x_2 ):\n", " \n", " return ( x_1 ** 3 + x_2 ) ** 2 + 2 * ( x_2 - x_1 - 4 ) ** 4" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x_1 x_2\n" ] } ], "source": [ "( x_1, x_2 ) = sym.symbols( 'x_1, x_2' )\n", "print( x_1, x_2 )" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x_1**6 + 2*x_1**4 + 32*x_1**3 + 192*x_1**2 + 6*x_1*x_2**5 + 512*x_1 + x_2**6 + x_2**4*(15*x_1**2 + 2) + x_2**3*(20*x_1**3 + 6*x_1) + x_2**2*(15*x_1**4 + 6*x_1**2 + 1) + x_2*(6*x_1**5 + 2*x_1**3) + 512\n" ] } ], "source": [ "print( f( x_1, x_2 ).subs( x_1, x_1 + x_2 ).expand().factor( x_2 ) )" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "( l, d_1, d_2 ) = sym.symbols( r'\\lambda, d_1, d_2' )" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(\\lambda*d_2 + x_2 + (\\lambda*d_1 + x_1)**3)**2 + 2*(-\\lambda*d_1 + \\lambda*d_2 - x_1 + x_2 - 4)**4\n" ] } ], "source": [ "print( f( x_1 + l * d_1, x_2 + l * d_2 ) )" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "import IPython\n", "def show_expr( expr ):\n", " IPython.display.display( IPython.display.Math( sym.latex( expr ) ) )" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\\left(\\lambda d_{2} + x_{2} + \\left(\\lambda d_{1} + x_{1}\\right)^{3}\\right)^{2} + 2 \\left(- \\lambda d_{1} + \\lambda d_{2} - x_{1} + x_{2} - 4\\right)^{4}$$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "show_expr( f( x_1 + l * d_1, x_2 + l * d_2 ) )" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$2$$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "g = f( x_1 + l * d_1, x_2 + l * d_2 ).subs( x_1, 0 )\\\n", " .subs( x_2, 0 )\\\n", " .subs( d_1, 1 )\\\n", " .subs( d_2, 1 )\n", "show_expr( g.simplify().diff( l, 2 ).subs( l, 0 ) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\min_{\\vect x \\in \\mathbb R^n} \\quad \\frac12\\| A \\vect x - \\vect b \\|^2\n", "$$\n", "$$\n", "A^TA\\vect x = A^T\\vect b\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Seja $A \\in \\mathbb R^{m \\times n}$, então\n", "$$\n", "A = U\\Sigma V^T,\n", "$$\n", "onde $U\\in \\mathbb R^{m \\times m}$ e $V \\in \\mathbb R^{n \\times n}$ são ortonormais e $\\sigma \\in \\mathbb R_+^{m \\times n}$ é diagonal." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[10. 10. 32. 5.47]\n", " [10. 10. 32. 5.47]\n", " [10. 10. 32. 5.47]]\n" ] } ], "source": [ "A = np.array( [ [ 10, 10, 32, 5.47 ], [ 10, 10, 32, 5.47 ], [ 10, 10, 32, 5.47 ] ] )\n", "print( A )" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-5.77350269e-01 -8.16496581e-01 2.22044605e-16]\n", " [-5.77350269e-01 4.08248290e-01 -7.07106781e-01]\n", " [-5.77350269e-01 4.08248290e-01 7.07106781e-01]]\n", "[6.13332104e+01 6.02391913e-15 0.00000000e+00]\n", "[[-0.28240015 -0.28240015 -0.90368049 -0.15447288]\n", " [ 0.95691309 -0.03004777 -0.28877024 -0.00511986]\n", " [ 0.05827428 -0.95025412 0.2935168 -0.08642399]\n", " [-0.03422854 -0.12792339 -0.11756362 0.98419653]]\n" ] } ], "source": [ "( U, S, V ) = np.linalg.svd( A, full_matrices = True )\n", "print( U )\n", "print( S )\n", "print( V )" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1.00000000e+00 1.23838084e-16 -4.67114000e-17 1.61659634e-17]\n", " [ 1.23838084e-16 1.00000000e+00 -8.65362637e-18 -5.60905669e-18]\n", " [-4.67114000e-17 -8.65362637e-18 1.00000000e+00 2.35693104e-18]\n", " [ 1.61659634e-17 -5.60905669e-18 2.35693104e-18 1.00000000e+00]]\n" ] } ], "source": [ "print( V @ V.T )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "A^TA = ( U\\Sigma V^T )^TU\\Sigma V^T = V\\Sigma^TU^T U\\Sigma V^T = V \\overline\\Sigma V^T\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "A^TA\\vect x = A^T\\vect b\n", "$$\n", "$$\n", "V \\overline\\Sigma V^T\\vect x = V\\Sigma^TU^T\\vect b\n", "$$\n", "$$\n", "V^T V\\overline\\Sigma V^T\\vect x = V^TV\\Sigma^TU^T\\vect b\n", "$$\n", "$$\n", "\\overline\\Sigma V^T\\vect x = \\Sigma^TU^T\\vect b\n", "$$" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1 0]\n", " [0 0]\n", " [0 0]]\n", "[[1 0 0]\n", " [0 0 0]]\n" ] } ], "source": [ "s1 = np.array( [ [ 1, 0 ], [0, 0], [ 0, 0 ] ] )\n", "print( s1 )\n", "print( s1.T )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Calcular $\\vect u = \\Sigma^T U^T \\vect b$\n", "2. Resolver $\\overline\\Sigma \\vect y = \\vect u$ (a solução é dada por componentes: $y_i = u_i / \\sigma_i^2$ ou $y_i = 0$ se $\\sigma_i = 0$)\n", "3. Calcular a solução de $V^T\\vect x = \\vect y$ ou seja, $\\vect x = V\\vect y$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "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.9" } }, "nbformat": 4, "nbformat_minor": 2 }