From 7f41cbb8b7fe46e4be3ea2894e6b09f673e62cbf Mon Sep 17 00:00:00 2001 From: "Ryan C. Cooper" Date: Mon, 18 Sep 2017 13:59:24 -0400 Subject: [PATCH] added root finding lectures --- ...06_roots and optimization-checkpoint.ipynb | 913 ----------------- .../06_roots and optimization.ipynb | 934 ------------------ 06_roots and optimization/bisect.m | 37 - 06_roots and optimization/falsepos.m | 39 - 06_roots and optimization/incsearch.m | 37 - 06_roots and optimization/lecture_06.md | 373 ------- 06_roots and optimization/octave-workspace | Bin 1332 -> 0 bytes 7 files changed, 2333 deletions(-) delete mode 100644 06_roots and optimization/.ipynb_checkpoints/06_roots and optimization-checkpoint.ipynb delete mode 100644 06_roots and optimization/06_roots and optimization.ipynb delete mode 100644 06_roots and optimization/bisect.m delete mode 100644 06_roots and optimization/falsepos.m delete mode 100644 06_roots and optimization/incsearch.m delete mode 100644 06_roots and optimization/lecture_06.md delete mode 100644 06_roots and optimization/octave-workspace diff --git a/06_roots and optimization/.ipynb_checkpoints/06_roots and optimization-checkpoint.ipynb b/06_roots and optimization/.ipynb_checkpoints/06_roots and optimization-checkpoint.ipynb deleted file mode 100644 index bcfbe7b..0000000 --- a/06_roots and optimization/.ipynb_checkpoints/06_roots and optimization-checkpoint.ipynb +++ /dev/null @@ -1,913 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 2, - "metadata": { - "collapsed": true, - "scrolled": true, - "slideshow": { - "slide_type": "skip" - } - }, - "outputs": [], - "source": [ - "%plot --format svg" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "# Roots and Optimization\n", - "## Bracketing ch. 5" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "source": [ - "Given a function, numerical or analytical, it's not always possible 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": 3, - "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": 5, - "metadata": { - "collapsed": false, - "scrolled": true, - "slideshow": { - "slide_type": "subslide" - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "ans = -0.056985\r\n" - ] - } - ], - "source": [ - "f_m(140)" - ] - }, - { - "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": 7, - "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/me3255_S2017/course_git/lecture_06/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": 9, - "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.42 143.94\n", - "\n" - ] - } - ], - "source": [ - "incsearch(f_m,50, 200,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": 10, - "metadata": { - "collapsed": false, - "scrolled": true, - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "ans = 0.020577\r\n" - ] - } - ], - "source": [ - "f_m(143.75)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## Bisect" - ] - }, - { - "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": 3, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "'bisect' is a function from the file /home/ryan/Documents/UConn/ME3255/course_git_F2017/06_roots and optimization/bisect.m\n", - "\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", - "\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 bisect" - ] - }, - { - "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 and optimization/06_roots and optimization.ipynb b/06_roots and optimization/06_roots and optimization.ipynb deleted file mode 100644 index 72dfcbf..0000000 --- a/06_roots and optimization/06_roots and optimization.ipynb +++ /dev/null @@ -1,934 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 5, - "metadata": { - "collapsed": true, - "scrolled": true, - "slideshow": { - "slide_type": "skip" - } - }, - "outputs": [], - "source": [ - "%plot --format svg" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "# Roots and Optimization\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": 6, - "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": 7, - "metadata": { - "collapsed": false, - "scrolled": true, - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "ans = 0.0053580\r\n" - ] - } - ], - "source": [ - "f_m(143)" - ] - }, - { - "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": 7, - "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/me3255_S2017/course_git/lecture_06/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": 9, - "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.42 143.94\n", - "\n" - ] - } - ], - "source": [ - "incsearch(f_m,50, 200,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 = 125\n", - "interval left f(x_l)= -4.6,f(x_r)= -0.4\n", - "interval right f(x_r)= -0.4,f(x_u)= 0.9\n", - "ans = -0.40860\n" - ] - } - ], - "source": [ - "x_l=50; x_u=200;\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": 3, - "metadata": { - "collapsed": false, - "slideshow": { - "slide_type": "subslide" - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "'bisect' is a function from the file /home/ryan/Documents/UConn/ME3255/course_git_F2017/06_roots and optimization/bisect.m\n", - "\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", - "\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 bisect" - ] - }, - { - "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 and optimization/bisect.m b/06_roots and optimization/bisect.m deleted file mode 100644 index c09ffbf..0000000 --- a/06_roots and optimization/bisect.m +++ /dev/null @@ -1,37 +0,0 @@ -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 and optimization/falsepos.m b/06_roots and optimization/falsepos.m deleted file mode 100644 index 0a3477c..0000000 --- a/06_roots and optimization/falsepos.m +++ /dev/null @@ -1,39 +0,0 @@ -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 and optimization/incsearch.m b/06_roots and optimization/incsearch.m deleted file mode 100644 index bd82554..0000000 --- a/06_roots and optimization/incsearch.m +++ /dev/null @@ -1,37 +0,0 @@ -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 and optimization/lecture_06.md b/06_roots and optimization/lecture_06.md deleted file mode 100644 index b3d2266..0000000 --- a/06_roots and optimization/lecture_06.md +++ /dev/null @@ -1,373 +0,0 @@ - - -```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 and optimization/octave-workspace b/06_roots and optimization/octave-workspace deleted file mode 100644 index 0918300ce10589a9faf7b968c7df71a640e71a70..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1332 zcmbu9OK4L;6o#*eMi5QGMG+Sr6lxksiXdIpcu*?VK9lCrq)Bfk_ZdiTHMuEOx=?(8 zsJKvYAu6c2sN$k4>Bj0t6hteys3NGiDC)vxSB*1s4oI=;!V4L==bt%`@4u5F!{?J$ zWKX1DTd!#v_v&KRwB|;!)D2$YwZO^4m*&s4zk4ioLpP>5x^?s?_dR4tPgYu=} z8Gd=fV=nj1inUdQpqqOW-tp>0qaI^>_sK}8+QI6_YCa1&EX=B`lZhf7EE;A$_gohG zrHMF;b~3SFXdNuVCjZsaI#Iw8zts)5o&6SnbaNL-)&EUZ0j));lYdYTDE04|34{b$ z<$kR^wbE=hWl2depb~4n@HEgi{o_e5oC(I3Kb}4Uhgu8%LQfxTlHouL9Acrh`l4*)r;*o_rBFoaw9C68kCDJc` z$q}D1*ewg(8zWA6P-Wdb@ya8IaDliLpehIZULbzs$l>8QDC|KUIULnNjk}=ZJll)J zwTN?_D1yQ}Q^fh+7l|{+IlDRW<~Z*SPTV=py~2rq3Fp^K