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",
+ ""
+ ]
+ },
+ {
+ "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": [
+ ""
+ ],
+ "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": [
+ ""
+ ],
+ "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?
+
+
+
+What cool personal projects are you working on?
+While we delve deeper into Matlab functions, could we review some of the basic logic
+operators it uses and command codes.
+
+I still dont know when these forms are technically due.
+
+ -by the following lecture
+
+I'm having trouble interfacing Atom with GitHub. Is there a simple tutorial for this?
+
+ -Mac? Seems there could be a bug that folks are working on
+
+What are the bear necessities of life?
+please go over how to "submit" the homeworks because it is still confusing
+
+Do you prefer Matlab or Octave?
+
+ -octave is my preference, but Matlab has some benefits
+
+Would you consider a country to be open-source?
+
+ -??
+
+Is there a way to download matlab for free?
+
+ -not legally
+
+how do you add files to current folder in matlab?
+
+ -you can do this either through a file browser or cli
+
+How should Homework 2 be submitted? By simply putting the function into the homework_1
+repository?
+
+ -yes
+
+How can we tell that these forms are being received?
+
+ -when you hit submit, the form says "form received"
+
+can you save scripted outputs from matlab/octave as an excel file?
+
+ -yes, easy way is open a file with a .csv extension then fprintf and separate everything with commas, harder way is to use the `xlswrite`
+
+
+Also, can you update your notes to show what happens when these things are run, as you do
+in class?"
+
+ -I always update the lecture notes after class so they should display what we did in class
+
+I have a little difficulty following along in class on my laptop when you have programs
+pre-written. Maybe if you posted those codes on Github so I could copy them when you
+switch to different desktops I would be able to follow along better.
+
+Kirk or Picard?
+
+ -Kirk
+
+Who is our TA?
+
+ -Peiyu Zhang peiyu.zhang@uconn.edu
+
+Can we download libraries of data like thermodynamic tables into matlab?
+
+-YES! [Matlab Steam Tables](http://bit.ly/2kZygu8)
+
+Will we use the Simulink addition to Matlab? I found it interesting and useful for
+evaluating ODEs in Linear systems.
+
+ -not in this class, everything in simulink has a matlab script/function that can be substituted, but many times its hidden by the gui. Here we want to look directly at our solvers
+
+
+# Roots and Optimization
+## Bracketing ch. 5
+
+When you are given a function, numerical or analytical, it's not always possible to solve directly for a given variable.
+
+Even for the freefall example we first explored,
+
+$v(t)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)$
+
+There is no way to solve for m in terms of the other variables.
+
+Instead, we can solve the problem by creating a new function f(m) where
+
+$f(m)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)-v(t)$.
+
+When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity)
+
+
+```octave
+setdefaults
+g=9.81; % acceleration due to gravity
+m=linspace(50, 200,100); % possible values for mass 50 to 200 kg
+c_d=0.25; % drag coefficient
+t=4; % at time = 4 seconds
+v=36; % speed must be 36 m/s
+f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m
+
+plot(m,f_m(m),m,zeros(length(m),1))
+axis([45 200 -5 1])
+```
+
+
+
+
+
+
+```octave
+f_m(145)
+```
+
+ ans = 0.045626
+
+
+Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0
+
+Better methods are the
+1. Bracketing methods
+2. Open methods
+
+Both need an initial guess.
+
+
+## Incremental method (Brute force)
+
+You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive.
+
+```matlab
+function xb = incsearch(func,xmin,xmax,ns)
+% incsearch: incremental search root locator
+% xb = incsearch(func,xmin,xmax,ns):
+% finds brackets of x that contain sign changes
+% of a function on an interval
+% input:
+% func = name of function
+% xmin, xmax = endpoints of interval
+% ns = number of subintervals (default = 50)
+% output:
+% xb(k,1) is the lower bound of the kth sign change
+% xb(k,2) is the upper bound of the kth sign change
+% If no brackets found, xb = [].
+if nargin < 3, error('at least 3 arguments required'), end
+if nargin < 4, ns = 50; end %if ns blank set to 50
+% Incremental search
+x = linspace(xmin,xmax,ns);
+f = func(x);
+nb = 0; xb = []; %xb is null unless sign change detected
+for k = 1:length(x)-1
+ if sign(f(k)) ~= sign(f(k+1)) %check for sign change
+ nb = nb + 1;
+ xb(nb,1) = x(k);
+ xb(nb,2) = x(k+1);
+ end
+end
+if isempty(xb) %display that no brackets were found
+ fprintf('no brackets found\n')
+ fprintf('check interval or increase ns\n')
+else
+ fprintf('number of brackets: %i\n',nb) %display number of brackets
+end
+```
+
+
+```octave
+help incsearch
+```
+
+ 'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/me3255_S2017/lecture_06/incsearch.m
+
+ incsearch: incremental search root locator
+ xb = incsearch(func,xmin,xmax,ns):
+ finds brackets of x that contain sign changes
+ of a function on an interval
+ input:
+ func = name of function
+ xmin, xmax = endpoints of interval
+ ns = number of subintervals (default = 50)
+ output:
+ xb(k,1) is the lower bound of the kth sign change
+ xb(k,2) is the upper bound of the kth sign change
+ If no brackets found, xb = [].
+
+
+ Additional help for built-in functions and operators is
+ available in the online version of the manual. Use the command
+ 'doc ' to search the manual index.
+
+ Help and information about Octave is also available on the WWW
+ at http://www.octave.org and via the help@octave.org
+ mailing list.
+
+
+
+```octave
+incsearch(f_m,50, 200,55)
+```
+
+ no brackets found
+ check interval or increase ns
+ ans = [](1x0)
+
+
+## Bisection method
+
+Divide interval in half until error is reduced to some level
+
+in previous example of freefall, choose x_l=50, x_u=200
+
+x_r = (50+200)/2 = 125
+
+f_m(125) = -0.408
+
+x_r= (125+200)/2 = 162.5
+
+f_m(162.5) = 0.3594
+
+x_r = (125+162.5)/2=143.75
+
+f_m(143.75)= 0.0206
+
+
+```octave
+f_m(143.75)
+```
+
+ ans = 0.020577
+
+
+Much better root locator, with 4 iterations, our function is already close to zero
+
+Automate this with a function:
+`bisect.m`
+
+```matlab
+function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
+% bisect: root location zeroes
+% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
+% uses bisection method to find the root of func
+% input:
+% func = name of function
+% xl, xu = lower and upper guesses
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by func
+% output:
+% root = real root
+% fx = function value at root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if test>0,error('no sign change'),end
+if nargin<4|isempty(es), es=0.0001;end
+if nargin<5|isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+ xrold = xr;
+ xr = (xl + xu)/2;
+ iter = iter + 1;
+ if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+ test = func(xl,varargin{:})*func(xr,varargin{:});
+ if test < 0
+ xu = xr;
+ elseif test > 0
+ xl = xr;
+ else
+ ea = 0;
+ end
+ if ea <= es | iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});
+```
+
+## False position (linear interpolation)
+
+Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner
+
+$ x_{r} = x_{u} - \frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$
+
+
+```octave
+xl=50; xu=200;
+xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu));
+
+plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0)
+```
+
+
+
+
+
+Much better root locator, with 4 iterations, our function is already close to zero
+
+Automate this with a function:
+`falsepos.m`
+
+```matlab
+function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)
+% falsepos: root location zeroes
+% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
+% uses false position method to find the root of func
+% input:
+% func = name of function
+% xl, xu = lower and upper guesses
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by func
+% output:
+% root = real root
+% fx = function value at root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if test>0,error('no sign change'),end
+if nargin<4|isempty(es), es=0.0001;end
+if nargin<5|isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+ xrold = xr;
+ % xr = (xl + xu)/2; % bisect method
+ xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method
+ iter = iter + 1;
+ if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+ test = func(xl,varargin{:})*func(xr,varargin{:});
+ if test < 0
+ xu = xr;
+ elseif test > 0
+ xl = xr;
+ else
+ ea = 0;
+ end
+ if ea <= es | iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});
+```
+
+
+```octave
+
+```
diff --git a/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 @@
+
\ 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 @@
+
\ No newline at end of file
diff --git a/lecture_07/.ipynb_checkpoints/lecture_07-checkpoint.ipynb b/lecture_07/.ipynb_checkpoints/lecture_07-checkpoint.ipynb
new file mode 100644
index 0000000..d7ff495
--- /dev/null
+++ b/lecture_07/.ipynb_checkpoints/lecture_07-checkpoint.ipynb
@@ -0,0 +1,1042 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "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": 20,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "error_approx = 1\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "f= @(x) exp(-x)-x;\n",
+ "df= @(x) -exp(-x)-1;\n",
+ "\n",
+ "x_i= 0;\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_r=x_i;\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "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": 22,
+ "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": 23,
+ "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": 7,
+ "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": 8,
+ "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": 9,
+ "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": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "error: 'plot_bool' undefined near line 12 column 6\n",
+ "error: called from\n",
+ " car_payments at line 12 column 3\n",
+ " mod_secant at line 22 column 8\n",
+ "error: 'Amt' undefined near line 1 column 14\n",
+ "error: evaluating argument list element number 1\n"
+ ]
+ }
+ ],
+ "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": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "error: 'Amt' undefined near line 1 column 1\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": 12,
+ "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": 13,
+ "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": 14,
+ "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": 15,
+ "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": [
+ ""
+ ],
+ "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": 16,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ea_ms =\n",
+ "\n",
+ " Columns 1 through 6:\n",
+ "\n",
+ " 2.3382e+03 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14\n",
+ "\n",
+ " Columns 7 through 12:\n",
+ "\n",
+ " 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14\n",
+ "\n",
+ " Columns 13 through 18:\n",
+ "\n",
+ " 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14\n",
+ "\n",
+ " Columns 19 and 20:\n",
+ "\n",
+ " 1.9171e-14 1.9171e-14\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "ea_ms"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "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": [
+ ""
+ ],
+ "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": 18,
+ "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": 19,
+ "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/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..d7ff495
--- /dev/null
+++ b/lecture_07/lecture_07.ipynb
@@ -0,0 +1,1042 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "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": 20,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "error_approx = 1\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "f= @(x) exp(-x)-x;\n",
+ "df= @(x) -exp(-x)-1;\n",
+ "\n",
+ "x_i= 0;\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_r=x_i;\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "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": 22,
+ "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": 23,
+ "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": 7,
+ "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": 8,
+ "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": 9,
+ "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": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "error: 'plot_bool' undefined near line 12 column 6\n",
+ "error: called from\n",
+ " car_payments at line 12 column 3\n",
+ " mod_secant at line 22 column 8\n",
+ "error: 'Amt' undefined near line 1 column 14\n",
+ "error: evaluating argument list element number 1\n"
+ ]
+ }
+ ],
+ "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": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "error: 'Amt' undefined near line 1 column 1\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": 12,
+ "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": 13,
+ "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": 14,
+ "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": 15,
+ "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": [
+ ""
+ ],
+ "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": 16,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ea_ms =\n",
+ "\n",
+ " Columns 1 through 6:\n",
+ "\n",
+ " 2.3382e+03 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14\n",
+ "\n",
+ " Columns 7 through 12:\n",
+ "\n",
+ " 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14\n",
+ "\n",
+ " Columns 13 through 18:\n",
+ "\n",
+ " 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14\n",
+ "\n",
+ " Columns 19 and 20:\n",
+ "\n",
+ " 1.9171e-14 1.9171e-14\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "ea_ms"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "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": [
+ ""
+ ],
+ "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": 18,
+ "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": 19,
+ "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..ffdc48d
--- /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;
+error_approx = abs((x_r-x_i)/x_r)
+x_r=x_i;
+
+```
+
+ 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.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)
+```
+
+ error: 'plot_bool' undefined near line 12 column 6
+ error: called from
+ car_payments at line 12 column 3
+ mod_secant at line 22 column 8
+ error: 'Amt' undefined near line 1 column 14
+ error: evaluating argument list element number 1
+
+
+
+```octave
+Amt*12*5
+```
+
+ error: 'Amt' undefined near line 1 column 1
+
+
+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
+
+
+
+
+
+
+
+```octave
+ea_ms
+```
+
+ ea_ms =
+
+ Columns 1 through 6:
+
+ 2.3382e+03 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14
+
+ Columns 7 through 12:
+
+ 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14
+
+ Columns 13 through 18:
+
+ 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14
+
+ Columns 19 and 20:
+
+ 1.9171e-14 1.9171e-14
+
+
+
+
+```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
+
+
+
+
+
+
+
+```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 0000000..79b130c
Binary files /dev/null and b/lecture_07/lecture_07.pdf differ
diff --git a/lecture_07/lecture_07_files/lecture_07_14_1.svg b/lecture_07/lecture_07_files/lecture_07_14_1.svg
new file mode 100644
index 0000000..9d54798
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_14_1.svg
@@ -0,0 +1,146 @@
+
\ No newline at end of file
diff --git a/lecture_07/lecture_07_files/lecture_07_22_1.svg b/lecture_07/lecture_07_files/lecture_07_22_1.svg
new file mode 100644
index 0000000..502c4c4
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_22_1.svg
@@ -0,0 +1,197 @@
+
\ No newline at end of file
diff --git a/lecture_07/lecture_07_files/lecture_07_24_1.svg b/lecture_07/lecture_07_files/lecture_07_24_1.svg
new file mode 100644
index 0000000..a8b614c
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_24_1.svg
@@ -0,0 +1,182 @@
+
\ No newline at end of file
diff --git a/lecture_07/mod_secant.m b/lecture_07/mod_secant.m
new file mode 100644
index 0000000..a3caa10
--- /dev/null
+++ b/lecture_07/mod_secant.m
@@ -0,0 +1,31 @@
+function [root,ea,iter]=mod_secant(func,dx,xr,es,maxit,varargin)
+% newtraph: Modified secant root location zeroes
+% [root,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...):
+% uses modified secant method to find the root of func
+% input:
+% func = name of function
+% dx = perturbation fraction
+% xr = initial guess
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by function
+% output:
+% root = real root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+if nargin<4 || isempty(es),es=0.0001;end
+if nargin<5 || isempty(maxit),maxit=50;end
+iter = 0;
+while (1)
+ xrold = xr;
+ dfunc=(func(xr+dx)-func(xr))./dx;
+ xr = xr - func(xr)/dfunc;
+ iter = iter + 1;
+ if xr ~= 0
+ ea = abs((xr - xrold)/xr) * 100;
+ end
+ if ea <= es || iter >= maxit, break, end
+end
+root = xr;
+
diff --git a/lecture_07/newtraph.m b/lecture_07/newtraph.m
new file mode 100644
index 0000000..94ca667
--- /dev/null
+++ b/lecture_07/newtraph.m
@@ -0,0 +1,29 @@
+function [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,varargin)
+% newtraph: Newton-Raphson root location zeroes
+% [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,p1,p2,...):
+% uses Newton-Raphson method to find the root of func
+% input:
+% func = name of function
+% dfunc = name of derivative of function
+% xr = initial guess
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by function
+% output:
+% root = real root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+if nargin<4 || isempty(es),es=0.0001;end
+if nargin<5 || isempty(maxit),maxit=50;end
+iter = 0;
+while (1)
+ xrold = xr;
+ xr = xr - func(xr)/dfunc(xr);
+ iter = iter + 1;
+ if xr ~= 0
+ ea = abs((xr - xrold)/xr) * 100;
+ end
+ if ea <= es || iter >= maxit, break, end
+end
+root = xr;
diff --git a/lecture_07/octave-workspace b/lecture_07/octave-workspace
new file mode 100644
index 0000000..a214fe9
Binary files /dev/null and b/lecture_07/octave-workspace differ