From 908bbf1e62e5b4fad1685c9e0921313fb305b061 Mon Sep 17 00:00:00 2001 From: "Ryan C. Cooper" Date: Thu, 2 Feb 2017 09:30:03 -0500 Subject: [PATCH 01/14] 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 From 2ce5ca32f0eeb211bb368079c412aabf8646d33b Mon Sep 17 00:00:00 2001 From: "Ryan C. Cooper" Date: Mon, 6 Feb 2017 23:06:31 -0500 Subject: [PATCH 02/14] updated lecture 07 --- lecture_07/bisect.m | 37 + lecture_07/car_payments.m | 17 + lecture_07/falsepos.m | 39 + lecture_07/fzerosimp.m | 41 + lecture_07/incsearch.m | 37 + lecture_07/lecture_07.ipynb | 1195 +++++++++++++++++ lecture_07/lecture_07.md | 388 ++++++ lecture_07/lecture_07.pdf | Bin 0 -> 68853 bytes .../lecture_07_files/lecture_07_14_1.svg | 146 ++ .../lecture_07_files/lecture_07_22_1.svg | 197 +++ .../lecture_07_files/lecture_07_24_1.svg | 182 +++ lecture_07/mod_secant.m | 31 + lecture_07/newtraph.m | 29 + 13 files changed, 2339 insertions(+) create mode 100644 lecture_07/bisect.m create mode 100644 lecture_07/car_payments.m create mode 100644 lecture_07/falsepos.m create mode 100644 lecture_07/fzerosimp.m create mode 100644 lecture_07/incsearch.m create mode 100644 lecture_07/lecture_07.ipynb create mode 100644 lecture_07/lecture_07.md create mode 100644 lecture_07/lecture_07.pdf create mode 100644 lecture_07/lecture_07_files/lecture_07_14_1.svg create mode 100644 lecture_07/lecture_07_files/lecture_07_22_1.svg create mode 100644 lecture_07/lecture_07_files/lecture_07_24_1.svg create mode 100644 lecture_07/mod_secant.m create mode 100644 lecture_07/newtraph.m diff --git a/lecture_07/bisect.m b/lecture_07/bisect.m new file mode 100644 index 0000000..9f696a0 --- /dev/null +++ b/lecture_07/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/lecture_07/car_payments.m b/lecture_07/car_payments.m new file mode 100644 index 0000000..9d5b2a5 --- /dev/null +++ b/lecture_07/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/lecture_07/falsepos.m b/lecture_07/falsepos.m new file mode 100644 index 0000000..d5575d5 --- /dev/null +++ b/lecture_07/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/lecture_07/fzerosimp.m b/lecture_07/fzerosimp.m new file mode 100644 index 0000000..05c7a9b --- /dev/null +++ b/lecture_07/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/lecture_07/incsearch.m b/lecture_07/incsearch.m new file mode 100644 index 0000000..bd82554 --- /dev/null +++ b/lecture_07/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_07/lecture_07.ipynb b/lecture_07/lecture_07.ipynb new file mode 100644 index 0000000..0c3c2ca --- /dev/null +++ b/lecture_07/lecture_07.ipynb @@ -0,0 +1,1195 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Roots: Open methods\n", + "## Newton-Raphson" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "error: 'c' undefined near line 1 column 1\n", + "error: 'x_r' undefined near line 1 column 21\n", + "error: evaluating argument list element number 1\n", + "error: 'x_r' undefined near line 1 column 5\n" + ] + } + ], + "source": [ + "f= @(x) exp(-x)-x;\n", + "df= @(x) -exp(-x)-1;\n", + "\n", + "x_i= 0;\n", + "c\n", + "error_approx = abs((x_r-x_i)/x_r)\n", + "x_i=x_r;\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.50000\n", + "error_approx = 1\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": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 0.56631\n", + "error_approx = 0.11709\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": 4, + "metadata": { + "collapsed": false + }, + "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": "markdown", + "metadata": {}, + "source": [ + "In the bungee jumper 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": 6, + "metadata": { + "collapsed": true + }, + "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": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 142.74\r\n" + ] + } + ], + "source": [ + "newtraph(f_m,df_m,140,0.00001)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Secant Methods" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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\n", + "\n", + "## Modified Secant method\n", + "\n", + "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": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 142.74\r\n" + ] + } + ], + "source": [ + "mod_secant(f_m,1e-6,50,0.00001)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Amt_numerical = 563.79\n", + "ans = -160.42\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t25000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tprinciple amount left ($)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\ttime (years)\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "Amt_numerical=mod_secant(@(A) car_payments(A,30000,0.05,5),1e-6,50,0.001)\n", + "car_payments(Amt,30000,0.05,5)" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 3.3968e+04\r\n" + ] + } + ], + "source": [ + "Amt*12*5" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "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": 27, + "metadata": { + "collapsed": false + }, + "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": {}, + "source": [ + "## Matlab's function\n", + "\n", + "Matlab and Octave combine bracketing and open methods in the `fzero` function. " + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": { + "collapsed": false + }, + "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": 40, + "metadata": { + "collapsed": false + }, + "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": {}, + "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": 84, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "warning: axis: omitting non-positive data in log plot\n", + "warning: called from\n", + " __line__ at line 120 column 16\n", + " line at line 56 column 8\n", + " __plt__>__plt2vv__ at line 500 column 10\n", + " __plt__>__plt2__ at line 246 column 14\n", + " __plt__ at line 133 column 15\n", + " semilogy at line 60 column 10\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t250\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t350\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t400\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "iterations = linspace(1,400,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f_m,df_m,200,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", + "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": 75, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ea_ms =\n", + "\n", + " Columns 1 through 7:\n", + "\n", + " 43.43883 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n", + "\n", + " Columns 8 through 14:\n", + "\n", + " 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n", + "\n", + " Columns 15 through 20:\n", + "\n", + " 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n", + "\n" + ] + } + ], + "source": [ + "ea_ms" + ] + }, + { + "cell_type": "code", + "execution_count": 101, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "warning: axis: omitting non-positive data in log plot\n", + "warning: called from\n", + " __line__ at line 120 column 16\n", + " line at line 56 column 8\n", + " __plt__>__plt2vv__ at line 500 column 10\n", + " __plt__>__plt2__ at line 246 column 14\n", + " __plt__ at line 133 column 15\n", + " semilogy at line 60 column 10\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n", + "warning: axis: omitting non-positive data in log plot\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "f= @(x) x^10-1;\n", + "df=@(x) 10*x^9;\n", + "iterations = linspace(1,50,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i));\n", + " [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i));\n", + " [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i));\n", + " [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i));\n", + "end\n", + " \n", + "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n", + "legend('newton-raphson','mod-secant','false point','bisection')" + ] + }, + { + "cell_type": "code", + "execution_count": 102, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ea_bs =\n", + "\n", + " Columns 1 through 6:\n", + "\n", + " 9.5357e+03 -4.7554e-01 -2.1114e-01 6.0163e-02 -2.4387e-03 6.1052e-04\n", + "\n", + " Columns 7 through 12:\n", + "\n", + " 2.2891e-04 -9.5367e-06 2.3842e-06 8.9407e-07 -2.2352e-07 9.3132e-09\n", + "\n", + " Columns 13 through 18:\n", + "\n", + " -2.3283e-09 -8.7311e-10 3.6380e-11 -9.0949e-12 -3.4106e-12 8.5265e-13\n", + "\n", + " Columns 19 and 20:\n", + "\n", + " -3.5527e-14 8.8818e-15\n", + "\n", + "ans = 16.208\n" + ] + } + ], + "source": [ + "ea_bs\n", + "newtraph(f,df,0.5,0,12)" + ] + }, + { + "cell_type": "code", + "execution_count": 93, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 1.9683e+23\r\n" + ] + } + ], + "source": [ + "df(300)" + ] + }, + { + "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_07/lecture_07.md b/lecture_07/lecture_07.md new file mode 100644 index 0000000..020699b --- /dev/null +++ b/lecture_07/lecture_07.md @@ -0,0 +1,388 @@ + + +```octave +%plot --format svg +``` + + +```octave +setdefaults +``` + +# Roots: Open methods +## Newton-Raphson + +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) + +$f'(x_{i})=\frac{f(x_{i})-0}{x_{i}-x_{i+1}}$ + +$x_{i+1}=x_{i}-\frac{f(x_{i})}{f'(x_{i})}$ + +Use Newton-Raphson to find solution when $e^{x}=x$ + + +```octave +f= @(x) exp(-x)-x; +df= @(x) -exp(-x)-1; + +x_i= 0; +c +error_approx = abs((x_r-x_i)/x_r) +x_i=x_r; + +``` + + error: 'c' undefined near line 1 column 1 + error: 'x_r' undefined near line 1 column 21 + error: evaluating argument list element number 1 + error: 'x_r' undefined near line 1 column 5 + + + +```octave +x_r = x_i-f(x_i)/df(x_i) +error_approx = abs((x_r-x_i)/x_r) +x_i=x_r; +``` + + x_r = 0.50000 + error_approx = 1 + + + +```octave +x_r = x_i-f(x_i)/df(x_i) +error_approx = abs((x_r-x_i)/x_r) +x_i=x_r; +``` + + x_r = 0.56631 + error_approx = 0.11709 + + + +```octave +x_r = x_i-f(x_i)/df(x_i) +error_approx = abs((x_r-x_i)/x_r) +x_i=x_r; +``` + + x_r = 0.56714 + error_approx = 0.0014673 + + +In the bungee jumper 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. + +$f(m)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)-v(t)$. + +to use the Newton-Raphson method, we need the derivative $\frac{df}{dm}$ + +$\frac{df}{dm}=\frac{1}{2}\sqrt{\frac{g}{mc_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)- +\frac{g}{2m}\mathrm{sech}^{2}(\sqrt{\frac{gc_{d}}{m}}t)$ + + +```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 +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; +``` + + +```octave +newtraph(f_m,df_m,140,0.00001) +``` + + ans = 142.74 + + +## Secant Methods + +Not always able to evaluate the derivative. Approximation of derivative: + +$f'(x_{i})=\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}$ + +$x_{i+1}=x_{i}-\frac{f(x_{i})}{f'(x_{i})}$ + +$x_{i+1}=x_{i}-\frac{f(x_{i})}{\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}}= + x_{i}-\frac{f(x_{i})(x_{i-1}-x_{i})}{f(x_{i-1})-f(x_{i})}$ + +What values should $x_{i}$ and $x_{i-1}$ take? + +To reduce arbitrary selection of variables, use the + +## Modified Secant method + +Change the x evaluations to a perturbation $\delta$. + +$x_{i+1}=x_{i}-\frac{f(x_{i})(\delta x_{i})}{f(x_{i}+\delta x_{i})-f(x_{i})}$ + + +```octave +mod_secant(f_m,1e-6,50,0.00001) +``` + + ans = 142.74 + + + +```octave +Amt_numerical=mod_secant(@(A) car_payments(A,30000,0.05,5),1e-6,50,0.001) +car_payments(Amt,30000,0.05,5) +``` + + Amt_numerical = 563.79 + ans = -160.42 + + + +![svg](lecture_07_files/lecture_07_14_1.svg) + + + +```octave +Amt*12*5 +``` + + ans = 3.3968e+04 + + +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. + + +```octave +% Amortization calculation +A = @(P,r,n) P*(r*(1+r)^n)./((1+r)^n-1); +Amt=A(30000,0.05/12,5*12) +``` + + Amt = 566.14 + + +## Matlab's function + +Matlab and Octave combine bracketing and open methods in the `fzero` function. + + +```octave +help fzero +``` + + 'fzero' is a function from the file /usr/share/octave/4.0.0/m/optimization/fzero.m + + -- Function File: fzero (FUN, X0) + -- Function File: fzero (FUN, X0, OPTIONS) + -- Function File: [X, FVAL, INFO, OUTPUT] = fzero (...) + Find a zero of a univariate function. + + FUN is a function handle, inline function, or string containing the + name of the function to evaluate. + + X0 should be a two-element vector specifying two points which + bracket a zero. In other words, there must be a change in sign of + the function between X0(1) and X0(2). More mathematically, the + following must hold + + sign (FUN(X0(1))) * sign (FUN(X0(2))) <= 0 + + If X0 is a single scalar then several nearby and distant values are + probed in an attempt to obtain a valid bracketing. If this is not + successful, the function fails. + + OPTIONS is a structure specifying additional options. Currently, + 'fzero' recognizes these options: "FunValCheck", "OutputFcn", + "TolX", "MaxIter", "MaxFunEvals". For a description of these + options, see *note optimset: XREFoptimset. + + On exit, the function returns X, the approximate zero point and + FVAL, the function value thereof. + + INFO is an exit flag that can have these values: + + * 1 The algorithm converged to a solution. + + * 0 Maximum number of iterations or function evaluations has + been reached. + + * -1 The algorithm has been terminated from user output + function. + + * -5 The algorithm may have converged to a singular point. + + OUTPUT is a structure containing runtime information about the + 'fzero' algorithm. Fields in the structure are: + + * iterations Number of iterations through loop. + + * nfev Number of function evaluations. + + * bracketx A two-element vector with the final bracketing of the + zero along the x-axis. + + * brackety A two-element vector with the final bracketing of the + zero along the y-axis. + + See also: optimset, fsolve. + + 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 +fzero(@(A) car_payments(A,30000,0.05,5,0),500) +``` + + ans = 563.79 + + +## Comparison of Solvers + +It's helpful to compare to the convergence of different routines to see how quickly you find a solution. + +Comparing the freefall example + + + +```octave +N=20; +iterations = linspace(1,400,N); +ea_nr=zeros(1,N); % appr error Newton-Raphson +ea_ms=zeros(1,N); % appr error Modified Secant +ea_fp=zeros(1,N); % appr error false point method +ea_bs=zeros(1,N); % appr error bisect method +for i=1:length(iterations) + [root_nr,ea_nr(i),iter_nr]=newtraph(f_m,df_m,200,0,iterations(i)); + [root_ms,ea_ms(i),iter_ms]=mod_secant(f_m,1e-6,300,0,iterations(i)); + [root_fp,ea_fp(i),iter_fp]=falsepos(f_m,1,300,0,iterations(i)); + [root_bs,ea_bs(i),iter_bs]=bisect(f_m,1,300,0,iterations(i)); +end + +semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs)) +legend('newton-raphson','mod-secant','false point','bisection') +``` + + warning: axis: omitting non-positive data in log plot + warning: called from + __line__ at line 120 column 16 + line at line 56 column 8 + __plt__>__plt2vv__ at line 500 column 10 + __plt__>__plt2__ at line 246 column 14 + __plt__ at line 133 column 15 + semilogy at line 60 column 10 + warning: axis: omitting non-positive data in log plot + warning: axis: omitting non-positive data in log plot + warning: axis: omitting non-positive data in log plot + + + +![svg](lecture_07_files/lecture_07_22_1.svg) + + + +```octave +ea_ms +``` + + ea_ms = + + Columns 1 through 7: + + 43.43883 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + + Columns 8 through 14: + + 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + + Columns 15 through 20: + + 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + + + + +```octave +N=20; +f= @(x) x^10-1; +df=@(x) 10*x^9; +iterations = linspace(1,50,N); +ea_nr=zeros(1,N); % appr error Newton-Raphson +ea_ms=zeros(1,N); % appr error Modified Secant +ea_fp=zeros(1,N); % appr error false point method +ea_bs=zeros(1,N); % appr error bisect method +for i=1:length(iterations) + [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i)); + [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i)); + [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i)); + [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i)); +end + +semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs)) +legend('newton-raphson','mod-secant','false point','bisection') +``` + + warning: axis: omitting non-positive data in log plot + warning: called from + __line__ at line 120 column 16 + line at line 56 column 8 + __plt__>__plt2vv__ at line 500 column 10 + __plt__>__plt2__ at line 246 column 14 + __plt__ at line 133 column 15 + semilogy at line 60 column 10 + warning: axis: omitting non-positive data in log plot + warning: axis: omitting non-positive data in log plot + warning: axis: omitting non-positive data in log plot + + + +![svg](lecture_07_files/lecture_07_24_1.svg) + + + +```octave +ea_bs +newtraph(f,df,0.5,0,12) +``` + + ea_bs = + + Columns 1 through 6: + + 9.5357e+03 -4.7554e-01 -2.1114e-01 6.0163e-02 -2.4387e-03 6.1052e-04 + + Columns 7 through 12: + + 2.2891e-04 -9.5367e-06 2.3842e-06 8.9407e-07 -2.2352e-07 9.3132e-09 + + Columns 13 through 18: + + -2.3283e-09 -8.7311e-10 3.6380e-11 -9.0949e-12 -3.4106e-12 8.5265e-13 + + Columns 19 and 20: + + -3.5527e-14 8.8818e-15 + + ans = 16.208 + + + +```octave +df(300) +``` + + ans = 1.9683e+23 + + + +```octave + +``` diff --git a/lecture_07/lecture_07.pdf b/lecture_07/lecture_07.pdf new file mode 100644 index 0000000000000000000000000000000000000000..4216be6e0844ae02657decc5b169f0f1569b11c0 GIT binary patch literal 68853 zcma&NbF3)A*Cl%HvG1{M+qP}nwr$(CZQHhO+kE$%%;f!E=9`!4KdN_ECp(?4%38Ho zts<5a5}~H2VTL3=y1RaZq-VgR!?QK8faKz$6*0GRGIpR9vC?-k7BV)pH8Q4^GPW^w zGQ*>1W})NZfpm0oFxIz*blccgleFCwf$Moyy;DdM7@Y3GC5J`=$td7%3cJiX=!>JC zK*5WL6d%6my}hErl?-nX+Z@lg1B~zFY}}pl;$)mWa}*+g_9aq4`q@9TLHu?=?)H`? zD1aj-=+F=u0vA$F0*@&QuR*~*boNlajzE{U_Xu3RT3p7}3Bn`H5*vtLkP#p!z31ts zHsRQrj`TH!I3`K4-BaC5GdMv0@;HvuBPy8O-55RCIn0OeAJa2TJ|x4}A|%c}*=4DrjO6XUyi?bE1&D$St@(r5!y2ZZr=(ZpCVT$6>8#_v zC*Cj>Ila|g-uC6-+qCW0N*1srRNF)Vo$D;o`LaeT$`CzWEtRek>BPQUUD?nI9}6xT zHMs#AflMIdH-!*DA4B`tMX@`L5F-7x$}-Jm&tm=@M;oKmjMP+tCfuRmG&@f?aMr>JtNWNZ zdu>aR8D#wV68=CMm1)*-)I$t)s{^hODP_JLh@oa7P4O4baYyp0S|T)x912!N1KxFW z2~g-McyfkuYl_7QmQ{O1psDRtnSIxI0+;~-sO>Ct_~UVcU(^3$;m;U7p@~QJ5ke@B zaCjk_2daH`#QenYyd5s96NOiTD@hIw3k^EW(qOMG|M}}8Fu^tIiXeXh5#CgW{I%9W zET<>NEeXn&A{HIyG8$Sxt^{>;bOAodXEp{Wjhs+Wcs@Lr4#c}q8Uu*iR~kqiiRorRSrX4A z5G*rs&}msTrAe|!9~Asjr66cp5h4ir6YpEGgl8Ogi7WTS@fT{oWc|Qsgx3VqDY&J= z!*YR8(Vw(w^5#=t#1h5m)jla_375eHZgr6(pD!d!2kO1=X(h%ixQ#j_f9het@MB$h zp^s}*aOI?Z6#o&rDU;j!+{Y6IuwGjkS4disGT7~Nx4Ciq5)J^czrI>KkCOc_00xk! z1YrK$Q#g1Tzk-Z#-yX(d?Dh{SX|kWBwNfXExjlivqMiVtgAi&>Ho6cy!wE2AX0U*A zVZjbtH=?J{68Wq2I^s$etP<%{K$)W(*Xe9IMgGUjx5Nx4F_~7hrvP(ohN4I&<Y1lhbITR_DglZl4_T(@arvcW}1Jm2KFQwtchC(orFeZH)dWMf{8Yi&W^@{%_cYk)E0PKiTCN)6({E~_uz?7vZ;6^$*aKV^T*rNCG|+5fGlA_TNC&+@NF$SJNtDC zwlEexlI+Vce*aFiP#x#{ou%U^MKmfQh-Cg6s4M3nI>zRdS9XWTs-`63e0&<23)w*+ zd_kkGQrG5JUlg(+2|L}I<#Iao>D|52XES#oEk8W_VO28{fG^zVOVhdC-AY&x$S!nb zR~7bXmN+hOKZ}r{AIaUzOiB#__WL*!gH@;k1hX6*sx-(+eyZ=gA%NfL8TS=!tL2(e zo(Kj1>|ERk%~&YIZe|eMz|NH!0gMD`A0|O7S43dI80g*|D6)K4tu<(g;g11dx76CC z3YWV;Q%uQ-pO|{ux%WJfqI433>9+I@6_}{Pq?myccrDnZfrJpk411_5uc{DsD`d#HzKXxXUKG5`!gDKCQA*1uB1C@k z4^$25rf>+Vi9@d@x5>}PZs!;KX4~Zz^>l7-gR33s^U-o_nsp7IXroi+Gy=g@$C%>JKPX$lFjzws_=@NuORoG^>e`hg=}O2QQ0=XbeJ%tq+`0&4P({_a`qI; zWkDB>OJi`S1LdLeduMmN&s?X2*a?kJL+P_s`}fgJO7JK0#`no8$o&vAj64wnQ9OD9 zf+?2-te>HtkeMM24n2L{g1ICVlDN}5k8LI{V&UXFpCS&a&ZfSWAUy&CWfm$X=bIt< z7%m|}dxQ?5HFz0UBuuz_l4zb-nb*$``3Lz&KSXCZwD{Eq#u4aKHgN*Dem)O@FjoOi z%N#f-3E6Wve{7s63BRVwPU`?m9T#1SqCYU-T4MwUkcv3hJT0wTfw5AYNTAyBk z9bL=0#k`2{OdN-&mD$|vgTH-2Y<#T!UGD;X6yzEaBH*h6kJ&cyA(b$0H~2kK%tZKLb|;_y_ncpmvoxGBjlO&)=R6k5HMdAiZCl4?2e^^I7e}lHM~gM zouJN!iX8AkxzEG`rr=6iZIQ2Rm$&>#0@%_hTZirHOwUz>)vSb9Znw>9<$3J+V}Y8d zXQJumu*c2ycqMMz#2gRH$xP$)1MxPEZi9s;aVQ&`Te;rtPE?;#Nm45b{x66Pl2S^b zd{q1|ATVHcq3_&;Zx)d~!kE8mf_vVOVvgMae+b{QdUuzVI?fe#AX|`ORH2?T_^+Sl zv3KX)@R%JQsi5($u4k;fR6z`6lPg<8g99@YKU>{gIXq|I-vY!fvEbZ-pi7o7i?HJ) z4Bj0R%)Tx^B=Uk|`o)QR>^U{|5j1dkHRoL$O+zMTKkQMr5%dtCN(!qZ)oXA8*!~W3 zSj7aA;YY(468_FzYb|CGjyg}rg|-RCXiB&FA^0ecNzk{*3;s?!|CtDHgn^mZwH2F1 z$T~BGxC{vJ;US3gkZ8;n(JR=+3R)WtN~t8nBn}esym_LP{hy}pv!xh|g?x`CNy(y3 zoi7Bd9la=}>+Ug=NT6m=3Izx|uf2iILNb^Jjc9ao<1CC)x#Bc?r`Zn8DK-%L&Zs1O zvI#*%Nc!NVZG{nk$d$IN6b5gNi3|ajvWUKBX(6y=_yG0~C2V}O&CWQuDe{JEE>AI6 z+06qwZ;vkYbbN|yD`=sn9nXMuBhq>xu(-RCrdMXbKJ&~B&x+7jVn_litG+PMvRbVB z#dRvr36^rnt~xh{UE>?%5xTcM`@|kSjm5+?I7p(HwD-9P0r66dp+yH}uF}RTS7eg3 z%V$p^asQ1fuAt&Gb7(vnczeT@&U)~-2?>A*#f{D;2MA~+U?*vWIvzB>hEcafDAmK@ zlwpj-0=O*H%~NN2s@q%!(tO(UMH&!++0+qlUB{wlGt{9d+}AH zY>Psmc}L^R_Tbv+J8eI2pQA(LB@nY0X&uR!x-9?Y~K36Klk$KbyX zH02g@0=Xs0sofq6=uHkjI4u6_uIpH3o{l^#yLRY*09Lj07ib)p_v2yIfCihP{J{DYkp-H^R%d^vxn+KCE?xWWM}HLh$y z(uqtgk`?YHs(O5OJ_H*9kXfmSauO%($!(W_Z63&6jy;giZiEmI(Gvk4;BtQ}X&}Ox z&AbW@58P4Rkfy1nbReIe^mKqJ&X5==IyMX>=UivtEu5`IbR6}LgdqL}i+$<4;bI4! zj4lx5sPQW12o5Id&Iq9Z$(~%61b=+RD{n!3`r3X>7&?||oQK>K4JRFd%N>8Pu#7VU zLkJi|yhRSRfeUhgZF>p!p_Xx3iGMEHDbzmH5wAl*x}gWx=LH{+L;t?6F;*{z-T4JR z>e~ASeoEcAc2oMJt8FlA4?Ev$)FDQirbD zYKl@)9x)*}HfrjacBD`|Qm`Ih!`<|P1$x-`$t z!^Bh9?&0;a`jSpOS7%YSXO*fM*?yCu!KD{HWn|IXr71p0OB#PD{pSSY<|)0$4W5!^{Mx34ak?r<`DntZ zY;5^JoyL!*oM60!hUxHT1F!@QYPct$ptr{w9{|khIokGbIIdBGSDv5f9u&)n(zAOl zTuhb_;^n{>AAY+zml&GOwH`FvgTgBF!tN0`GqQI9;q;H9jPM3yyRaT$t+?VIZ|^zB zd`p~07VSn+Smd`%mdSp4>5Dnj^#I+iv**E+Xo+miO${qDOHWOnM6$-k0xZ{vhWHjB zS$QkJ>3c$8`JLv-dJ?3qI0xa;lHuwgCA`ZvTc^U&xr`LTo;29oaq^a1$4joS_Dx%5 znRbyz=bTY_?zLgr336WI(k+_wO$y4E;14hu_3F=bChru#760Jh)eN>lR{};chq&uKnbiFF6%-1RuIX5EyPM4{W>ojbeHVC+d8ty3UN8qt# z%7;Id1`=jT1zRqX29}k+b2D9P8y9s6%6sen)9W}7;LovfsU+RtNKc@aMDRE)VKUOF zf`n7t(x|0e2xs+gww4{{W*?gI>QEG zGV`4j`klJGl$pGyknvek_-d|@`ZMErOQ&Hx14$ke2SY8J9m4u+yxU;9ai1yD#xI(+ z!yd0pw)Cjt0nDTM5ax%K0NR?7PRlKztqbgHM9sq@F!X1Tu^7{9 zqWsQ*SJ%n|q(76WukpKHtdm~Z29CX3HVe_eTs31k48D-jt;YHL#M@N$5%NA-t8YFf zg%t&Og^$iabwfQ>58EC+LD-KEoEddf2R!ZUI>z>$-;)H|+#fJX!-=gT<{&L{OqMj| zpYAg(o7*Ql%8FYY7PV7K62D!_bZKT@#!aZKjBE5STuaKSyH<4$2}hAd=^AUa0QZkbsh;m(+T|rKr{TYd5KpVrPRkZt+M*xvzn>BrzjCB> z;&@8$RqMr@`UDjfVwrAq4^xnQh?-3I!3P!4cQ(*AoHT z2Lkfb?UljP&nvIN_RZr zrT-H@3F#;a9cDiS47l6?(HF+%52@Zh9rklVG1%kx z4M3|_q)bXX%bS5$s-O~pOHp4ryWoPiF%e=#LtZa*s<(&Pe?P>)dViHQ2QB{Giwm%a zY4N(o{t|=~6DfXXm|k94nXmD3SKqImFka5V&QQaiVPNb#_iKAAI->fE3~CG`B@veD zoNKVutk+Ba&7n@fJJ{=ecQ;iIhU@*^{i>6DXVBx+lEZ7+ct7d;9`n2RoYO|^Eo*(9 zGelG%?I)4YtI+{S66Oc&n}*dSIQJFhEnU=a{g}6T9--8+!O674%?vwC9o+56JHwpcqpwmz3PVSk( zJ&sae$_&d323Lv@DsuOkO-M%u$#0$x8vyEEbtMvJOmZwmD=KwZ*4-QlEMYl&HmAgH zbC49^wqhojNrnU)qgYI(TOwg(&yr$mz_?>**m-d=yf80a$T_gPgs~NY5ih;)&PWoTSUlb#d_Zt+uWFnm157Uwg`%C{$-_lM8%E2K1EAVD1ep(#F3nK9a z)Q4_TMK`44HW(lyIaep+0@i(rkQ2YkrYi&E8v2K$U`t$QI66ixFsowTZqWk7W++^4jBDx~?svxPv!N=}#R~}Ugy)gr*U$HP zZ(jX<^aP7-<7=XZ8?&Y7YikH7HW`QZEg{`EuxD+hIoId*5xa-Du;gPEvxS<9oLW-K z;>m+|tJ^KqWb?>%h@*k>4^`b1?jV)fAE8NbzC;=_dMy(@6T&v*ZnGAO+$c*rcsfb& zG<6^slG0{w!DXh+e*t%3U62wb4(T(_YSP9iu(uu)QXgR*oJt@~)^jk?3MyFr{T0q5 zcglEG8ns@e9dTnE0}K@` z?R7%*v%h&)@UKw0F~;Y~&~PLvnc!H+3*eqwf^4u;uOQ1gAG3aO^Y~l=L-PGgRAix04vMJG9p;bHA