diff --git a/HW3/README.md b/HW3/README.md
new file mode 100644
index 0000000..45214a4
--- /dev/null
+++ b/HW3/README.md
@@ -0,0 +1,73 @@
+# Homework #3
+## due 2/15/17 by 11:59pm
+
+
+1. Create a new github repository called 'roots_and_optimization'.
+
+ a. Add rcc02007 and pez16103 as collaborators.
+
+ b. Clone the repository to your computer.
+
+ c. Copy your `projectile.m` function into the 'roots_and_optimization' folder.
+ *Disable the plotting routine for the solvers*
+
+ d. Use the four solvers `falsepos.m`, `bisect.m`, `newtraph.m` and `mod_secant.m`
+ to solve for the angle needed to reach h=1.72 m, with an initial speed of 15 m/s.
+
+ e. The `newtraph.m` function needs a derivative, calculate the derivative of your
+ function with respect to theta, `dprojectile_dtheta.m`. This function should
+ output d(h)/d(theta).
+
+
+ f. In your `README.md` file, document the following under the heading `#
+ Homework #3`:
+
+ i. Compare the number of iterations that each function needed to reach an
+ accuracy of 0.00001%. Include a table in your README.md with:
+
+ ```
+ | solver | initial guess(es) | ea | number of iterations|
+ | --- | --- | --- | --- |
+ |falsepos | | | |
+ |bisect | | | |
+ |newtraph | | | |
+ |mod_secant | | | |
+ ```
+
+ ii. Compare the convergence of the 4 methods. Plot the approximate error vs the
+ number of iterations that the solver has calculated. Save the plot as
+ `convergence.png` and display the plot in your `README.md` with:
+
+ ``
+
+ iii. In the `README.md` provide a description of the files used to create the
+ table and the convergence plot.
+
+2. The Newton-Raphson method and the modified secant method do not always converge to a
+solution. One simple example is the function f(x) = x*exp(-x^2). The root is at 0, but
+using the numerical solvers, `newtraph.m` and `mod_secant.m`, there are certain initial
+guesses that do not converge.
+
+ a. Calculate the first 5 iterations for the Newton-Raphson method with an initial
+ guess of x_i=2 for f(x)=x*exp(-x^2).
+
+ b. Add the results to a table in the `README.md` with:
+
+ ```
+ ### divergence of Newton-Raphson method
+
+ | iteration | x_i | approx error |
+ | --- | --- | --- |
+ | 0 | 2 | n/a |
+ | 1 | | |
+ | 2 | | |
+ | 3 | | |
+ | 4 | | |
+ | 5 | | |
+ ```
+
+ c. Repeat steps a-b for an initial guess of 0.2. (But change the heading from
+ 'divergence' to 'convergence')
+
+3. Commit your changes to your repository. Sync your local repository with github. Then
+copy and paste the "clone URL" into the following Google Form [Homework 3](https://goo.gl/forms/UJBGwp0fQcSxImkq2)
diff --git a/README.md b/README.md
index 8059993..73ca14a 100644
--- a/README.md
+++ b/README.md
@@ -74,8 +74,8 @@ general, I will not post homework solutions.
|3|1/31||Consistent Coding habits|
| |2/2|5|Root Finding|
|4|2/7|6|Root Finding con’d|
-| |2/9|7|Optimization|
-|5|2/14||Intro to Linear Algebra|
+| |2/9|7| **Snow Day**|
+|5|2/14|| Optimization |
| |2/16|8|Linear Algebra|
|6|2/21|9|Linear systems: Gauss elimination|
| |2/23|10|Linear Systems: LU factorization|
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/.lecture_07.md.swp b/lecture_07/.lecture_07.md.swp
new file mode 100644
index 0000000..0998c86
Binary files /dev/null and b/lecture_07/.lecture_07.md.swp differ
diff --git a/lecture_07/.newtraph.m.swp b/lecture_07/.newtraph.m.swp
new file mode 100644
index 0000000..9759dd2
Binary files /dev/null and b/lecture_07/.newtraph.m.swp differ
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..08d29c4
--- /dev/null
+++ b/lecture_07/lecture_07.ipynb
@@ -0,0 +1,1458 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "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": 2,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 0.50000\n",
+ "error_approx = 1\n"
+ ]
+ }
+ ],
+ "source": [
+ "f= @(x) exp(-x)-x;\n",
+ "df= @(x) -exp(-x)-1;\n",
+ "\n",
+ "x_i= 0;\n",
+ "x_r = x_i-f(x_i)/df(x_i)\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 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": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 0.56714\n",
+ "error_approx = 2.2106e-07\n"
+ ]
+ }
+ ],
+ "source": [
+ "x_r = x_i-f(x_i)/df(x_i)\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;"
+ ]
+ },
+ {
+ "cell_type": "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": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "root = 142.74\n",
+ "ea = 8.0930e-06\n",
+ "iter = 48\n"
+ ]
+ }
+ ],
+ "source": [
+ "[root,ea,iter]=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": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "root = 142.74\n",
+ "ea = 3.0615e-07\n",
+ "iter = 7\n"
+ ]
+ }
+ ],
+ "source": [
+ "[root,ea,iter]=mod_secant(f_m,1,50,0.00001)"
+ ]
+ },
+ {
+ "cell_type": "raw",
+ "metadata": {},
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 1.1185e+04\r\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "car_payments(400,30000,0.05,5,1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Amt_numerical = 5467.0\n",
+ "ans = 3.9755e-04\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "Amt_numerical=mod_secant(@(A) car_payments(A,700000,0.0875,30,0),1e-6,50,0.001)\n",
+ "car_payments(Amt_numerical,700000,0.0875,30,1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 1.9681e+06\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "Amt_numerical*12*30"
+ ]
+ },
+ {
+ "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": 19,
+ "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": 20,
+ "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": 23,
+ "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,300,0,iterations(i));\n",
+ " [root_ms,ea_ms(i),iter_ms]=mod_secant(f_m,1e-6,300,0,iterations(i));\n",
+ " [root_fp,ea_fp(i),iter_fp]=falsepos(f_m,1,300,0,iterations(i));\n",
+ " [root_bs,ea_bs(i),iter_bs]=bisect(f_m,1,300,0,iterations(i));\n",
+ "end\n",
+ "\n",
+ "setdefaults\n",
+ "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n",
+ "legend('newton-raphson','mod-secant','false point','bisection')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ea_nr =\n",
+ "\n",
+ " Columns 1 through 8:\n",
+ "\n",
+ " 6.36591 0.06436 0.00052 0.00000 0.00000 0.00000 0.00000 0.00000\n",
+ "\n",
+ " Columns 9 through 16:\n",
+ "\n",
+ " 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n",
+ "\n",
+ " Columns 17 through 20:\n",
+ "\n",
+ " 0.00000 0.00000 0.00000 0.00000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "ea_nr"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: called from\n",
+ " __line__ at line 120 column 16\n",
+ " line at line 56 column 8\n",
+ " __plt__>__plt2vv__ at line 500 column 10\n",
+ " __plt__>__plt2__ at line 246 column 14\n",
+ " __plt__ at line 133 column 15\n",
+ " semilogy at line 60 column 10\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=20;\n",
+ "f= @(x) x^10-1;\n",
+ "df=@(x) 10*x^9;\n",
+ "iterations = linspace(1,50,N);\n",
+ "ea_nr=zeros(1,N); % appr error Newton-Raphson\n",
+ "ea_ms=zeros(1,N); % appr error Modified Secant\n",
+ "ea_fp=zeros(1,N); % appr error false point method\n",
+ "ea_bs=zeros(1,N); % appr error bisect method\n",
+ "for i=1:length(iterations)\n",
+ " [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i));\n",
+ " [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i));\n",
+ " [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i));\n",
+ " [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i));\n",
+ "end\n",
+ " \n",
+ "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n",
+ "legend('newton-raphson','mod-secant','false point','bisection')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ea_nr =\n",
+ "\n",
+ " Columns 1 through 7:\n",
+ "\n",
+ " 99.03195 11.11111 11.11111 11.11111 11.11111 11.11111 11.11111\n",
+ "\n",
+ " Columns 8 through 14:\n",
+ "\n",
+ " 11.11111 11.11111 11.11111 11.11109 11.11052 11.10624 10.99684\n",
+ "\n",
+ " Columns 15 through 20:\n",
+ "\n",
+ " 8.76956 2.12993 0.00000 0.00000 0.00000 0.00000\n",
+ "\n",
+ "ans = 16.208\n"
+ ]
+ }
+ ],
+ "source": [
+ "ea_nr\n",
+ "newtraph(f,df,0.5,0,12)"
+ ]
+ },
+ {
+ "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": 28,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "f =\n",
+ "\n",
+ "@(x) tan (x) - (x - 1) .^ 2\n",
+ "\n",
+ "ans = 0.37375\n"
+ ]
+ }
+ ],
+ "source": [
+ "% our class function\n",
+ "f= @(x) tan(x)-(x-1).^2\n",
+ "mod_secant(f,1e-3,1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = -3.5577e-13\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "f(ans)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 30,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 0.39218\n",
+ "ans = 0.39219\n"
+ ]
+ }
+ ],
+ "source": [
+ "tan(0.37375)\n",
+ "(0.37375-1)^2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " Columns 1 through 8:\n",
+ "\n",
+ " -1.0000 1.5574 -3.1850 -4.1425 -7.8422 -19.3805 -25.2910 -35.1286\n",
+ "\n",
+ " Columns 9 through 11:\n",
+ "\n",
+ " -55.7997 -64.4523 -80.3516\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "f([0:10])"
+ ]
+ },
+ {
+ "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..3aacefb
--- /dev/null
+++ b/lecture_07/lecture_07.md
@@ -0,0 +1,448 @@
+
+
+```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;
+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
+
+
+
+```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 = 2.2106e-07
+
+
+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
+[root,ea,iter]=newtraph(f_m,df_m,140,0.00001)
+```
+
+ root = 142.74
+ ea = 8.0930e-06
+ iter = 48
+
+
+## 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
+[root,ea,iter]=mod_secant(f_m,1,50,0.00001)
+```
+
+ root = 142.74
+ ea = 3.0615e-07
+ iter = 7
+
+
+
+```octave
+car_payments(400,30000,0.05,5,1)
+```
+
+ ans = 1.1185e+04
+
+
+
+
+
+
+
+```octave
+Amt_numerical=mod_secant(@(A) car_payments(A,700000,0.0875,30,0),1e-6,50,0.001)
+car_payments(Amt_numerical,700000,0.0875,30,1)
+```
+
+ Amt_numerical = 5467.0
+ ans = 3.9755e-04
+
+
+
+
+
+
+
+```octave
+Amt_numerical*12*30
+```
+
+ ans = 1.9681e+06
+
+
+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,300,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
+
+setdefaults
+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_nr
+```
+
+ ea_nr =
+
+ Columns 1 through 8:
+
+ 6.36591 0.06436 0.00052 0.00000 0.00000 0.00000 0.00000 0.00000
+
+ Columns 9 through 16:
+
+ 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+
+ Columns 17 through 20:
+
+ 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
+
+
+
+
+
+
+
+```octave
+ea_nr
+newtraph(f,df,0.5,0,12)
+```
+
+ ea_nr =
+
+ Columns 1 through 7:
+
+ 99.03195 11.11111 11.11111 11.11111 11.11111 11.11111 11.11111
+
+ Columns 8 through 14:
+
+ 11.11111 11.11111 11.11111 11.11109 11.11052 11.10624 10.99684
+
+ Columns 15 through 20:
+
+ 8.76956 2.12993 0.00000 0.00000 0.00000 0.00000
+
+ ans = 16.208
+
+
+
+```octave
+df(300)
+```
+
+ ans = 1.9683e+23
+
+
+
+```octave
+% our class function
+f= @(x) tan(x)-(x-1).^2
+mod_secant(f,1e-3,1)
+```
+
+ f =
+
+ @(x) tan (x) - (x - 1) .^ 2
+
+ ans = 0.37375
+
+
+
+```octave
+f(ans)
+```
+
+ ans = -3.5577e-13
+
+
+
+```octave
+tan(0.37375)
+(0.37375-1)^2
+```
+
+ ans = 0.39218
+ ans = 0.39219
+
+
+
+```octave
+f([0:10])
+```
+
+ ans =
+
+ Columns 1 through 8:
+
+ -1.0000 1.5574 -3.1850 -4.1425 -7.8422 -19.3805 -25.2910 -35.1286
+
+ Columns 9 through 11:
+
+ -55.7997 -64.4523 -80.3516
+
+
+
+
+```octave
+
+```
diff --git a/lecture_07/lecture_07.pdf b/lecture_07/lecture_07.pdf
new file mode 100644
index 0000000..70d8fc2
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_15_1.svg b/lecture_07/lecture_07_files/lecture_07_15_1.svg
new file mode 100644
index 0000000..903c445
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_15_1.svg
@@ -0,0 +1,131 @@
+
\ No newline at end of file
diff --git a/lecture_07/lecture_07_files/lecture_07_16_1.svg b/lecture_07/lecture_07_files/lecture_07_16_1.svg
new file mode 100644
index 0000000..a3a2125
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_16_1.svg
@@ -0,0 +1,151 @@
+
\ 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..b511b10
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_24_1.svg
@@ -0,0 +1,197 @@
+
\ No newline at end of file
diff --git a/lecture_07/lecture_07_files/lecture_07_26_1.svg b/lecture_07/lecture_07_files/lecture_07_26_1.svg
new file mode 100644
index 0000000..a8b614c
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_26_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
diff --git a/lecture_08/.ipynb_checkpoints/lecture_08-checkpoint.ipynb b/lecture_08/.ipynb_checkpoints/lecture_08-checkpoint.ipynb
new file mode 100644
index 0000000..2fd6442
--- /dev/null
+++ b/lecture_08/.ipynb_checkpoints/lecture_08-checkpoint.ipynb
@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_08/Auchain_model.gif b/lecture_08/Auchain_model.gif
new file mode 100644
index 0000000..2f8c78c
Binary files /dev/null and b/lecture_08/Auchain_model.gif differ
diff --git a/lecture_08/Auchain_model.png b/lecture_08/Auchain_model.png
new file mode 100644
index 0000000..dce9cd7
Binary files /dev/null and b/lecture_08/Auchain_model.png differ
diff --git a/lecture_08/au_chain.jpg b/lecture_08/au_chain.jpg
new file mode 100644
index 0000000..492fee8
Binary files /dev/null and b/lecture_08/au_chain.jpg differ
diff --git a/lecture_08/goldenratio.png b/lecture_08/goldenratio.png
new file mode 100644
index 0000000..9493c8c
Binary files /dev/null and b/lecture_08/goldenratio.png differ
diff --git a/lecture_08/goldmin.m b/lecture_08/goldmin.m
new file mode 100644
index 0000000..fb4ce0b
--- /dev/null
+++ b/lecture_08/goldmin.m
@@ -0,0 +1,36 @@
+function [x,fx,ea,iter]=goldmin(f,xl,xu,es,maxit,varargin)
+% goldmin: minimization golden section search
+% [x,fx,ea,iter]=goldmin(f,xl,xu,es,maxit,p1,p2,...):
+% uses golden section search to find the minimum of f
+% input:
+% f = 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 f
+% output:
+% x = location of minimum
+% fx = minimum function value
+% 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
+phi=(1+sqrt(5))/2;
+iter=0;
+while(1)
+ d = (phi-1)*(xu - xl);
+ x1 = xl + d;
+ x2 = xu - d;
+ if f(x1,varargin{:}) < f(x2,varargin{:})
+ xopt = x1;
+ xl = x2;
+ else
+ xopt = x2;
+ xu = x1;
+ end
+ iter=iter+1;
+ if xopt~=0, ea = (2 - phi) * abs((xu - xl) / xopt) * 100;end
+ if ea <= es | iter >= maxit,break,end
+end
+x=xopt;fx=f(xopt,varargin{:});
diff --git a/lecture_08/lecture_08.ipynb b/lecture_08/lecture_08.ipynb
new file mode 100644
index 0000000..32b973c
--- /dev/null
+++ b/lecture_08/lecture_08.ipynb
@@ -0,0 +1,2762 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "function h = parabola(x)\n",
+ " y=1;\n",
+ " h=sum(x.^2+y.^2-1);\n",
+ "end\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 0\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "fzero(@(x) parabola(x),1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Optimization\n",
+ "\n",
+ "Many problems involve finding a minimum or maximum based on given constraints. Engineers and scientists typically use energy balance equations to find the conditions of minimum energy, but this value may never go to 0 and the actual value of energy may not be of interest. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The Lennard-Jones potential is commonly used to model interatomic bonding. \n",
+ "\n",
+ "$E_{LJ}(x)=4\\epsilon \\left(\\left(\\frac{\\sigma}{x}\\right)^{12}-\\left(\\frac{\\sigma}{x}\\right)^{6}\\right)$\n",
+ "\n",
+ "Considering a 1-D gold chain, we can calculate the bond length, $x_{b}$, with no force applied to the chain and even for tension, F. This will allow us to calculate the nonlinear spring constant of a 1-D gold chain. \n",
+ "\n",
+ "\n",
+ "\n",
+ "Computational Tools to Study and Predict the Long-Term Stability of Nanowires.\n",
+ "By Martin E. Zoloff Michoff, Patricio Vélez, Sergio A. Dassie and Ezequiel P. M. Leiva \n",
+ "\n",
+ "\n",
+ "\n",
+ "[Single atom gold chain mechanics](http://www.uam.es/personal_pdi/ciencias/agrait/)\n",
+ "\n",
+ "### First, let's find the minimum energy $\\min(E_{LJ}(x))$\n",
+ "\n",
+ "## Brute force"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "setdefaults\n",
+ "epsilon = 0.039; % kcal/mol\n",
+ "sigma = 2.934; % Angstrom\n",
+ "x=linspace(2.8,6,200); % bond length in Angstrom\n",
+ "\n",
+ "Ex = lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "[Emin,imin]=min(Ex);\n",
+ "\n",
+ "plot(x,Ex,x(imin),Emin,'o')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 3.2824\n",
+ "ans = 3.2985\n",
+ "ans = 3.3146\n"
+ ]
+ }
+ ],
+ "source": [
+ "x(imin-1)\n",
+ "x(imin)\n",
+ "x(imin+1)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Golden Search Algorithm\n",
+ "\n",
+ "We can't just look for a sign change for the problem (unless we can take a derivative) so we need a new approach to determine whether we have a maximum between the two bounds.\n",
+ "\n",
+ "Rather than using the midpoint of initial bounds, here the problem is more difficult. We need to compare the values of 4 function evaluations. The golden search uses the golden ratio to determine two interior points. \n",
+ "\n",
+ "\n",
+ "\n",
+ "Start with bounds of 2.5 and 6 Angstrom. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "current_min = -0.019959\r\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% define Au atomic potential\n",
+ "epsilon = 0.039; % kcal/mol\n",
+ "sigma = 2.934; % Angstrom\n",
+ "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "% calculate golden ratio\n",
+ "phi = 1/2+sqrt(5)/2;\n",
+ "% set initial limits\n",
+ "x_l=2.8; \n",
+ "x_u=6; \n",
+ "\n",
+ "% Iteration #1\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\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",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\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 \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_6a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Iteration #2\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\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",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\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 \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_6a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Iteration #3\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\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",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\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 \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_6a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Iteration #3\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\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",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\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",
+ "\tgnuplot_plot_3a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_4a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% define Au atomic potential\n",
+ "epsilon = 0.039; % kcal/mol\n",
+ "sigma = 2.934; % Angstrom\n",
+ "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "% set initial limits\n",
+ "x_l=2.8; \n",
+ "x_u=3.5; \n",
+ "\n",
+ "% Iteration #1\n",
+ "x1=x_l;\n",
+ "x2=mean([x_l,x_u]);\n",
+ "x3=x_u;\n",
+ "\n",
+ "% evaluate Au_x(x1), Au_x(x2) and Au_x(x3)\n",
+ " \n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "f3=Au_x(x3);\n",
+ "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n",
+ "x_fit = linspace(x1,x3,20);\n",
+ "y_fit = polyval(p,x_fit);\n",
+ "\n",
+ "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n",
+ "hold on\n",
+ "if f2x2\n",
+ " plot(x4,f4,'*',[x1,x2],[f1,f2])\n",
+ " x1=x2;\n",
+ " f1=f2;\n",
+ " else\n",
+ " plot(x4,f4,'*',[x3,x2],[f3,f2])\n",
+ " x3=x2;\n",
+ " f3=f2;\n",
+ " end\n",
+ " x2=x4; f2=f4;\n",
+ "else\n",
+ " error('no minimum in bracket')\n",
+ "end\n",
+ "hold off"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n",
+ "x_fit = linspace(x1,x3,20);\n",
+ "y_fit = polyval(p,x_fit);\n",
+ "\n",
+ "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n",
+ "hold on\n",
+ "if f2x2\n",
+ " plot(x4,f4,'*',[x1,x2],[f1,f2])\n",
+ " x1=x2;\n",
+ " f1=f2;\n",
+ " else\n",
+ " plot(x4,f4,'*',[x3,x2],[f3,f2])\n",
+ " x3=x2;\n",
+ " f3=f2;\n",
+ " end\n",
+ " x2=x4; f2=f4;\n",
+ "else\n",
+ " error('no minimum in bracket')\n",
+ "end\n",
+ "hold off\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n",
+ "x_fit = linspace(x1,x3,20);\n",
+ "y_fit = polyval(p,x_fit);\n",
+ "\n",
+ "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n",
+ "hold on\n",
+ "if f2x2\n",
+ " plot(x4,f4,'*',[x1,x2],[f1,f2])\n",
+ " x1=x2;\n",
+ " f1=f2;\n",
+ " else\n",
+ " plot(x4,f4,'*',[x3,x2],[f3,f2])\n",
+ " x3=x2;\n",
+ " f3=f2;\n",
+ " end\n",
+ " x2=x4; f2=f4;\n",
+ "else\n",
+ " error('no minimum in bracket')\n",
+ "end\n",
+ "hold off\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Parabolic interpolation does not converge in many scenarios even though it it a bracketing method. Instead, functions like `fminbnd` in Matlab and Octave use a combination of the two (Golden Ratio and Parabolic)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Using the solutions to minimization for the nonlinear spring constant\n",
+ "\n",
+ "Now, we have two routines to find minimums of a univariate function (Golden Ratio and Parabolic). Let's use these to solve for the minimum energy given a range of applied forces to the single atom gold chain\n",
+ "\n",
+ "$E_{total}(\\Delta x) = E_{LJ}(x_{min}+\\Delta x) - F \\cdot \\Delta x$\n",
+ "\n",
+ "$1 aJ = 10^{-18} J = 1~nN* 1~nm = 10^{-9}N * 10^{-9} m$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "warning: Matlab-style short-circuit operation performed for operator |\n",
+ "warning: called from\n",
+ " goldmin at line 17 column 1\n",
+ "warning: Matlab-style short-circuit operation performed for operator |\n",
+ "warning: Matlab-style short-circuit operation performed for operator |\n",
+ "xmin = 0.32933\n",
+ "Emin = -2.7096e-04\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "epsilon = 0.039; % kcal/mol\n",
+ "epsilon = epsilon*6.9477e-21; % J/atom\n",
+ "epsilon = epsilon*1e18; % aJ/J\n",
+ "% final units for epsilon are aJ\n",
+ "\n",
+ "sigma = 2.934; % Angstrom\n",
+ "sigma = sigma*0.10; % nm/Angstrom\n",
+ "x=linspace(2.8,6,200)*0.10; % bond length in um\n",
+ "\n",
+ "Ex = lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "%[Emin,imin]=min(Ex);\n",
+ "\n",
+ "[xmin,Emin] = goldmin(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)\n",
+ "\n",
+ "plot(x,Ex,xmin,Emin,'o')\n",
+ "ylabel('Lennard Jones Potential (aJ/atom)')\n",
+ "xlabel('bond length (nm)')\n",
+ "\n",
+ "Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=50;\n",
+ "warning('off')\n",
+ "dx = zeros(1,N); % [in nm]\n",
+ "F_applied=linspace(0,0.0022,N); % [in nN]\n",
+ "for i=1:N\n",
+ " optmin=goldmin(@(dx) Etotal(dx,F_applied(i)),-0.001,0.035);\n",
+ " dx(i)=optmin;\n",
+ "end\n",
+ "\n",
+ "plot(dx,F_applied)\n",
+ "xlabel('dx (nm)')\n",
+ "ylabel('Force (nN)')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For this function, it is possible to take a derivative and compare the analytical result:\n",
+ "\n",
+ "dE/dx = 0 = d(E_LJ)/dx - F\n",
+ "\n",
+ "F= d(E_LJ)/dx"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 32,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "F =\n",
+ "\n",
+ "@(dx) 4 * epsilon * 6 * (sigma ^ 6 ./ (xmin + dx) .^ 7 - 2 * sigma ^ 12 ./ (xmin + dx) .^ 13)\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "dx_full=linspace(0,0.06,50);\n",
+ "F= @(dx) 4*epsilon*6*(sigma^6./(xmin+dx).^7-2*sigma^12./(xmin+dx).^13)\n",
+ "plot(dx_full,F(dx_full),dx,F_applied)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "## Curve-fitting\n",
+ "Another example is minimizing error in your approximation of a function. If you have data (now we have Force-displacement data) we can fit this to a function, such as:\n",
+ "\n",
+ "$F(x) = K_{1}\\Delta x + \\frac{1}{2} K_{2}(\\Delta x)^{2}$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 33,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "function SSE = sse_of_parabola(K,xdata,ydata)\n",
+ " % calculate the sum of squares error for a parabola given a function, func, and xdata and ydata\n",
+ " % output is SSE=sum of squares error\n",
+ " K1=K(1);\n",
+ " K2=K(2);\n",
+ " y_function = K1*xdata+1/2*K2*xdata.^2;\n",
+ " SSE = sum((ydata-y_function).^2);\n",
+ "end\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "The nonlinear spring constants are K1=0.16 nN/nm and K2=-5.98 nN/nm^2\n",
+ "The mininum sum of squares error = 7.35e-08\n"
+ ]
+ }
+ ],
+ "source": [
+ "[K,SSE_min]=fminsearch(@(K) sse_of_parabola(K,dx,F_applied),[1,1]);\n",
+ "fprintf('\\nThe nonlinear spring constants are K1=%1.2f nN/nm and K2=%1.2f nN/nm^2\\n',K)\n",
+ "fprintf('The mininum sum of squares error = %1.2e',SSE_min)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 35,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2)"
+ ]
+ },
+ {
+ "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_08/lecture_08.md b/lecture_08/lecture_08.md
new file mode 100644
index 0000000..79edc0a
--- /dev/null
+++ b/lecture_08/lecture_08.md
@@ -0,0 +1,626 @@
+
+
+```octave
+%plot --format svg
+```
+
+# Optimization
+
+Many problems involve finding a minimum or maximum based on given constraints. Engineers and scientists typically use energy balance equations to find the conditions of minimum energy, but this value may never go to 0 and the actual value of energy may not be of interest.
+
+The Lennard-Jones potential is commonly used to model interatomic bonding.
+
+$E_{LJ}(x)=4\epsilon \left(\left(\frac{\sigma}{x}\right)^{12}-\left(\frac{\sigma}{x}\right)^{6}\right)$
+
+Considering a 1-D gold chain, we can calculate the bond length, $x_{b}$, with no force applied to the chain and even for tension, F. This will allow us to calculate the nonlinear spring constant of a 1-D gold chain.
+
+
+
+Computational Tools to Study and Predict the Long-Term Stability of Nanowires.
+By Martin E. Zoloff Michoff, Patricio Vélez, Sergio A. Dassie and Ezequiel P. M. Leiva
+
+
+
+[Single atom gold chain mechanics](http://www.uam.es/personal_pdi/ciencias/agrait/)
+
+### First, let's find the minimum energy $\min(E_{LJ}(x))$
+
+## Brute force
+
+
+```octave
+setdefaults
+epsilon = 0.039; % kcal/mol
+sigma = 2.934; % Angstrom
+x=linspace(2.8,6,200); % bond length in Angstrom
+
+Ex = lennard_jones(x,sigma,epsilon);
+
+[Emin,imin]=min(Ex);
+
+plot(x,Ex,x(imin),Emin,'o')
+```
+
+
+
+
+
+
+```octave
+x(imin-1)
+x(imin)
+x(imin+1)
+```
+
+ ans = 3.2824
+ ans = 3.2985
+ ans = 3.3146
+
+
+## Golden Search Algorithm
+
+We can't just look for a sign change for the problem (unless we can take a derivative) so we need a new approach to determine whether we have a maximum between the two bounds.
+
+Rather than using the midpoint of initial bounds, here the problem is more difficult. We need to compare the values of 4 function evaluations. The golden search uses the golden ratio to determine two interior points.
+
+
+
+Start with bounds of 2.5 and 6 Angstrom.
+
+
+```octave
+% define Au atomic potential
+epsilon = 0.039; % kcal/mol
+sigma = 2.934; % Angstrom
+Au_x= @(x) lennard_jones(x,sigma,epsilon);
+
+% calculate golden ratio
+phi = 1/2+sqrt(5)/2;
+% set initial limits
+x_l=2.8;
+x_u=6;
+
+% Iteration #1
+d=(phi-1)*(x_u-x_l);
+
+x1=x_l+d; % define point 1
+x2=x_u-d; % define point 2
+
+
+% evaluate Au_x(x1) and Au_x(x2)
+
+f1=Au_x(x1);
+f2=Au_x(x2);
+plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')
+hold on;
+
+if f2x2
+ plot(x4,f4,'*',[x1,x2],[f1,f2])
+ x1=x2;
+ f1=f2;
+ else
+ plot(x4,f4,'*',[x3,x2],[f3,f2])
+ x3=x2;
+ f3=f2;
+ end
+ x2=x4; f2=f4;
+else
+ error('no minimum in bracket')
+end
+hold off
+```
+
+
+
+
+
+
+```octave
+p = polyfit([x1,x2,x3],[f1,f2,f3],2);
+x_fit = linspace(x1,x3,20);
+y_fit = polyval(p,x_fit);
+
+plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')
+hold on
+if f2x2
+ plot(x4,f4,'*',[x1,x2],[f1,f2])
+ x1=x2;
+ f1=f2;
+ else
+ plot(x4,f4,'*',[x3,x2],[f3,f2])
+ x3=x2;
+ f3=f2;
+ end
+ x2=x4; f2=f4;
+else
+ error('no minimum in bracket')
+end
+hold off
+
+```
+
+
+
+
+
+
+```octave
+p = polyfit([x1,x2,x3],[f1,f2,f3],2);
+x_fit = linspace(x1,x3,20);
+y_fit = polyval(p,x_fit);
+
+plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')
+hold on
+if f2x2
+ plot(x4,f4,'*',[x1,x2],[f1,f2])
+ x1=x2;
+ f1=f2;
+ else
+ plot(x4,f4,'*',[x3,x2],[f3,f2])
+ x3=x2;
+ f3=f2;
+ end
+ x2=x4; f2=f4;
+else
+ error('no minimum in bracket')
+end
+hold off
+
+```
+
+
+
+
+
+Parabolic interpolation does not converge in many scenarios even though it it a bracketing method. Instead, functions like `fminbnd` in Matlab and Octave use a combination of the two (Golden Ratio and Parabolic)
+
+## Using the solutions to minimization for the nonlinear spring constant
+
+Now, we have two routines to find minimums of a univariate function (Golden Ratio and Parabolic). Let's use these to solve for the minimum energy given a range of applied forces to the single atom gold chain
+
+$E_{total}(\Delta x) = E_{LJ}(x_{min}+\Delta x) - F \cdot \Delta x$
+
+
+```octave
+epsilon = 0.039; % kcal/mol
+epsilon = epsilon*6.9477e-21; % J/atom
+epsilon = epsilon*1e18; % aJ/J
+% final units for epsilon are aJ
+
+sigma = 2.934; % Angstrom
+sigma = sigma*0.10; % nm/Angstrom
+x=linspace(2.8,6,200)*0.10; % bond length in um
+
+Ex = lennard_jones(x,sigma,epsilon);
+
+%[Emin,imin]=min(Ex);
+
+[xmin,Emin] = fminbnd(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)
+
+plot(x,Ex,xmin,Emin,'o')
+ylabel('Lennard Jones Potential (aJ/atom)')
+xlabel('bond length (nm)')
+
+Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;
+
+```
+
+ xmin = 0.32933
+ Emin = -2.7096e-04
+
+
+
+
+
+
+
+```octave
+N=50;
+dx = zeros(1,N); % [in nm]
+F_applied=linspace(0,0.0022,N); % [in nN]
+for i=1:N
+ optmin=goldmin(@(dx) Etotal(dx,F_applied(i)),-0.001,0.035);
+ dx(i)=optmin;
+end
+
+plot(dx,F_applied)
+xlabel('dx (nm)')
+ylabel('Force (nN)')
+```
+
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: called from
+ goldmin at line 17 column 1
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+ warning: Matlab-style short-circuit operation performed for operator |
+
+
+
+
+
+
+For this function, it is possible to take a derivative and compare the analytical result:
+
+
+```octave
+dx_full=linspace(0,0.06,50);
+F= @(dx) 4*epsilon*6*(sigma^6./(xmin+dx).^7-2*sigma^12./(xmin+dx).^13)
+plot(dx_full,F(dx_full),dx,F_applied)
+```
+
+ F =
+
+ @(dx) 4 * epsilon * 6 * (sigma ^ 6 ./ (xmin + dx) .^ 7 - 2 * sigma ^ 12 ./ (xmin + dx) .^ 13)
+
+
+
+
+
+
+
+## Curve-fitting
+Another example is minimizing error in your approximation of a function. If you have data (now we have Force-displacement data) we can fit this to a function, such as:
+
+$F(x) = K_{1}\Delta x + \frac{1}{2} K_{2}(\Delta x)^{2}$
+
+
+
+```octave
+function SSE = sse_of_parabola(K,xdata,ydata)
+ % calculate the sum of squares error for a parabola given a function, func, and xdata and ydata
+ % output is SSE=sum of squares error
+ K1=K(1);
+ K2=K(2);
+ y_function = K1*xdata+1/2*K2*xdata.^2;
+ SSE = sum((ydata-y_function).^2);
+end
+
+```
+
+
+```octave
+[K,SSE_min]=fminsearch(@(K) sse_of_parabola(K,dx,F_applied),[1,1]);
+fprintf('\nThe nonlinear spring constants are K1=%1.2f nN/nm and K2=%1.2f nN/nm^2\n',K)
+fprintf('The mininum sum of squares error = %1.2e',SSE_min)
+```
+
+
+ The nonlinear spring constants are K1=0.16 nN/nm and K2=-5.98 nN/nm^2
+ The mininum sum of squares error = 7.35e-08
+
+
+
+```octave
+plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2)
+```
+
+
+
+
diff --git a/lecture_08/lecture_08.pdf b/lecture_08/lecture_08.pdf
new file mode 100644
index 0000000..c251694
Binary files /dev/null and b/lecture_08/lecture_08.pdf differ
diff --git a/lecture_08/lecture_08_files/lecture_08_11_0.svg b/lecture_08/lecture_08_files/lecture_08_11_0.svg
new file mode 100644
index 0000000..c61f790
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_11_0.svg
@@ -0,0 +1,168 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_12_0.svg b/lecture_08/lecture_08_files/lecture_08_12_0.svg
new file mode 100644
index 0000000..b641032
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_12_0.svg
@@ -0,0 +1,163 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_13_0.svg b/lecture_08/lecture_08_files/lecture_08_13_0.svg
new file mode 100644
index 0000000..3e467e5
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_13_0.svg
@@ -0,0 +1,163 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_16_1.svg b/lecture_08/lecture_08_files/lecture_08_16_1.svg
new file mode 100644
index 0000000..ff27c29
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_16_1.svg
@@ -0,0 +1,142 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_17_1.svg b/lecture_08/lecture_08_files/lecture_08_17_1.svg
new file mode 100644
index 0000000..1f3103b
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_17_1.svg
@@ -0,0 +1,136 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_19_1.svg b/lecture_08/lecture_08_files/lecture_08_19_1.svg
new file mode 100644
index 0000000..2a891df
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_19_1.svg
@@ -0,0 +1,135 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_23_0.svg b/lecture_08/lecture_08_files/lecture_08_23_0.svg
new file mode 100644
index 0000000..8fa3bb1
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_23_0.svg
@@ -0,0 +1,196 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_3_0.svg b/lecture_08/lecture_08_files/lecture_08_3_0.svg
new file mode 100644
index 0000000..2c38150
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_3_0.svg
@@ -0,0 +1,147 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_6_1.svg b/lecture_08/lecture_08_files/lecture_08_6_1.svg
new file mode 100644
index 0000000..ecf7a18
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_6_1.svg
@@ -0,0 +1,169 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_7_1.svg b/lecture_08/lecture_08_files/lecture_08_7_1.svg
new file mode 100644
index 0000000..e63011c
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_7_1.svg
@@ -0,0 +1,169 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_8_1.svg b/lecture_08/lecture_08_files/lecture_08_8_1.svg
new file mode 100644
index 0000000..873b04d
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_8_1.svg
@@ -0,0 +1,169 @@
+
\ No newline at end of file
diff --git a/lecture_08/lecture_08_files/lecture_08_9_1.svg b/lecture_08/lecture_08_files/lecture_08_9_1.svg
new file mode 100644
index 0000000..28fd681
--- /dev/null
+++ b/lecture_08/lecture_08_files/lecture_08_9_1.svg
@@ -0,0 +1,169 @@
+
\ No newline at end of file
diff --git a/lecture_08/lennard_jones.m b/lecture_08/lennard_jones.m
new file mode 100644
index 0000000..d18a6ad
--- /dev/null
+++ b/lecture_08/lennard_jones.m
@@ -0,0 +1,4 @@
+function E_LJ =lennard_jones(x,sigma,epsilon)
+ E_LJ = 4*epsilon*((sigma./x).^12-(sigma./x).^6);
+end
+
diff --git a/lecture_08/octave-workspace b/lecture_08/octave-workspace
new file mode 100644
index 0000000..8c437bb
Binary files /dev/null and b/lecture_08/octave-workspace differ