From 908bbf1e62e5b4fad1685c9e0921313fb305b061 Mon Sep 17 00:00:00 2001 From: "Ryan C. Cooper" Date: Thu, 2 Feb 2017 09:30:03 -0500 Subject: [PATCH] updated lecture 06 --- lecture_06/bisect.m | 37 + lecture_06/falsepos.m | 39 + lecture_06/incsearch.m | 37 + lecture_06/lecture_06.ipynb | 834 ++++++++++++++++++ lecture_06/lecture_06.md | 373 ++++++++ .../lecture_06_files/lecture_06_18_0.svg | 144 +++ .../lecture_06_files/lecture_06_5_0.svg | 145 +++ 7 files changed, 1609 insertions(+) create mode 100644 lecture_06/bisect.m create mode 100644 lecture_06/falsepos.m create mode 100644 lecture_06/incsearch.m create mode 100644 lecture_06/lecture_06.ipynb create mode 100644 lecture_06/lecture_06.md create mode 100644 lecture_06/lecture_06_files/lecture_06_18_0.svg create mode 100644 lecture_06/lecture_06_files/lecture_06_5_0.svg diff --git a/lecture_06/bisect.m b/lecture_06/bisect.m new file mode 100644 index 0000000..c09ffbf --- /dev/null +++ b/lecture_06/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/lecture_06/falsepos.m b/lecture_06/falsepos.m new file mode 100644 index 0000000..0a3477c --- /dev/null +++ b/lecture_06/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/lecture_06/incsearch.m b/lecture_06/incsearch.m new file mode 100644 index 0000000..bd82554 --- /dev/null +++ b/lecture_06/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/lecture_06/lecture_06.ipynb b/lecture_06/lecture_06.ipynb new file mode 100644 index 0000000..f44f84c --- /dev/null +++ b/lecture_06/lecture_06.ipynb @@ -0,0 +1,834 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "my_caller.m\n", + "```matlab\n", + "function [vx,vy] = my_caller(max_time)\n", + " N=100;\n", + " t=linspace(0,max_time,N);\n", + " [x,y]=my_function(max_time);\n", + " vx=diff(x)./diff(t);\n", + " vy=diff(y)./diff(t);\n", + "end\n", + "```\n", + "\n", + "my_function.m\n", + "```matlab\n", + "function [x,y] = my_function(max_time)\n", + " N=100;\n", + " t=linspace(0,max_time,N);\n", + " x=t.^2;\n", + " y=2*t;\n", + "end\n", + "```\n", + "\n", + "In order to use `my_caller.m` where does `my_function.m` need to be saved?\n", + "![responses](q1.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "What cool personal projects are you working on?\n", + "While we delve deeper into Matlab functions, could we review some of the basic logic\n", + "operators it uses and command codes. \n", + "\n", + "I still dont know when these forms are technically due. \n", + " \n", + " -by the following lecture\n", + "\n", + "I'm having trouble interfacing Atom with GitHub. Is there a simple tutorial for this?\n", + " \n", + " -Mac? Seems there could be a bug that folks are working on\n", + "\n", + "What are the bear necessities of life? \n", + "please go over how to \"submit\" the homeworks because it is still confusing\n", + "\n", + "Do you prefer Matlab or Octave?\n", + " \n", + " -octave is my preference, but Matlab has some benefits\n", + "\n", + "Would you consider a country to be open-source?\n", + " \n", + " -??\n", + "\n", + "Is there a way to download matlab for free?\n", + " \n", + " -not legally\n", + "\n", + "how do you add files to current folder in matlab?\n", + " \n", + " -you can do this either through a file browser or cli\n", + "\n", + "How should Homework 2 be submitted? By simply putting the function into the homework_1\n", + "repository?\n", + " \n", + " -yes\n", + " \n", + "How can we tell that these forms are being received?\n", + " \n", + " -when you hit submit, the form says \"form received\"\n", + " \n", + "can you save scripted outputs from matlab/octave as an excel file?\n", + " \n", + " -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`\n", + " \n", + "\n", + "Also, can you update your notes to show what happens when these things are run, as you do\n", + "in class?\"\n", + " \n", + " -I always update the lecture notes after class so they should display what we did in class\n", + " \n", + "I have a little difficulty following along in class on my laptop when you have programs\n", + "pre-written. Maybe if you posted those codes on Github so I could copy them when you\n", + "switch to different desktops I would be able to follow along better.\n", + "\n", + "Kirk or Picard?\n", + " \n", + " -Kirk\n", + " \n", + "Who is our TA?\n", + " \n", + " -Peiyu Zhang peiyu.zhang@uconn.edu\n", + "\n", + "Can we download libraries of data like thermodynamic tables into matlab?\n", + " \n", + "-YES! [Matlab Steam Tables](http://bit.ly/2kZygu8)\n", + "\n", + "Will we use the Simulink addition to Matlab? I found it interesting and useful for\n", + "evaluating ODEs in Linear systems.\n", + " \n", + " -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\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Roots and Optimization\n", + "## Bracketing ch. 5" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When you are given a function, numerical or analytical, it's not always possible to solve directly for a given variable. \n", + "\n", + "Even for the freefall example we first explored, \n", + "\n", + "$v(t)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)$\n", + "\n", + "There is no way to solve for m in terms of the other variables. \n", + "\n", + "Instead, we can 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": 4, + "metadata": { + "collapsed": false + }, + "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": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 0.045626\r\n" + ] + } + ], + "source": [ + "f_m(145)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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": {}, + "source": [ + "## Incremental method (Brute force)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. \n", + "\n", + "```matlab\n", + "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\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/me3255_S2017/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": 24, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "no brackets found\n", + "check interval or increase ns\n", + "ans = [](1x0)\n" + ] + } + ], + "source": [ + "incsearch(f_m,50, 200,55)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "## Bisection method" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 0.020577\r\n" + ] + } + ], + "source": [ + "f_m(143.75)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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`\n", + "\n", + "```matlab\n", + "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{:});\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## False position (linear interpolation)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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 + }, + "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", + "xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu));\n", + "\n", + "plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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`\n", + "\n", + "```matlab\n", + "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{:});\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "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/lecture_06/lecture_06.md b/lecture_06/lecture_06.md new file mode 100644 index 0000000..b3d2266 --- /dev/null +++ b/lecture_06/lecture_06.md @@ -0,0 +1,373 @@ + + +```octave +%plot --format svg +``` + +my_caller.m +```matlab +function [vx,vy] = my_caller(max_time) + N=100; + t=linspace(0,max_time,N); + [x,y]=my_function(max_time); + vx=diff(x)./diff(t); + vy=diff(y)./diff(t); +end +``` + +my_function.m +```matlab +function [x,y] = my_function(max_time) + N=100; + t=linspace(0,max_time,N); + x=t.^2; + y=2*t; +end +``` + +In order to use `my_caller.m` where does `my_function.m` need to be saved? +![responses](q1.png) + + +What cool personal projects are you working on? +While we delve deeper into Matlab functions, could we review some of the basic logic +operators it uses and command codes. + +I still dont know when these forms are technically due. + + -by the following lecture + +I'm having trouble interfacing Atom with GitHub. Is there a simple tutorial for this? + + -Mac? Seems there could be a bug that folks are working on + +What are the bear necessities of life? +please go over how to "submit" the homeworks because it is still confusing + +Do you prefer Matlab or Octave? + + -octave is my preference, but Matlab has some benefits + +Would you consider a country to be open-source? + + -?? + +Is there a way to download matlab for free? + + -not legally + +how do you add files to current folder in matlab? + + -you can do this either through a file browser or cli + +How should Homework 2 be submitted? By simply putting the function into the homework_1 +repository? + + -yes + +How can we tell that these forms are being received? + + -when you hit submit, the form says "form received" + +can you save scripted outputs from matlab/octave as an excel file? + + -yes, easy way is open a file with a .csv extension then fprintf and separate everything with commas, harder way is to use the `xlswrite` + + +Also, can you update your notes to show what happens when these things are run, as you do +in class?" + + -I always update the lecture notes after class so they should display what we did in class + +I have a little difficulty following along in class on my laptop when you have programs +pre-written. Maybe if you posted those codes on Github so I could copy them when you +switch to different desktops I would be able to follow along better. + +Kirk or Picard? + + -Kirk + +Who is our TA? + + -Peiyu Zhang peiyu.zhang@uconn.edu + +Can we download libraries of data like thermodynamic tables into matlab? + +-YES! [Matlab Steam Tables](http://bit.ly/2kZygu8) + +Will we use the Simulink addition to Matlab? I found it interesting and useful for +evaluating ODEs in Linear systems. + + -not in this class, everything in simulink has a matlab script/function that can be substituted, but many times its hidden by the gui. Here we want to look directly at our solvers + + +# Roots and Optimization +## Bracketing ch. 5 + +When you are given a function, numerical or analytical, it's not always possible to solve directly for a given variable. + +Even for the freefall example we first explored, + +$v(t)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)$ + +There is no way to solve for m in terms of the other variables. + +Instead, we can solve the problem by creating a new function f(m) where + +$f(m)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)-v(t)$. + +When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity) + + +```octave +setdefaults +g=9.81; % acceleration due to gravity +m=linspace(50, 200,100); % possible values for mass 50 to 200 kg +c_d=0.25; % drag coefficient +t=4; % at time = 4 seconds +v=36; % speed must be 36 m/s +f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m + +plot(m,f_m(m),m,zeros(length(m),1)) +axis([45 200 -5 1]) +``` + + +![svg](lecture_06_files/lecture_06_5_0.svg) + + + +```octave +f_m(145) +``` + + ans = 0.045626 + + +Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0 + +Better methods are the +1. Bracketing methods +2. Open methods + +Both need an initial guess. + + +## Incremental method (Brute force) + +You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. + +```matlab +function xb = incsearch(func,xmin,xmax,ns) +% incsearch: incremental search root locator +% xb = incsearch(func,xmin,xmax,ns): +% finds brackets of x that contain sign changes +% of a function on an interval +% input: +% func = name of function +% xmin, xmax = endpoints of interval +% ns = number of subintervals (default = 50) +% output: +% xb(k,1) is the lower bound of the kth sign change +% xb(k,2) is the upper bound of the kth sign change +% If no brackets found, xb = []. +if nargin < 3, error('at least 3 arguments required'), end +if nargin < 4, ns = 50; end %if ns blank set to 50 +% Incremental search +x = linspace(xmin,xmax,ns); +f = func(x); +nb = 0; xb = []; %xb is null unless sign change detected +for k = 1:length(x)-1 + if sign(f(k)) ~= sign(f(k+1)) %check for sign change + nb = nb + 1; + xb(nb,1) = x(k); + xb(nb,2) = x(k+1); + end +end +if isempty(xb) %display that no brackets were found + fprintf('no brackets found\n') + fprintf('check interval or increase ns\n') +else + fprintf('number of brackets: %i\n',nb) %display number of brackets +end +``` + + +```octave +help incsearch +``` + + 'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/me3255_S2017/lecture_06/incsearch.m + + incsearch: incremental search root locator + xb = incsearch(func,xmin,xmax,ns): + finds brackets of x that contain sign changes + of a function on an interval + input: + func = name of function + xmin, xmax = endpoints of interval + ns = number of subintervals (default = 50) + output: + xb(k,1) is the lower bound of the kth sign change + xb(k,2) is the upper bound of the kth sign change + If no brackets found, xb = []. + + + Additional help for built-in functions and operators is + available in the online version of the manual. Use the command + 'doc ' to search the manual index. + + Help and information about Octave is also available on the WWW + at http://www.octave.org and via the help@octave.org + mailing list. + + + +```octave +incsearch(f_m,50, 200,55) +``` + + no brackets found + check interval or increase ns + ans = [](1x0) + + +## Bisection method + +Divide interval in half until error is reduced to some level + +in previous example of freefall, choose x_l=50, x_u=200 + +x_r = (50+200)/2 = 125 + +f_m(125) = -0.408 + +x_r= (125+200)/2 = 162.5 + +f_m(162.5) = 0.3594 + +x_r = (125+162.5)/2=143.75 + +f_m(143.75)= 0.0206 + + +```octave +f_m(143.75) +``` + + ans = 0.020577 + + +Much better root locator, with 4 iterations, our function is already close to zero + +Automate this with a function: +`bisect.m` + +```matlab +function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + xr = (xl + xu)/2; + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); +``` + +## False position (linear interpolation) + +Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner + +$ x_{r} = x_{u} - \frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$ + + +```octave +xl=50; xu=200; +xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); + +plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0) +``` + + +![svg](lecture_06_files/lecture_06_18_0.svg) + + +Much better root locator, with 4 iterations, our function is already close to zero + +Automate this with a function: +`falsepos.m` + +```matlab +function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin) +% falsepos: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses false position method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + % xr = (xl + xu)/2; % bisect method + xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); +``` + + +```octave + +``` diff --git a/lecture_06/lecture_06_files/lecture_06_18_0.svg b/lecture_06/lecture_06_files/lecture_06_18_0.svg new file mode 100644 index 0000000..9652b37 --- /dev/null +++ b/lecture_06/lecture_06_files/lecture_06_18_0.svg @@ -0,0 +1,144 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -5 + + + + + -4 + + + + + -3 + + + + + -2 + + + + + -1 + + + + + 0 + + + + + 1 + + + + + 0 + + + + + 50 + + + + + 100 + + + + + 150 + + + + + 200 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + gnuplot_plot_3a + + + + + + gnuplot_plot_4a + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_06/lecture_06_files/lecture_06_5_0.svg b/lecture_06/lecture_06_files/lecture_06_5_0.svg new file mode 100644 index 0000000..3ecaecd --- /dev/null +++ b/lecture_06/lecture_06_files/lecture_06_5_0.svg @@ -0,0 +1,145 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -5 + + + + + -4 + + + + + -3 + + + + + -2 + + + + + -1 + + + + + 0 + + + + + 1 + + + + + 60 + + + + + 80 + + + + + 100 + + + + + 120 + + + + + 140 + + + + + 160 + + + + + 180 + + + + + 200 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + + + + + + + + \ No newline at end of file