From d7c884f2853fbb50319e5e6b8022cdeda57ae887 Mon Sep 17 00:00:00 2001 From: "Ryan C. Cooper" Date: Mon, 18 Sep 2017 13:57:54 -0400 Subject: [PATCH] added root finding lectures --- .../06_roots-1-checkpoint.ipynb | 910 ++++++++++ 06_roots-1/06_roots-1.ipynb | 910 ++++++++++ 06_roots-1/bisect.m | 37 + 06_roots-1/falsepos.m | 39 + 06_roots-1/incsearch.m | 37 + 06_roots-1/lecture_06.md | 373 ++++ 06_roots-1/octave-workspace | Bin 0 -> 1373 bytes ..._roots_and_optimization-2-checkpoint.ipynb | 1525 +++++++++++++++++ 07_roots-2/.lecture_07.md.swp | Bin 0 -> 16384 bytes 07_roots-2/.newtraph.m.swp | Bin 0 -> 12288 bytes 07_roots-2/07_roots_and_optimization-2.ipynb | 1525 +++++++++++++++++ 07_roots-2/bisect.m | 37 + 07_roots-2/car_payments.m | 17 + 07_roots-2/falsepos.m | 39 + 07_roots-2/fzerosimp.m | 41 + 07_roots-2/incsearch.m | 37 + 07_roots-2/mod_secant.m | 31 + 07_roots-2/newtraph.m | 29 + 07_roots-2/octave-workspace | Bin 0 -> 1948 bytes 19 files changed, 5587 insertions(+) create mode 100644 06_roots-1/.ipynb_checkpoints/06_roots-1-checkpoint.ipynb create mode 100644 06_roots-1/06_roots-1.ipynb create mode 100644 06_roots-1/bisect.m create mode 100644 06_roots-1/falsepos.m create mode 100644 06_roots-1/incsearch.m create mode 100644 06_roots-1/lecture_06.md create mode 100644 06_roots-1/octave-workspace create mode 100644 07_roots-2/.ipynb_checkpoints/07_roots_and_optimization-2-checkpoint.ipynb create mode 100644 07_roots-2/.lecture_07.md.swp create mode 100644 07_roots-2/.newtraph.m.swp create mode 100644 07_roots-2/07_roots_and_optimization-2.ipynb create mode 100644 07_roots-2/bisect.m create mode 100644 07_roots-2/car_payments.m create mode 100644 07_roots-2/falsepos.m create mode 100644 07_roots-2/fzerosimp.m create mode 100644 07_roots-2/incsearch.m create mode 100644 07_roots-2/mod_secant.m create mode 100644 07_roots-2/newtraph.m create mode 100644 07_roots-2/octave-workspace diff --git a/06_roots-1/.ipynb_checkpoints/06_roots-1-checkpoint.ipynb b/06_roots-1/.ipynb_checkpoints/06_roots-1-checkpoint.ipynb new file mode 100644 index 0000000..e95edef --- /dev/null +++ b/06_roots-1/.ipynb_checkpoints/06_roots-1-checkpoint.ipynb @@ -0,0 +1,910 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true, + "scrolled": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Roots of Nonlinear functions\n", + "## Bracketing ch. 5" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## Not always possible to solve for a given variable. \n", + "\n", + "### Freefall example:\n", + "If an object, with drag coefficient of 0.25 kg/m reaches a velocity of 36 m/s after 4 seconds of freefalling, what is its mass?\n", + "\n", + "$v(t)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)$\n", + "\n", + "Cannot solve for m \n", + "\n", + "Instead, solve the problem by creating a new function f(m) where\n", + "\n", + "$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$. \n", + "\n", + "When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t140\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t160\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t180\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "setdefaults\n", + "g=9.81; % acceleration due to gravity\n", + "m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n", + "c_d=0.25; % drag coefficient\n", + "t=4; % at time = 4 seconds\n", + "v=36; % speed must be 36 m/s\n", + "f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n", + "\n", + "plot(m,f_m(m),m,zeros(length(m),1))\n", + "axis([45 200 -5 1])" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 0.045626\r\n" + ] + } + ], + "source": [ + "f_m(145)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0\n", + "\n", + "Better methods are the \n", + "1. Bracketing methods\n", + "2. Open methods\n", + "\n", + "Both need an initial guess. \n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Incremental method (Brute force)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function xb = incsearch(func,xmin,xmax,ns)\n", + "% incsearch: incremental search root locator\n", + "% xb = incsearch(func,xmin,xmax,ns):\n", + "% finds brackets of x that contain sign changes\n", + "% of a function on an interval\n", + "% input:\n", + "% func = name of function\n", + "% xmin, xmax = endpoints of interval\n", + "% ns = number of subintervals (default = 50)\n", + "% output:\n", + "% xb(k,1) is the lower bound of the kth sign change\n", + "% xb(k,2) is the upper bound of the kth sign change\n", + "% If no brackets found, xb = [].\n", + "if nargin < 3, error('at least 3 arguments required'), end\n", + "if nargin < 4, ns = 50; end %if ns blank set to 50\n", + "% Incremental search\n", + "x = linspace(xmin,xmax,ns);\n", + "f = func(x);\n", + "nb = 0; xb = []; %xb is null unless sign change detected\n", + "for k = 1:length(x)-1\n", + " if sign(f(k)) ~= sign(f(k+1)) %check for sign change\n", + " nb = nb + 1;\n", + " xb(nb,1) = x(k);\n", + " xb(nb,2) = x(k+1);\n", + " end\n", + "end\n", + "if isempty(xb) %display that no brackets were found\n", + " fprintf('no brackets found\\n')\n", + " fprintf('check interval or increase ns\\n')\n", + "else\n", + " fprintf('number of brackets: %i\\n',nb) %display number of brackets\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/course_git_F2017/06_roots and optimization/incsearch.m\n", + "\n", + " incsearch: incremental search root locator\n", + " xb = incsearch(func,xmin,xmax,ns):\n", + " finds brackets of x that contain sign changes\n", + " of a function on an interval\n", + " input:\n", + " func = name of function\n", + " xmin, xmax = endpoints of interval\n", + " ns = number of subintervals (default = 50)\n", + " output:\n", + " xb(k,1) is the lower bound of the kth sign change\n", + " xb(k,2) is the upper bound of the kth sign change\n", + " If no brackets found, xb = [].\n", + "\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help incsearch" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "number of brackets: 1\n", + "ans =\n", + "\n", + " 142.73 142.83\n", + "\n" + ] + } + ], + "source": [ + "incsearch(f_m,140, 150,100)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Bisection method" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Divide interval in half until error is reduced to some level\n", + "\n", + "in previous example of freefall, choose x_l=50, x_u=200\n", + "\n", + "x_r = (50+200)/2 = 125\n", + "\n", + "f_m(125) = -0.408\n", + "\n", + "x_r= (125+200)/2 = 162.5\n", + "\n", + "f_m(162.5) = 0.3594\n", + "\n", + "x_r = (125+162.5)/2=143.75\n", + "\n", + "f_m(143.75)= 0.0206" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 134.38\n", + "interval left f(x_l)= -0.4,f(x_r)= -0.2\n", + "interval right f(x_r)= -0.2,f(x_u)= 0.0\n", + "ans = -0.18060\n" + ] + } + ], + "source": [ + "x_l=125; x_u=143.75;\n", + "x_r=(x_l+x_u)/2\n", + "fprintf('interval left f(x_l)= %1.1f,f(x_r)= %1.1f\\n',f_m(x_l),f_m(x_r))\n", + "fprintf('interval right f(x_r)= %1.1f,f(x_u)= %1.1f\\n',f_m(x_r),f_m(x_u))\n", + "f_m(x_r)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Bisect Function" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Much better root locator, with 4 iterations, our function is already close to zero\n", + "\n", + "Automate this with a function:\n", + "`bisect.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)\n", + "% bisect: root location zeroes\n", + "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + "% uses bisection method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% xl, xu = lower and upper guesses\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by func\n", + "% output:\n", + "% root = real root\n", + "% fx = function value at root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "test = func(xl,varargin{:})*func(xu,varargin{:});\n", + "if test>0,error('no sign change'),end\n", + "if nargin<4|isempty(es), es=0.0001;end\n", + "if nargin<5|isempty(maxit), maxit=50;end\n", + "iter = 0; xr = xl; ea = 100;\n", + "while (1)\n", + " xrold = xr;\n", + " xr = (xl + xu)/2;\n", + " iter = iter + 1;\n", + " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n", + " test = func(xl,varargin{:})*func(xr,varargin{:});\n", + " if test < 0\n", + " xu = xr;\n", + " elseif test > 0\n", + " xl = xr;\n", + " else\n", + " ea = 0;\n", + " end\n", + " if ea <= es | iter >= maxit,break,end\n", + "end\n", + "root = xr; fx = func(xr, varargin{:});" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mass_at_36ms = 142.74\r\n" + ] + } + ], + "source": [ + "Mass_at_36ms=bisect(f_m,50,200)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## False position (linear interpolation)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner\n", + "\n", + "$ x_{r} = x_{u} - \\frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t\t \n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%xl=50; xu=200; \n", + "xu=xr;\n", + "xl=xl;\n", + "xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu));\n", + "plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## False Position" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Much better root locator, with 4 iterations, our function is already close to zero\n", + "\n", + "Automate this with a function:\n", + "`falsepos.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)\n", + "% falsepos: root location zeroes\n", + "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + "% uses false position method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% xl, xu = lower and upper guesses\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by func\n", + "% output:\n", + "% root = real root\n", + "% fx = function value at root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "test = func(xl,varargin{:})*func(xu,varargin{:});\n", + "if test>0,error('no sign change'),end\n", + "if nargin<4|isempty(es), es=0.0001;end\n", + "if nargin<5|isempty(maxit), maxit=50;end\n", + "iter = 0; xr = xl; ea = 100;\n", + "while (1)\n", + " xrold = xr;\n", + " % xr = (xl + xu)/2; % bisect method\n", + " xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method\n", + " iter = iter + 1;\n", + " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n", + " test = func(xl,varargin{:})*func(xr,varargin{:});\n", + " if test < 0\n", + " xu = xr;\n", + " elseif test > 0\n", + " xl = xr;\n", + " else\n", + " ea = 0;\n", + " end\n", + " if ea <= es | iter >= maxit,break,end\n", + "end\n", + "root = xr; fx = func(xr, varargin{:});" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/06_roots-1/06_roots-1.ipynb b/06_roots-1/06_roots-1.ipynb new file mode 100644 index 0000000..e95edef --- /dev/null +++ b/06_roots-1/06_roots-1.ipynb @@ -0,0 +1,910 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true, + "scrolled": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Roots of Nonlinear functions\n", + "## Bracketing ch. 5" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## Not always possible to solve for a given variable. \n", + "\n", + "### Freefall example:\n", + "If an object, with drag coefficient of 0.25 kg/m reaches a velocity of 36 m/s after 4 seconds of freefalling, what is its mass?\n", + "\n", + "$v(t)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)$\n", + "\n", + "Cannot solve for m \n", + "\n", + "Instead, solve the problem by creating a new function f(m) where\n", + "\n", + "$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$. \n", + "\n", + "When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t140\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t160\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t180\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "setdefaults\n", + "g=9.81; % acceleration due to gravity\n", + "m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n", + "c_d=0.25; % drag coefficient\n", + "t=4; % at time = 4 seconds\n", + "v=36; % speed must be 36 m/s\n", + "f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n", + "\n", + "plot(m,f_m(m),m,zeros(length(m),1))\n", + "axis([45 200 -5 1])" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 0.045626\r\n" + ] + } + ], + "source": [ + "f_m(145)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0\n", + "\n", + "Better methods are the \n", + "1. Bracketing methods\n", + "2. Open methods\n", + "\n", + "Both need an initial guess. \n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Incremental method (Brute force)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function xb = incsearch(func,xmin,xmax,ns)\n", + "% incsearch: incremental search root locator\n", + "% xb = incsearch(func,xmin,xmax,ns):\n", + "% finds brackets of x that contain sign changes\n", + "% of a function on an interval\n", + "% input:\n", + "% func = name of function\n", + "% xmin, xmax = endpoints of interval\n", + "% ns = number of subintervals (default = 50)\n", + "% output:\n", + "% xb(k,1) is the lower bound of the kth sign change\n", + "% xb(k,2) is the upper bound of the kth sign change\n", + "% If no brackets found, xb = [].\n", + "if nargin < 3, error('at least 3 arguments required'), end\n", + "if nargin < 4, ns = 50; end %if ns blank set to 50\n", + "% Incremental search\n", + "x = linspace(xmin,xmax,ns);\n", + "f = func(x);\n", + "nb = 0; xb = []; %xb is null unless sign change detected\n", + "for k = 1:length(x)-1\n", + " if sign(f(k)) ~= sign(f(k+1)) %check for sign change\n", + " nb = nb + 1;\n", + " xb(nb,1) = x(k);\n", + " xb(nb,2) = x(k+1);\n", + " end\n", + "end\n", + "if isempty(xb) %display that no brackets were found\n", + " fprintf('no brackets found\\n')\n", + " fprintf('check interval or increase ns\\n')\n", + "else\n", + " fprintf('number of brackets: %i\\n',nb) %display number of brackets\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/course_git_F2017/06_roots and optimization/incsearch.m\n", + "\n", + " incsearch: incremental search root locator\n", + " xb = incsearch(func,xmin,xmax,ns):\n", + " finds brackets of x that contain sign changes\n", + " of a function on an interval\n", + " input:\n", + " func = name of function\n", + " xmin, xmax = endpoints of interval\n", + " ns = number of subintervals (default = 50)\n", + " output:\n", + " xb(k,1) is the lower bound of the kth sign change\n", + " xb(k,2) is the upper bound of the kth sign change\n", + " If no brackets found, xb = [].\n", + "\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help incsearch" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "number of brackets: 1\n", + "ans =\n", + "\n", + " 142.73 142.83\n", + "\n" + ] + } + ], + "source": [ + "incsearch(f_m,140, 150,100)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Bisection method" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Divide interval in half until error is reduced to some level\n", + "\n", + "in previous example of freefall, choose x_l=50, x_u=200\n", + "\n", + "x_r = (50+200)/2 = 125\n", + "\n", + "f_m(125) = -0.408\n", + "\n", + "x_r= (125+200)/2 = 162.5\n", + "\n", + "f_m(162.5) = 0.3594\n", + "\n", + "x_r = (125+162.5)/2=143.75\n", + "\n", + "f_m(143.75)= 0.0206" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 134.38\n", + "interval left f(x_l)= -0.4,f(x_r)= -0.2\n", + "interval right f(x_r)= -0.2,f(x_u)= 0.0\n", + "ans = -0.18060\n" + ] + } + ], + "source": [ + "x_l=125; x_u=143.75;\n", + "x_r=(x_l+x_u)/2\n", + "fprintf('interval left f(x_l)= %1.1f,f(x_r)= %1.1f\\n',f_m(x_l),f_m(x_r))\n", + "fprintf('interval right f(x_r)= %1.1f,f(x_u)= %1.1f\\n',f_m(x_r),f_m(x_u))\n", + "f_m(x_r)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Bisect Function" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Much better root locator, with 4 iterations, our function is already close to zero\n", + "\n", + "Automate this with a function:\n", + "`bisect.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)\n", + "% bisect: root location zeroes\n", + "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + "% uses bisection method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% xl, xu = lower and upper guesses\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by func\n", + "% output:\n", + "% root = real root\n", + "% fx = function value at root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "test = func(xl,varargin{:})*func(xu,varargin{:});\n", + "if test>0,error('no sign change'),end\n", + "if nargin<4|isempty(es), es=0.0001;end\n", + "if nargin<5|isempty(maxit), maxit=50;end\n", + "iter = 0; xr = xl; ea = 100;\n", + "while (1)\n", + " xrold = xr;\n", + " xr = (xl + xu)/2;\n", + " iter = iter + 1;\n", + " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n", + " test = func(xl,varargin{:})*func(xr,varargin{:});\n", + " if test < 0\n", + " xu = xr;\n", + " elseif test > 0\n", + " xl = xr;\n", + " else\n", + " ea = 0;\n", + " end\n", + " if ea <= es | iter >= maxit,break,end\n", + "end\n", + "root = xr; fx = func(xr, varargin{:});" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mass_at_36ms = 142.74\r\n" + ] + } + ], + "source": [ + "Mass_at_36ms=bisect(f_m,50,200)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## False position (linear interpolation)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner\n", + "\n", + "$ x_{r} = x_{u} - \\frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t\t \n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%xl=50; xu=200; \n", + "xu=xr;\n", + "xl=xl;\n", + "xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu));\n", + "plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## False Position" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Much better root locator, with 4 iterations, our function is already close to zero\n", + "\n", + "Automate this with a function:\n", + "`falsepos.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)\n", + "% falsepos: root location zeroes\n", + "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + "% uses false position method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% xl, xu = lower and upper guesses\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by func\n", + "% output:\n", + "% root = real root\n", + "% fx = function value at root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "test = func(xl,varargin{:})*func(xu,varargin{:});\n", + "if test>0,error('no sign change'),end\n", + "if nargin<4|isempty(es), es=0.0001;end\n", + "if nargin<5|isempty(maxit), maxit=50;end\n", + "iter = 0; xr = xl; ea = 100;\n", + "while (1)\n", + " xrold = xr;\n", + " % xr = (xl + xu)/2; % bisect method\n", + " xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method\n", + " iter = iter + 1;\n", + " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n", + " test = func(xl,varargin{:})*func(xr,varargin{:});\n", + " if test < 0\n", + " xu = xr;\n", + " elseif test > 0\n", + " xl = xr;\n", + " else\n", + " ea = 0;\n", + " end\n", + " if ea <= es | iter >= maxit,break,end\n", + "end\n", + "root = xr; fx = func(xr, varargin{:});" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/06_roots-1/bisect.m b/06_roots-1/bisect.m new file mode 100644 index 0000000..c09ffbf --- /dev/null +++ b/06_roots-1/bisect.m @@ -0,0 +1,37 @@ +function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + xr = (xl + xu)/2; + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); \ No newline at end of file diff --git a/06_roots-1/falsepos.m b/06_roots-1/falsepos.m new file mode 100644 index 0000000..0a3477c --- /dev/null +++ b/06_roots-1/falsepos.m @@ -0,0 +1,39 @@ +function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + xr = (xl + xu)/2; + % xr = (xl + xu)/2; % bisect method + xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); diff --git a/06_roots-1/incsearch.m b/06_roots-1/incsearch.m new file mode 100644 index 0000000..bd82554 --- /dev/null +++ b/06_roots-1/incsearch.m @@ -0,0 +1,37 @@ +function xb = incsearch(func,xmin,xmax,ns) +% incsearch: incremental search root locator +% xb = incsearch(func,xmin,xmax,ns): +% finds brackets of x that contain sign changes +% of a function on an interval +% input: +% func = name of function +% xmin, xmax = endpoints of interval +% ns = number of subintervals (default = 50) +% output: +% xb(k,1) is the lower bound of the kth sign change +% xb(k,2) is the upper bound of the kth sign change +% If no brackets found, xb = []. +if nargin < 3, error('at least 3 arguments required'), end +if nargin < 4, ns = 50; end %if ns blank set to 50 +% Incremental search +x = linspace(xmin,xmax,ns); +f = func(x); +nb = 0; xb = []; %xb is null unless sign change detected +%for k = 1:length(x)-1 +% if sign(f(k)) ~= sign(f(k+1)) %check for sign change +% nb = nb + 1; +% xb(nb,1) = x(k); +% xb(nb,2) = x(k+1); +% end +%end +sign_change = diff(sign(f)); +[~,i_change] = find(sign_change~=0); +nb=length(i_change); +xb=[x(i_change)',x(i_change+1)']; + +if isempty(xb) %display that no brackets were found + fprintf('no brackets found\n') + fprintf('check interval or increase ns\n') +else + fprintf('number of brackets: %i\n',nb) %display number of brackets +end diff --git a/06_roots-1/lecture_06.md b/06_roots-1/lecture_06.md new file mode 100644 index 0000000..b3d2266 --- /dev/null +++ b/06_roots-1/lecture_06.md @@ -0,0 +1,373 @@ + + +```octave +%plot --format svg +``` + +my_caller.m +```matlab +function [vx,vy] = my_caller(max_time) + N=100; + t=linspace(0,max_time,N); + [x,y]=my_function(max_time); + vx=diff(x)./diff(t); + vy=diff(y)./diff(t); +end +``` + +my_function.m +```matlab +function [x,y] = my_function(max_time) + N=100; + t=linspace(0,max_time,N); + x=t.^2; + y=2*t; +end +``` + +In order to use `my_caller.m` where does `my_function.m` need to be saved? +![responses](q1.png) + + +What cool personal projects are you working on? +While we delve deeper into Matlab functions, could we review some of the basic logic +operators it uses and command codes. + +I still dont know when these forms are technically due. + + -by the following lecture + +I'm having trouble interfacing Atom with GitHub. Is there a simple tutorial for this? + + -Mac? Seems there could be a bug that folks are working on + +What are the bear necessities of life? +please go over how to "submit" the homeworks because it is still confusing + +Do you prefer Matlab or Octave? + + -octave is my preference, but Matlab has some benefits + +Would you consider a country to be open-source? + + -?? + +Is there a way to download matlab for free? + + -not legally + +how do you add files to current folder in matlab? + + -you can do this either through a file browser or cli + +How should Homework 2 be submitted? By simply putting the function into the homework_1 +repository? + + -yes + +How can we tell that these forms are being received? + + -when you hit submit, the form says "form received" + +can you save scripted outputs from matlab/octave as an excel file? + + -yes, easy way is open a file with a .csv extension then fprintf and separate everything with commas, harder way is to use the `xlswrite` + + +Also, can you update your notes to show what happens when these things are run, as you do +in class?" + + -I always update the lecture notes after class so they should display what we did in class + +I have a little difficulty following along in class on my laptop when you have programs +pre-written. Maybe if you posted those codes on Github so I could copy them when you +switch to different desktops I would be able to follow along better. + +Kirk or Picard? + + -Kirk + +Who is our TA? + + -Peiyu Zhang peiyu.zhang@uconn.edu + +Can we download libraries of data like thermodynamic tables into matlab? + +-YES! [Matlab Steam Tables](http://bit.ly/2kZygu8) + +Will we use the Simulink addition to Matlab? I found it interesting and useful for +evaluating ODEs in Linear systems. + + -not in this class, everything in simulink has a matlab script/function that can be substituted, but many times its hidden by the gui. Here we want to look directly at our solvers + + +# Roots and Optimization +## Bracketing ch. 5 + +When you are given a function, numerical or analytical, it's not always possible to solve directly for a given variable. + +Even for the freefall example we first explored, + +$v(t)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)$ + +There is no way to solve for m in terms of the other variables. + +Instead, we can solve the problem by creating a new function f(m) where + +$f(m)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)-v(t)$. + +When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity) + + +```octave +setdefaults +g=9.81; % acceleration due to gravity +m=linspace(50, 200,100); % possible values for mass 50 to 200 kg +c_d=0.25; % drag coefficient +t=4; % at time = 4 seconds +v=36; % speed must be 36 m/s +f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m + +plot(m,f_m(m),m,zeros(length(m),1)) +axis([45 200 -5 1]) +``` + + +![svg](lecture_06_files/lecture_06_5_0.svg) + + + +```octave +f_m(145) +``` + + ans = 0.045626 + + +Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0 + +Better methods are the +1. Bracketing methods +2. Open methods + +Both need an initial guess. + + +## Incremental method (Brute force) + +You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. + +```matlab +function xb = incsearch(func,xmin,xmax,ns) +% incsearch: incremental search root locator +% xb = incsearch(func,xmin,xmax,ns): +% finds brackets of x that contain sign changes +% of a function on an interval +% input: +% func = name of function +% xmin, xmax = endpoints of interval +% ns = number of subintervals (default = 50) +% output: +% xb(k,1) is the lower bound of the kth sign change +% xb(k,2) is the upper bound of the kth sign change +% If no brackets found, xb = []. +if nargin < 3, error('at least 3 arguments required'), end +if nargin < 4, ns = 50; end %if ns blank set to 50 +% Incremental search +x = linspace(xmin,xmax,ns); +f = func(x); +nb = 0; xb = []; %xb is null unless sign change detected +for k = 1:length(x)-1 + if sign(f(k)) ~= sign(f(k+1)) %check for sign change + nb = nb + 1; + xb(nb,1) = x(k); + xb(nb,2) = x(k+1); + end +end +if isempty(xb) %display that no brackets were found + fprintf('no brackets found\n') + fprintf('check interval or increase ns\n') +else + fprintf('number of brackets: %i\n',nb) %display number of brackets +end +``` + + +```octave +help incsearch +``` + + 'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/me3255_S2017/lecture_06/incsearch.m + + incsearch: incremental search root locator + xb = incsearch(func,xmin,xmax,ns): + finds brackets of x that contain sign changes + of a function on an interval + input: + func = name of function + xmin, xmax = endpoints of interval + ns = number of subintervals (default = 50) + output: + xb(k,1) is the lower bound of the kth sign change + xb(k,2) is the upper bound of the kth sign change + If no brackets found, xb = []. + + + Additional help for built-in functions and operators is + available in the online version of the manual. Use the command + 'doc ' to search the manual index. + + Help and information about Octave is also available on the WWW + at http://www.octave.org and via the help@octave.org + mailing list. + + + +```octave +incsearch(f_m,50, 200,55) +``` + + no brackets found + check interval or increase ns + ans = [](1x0) + + +## Bisection method + +Divide interval in half until error is reduced to some level + +in previous example of freefall, choose x_l=50, x_u=200 + +x_r = (50+200)/2 = 125 + +f_m(125) = -0.408 + +x_r= (125+200)/2 = 162.5 + +f_m(162.5) = 0.3594 + +x_r = (125+162.5)/2=143.75 + +f_m(143.75)= 0.0206 + + +```octave +f_m(143.75) +``` + + ans = 0.020577 + + +Much better root locator, with 4 iterations, our function is already close to zero + +Automate this with a function: +`bisect.m` + +```matlab +function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + xr = (xl + xu)/2; + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); +``` + +## False position (linear interpolation) + +Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner + +$ x_{r} = x_{u} - \frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$ + + +```octave +xl=50; xu=200; +xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); + +plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0) +``` + + +![svg](lecture_06_files/lecture_06_18_0.svg) + + +Much better root locator, with 4 iterations, our function is already close to zero + +Automate this with a function: +`falsepos.m` + +```matlab +function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin) +% falsepos: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses false position method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + % xr = (xl + xu)/2; % bisect method + xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); +``` + + +```octave + +``` diff --git a/06_roots-1/octave-workspace b/06_roots-1/octave-workspace new file mode 100644 index 0000000000000000000000000000000000000000..05454dadf909ce344f6be3b8fe0223c066b434cc GIT binary patch literal 1373 zcmbu9Pe>F|9LFc34EnP|2Z;_o3T$`5l_>TQw*4p!OxJbQRd?NW=I#6ycE;Kr#llNe zhbZV!&_ODw=wQ*oQ+TQ_Q50nmArVw`km}IUliPdm_k$44O9L~^{J!t~{(QbOlQBa+ zX}2D1?bdc^nwI9Zn$AQ0P|~l7M@_456_=XHE4zRRy$0 zrB41sJ+9P^GZP32+UH?)Y-*)mugj8>poA`tPI%tXJpKDwCtQrAmcLv$1t*(I;X+3j zY?I-JZg|x)Q<;0&4ObKMOLwoNp(nX$T~GHwqg)=z!21rlA%p8Wd*R`%JbH2c_AH3a zB10DUy~%<~ljy_!*ZM#d5taHdPqrV_^0|Jvkrst^4}fS!6c~WPoH#xVfGVHNfok(+ z4rCXpa;<>9MLW8{fP z9`1=OOZy7Mr2w`_zxcgCd`94aEbw52I2A#a^@_x+h#cZ2;#LA*4)&`={6>+(<55u9 zBRX<8t%DkOLC1Oa7{t}Uxy~7&@Xi`I-zS4ObDXoC6K{_5-r~faEp_?K~hy-Yro zk;Bt6`B6q5%p_k-mc7B4*BdNhc1`=c9HiZmwb1T!%mO<_mD&2(R<*b4{tqsAIARc i+mwHWDZL`DtG{~{0A2m)ru6Fi`fG#qiS_km4u1fOPS=V6 literal 0 HcmV?d00001 diff --git a/07_roots-2/.ipynb_checkpoints/07_roots_and_optimization-2-checkpoint.ipynb b/07_roots-2/.ipynb_checkpoints/07_roots_and_optimization-2-checkpoint.ipynb new file mode 100644 index 0000000..4a2ba51 --- /dev/null +++ b/07_roots-2/.ipynb_checkpoints/07_roots_and_optimization-2-checkpoint.ipynb @@ -0,0 +1,1525 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Roots: Open methods\n", + "## Newton-Raphson" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "First-order approximation for the location of the root (i.e. assume the slope at the given point is constant, what is the solution when f(x)=0)\n", + "\n", + "$f'(x_{i})=\\frac{f(x_{i})-0}{x_{i}-x_{i+1}}$\n", + "\n", + "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n", + "\n", + "Use Newton-Raphson to find solution when $e^{-x}=x$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "Find x when $e^{-x}=x$" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.50000\n", + "error_approx = 1\n" + ] + } + ], + "source": [ + "f= @(x) exp(-x)-x;\n", + "df= @(x) -exp(-x)-1;\n", + "\n", + "x_i= 0;\n", + "x_r = x_i-f(x_i)/df(x_i)\n", + "error_approx = abs((x_r-x_i)/x_r)\n", + "x_i=x_r;\n" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.56714\n", + "error_approx = 0.0014673\n" + ] + } + ], + "source": [ + "x_r = x_i-f(x_i)/df(x_i)\n", + "error_approx = abs((x_r-x_i)/x_r)\n", + "x_i=x_r;" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.56714\n", + "error_approx = 2.2106e-07\n" + ] + } + ], + "source": [ + "x_r = x_i-f(x_i)/df(x_i)\n", + "error_approx = abs((x_r-x_i)/x_r)\n", + "x_i=x_r;" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.56714\n", + "error_approx = 5.0897e-15\n" + ] + } + ], + "source": [ + "x_r = x_i-f(x_i)/df(x_i)\n", + "error_approx = abs((x_r-x_i)/x_r)\n", + "x_i=x_r;" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Create function `newtraph.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,varargin)\n", + "% newtraph: Newton-Raphson root location zeroes\n", + "% [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,p1,p2,...):\n", + "% uses Newton-Raphson method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% dfunc = name of derivative of function\n", + "% xr = initial guess\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by function\n", + "% output:\n", + "% root = real root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "if nargin<4 || isempty(es),es=0.0001;end\n", + "if nargin<5 || isempty(maxit),maxit=50;end\n", + "iter = 0;\n", + "while (1)\n", + " xrold = xr;\n", + " xr = xr - func(xr)/dfunc(xr);\n", + " iter = iter + 1;\n", + " if xr ~= 0 \n", + " ea = abs((xr - xrold)/xr) * 100;\n", + " end\n", + " if ea <= es || iter >= maxit, break, end\n", + "end\n", + "root = xr;" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "In the freefall example, we created a function f(m) that when f(m)=0, then the mass had been chosen such that at t=4 s, the velocity is 36 m/s. \n", + "\n", + "$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$.\n", + "\n", + "to use the Newton-Raphson method, we need the derivative $\\frac{df}{dm}$\n", + "\n", + "$\\frac{df}{dm}=\\frac{1}{2}\\sqrt{\\frac{g}{mc_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-\n", + "\\frac{g}{2m}\\mathrm{sech}^{2}(\\sqrt{\\frac{gc_{d}}{m}}t)$" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "setdefaults\n", + "g=9.81; % acceleration due to gravity\n", + "m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n", + "c_d=0.25; % drag coefficient\n", + "t=4; % at time = 4 seconds\n", + "v=36; % speed must be 36 m/s\n", + "f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n", + "df_m = @(m) 1/2*sqrt(g./m/c_d).*tanh(sqrt(g*c_d./m)*t)-g/2./m*sech(sqrt(g*c_d./m)*t).^2;" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "root = 142.74\n", + "ea = 1.8806e-04\n", + "iter = 50\n" + ] + } + ], + "source": [ + "[root,ea,iter]=newtraph(f_m,df_m,50,0.0001)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Secant Methods" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Not always able to evaluate the derivative. Approximation of derivative:\n", + "\n", + "$f'(x_{i})=\\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}$\n", + "\n", + "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n", + "\n", + "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{\\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}}=\n", + " x_{i}-\\frac{f(x_{i})(x_{i-1}-x_{i})}{f(x_{i-1})-f(x_{i})}$\n", + " \n", + "What values should $x_{i}$ and $x_{i-1}$ take?\n", + "\n", + "To reduce arbitrary selection of variables, use the" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Modified Secant method" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Change the x evaluations to a perturbation $\\delta$. \n", + "\n", + "$x_{i+1}=x_{i}-\\frac{f(x_{i})(\\delta x_{i})}{f(x_{i}+\\delta x_{i})-f(x_{i})}$" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "root = 142.74\n", + "ea = 3.0615e-07\n", + "iter = 7\n" + ] + } + ], + "source": [ + "[root,ea,iter]=mod_secant(f_m,1,50,0.00001)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 1.1185e+04\r\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t25000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tprinciple amount left ($)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\ttime (years)\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "car_payments(400,30000,0.05,5,1)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Amt_numerical = 5467.0\n", + "ans = 3.9755e-04\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t400000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t500000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t600000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t700000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t25\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tprinciple amount left ($)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\ttime (years)\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "Amt_numerical=mod_secant(@(A) car_payments(A,700000,0.0875,30,0),1e-6,50,0.001)\n", + "car_payments(Amt_numerical,700000,0.0875,30,1)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 1.9681e+06\r\n" + ] + } + ], + "source": [ + "Amt_numerical*12*30" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Amortization\n", + "\n", + "Amortization calculation makes the same calculation for the monthly payment amount, A, paying off the principle amount, P, over n pay periods with monthly interest rate, r. " + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Amt = 566.14\r\n" + ] + } + ], + "source": [ + "% Amortization calculation\n", + "A = @(P,r,n) P*(r*(1+r)^n)./((1+r)^n-1);\n", + "Amt=A(30000,0.05/12,5*12)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Matlab's function\n", + "\n", + "Matlab and Octave combine bracketing and open methods in the `fzero` function. " + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'fzero' is a function from the file /usr/share/octave/4.0.0/m/optimization/fzero.m\n", + "\n", + " -- Function File: fzero (FUN, X0)\n", + " -- Function File: fzero (FUN, X0, OPTIONS)\n", + " -- Function File: [X, FVAL, INFO, OUTPUT] = fzero (...)\n", + " Find a zero of a univariate function.\n", + "\n", + " FUN is a function handle, inline function, or string containing the\n", + " name of the function to evaluate.\n", + "\n", + " X0 should be a two-element vector specifying two points which\n", + " bracket a zero. In other words, there must be a change in sign of\n", + " the function between X0(1) and X0(2). More mathematically, the\n", + " following must hold\n", + "\n", + " sign (FUN(X0(1))) * sign (FUN(X0(2))) <= 0\n", + "\n", + " If X0 is a single scalar then several nearby and distant values are\n", + " probed in an attempt to obtain a valid bracketing. If this is not\n", + " successful, the function fails.\n", + "\n", + " OPTIONS is a structure specifying additional options. Currently,\n", + " 'fzero' recognizes these options: \"FunValCheck\", \"OutputFcn\",\n", + " \"TolX\", \"MaxIter\", \"MaxFunEvals\". For a description of these\n", + " options, see *note optimset: XREFoptimset.\n", + "\n", + " On exit, the function returns X, the approximate zero point and\n", + " FVAL, the function value thereof.\n", + "\n", + " INFO is an exit flag that can have these values:\n", + "\n", + " * 1 The algorithm converged to a solution.\n", + "\n", + " * 0 Maximum number of iterations or function evaluations has\n", + " been reached.\n", + "\n", + " * -1 The algorithm has been terminated from user output\n", + " function.\n", + "\n", + " * -5 The algorithm may have converged to a singular point.\n", + "\n", + " OUTPUT is a structure containing runtime information about the\n", + " 'fzero' algorithm. Fields in the structure are:\n", + "\n", + " * iterations Number of iterations through loop.\n", + "\n", + " * nfev Number of function evaluations.\n", + "\n", + " * bracketx A two-element vector with the final bracketing of the\n", + " zero along the x-axis.\n", + "\n", + " * brackety A two-element vector with the final bracketing of the\n", + " zero along the y-axis.\n", + "\n", + " See also: optimset, fsolve.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help fzero" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 563.79\r\n" + ] + } + ], + "source": [ + "fzero(@(A) car_payments(A,30000,0.05,5,0),500)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Comparison of Solvers\n", + "\n", + "It's helpful to compare to the convergence of different routines to see how quickly you find a solution. \n", + "\n", + "Comparing the freefall example\n" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "warning: axis: omitting non-positive data in log plot\n", + "warning: called from\n", + " __line__ at line 120 column 16\n", + " line at line 56 column 8\n", + " __plt__>__plt2vv__ at line 500 column 10\n", + " __plt__>__plt2__ at line 246 column 14\n", + " __plt__ at line 133 column 15\n", + " semilogy at line 60 column 10\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t250\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t350\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t400\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "iterations = linspace(1,400,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f_m,df_m,300,0,iterations(i));\n", + " [root_ms,ea_ms(i),iter_ms]=mod_secant(f_m,1e-6,300,0,iterations(i));\n", + " [root_fp,ea_fp(i),iter_fp]=falsepos(f_m,1,300,0,iterations(i));\n", + " [root_bs,ea_bs(i),iter_bs]=bisect(f_m,1,300,0,iterations(i));\n", + "end\n", + "\n", + "setdefaults\n", + "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n", + "legend('newton-raphson','mod-secant','false point','bisection')" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ea_nr =\n", + "\n", + " Columns 1 through 8:\n", + "\n", + " 6.36591 0.06436 0.00052 0.00000 0.00000 0.00000 0.00000 0.00000\n", + "\n", + " Columns 9 through 16:\n", + "\n", + " 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n", + "\n", + " Columns 17 through 20:\n", + "\n", + " 0.00000 0.00000 0.00000 0.00000\n", + "\n" + ] + } + ], + "source": [ + "ea_nr" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "warning: axis: omitting non-positive data in log plot\n", + "warning: called from\n", + " __line__ at line 120 column 16\n", + " line at line 56 column 8\n", + " __plt__>__plt2vv__ at line 500 column 10\n", + " __plt__>__plt2__ at line 246 column 14\n", + " __plt__ at line 133 column 15\n", + " semilogy at line 60 column 10\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "f= @(x) x^10-1;\n", + "df=@(x) 10*x^9;\n", + "iterations = linspace(1,50,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i));\n", + " [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i));\n", + " [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i));\n", + " [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i));\n", + "end\n", + " \n", + "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n", + "legend('newton-raphson','mod-secant','false point','bisection')" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ea_nr =\n", + "\n", + " Columns 1 through 7:\n", + "\n", + " 99.03195 11.11111 11.11111 11.11111 11.11111 11.11111 11.11111\n", + "\n", + " Columns 8 through 14:\n", + "\n", + " 11.11111 11.11111 11.11111 11.11109 11.11052 11.10624 10.99684\n", + "\n", + " Columns 15 through 20:\n", + "\n", + " 8.76956 2.12993 0.00000 0.00000 0.00000 0.00000\n", + "\n", + "ans = 16.208\n" + ] + } + ], + "source": [ + "ea_nr\n", + "newtraph(f,df,0.5,0,12)" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/07_roots-2/.lecture_07.md.swp b/07_roots-2/.lecture_07.md.swp new file mode 100644 index 0000000000000000000000000000000000000000..0998c8680bc72b8c2871c66135b91011671ee56b GIT binary patch literal 16384 zcmeI3e~cs7UBD+?nl?=?jS_@FBJk4e>Fhmwc4q9c*LO0xO`1vtE<_0jq`6yX)*jnK zJTvahc-NP$8{nT5HIbsIkhVgV8dVi)fj0f4O>$`pO+zY#N(%_3Ex)QH5(QdR8d6n7 z4WIXB{4=|`4M;>FkMyx;-n^gR_kG{*`@Z#V8Nc$5Q=+PuHLfc)tvmeCuN_)@`$bn@ zrfJbRCy>b(y4r5!dl$;Y`}fES*WMO32fiC5@!IJd!ys5YbyLZ*?KR((Uya)=!>q0~ z!$B0gjgFVB_1tDMh}?!zow?Qh)&bgEW?4#LDS=rD9M(#fX~>l8-g>RL`l{2jFPHaA z2`nYBl)zE~O9?C`u#~`30!s-jCGh`T0!janb`N9xv22{5%DyjJ_;B}CI!ynSG zqTTPv0|5U1KaT%jev7942wsG5!r#H)!Z+Zn@HG4h{4snSJ^&BF8{rc8+7+7iEc_w- z0SsXWwxI{N!`t9SI00`34ZeR^(_VsqhwsAQ!R;54mZa4+M2)_V74=14j zIXDVez+t!?HaToM3uoX~sKR0R9yLhH*zI|t6=Dcy=Gbk zd)2h^+J+EyaWZ#HAu0rUam#8hTFsfqqWor%*Vl5{ZPCo1(2j@=Ka3LZZYS}=Ks23R zbI?=w8cl@F#MyRfeZNe?ckXgyk#t=VJH9(xUbMqVz3{^z>GsZvzH?6gTt+y4I0%x0 z*eJ*(FX)J{-BvmKkry<*e$Socxuqb&Z8s8u%qseB-%9sn;s&rhCkgYbX<3BjC3?k;ixu*lNj`8?AQj9r1SOY;O9KokrX1 zx$(>ocN%7;Vd~_{%ej%QOqwg{wW@8qMWd|JnWbrCvO{*cQq}qGg6v7ti5i();@n)$ zLRHSIn%MQKT__m^lX;!3tL}%bM(j47AjzG~ZR8RDY-O*kz{nR&w^%9Is#b&g=)z=* z4>V|dzd@V2X|CJIxSYv&t}|r@GS$<9Eyps~wiKYpJWBagmDpS-oxGB6R7@$&G&u^x zM7}V~mR{wzUeu!9ArlL(Q;;>BofS*F;TO#Dc#XvouM8Wz-Z67jm*1n}nRetfciTDj zl+T@MxxK^@>C%^(B<+VMxC2mu*DSS8ShAnJzC?dBtXtExnEiZ|j=$wdMsg-iAhHbI!M4q#S zW){RCR!wW#uXde;7rg=NGVX?hUQ4Kfj#g07?7nF7l(2NKSF%+^A1U#b$y5EnZBXoHRLJS=IOk#^g(v4Ho33>B zS{R0HCTiQ&#YVp$g+tFz=M9lhUmjO;mrdCzHDz&T3&yHkUZ2zDb=|0y40O4ok`JKE zfxE-B_q#c%$t|VM<#hGFLYw7EkU3Z#Cv$#Ym}}N*{H`d;b@VlVt=VYh_0_})y16uQ zl@wm&SCf3Pvu1I-iY?jWt-fiU&}Qn*ieI40?Gxe%`W6P~u=cTN4}w%_sYu(dmn!mg z+;`D=e}L`Yaz&{k{Iyt1>Sa|l5s8OGsf#iy-wcCRtkL7Tp+256UV|f7Q0EyZ99Wr5=+}Nv7V{^bw<0EKRSc9bR&OP%5${2YmIVe zq|mZ;&}lYyTcgpgKN=TLg*=mn=TmC51W0QJDbA`yCmK&W>lQ{Qsw#Fjcs&H%# z#LehW7!+?sg<=Mi8kU0CaYZ0yFa7pZp~ZY()B;13>T_-Rqum2`WIM#lY>jNKCDBPjrjt7NyTPo|-?TC{eO;%8;*iqp zhB069pxI3eQ4T#4m_<3lZMPRTnLOdenU+oV2yPUGQNvM-JoPu0&JPPK6Mr*oD4%5L zcBN`&{$|+la2e4FZK1TY73Z){5!TYnuW^&t(&9WTb4hATdJJz5@3&BwU_)n^>=a6w zaem#pVVYH=w%CbE$xMCk|B5Wkh34hPC~Y1A|Ig<_)**}&8s<~l&S72@-C;jh9OjEA zrY}nyCk;=`E&S8?-Gx=g`fGc*+BobDlnqfjyy9-|7Kdy$#PFN&PPiR@7A}E* zAVz-@?uSiShpXUO_OBm?7BpcC9FY6k5?EluQMeXlO#ahwC0qvIC9eKgcnBVZ9f;r> z_(}Lr;_L6jpTX1cMR*Fn0FT0_;bC|XVz>!@0-h)C{w(|++zwa4UlMP>A0B{v;Rw7y zEd5RRNB9`r4{f*!ejHvPj{ZD+5E^g<{)L$NpW$2ZPjDU{hjE``{k9 z8}5WV;8M64zD1n;JUj^ZgN&(ra0CQg0f*soxC}0WA6(3t0lWy`fUm>X;92-GJO&?u z`{2EBFWdt%-hL;Xh8ti5WW0Sn*r3BIgJQM$0apYLg_#%nrm&`&+Jg8+);<%Mw`Is(lJV}?kH{|L zS#m>0z}{LrBKRL6(d@|$>3IJmYnf~joLVa zVWW7em{GB$hMI<6;{+!?PSHpZlXll9&6?%;W@|OwD48{>%?IqkeqGdLgX&|_>{QRI zGp)K#BVe1B>cBK|BW~k(&Mpv;a5Sef&Cl-vPMfU7)21CK3OG7DE}WrI-#_E6uWRZ_gC4v4*egh3en9pRAKCN4kpA z8fBK5dQ59tPih$_a^{S&M|9JULSGGF^5Qj3om3gE=XTtnl{+exi;0#e1B7))3rDdy zMZ!%sWjwT<9-C6NUE^uXMv1B)puDeI_ouDBQie|AJK_t5oM3Zg=oHp*Y&KIkYTW_y!9`Ct9 zC+W@zUzvlgxPG@AF%YwGvpTPk`*xw?L6M4Rd?2WVuz%*6d7mlaze@EB2al<7+P{H? zjr??&gwB)G7j@psKx8zf+@@))4mWGlKx(*Y8ZxlDxo#13q1N-BPR%=0KXpN`a;Gw< zB_^(@I`IjMuIN=NI%tE_e0wxKv+=6T_yx*T^isvHNgqfbDrL^PvH>bz_Zs`PN6S-( zG`sOq+B$Yf>S12o2XQ>*g%kh(p7mn5;`Jyg+Z`LdiBDNNp4YNxkSgObSOFmcuf#QOc{}qEf}!P4Ef=pLsFzcf}llXoWzq&n7m_; zq%8-5Gk0m`4Ba~iom-$&ySM*GhIY^Qj*@MtgPSqMBkWwMo-9P03x9qnggubjJwELkD_dlzIl+GD0!FvZs40}LtLV`!qjhxe&T@Y+ zhzjL~t%v=%w;GL1V=J|_oz+pUMs;J=R=gT%wO?CV?sSX)kc?fJ027!auqc-Mz1Tgx zeCad#_`;((J5QMa6JP>NfC(@GCcp%k02BB>6R67t@d6EgFm3YUbX=Glr&Au7025#W zOn?b60Vco%m;e)C0!)AjFoE}wfJ}w>72o>5{0PP4`~U3!|1al-_yKwVeF`l>f1eZL z59l>?2yH^=p|jANMInBNeu8$Po6slFB6J4&2l4-ceuI97UP9kPPob|Nm;XM*7ZYFt zOn?b60Vco%m;e)C0`Cn0YfMcEP3*c*Iu%6On3SZAgj7Y(pHo&VOE(fKp-@;&oT@vFpl7 zS5*a}s*<=H$8m4YaXD1h8r!)lYs!_ZYFegDmrYH`*&dE2tG;V8tJ13ogDh9*s88J# zLo{s0zGGQuMpxe6f-6gFCRrhCg%~-^dy0HI>RbuD-*l0{J|}mNho-5Ke$AW9URd_S zZnx`}Ax$%v0G_3^vQR)>QB$duh6kRlv-;OO6ujcbqA5t`x!IRPpmkW=GO3fQv21d; zM+R|Q4K_wXsj6JF|El=dy#+b#z%Os6Q5p){ZdBC*@Jozn8ihJOsj}H#O8_&s{zuii z$qCOBHop)p3G`XfeVjTs`ubNmcwQMzK(J#{stt`Zjd2H<-KL)=0E%ziJXA7teSe$G zPOW!bj?kAFPAv9DyN$vl+<{9v^In_J>yzYW(^b literal 0 HcmV?d00001 diff --git a/07_roots-2/07_roots_and_optimization-2.ipynb b/07_roots-2/07_roots_and_optimization-2.ipynb new file mode 100644 index 0000000..4a2ba51 --- /dev/null +++ b/07_roots-2/07_roots_and_optimization-2.ipynb @@ -0,0 +1,1525 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Roots: Open methods\n", + "## Newton-Raphson" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "First-order approximation for the location of the root (i.e. assume the slope at the given point is constant, what is the solution when f(x)=0)\n", + "\n", + "$f'(x_{i})=\\frac{f(x_{i})-0}{x_{i}-x_{i+1}}$\n", + "\n", + "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n", + "\n", + "Use Newton-Raphson to find solution when $e^{-x}=x$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "Find x when $e^{-x}=x$" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.50000\n", + "error_approx = 1\n" + ] + } + ], + "source": [ + "f= @(x) exp(-x)-x;\n", + "df= @(x) -exp(-x)-1;\n", + "\n", + "x_i= 0;\n", + "x_r = x_i-f(x_i)/df(x_i)\n", + "error_approx = abs((x_r-x_i)/x_r)\n", + "x_i=x_r;\n" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.56714\n", + "error_approx = 0.0014673\n" + ] + } + ], + "source": [ + "x_r = x_i-f(x_i)/df(x_i)\n", + "error_approx = abs((x_r-x_i)/x_r)\n", + "x_i=x_r;" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.56714\n", + "error_approx = 2.2106e-07\n" + ] + } + ], + "source": [ + "x_r = x_i-f(x_i)/df(x_i)\n", + "error_approx = abs((x_r-x_i)/x_r)\n", + "x_i=x_r;" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.56714\n", + "error_approx = 5.0897e-15\n" + ] + } + ], + "source": [ + "x_r = x_i-f(x_i)/df(x_i)\n", + "error_approx = abs((x_r-x_i)/x_r)\n", + "x_i=x_r;" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Create function `newtraph.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,varargin)\n", + "% newtraph: Newton-Raphson root location zeroes\n", + "% [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,p1,p2,...):\n", + "% uses Newton-Raphson method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% dfunc = name of derivative of function\n", + "% xr = initial guess\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by function\n", + "% output:\n", + "% root = real root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "if nargin<4 || isempty(es),es=0.0001;end\n", + "if nargin<5 || isempty(maxit),maxit=50;end\n", + "iter = 0;\n", + "while (1)\n", + " xrold = xr;\n", + " xr = xr - func(xr)/dfunc(xr);\n", + " iter = iter + 1;\n", + " if xr ~= 0 \n", + " ea = abs((xr - xrold)/xr) * 100;\n", + " end\n", + " if ea <= es || iter >= maxit, break, end\n", + "end\n", + "root = xr;" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "In the freefall example, we created a function f(m) that when f(m)=0, then the mass had been chosen such that at t=4 s, the velocity is 36 m/s. \n", + "\n", + "$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$.\n", + "\n", + "to use the Newton-Raphson method, we need the derivative $\\frac{df}{dm}$\n", + "\n", + "$\\frac{df}{dm}=\\frac{1}{2}\\sqrt{\\frac{g}{mc_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-\n", + "\\frac{g}{2m}\\mathrm{sech}^{2}(\\sqrt{\\frac{gc_{d}}{m}}t)$" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "setdefaults\n", + "g=9.81; % acceleration due to gravity\n", + "m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n", + "c_d=0.25; % drag coefficient\n", + "t=4; % at time = 4 seconds\n", + "v=36; % speed must be 36 m/s\n", + "f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n", + "df_m = @(m) 1/2*sqrt(g./m/c_d).*tanh(sqrt(g*c_d./m)*t)-g/2./m*sech(sqrt(g*c_d./m)*t).^2;" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "root = 142.74\n", + "ea = 1.8806e-04\n", + "iter = 50\n" + ] + } + ], + "source": [ + "[root,ea,iter]=newtraph(f_m,df_m,50,0.0001)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Secant Methods" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Not always able to evaluate the derivative. Approximation of derivative:\n", + "\n", + "$f'(x_{i})=\\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}$\n", + "\n", + "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n", + "\n", + "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{\\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}}=\n", + " x_{i}-\\frac{f(x_{i})(x_{i-1}-x_{i})}{f(x_{i-1})-f(x_{i})}$\n", + " \n", + "What values should $x_{i}$ and $x_{i-1}$ take?\n", + "\n", + "To reduce arbitrary selection of variables, use the" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Modified Secant method" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Change the x evaluations to a perturbation $\\delta$. \n", + "\n", + "$x_{i+1}=x_{i}-\\frac{f(x_{i})(\\delta x_{i})}{f(x_{i}+\\delta x_{i})-f(x_{i})}$" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "root = 142.74\n", + "ea = 3.0615e-07\n", + "iter = 7\n" + ] + } + ], + "source": [ + "[root,ea,iter]=mod_secant(f_m,1,50,0.00001)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 1.1185e+04\r\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t25000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tprinciple amount left ($)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\ttime (years)\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "car_payments(400,30000,0.05,5,1)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Amt_numerical = 5467.0\n", + "ans = 3.9755e-04\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t400000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t500000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t600000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t700000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t25\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tprinciple amount left ($)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\ttime (years)\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "Amt_numerical=mod_secant(@(A) car_payments(A,700000,0.0875,30,0),1e-6,50,0.001)\n", + "car_payments(Amt_numerical,700000,0.0875,30,1)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 1.9681e+06\r\n" + ] + } + ], + "source": [ + "Amt_numerical*12*30" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Amortization\n", + "\n", + "Amortization calculation makes the same calculation for the monthly payment amount, A, paying off the principle amount, P, over n pay periods with monthly interest rate, r. " + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Amt = 566.14\r\n" + ] + } + ], + "source": [ + "% Amortization calculation\n", + "A = @(P,r,n) P*(r*(1+r)^n)./((1+r)^n-1);\n", + "Amt=A(30000,0.05/12,5*12)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Matlab's function\n", + "\n", + "Matlab and Octave combine bracketing and open methods in the `fzero` function. " + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'fzero' is a function from the file /usr/share/octave/4.0.0/m/optimization/fzero.m\n", + "\n", + " -- Function File: fzero (FUN, X0)\n", + " -- Function File: fzero (FUN, X0, OPTIONS)\n", + " -- Function File: [X, FVAL, INFO, OUTPUT] = fzero (...)\n", + " Find a zero of a univariate function.\n", + "\n", + " FUN is a function handle, inline function, or string containing the\n", + " name of the function to evaluate.\n", + "\n", + " X0 should be a two-element vector specifying two points which\n", + " bracket a zero. In other words, there must be a change in sign of\n", + " the function between X0(1) and X0(2). More mathematically, the\n", + " following must hold\n", + "\n", + " sign (FUN(X0(1))) * sign (FUN(X0(2))) <= 0\n", + "\n", + " If X0 is a single scalar then several nearby and distant values are\n", + " probed in an attempt to obtain a valid bracketing. If this is not\n", + " successful, the function fails.\n", + "\n", + " OPTIONS is a structure specifying additional options. Currently,\n", + " 'fzero' recognizes these options: \"FunValCheck\", \"OutputFcn\",\n", + " \"TolX\", \"MaxIter\", \"MaxFunEvals\". For a description of these\n", + " options, see *note optimset: XREFoptimset.\n", + "\n", + " On exit, the function returns X, the approximate zero point and\n", + " FVAL, the function value thereof.\n", + "\n", + " INFO is an exit flag that can have these values:\n", + "\n", + " * 1 The algorithm converged to a solution.\n", + "\n", + " * 0 Maximum number of iterations or function evaluations has\n", + " been reached.\n", + "\n", + " * -1 The algorithm has been terminated from user output\n", + " function.\n", + "\n", + " * -5 The algorithm may have converged to a singular point.\n", + "\n", + " OUTPUT is a structure containing runtime information about the\n", + " 'fzero' algorithm. Fields in the structure are:\n", + "\n", + " * iterations Number of iterations through loop.\n", + "\n", + " * nfev Number of function evaluations.\n", + "\n", + " * bracketx A two-element vector with the final bracketing of the\n", + " zero along the x-axis.\n", + "\n", + " * brackety A two-element vector with the final bracketing of the\n", + " zero along the y-axis.\n", + "\n", + " See also: optimset, fsolve.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help fzero" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 563.79\r\n" + ] + } + ], + "source": [ + "fzero(@(A) car_payments(A,30000,0.05,5,0),500)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Comparison of Solvers\n", + "\n", + "It's helpful to compare to the convergence of different routines to see how quickly you find a solution. \n", + "\n", + "Comparing the freefall example\n" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "warning: axis: omitting non-positive data in log plot\n", + "warning: called from\n", + " __line__ at line 120 column 16\n", + " line at line 56 column 8\n", + " __plt__>__plt2vv__ at line 500 column 10\n", + " __plt__>__plt2__ at line 246 column 14\n", + " __plt__ at line 133 column 15\n", + " semilogy at line 60 column 10\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t250\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t350\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t400\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "iterations = linspace(1,400,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f_m,df_m,300,0,iterations(i));\n", + " [root_ms,ea_ms(i),iter_ms]=mod_secant(f_m,1e-6,300,0,iterations(i));\n", + " [root_fp,ea_fp(i),iter_fp]=falsepos(f_m,1,300,0,iterations(i));\n", + " [root_bs,ea_bs(i),iter_bs]=bisect(f_m,1,300,0,iterations(i));\n", + "end\n", + "\n", + "setdefaults\n", + "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n", + "legend('newton-raphson','mod-secant','false point','bisection')" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ea_nr =\n", + "\n", + " Columns 1 through 8:\n", + "\n", + " 6.36591 0.06436 0.00052 0.00000 0.00000 0.00000 0.00000 0.00000\n", + "\n", + " Columns 9 through 16:\n", + "\n", + " 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n", + "\n", + " Columns 17 through 20:\n", + "\n", + " 0.00000 0.00000 0.00000 0.00000\n", + "\n" + ] + } + ], + "source": [ + "ea_nr" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "warning: axis: omitting non-positive data in log plot\n", + "warning: called from\n", + " __line__ at line 120 column 16\n", + " line at line 56 column 8\n", + " __plt__>__plt2vv__ at line 500 column 10\n", + " __plt__>__plt2__ at line 246 column 14\n", + " __plt__ at line 133 column 15\n", + " semilogy at line 60 column 10\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "f= @(x) x^10-1;\n", + "df=@(x) 10*x^9;\n", + "iterations = linspace(1,50,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i));\n", + " [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i));\n", + " [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i));\n", + " [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i));\n", + "end\n", + " \n", + "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n", + "legend('newton-raphson','mod-secant','false point','bisection')" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ea_nr =\n", + "\n", + " Columns 1 through 7:\n", + "\n", + " 99.03195 11.11111 11.11111 11.11111 11.11111 11.11111 11.11111\n", + "\n", + " Columns 8 through 14:\n", + "\n", + " 11.11111 11.11111 11.11111 11.11109 11.11052 11.10624 10.99684\n", + "\n", + " Columns 15 through 20:\n", + "\n", + " 8.76956 2.12993 0.00000 0.00000 0.00000 0.00000\n", + "\n", + "ans = 16.208\n" + ] + } + ], + "source": [ + "ea_nr\n", + "newtraph(f,df,0.5,0,12)" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/07_roots-2/bisect.m b/07_roots-2/bisect.m new file mode 100644 index 0000000..9f696a0 --- /dev/null +++ b/07_roots-2/bisect.m @@ -0,0 +1,37 @@ +function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4||isempty(es), es=0.0001;end +if nargin<5||isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + xr = (xl + xu)/2; + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es || iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); diff --git a/07_roots-2/car_payments.m b/07_roots-2/car_payments.m new file mode 100644 index 0000000..9d5b2a5 --- /dev/null +++ b/07_roots-2/car_payments.m @@ -0,0 +1,17 @@ +function amount_left = car_payments(monthly_payment,price,apr,no_of_years,plot_bool) + interest_per_month = apr/12; + number_of_months = no_of_years*12; + principle=price; + P_vector=zeros(1,number_of_months); + for i = 1:number_of_months + principle=principle-monthly_payment; + principle=(1+interest_per_month)*principle; + P_vector(i)=principle; + end + amount_left=principle; + if plot_bool + plot([1:number_of_months]/12, P_vector) + xlabel('time (years)') + ylabel('principle amount left ($)') + end +end diff --git a/07_roots-2/falsepos.m b/07_roots-2/falsepos.m new file mode 100644 index 0000000..d5575d5 --- /dev/null +++ b/07_roots-2/falsepos.m @@ -0,0 +1,39 @@ +function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4||isempty(es), es=0.0001;end +if nargin<5||isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + xr = (xl + xu)/2; + % xr = (xl + xu)/2; % bisect method + xr=xu - (func(xu)*(xl-xu))/(func(xl)-func(xu)); % false position method + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es || iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); diff --git a/07_roots-2/fzerosimp.m b/07_roots-2/fzerosimp.m new file mode 100644 index 0000000..05c7a9b --- /dev/null +++ b/07_roots-2/fzerosimp.m @@ -0,0 +1,41 @@ +function b = fzerosimp(xl,xu) +a = xl; b = xu; fa = f(a); fb = f(b); +c = a; fc = fa; d = b - c; e = d; +while (1) + if fb == 0, break, end + if sign(fa) == sign(fb) %If needed, rearrange points + a = c; fa = fc; d = b - c; e = d; + end + if abs(fa) < abs(fb) + c = b; b = a; a = c; + fc = fb; fb = fa; fa = fc; + end + m = 0.5*(a - b); %Termination test and possible exit + tol = 2 * eps * max(abs(b), 1); + if abs(m) <= tol | fb == 0. + break + end + %Choose open methods or bisection + if abs(e) >= tol & abs(fc) > abs(fb) + s = fb/fc; + if a == c %Secant method + p = 2*m*s; + q = 1 - s; + else %Inverse quadratic interpolation + q = fc/fa; r = fb/fa; + p = s * (2*m*q * (q - r) - (b - c)*(r - 1)); + q = (q - 1)*(r - 1)*(s - 1); + end + if p > 0, q = -q; else p = -p; end; + if 2*p < 3*m*q - abs(tol*q) & p < abs(0.5*e*q) + e = d; d = p/q; + else + d = m; e = m; + end + else %Bisection + d = m; e = m; + end + c = b; fc = fb; + if abs(d) > tol, b=b+d; else b=b-sign(b-a)*tol; end + fb = f(b); +end \ No newline at end of file diff --git a/07_roots-2/incsearch.m b/07_roots-2/incsearch.m new file mode 100644 index 0000000..bd82554 --- /dev/null +++ b/07_roots-2/incsearch.m @@ -0,0 +1,37 @@ +function xb = incsearch(func,xmin,xmax,ns) +% incsearch: incremental search root locator +% xb = incsearch(func,xmin,xmax,ns): +% finds brackets of x that contain sign changes +% of a function on an interval +% input: +% func = name of function +% xmin, xmax = endpoints of interval +% ns = number of subintervals (default = 50) +% output: +% xb(k,1) is the lower bound of the kth sign change +% xb(k,2) is the upper bound of the kth sign change +% If no brackets found, xb = []. +if nargin < 3, error('at least 3 arguments required'), end +if nargin < 4, ns = 50; end %if ns blank set to 50 +% Incremental search +x = linspace(xmin,xmax,ns); +f = func(x); +nb = 0; xb = []; %xb is null unless sign change detected +%for k = 1:length(x)-1 +% if sign(f(k)) ~= sign(f(k+1)) %check for sign change +% nb = nb + 1; +% xb(nb,1) = x(k); +% xb(nb,2) = x(k+1); +% end +%end +sign_change = diff(sign(f)); +[~,i_change] = find(sign_change~=0); +nb=length(i_change); +xb=[x(i_change)',x(i_change+1)']; + +if isempty(xb) %display that no brackets were found + fprintf('no brackets found\n') + fprintf('check interval or increase ns\n') +else + fprintf('number of brackets: %i\n',nb) %display number of brackets +end diff --git a/07_roots-2/mod_secant.m b/07_roots-2/mod_secant.m new file mode 100644 index 0000000..a3caa10 --- /dev/null +++ b/07_roots-2/mod_secant.m @@ -0,0 +1,31 @@ +function [root,ea,iter]=mod_secant(func,dx,xr,es,maxit,varargin) +% newtraph: Modified secant root location zeroes +% [root,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...): +% uses modified secant method to find the root of func +% input: +% func = name of function +% dx = perturbation fraction +% xr = initial guess +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by function +% output: +% root = real root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +if nargin<4 || isempty(es),es=0.0001;end +if nargin<5 || isempty(maxit),maxit=50;end +iter = 0; +while (1) + xrold = xr; + dfunc=(func(xr+dx)-func(xr))./dx; + xr = xr - func(xr)/dfunc; + iter = iter + 1; + if xr ~= 0 + ea = abs((xr - xrold)/xr) * 100; + end + if ea <= es || iter >= maxit, break, end +end +root = xr; + diff --git a/07_roots-2/newtraph.m b/07_roots-2/newtraph.m new file mode 100644 index 0000000..94ca667 --- /dev/null +++ b/07_roots-2/newtraph.m @@ -0,0 +1,29 @@ +function [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,varargin) +% newtraph: Newton-Raphson root location zeroes +% [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,p1,p2,...): +% uses Newton-Raphson method to find the root of func +% input: +% func = name of function +% dfunc = name of derivative of function +% xr = initial guess +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by function +% output: +% root = real root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +if nargin<4 || isempty(es),es=0.0001;end +if nargin<5 || isempty(maxit),maxit=50;end +iter = 0; +while (1) + xrold = xr; + xr = xr - func(xr)/dfunc(xr); + iter = iter + 1; + if xr ~= 0 + ea = abs((xr - xrold)/xr) * 100; + end + if ea <= es || iter >= maxit, break, end +end +root = xr; diff --git a/07_roots-2/octave-workspace b/07_roots-2/octave-workspace new file mode 100644 index 0000000000000000000000000000000000000000..70eb62be42b668ffeb3338f31f05f120540b94dd GIT binary patch literal 1948 zcmcJPTWAzl7{|x9HkGRJ@(`pCIVi+Ua983b^dW9ODhkoWxG{H=%$}V&GrM$mb~8H( z+v1CO0Z}U0(ms^7pyGp-3O&s7XLir| zuD}0$)8j;PO2u2_UDk(|Wyzq7t7TO`!o^Yf{7ki4U4Xwozas2az1am5TD<5G5HUWR0W*Ub(p$o2G} zYCS%@$-s$n#Xw|b;O5n0+=6x)6oan{#mRDV1^tgznnheyN+K4gU!s+YbUnL(56tq4 zFjW1^VZ_Es+KCmKMXOjZHi#Bco(Llm^F?C4DBu>>+$<6;B9cK?#AXj$XkVK&qf11| zeeJKkFo<&AWStln8{P(lRPynji>vX_)6~VE*8rVHHK%&RS_TAD$(rN8&TagB^O;x& z!Bk;b3~gB|g~dvZ;orT7g9XoCv{ROM*JL?bYfWB~q(W`ZOOZ`mbQo0Rf7tU9V(NX8 zB_dOJ&s9ZlSYK#WUauLwow-Pb+PZcII<>obcWFWj5+|@SJAgDgW?DMo_qxV?&#!dC z!TRL%g9AHYd&5|CYuirvgtqIt;LoNX^B1ml!Lh`pCnt|~!>+b_?zi2$U>TX$_rUFK zuu&xvO%|d%jSMor483K_w{ueG~Tmm?(cruaF*O> zyEMEd_xqDH+$HxtBMtuy_ity64jJZfIb-z5Fc0AvT^#1J-!b|)%;%9~baI)~R@dm| zGOyFF(amLU3zgANF+Wck9TjuDs*IkBdB!}WtH&IUct&53`Tpw}oqgu8(>Hqi%=?^g zboZIVlC04`%N&BNc@MIz!;P$YA96q=78Z-OU&cS2wd+pYm2dfOzInRoyHc@pxXzva Tq3^#rzo_X^2cOrD_w4z9m{@J! literal 0 HcmV?d00001