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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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?
+
+
+
+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])
+```
+
+
+
+
+
+
+```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)
+```
+
+
+
+
+
+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 0000000..05454da
Binary files /dev/null and b/06_roots-1/octave-workspace differ
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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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 0000000..0998c86
Binary files /dev/null and b/07_roots-2/.lecture_07.md.swp differ
diff --git a/07_roots-2/.newtraph.m.swp b/07_roots-2/.newtraph.m.swp
new file mode 100644
index 0000000..9759dd2
Binary files /dev/null and b/07_roots-2/.newtraph.m.swp differ
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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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 0000000..70eb62b
Binary files /dev/null and b/07_roots-2/octave-workspace differ