diff --git a/lecture_07/bisect.m b/lecture_07/bisect.m
new file mode 100644
index 0000000..9f696a0
--- /dev/null
+++ b/lecture_07/bisect.m
@@ -0,0 +1,37 @@
+function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
+% bisect: root location zeroes
+% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
+% uses bisection method to find the root of func
+% input:
+% func = name of function
+% xl, xu = lower and upper guesses
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by func
+% output:
+% root = real root
+% fx = function value at root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if test>0,error('no sign change'),end
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+ xrold = xr;
+ xr = (xl + xu)/2;
+ iter = iter + 1;
+ if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+ test = func(xl,varargin{:})*func(xr,varargin{:});
+ if test < 0
+ xu = xr;
+ elseif test > 0
+ xl = xr;
+ else
+ ea = 0;
+ end
+ if ea <= es || iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});
diff --git a/lecture_07/car_payments.m b/lecture_07/car_payments.m
new file mode 100644
index 0000000..9d5b2a5
--- /dev/null
+++ b/lecture_07/car_payments.m
@@ -0,0 +1,17 @@
+function amount_left = car_payments(monthly_payment,price,apr,no_of_years,plot_bool)
+ interest_per_month = apr/12;
+ number_of_months = no_of_years*12;
+ principle=price;
+ P_vector=zeros(1,number_of_months);
+ for i = 1:number_of_months
+ principle=principle-monthly_payment;
+ principle=(1+interest_per_month)*principle;
+ P_vector(i)=principle;
+ end
+ amount_left=principle;
+ if plot_bool
+ plot([1:number_of_months]/12, P_vector)
+ xlabel('time (years)')
+ ylabel('principle amount left ($)')
+ end
+end
diff --git a/lecture_07/falsepos.m b/lecture_07/falsepos.m
new file mode 100644
index 0000000..d5575d5
--- /dev/null
+++ b/lecture_07/falsepos.m
@@ -0,0 +1,39 @@
+function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)
+% bisect: root location zeroes
+% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
+% uses bisection method to find the root of func
+% input:
+% func = name of function
+% xl, xu = lower and upper guesses
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by func
+% output:
+% root = real root
+% fx = function value at root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if test>0,error('no sign change'),end
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+ xrold = xr;
+ xr = (xl + xu)/2;
+ % xr = (xl + xu)/2; % bisect method
+ xr=xu - (func(xu)*(xl-xu))/(func(xl)-func(xu)); % false position method
+ iter = iter + 1;
+ if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+ test = func(xl,varargin{:})*func(xr,varargin{:});
+ if test < 0
+ xu = xr;
+ elseif test > 0
+ xl = xr;
+ else
+ ea = 0;
+ end
+ if ea <= es || iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});
diff --git a/lecture_07/fzerosimp.m b/lecture_07/fzerosimp.m
new file mode 100644
index 0000000..05c7a9b
--- /dev/null
+++ b/lecture_07/fzerosimp.m
@@ -0,0 +1,41 @@
+function b = fzerosimp(xl,xu)
+a = xl; b = xu; fa = f(a); fb = f(b);
+c = a; fc = fa; d = b - c; e = d;
+while (1)
+ if fb == 0, break, end
+ if sign(fa) == sign(fb) %If needed, rearrange points
+ a = c; fa = fc; d = b - c; e = d;
+ end
+ if abs(fa) < abs(fb)
+ c = b; b = a; a = c;
+ fc = fb; fb = fa; fa = fc;
+ end
+ m = 0.5*(a - b); %Termination test and possible exit
+ tol = 2 * eps * max(abs(b), 1);
+ if abs(m) <= tol | fb == 0.
+ break
+ end
+ %Choose open methods or bisection
+ if abs(e) >= tol & abs(fc) > abs(fb)
+ s = fb/fc;
+ if a == c %Secant method
+ p = 2*m*s;
+ q = 1 - s;
+ else %Inverse quadratic interpolation
+ q = fc/fa; r = fb/fa;
+ p = s * (2*m*q * (q - r) - (b - c)*(r - 1));
+ q = (q - 1)*(r - 1)*(s - 1);
+ end
+ if p > 0, q = -q; else p = -p; end;
+ if 2*p < 3*m*q - abs(tol*q) & p < abs(0.5*e*q)
+ e = d; d = p/q;
+ else
+ d = m; e = m;
+ end
+ else %Bisection
+ d = m; e = m;
+ end
+ c = b; fc = fb;
+ if abs(d) > tol, b=b+d; else b=b-sign(b-a)*tol; end
+ fb = f(b);
+end
\ No newline at end of file
diff --git a/lecture_07/incsearch.m b/lecture_07/incsearch.m
new file mode 100644
index 0000000..bd82554
--- /dev/null
+++ b/lecture_07/incsearch.m
@@ -0,0 +1,37 @@
+function xb = incsearch(func,xmin,xmax,ns)
+% incsearch: incremental search root locator
+% xb = incsearch(func,xmin,xmax,ns):
+% finds brackets of x that contain sign changes
+% of a function on an interval
+% input:
+% func = name of function
+% xmin, xmax = endpoints of interval
+% ns = number of subintervals (default = 50)
+% output:
+% xb(k,1) is the lower bound of the kth sign change
+% xb(k,2) is the upper bound of the kth sign change
+% If no brackets found, xb = [].
+if nargin < 3, error('at least 3 arguments required'), end
+if nargin < 4, ns = 50; end %if ns blank set to 50
+% Incremental search
+x = linspace(xmin,xmax,ns);
+f = func(x);
+nb = 0; xb = []; %xb is null unless sign change detected
+%for k = 1:length(x)-1
+% if sign(f(k)) ~= sign(f(k+1)) %check for sign change
+% nb = nb + 1;
+% xb(nb,1) = x(k);
+% xb(nb,2) = x(k+1);
+% end
+%end
+sign_change = diff(sign(f));
+[~,i_change] = find(sign_change~=0);
+nb=length(i_change);
+xb=[x(i_change)',x(i_change+1)'];
+
+if isempty(xb) %display that no brackets were found
+ fprintf('no brackets found\n')
+ fprintf('check interval or increase ns\n')
+else
+ fprintf('number of brackets: %i\n',nb) %display number of brackets
+end
diff --git a/lecture_07/lecture_07.ipynb b/lecture_07/lecture_07.ipynb
new file mode 100644
index 0000000..0c3c2ca
--- /dev/null
+++ b/lecture_07/lecture_07.ipynb
@@ -0,0 +1,1195 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 30,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Roots: Open methods\n",
+ "## Newton-Raphson"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First-order approximation for the location of the root (i.e. assume the slope at the given point is constant, what is the solution when f(x)=0)\n",
+ "\n",
+ "$f'(x_{i})=\\frac{f(x_{i})-0}{x_{i}-x_{i+1}}$\n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n",
+ "\n",
+ "Use Newton-Raphson to find solution when $e^{x}=x$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "error: 'c' undefined near line 1 column 1\n",
+ "error: 'x_r' undefined near line 1 column 21\n",
+ "error: evaluating argument list element number 1\n",
+ "error: 'x_r' undefined near line 1 column 5\n"
+ ]
+ }
+ ],
+ "source": [
+ "f= @(x) exp(-x)-x;\n",
+ "df= @(x) -exp(-x)-1;\n",
+ "\n",
+ "x_i= 0;\n",
+ "c\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 0.50000\n",
+ "error_approx = 1\n"
+ ]
+ }
+ ],
+ "source": [
+ "x_r = x_i-f(x_i)/df(x_i)\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 0.56631\n",
+ "error_approx = 0.11709\n"
+ ]
+ }
+ ],
+ "source": [
+ "x_r = x_i-f(x_i)/df(x_i)\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 0.56714\n",
+ "error_approx = 0.0014673\n"
+ ]
+ }
+ ],
+ "source": [
+ "x_r = x_i-f(x_i)/df(x_i)\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In the bungee jumper example, we created a function f(m) that when f(m)=0, then the mass had been chosen such that at t=4 s, the velocity is 36 m/s. \n",
+ "\n",
+ "$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$.\n",
+ "\n",
+ "to use the Newton-Raphson method, we need the derivative $\\frac{df}{dm}$\n",
+ "\n",
+ "$\\frac{df}{dm}=\\frac{1}{2}\\sqrt{\\frac{g}{mc_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-\n",
+ "\\frac{g}{2m}\\mathrm{sech}^{2}(\\sqrt{\\frac{gc_{d}}{m}}t)$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults\n",
+ "g=9.81; % acceleration due to gravity\n",
+ "m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n",
+ "c_d=0.25; % drag coefficient\n",
+ "t=4; % at time = 4 seconds\n",
+ "v=36; % speed must be 36 m/s\n",
+ "f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n",
+ "df_m = @(m) 1/2*sqrt(g./m/c_d).*tanh(sqrt(g*c_d./m)*t)-g/2./m*sech(sqrt(g*c_d./m)*t).^2;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 142.74\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "newtraph(f_m,df_m,140,0.00001)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Secant Methods"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Not always able to evaluate the derivative. Approximation of derivative:\n",
+ "\n",
+ "$f'(x_{i})=\\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}$\n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{\\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}}=\n",
+ " x_{i}-\\frac{f(x_{i})(x_{i-1}-x_{i})}{f(x_{i-1})-f(x_{i})}$\n",
+ " \n",
+ "What values should $x_{i}$ and $x_{i-1}$ take?\n",
+ "\n",
+ "To reduce arbitrary selection of variables, use the\n",
+ "\n",
+ "## Modified Secant method\n",
+ "\n",
+ "Change the x evaluations to a perturbation $\\delta$. \n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})(\\delta x_{i})}{f(x_{i}+\\delta x_{i})-f(x_{i})}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 142.74\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "mod_secant(f_m,1e-6,50,0.00001)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 32,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Amt_numerical = 563.79\n",
+ "ans = -160.42\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "Amt_numerical=mod_secant(@(A) car_payments(A,30000,0.05,5),1e-6,50,0.001)\n",
+ "car_payments(Amt,30000,0.05,5)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 3.3968e+04\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "Amt*12*5"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Amortization calculation makes the same calculation for the monthly payment amount, A, paying off the principle amount, P, over n pay periods with monthly interest rate, r. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Amt = 566.14\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "% Amortization calculation\n",
+ "A = @(P,r,n) P*(r*(1+r)^n)./((1+r)^n-1);\n",
+ "Amt=A(30000,0.05/12,5*12)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Matlab's function\n",
+ "\n",
+ "Matlab and Octave combine bracketing and open methods in the `fzero` function. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 33,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'fzero' is a function from the file /usr/share/octave/4.0.0/m/optimization/fzero.m\n",
+ "\n",
+ " -- Function File: fzero (FUN, X0)\n",
+ " -- Function File: fzero (FUN, X0, OPTIONS)\n",
+ " -- Function File: [X, FVAL, INFO, OUTPUT] = fzero (...)\n",
+ " Find a zero of a univariate function.\n",
+ "\n",
+ " FUN is a function handle, inline function, or string containing the\n",
+ " name of the function to evaluate.\n",
+ "\n",
+ " X0 should be a two-element vector specifying two points which\n",
+ " bracket a zero. In other words, there must be a change in sign of\n",
+ " the function between X0(1) and X0(2). More mathematically, the\n",
+ " following must hold\n",
+ "\n",
+ " sign (FUN(X0(1))) * sign (FUN(X0(2))) <= 0\n",
+ "\n",
+ " If X0 is a single scalar then several nearby and distant values are\n",
+ " probed in an attempt to obtain a valid bracketing. If this is not\n",
+ " successful, the function fails.\n",
+ "\n",
+ " OPTIONS is a structure specifying additional options. Currently,\n",
+ " 'fzero' recognizes these options: \"FunValCheck\", \"OutputFcn\",\n",
+ " \"TolX\", \"MaxIter\", \"MaxFunEvals\". For a description of these\n",
+ " options, see *note optimset: XREFoptimset.\n",
+ "\n",
+ " On exit, the function returns X, the approximate zero point and\n",
+ " FVAL, the function value thereof.\n",
+ "\n",
+ " INFO is an exit flag that can have these values:\n",
+ "\n",
+ " * 1 The algorithm converged to a solution.\n",
+ "\n",
+ " * 0 Maximum number of iterations or function evaluations has\n",
+ " been reached.\n",
+ "\n",
+ " * -1 The algorithm has been terminated from user output\n",
+ " function.\n",
+ "\n",
+ " * -5 The algorithm may have converged to a singular point.\n",
+ "\n",
+ " OUTPUT is a structure containing runtime information about the\n",
+ " 'fzero' algorithm. Fields in the structure are:\n",
+ "\n",
+ " * iterations Number of iterations through loop.\n",
+ "\n",
+ " * nfev Number of function evaluations.\n",
+ "\n",
+ " * bracketx A two-element vector with the final bracketing of the\n",
+ " zero along the x-axis.\n",
+ "\n",
+ " * brackety A two-element vector with the final bracketing of the\n",
+ " zero along the y-axis.\n",
+ "\n",
+ " See also: optimset, fsolve.\n",
+ "\n",
+ "Additional help for built-in functions and operators is\n",
+ "available in the online version of the manual. Use the command\n",
+ "'doc ' to search the manual index.\n",
+ "\n",
+ "Help and information about Octave is also available on the WWW\n",
+ "at http://www.octave.org and via the help@octave.org\n",
+ "mailing list.\n"
+ ]
+ }
+ ],
+ "source": [
+ "help fzero"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 40,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 563.79\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "fzero(@(A) car_payments(A,30000,0.05,5,0),500)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Comparison of Solvers\n",
+ "\n",
+ "It's helpful to compare to the convergence of different routines to see how quickly you find a solution. \n",
+ "\n",
+ "Comparing the freefall example\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 84,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: called from\n",
+ " __line__ at line 120 column 16\n",
+ " line at line 56 column 8\n",
+ " __plt__>__plt2vv__ at line 500 column 10\n",
+ " __plt__>__plt2__ at line 246 column 14\n",
+ " __plt__ at line 133 column 15\n",
+ " semilogy at line 60 column 10\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=20;\n",
+ "iterations = linspace(1,400,N);\n",
+ "ea_nr=zeros(1,N); % appr error Newton-Raphson\n",
+ "ea_ms=zeros(1,N); % appr error Modified Secant\n",
+ "ea_fp=zeros(1,N); % appr error false point method\n",
+ "ea_bs=zeros(1,N); % appr error bisect method\n",
+ "for i=1:length(iterations)\n",
+ " [root_nr,ea_nr(i),iter_nr]=newtraph(f_m,df_m,200,0,iterations(i));\n",
+ " [root_ms,ea_ms(i),iter_ms]=mod_secant(f_m,1e-6,300,0,iterations(i));\n",
+ " [root_fp,ea_fp(i),iter_fp]=falsepos(f_m,1,300,0,iterations(i));\n",
+ " [root_bs,ea_bs(i),iter_bs]=bisect(f_m,1,300,0,iterations(i));\n",
+ "end\n",
+ " \n",
+ "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n",
+ "legend('newton-raphson','mod-secant','false point','bisection')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 75,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ea_ms =\n",
+ "\n",
+ " Columns 1 through 7:\n",
+ "\n",
+ " 43.43883 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n",
+ "\n",
+ " Columns 8 through 14:\n",
+ "\n",
+ " 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n",
+ "\n",
+ " Columns 15 through 20:\n",
+ "\n",
+ " 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "ea_ms"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 101,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: called from\n",
+ " __line__ at line 120 column 16\n",
+ " line at line 56 column 8\n",
+ " __plt__>__plt2vv__ at line 500 column 10\n",
+ " __plt__>__plt2__ at line 246 column 14\n",
+ " __plt__ at line 133 column 15\n",
+ " semilogy at line 60 column 10\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=20;\n",
+ "f= @(x) x^10-1;\n",
+ "df=@(x) 10*x^9;\n",
+ "iterations = linspace(1,50,N);\n",
+ "ea_nr=zeros(1,N); % appr error Newton-Raphson\n",
+ "ea_ms=zeros(1,N); % appr error Modified Secant\n",
+ "ea_fp=zeros(1,N); % appr error false point method\n",
+ "ea_bs=zeros(1,N); % appr error bisect method\n",
+ "for i=1:length(iterations)\n",
+ " [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i));\n",
+ " [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i));\n",
+ " [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i));\n",
+ " [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i));\n",
+ "end\n",
+ " \n",
+ "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n",
+ "legend('newton-raphson','mod-secant','false point','bisection')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 102,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ea_bs =\n",
+ "\n",
+ " Columns 1 through 6:\n",
+ "\n",
+ " 9.5357e+03 -4.7554e-01 -2.1114e-01 6.0163e-02 -2.4387e-03 6.1052e-04\n",
+ "\n",
+ " Columns 7 through 12:\n",
+ "\n",
+ " 2.2891e-04 -9.5367e-06 2.3842e-06 8.9407e-07 -2.2352e-07 9.3132e-09\n",
+ "\n",
+ " Columns 13 through 18:\n",
+ "\n",
+ " -2.3283e-09 -8.7311e-10 3.6380e-11 -9.0949e-12 -3.4106e-12 8.5265e-13\n",
+ "\n",
+ " Columns 19 and 20:\n",
+ "\n",
+ " -3.5527e-14 8.8818e-15\n",
+ "\n",
+ "ans = 16.208\n"
+ ]
+ }
+ ],
+ "source": [
+ "ea_bs\n",
+ "newtraph(f,df,0.5,0,12)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 93,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 1.9683e+23\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "df(300)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_07/lecture_07.md b/lecture_07/lecture_07.md
new file mode 100644
index 0000000..020699b
--- /dev/null
+++ b/lecture_07/lecture_07.md
@@ -0,0 +1,388 @@
+
+
+```octave
+%plot --format svg
+```
+
+
+```octave
+setdefaults
+```
+
+# Roots: Open methods
+## Newton-Raphson
+
+First-order approximation for the location of the root (i.e. assume the slope at the given point is constant, what is the solution when f(x)=0)
+
+$f'(x_{i})=\frac{f(x_{i})-0}{x_{i}-x_{i+1}}$
+
+$x_{i+1}=x_{i}-\frac{f(x_{i})}{f'(x_{i})}$
+
+Use Newton-Raphson to find solution when $e^{x}=x$
+
+
+```octave
+f= @(x) exp(-x)-x;
+df= @(x) -exp(-x)-1;
+
+x_i= 0;
+c
+error_approx = abs((x_r-x_i)/x_r)
+x_i=x_r;
+
+```
+
+ error: 'c' undefined near line 1 column 1
+ error: 'x_r' undefined near line 1 column 21
+ error: evaluating argument list element number 1
+ error: 'x_r' undefined near line 1 column 5
+
+
+
+```octave
+x_r = x_i-f(x_i)/df(x_i)
+error_approx = abs((x_r-x_i)/x_r)
+x_i=x_r;
+```
+
+ x_r = 0.50000
+ error_approx = 1
+
+
+
+```octave
+x_r = x_i-f(x_i)/df(x_i)
+error_approx = abs((x_r-x_i)/x_r)
+x_i=x_r;
+```
+
+ x_r = 0.56631
+ error_approx = 0.11709
+
+
+
+```octave
+x_r = x_i-f(x_i)/df(x_i)
+error_approx = abs((x_r-x_i)/x_r)
+x_i=x_r;
+```
+
+ x_r = 0.56714
+ error_approx = 0.0014673
+
+
+In the bungee jumper example, we created a function f(m) that when f(m)=0, then the mass had been chosen such that at t=4 s, the velocity is 36 m/s.
+
+$f(m)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)-v(t)$.
+
+to use the Newton-Raphson method, we need the derivative $\frac{df}{dm}$
+
+$\frac{df}{dm}=\frac{1}{2}\sqrt{\frac{g}{mc_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)-
+\frac{g}{2m}\mathrm{sech}^{2}(\sqrt{\frac{gc_{d}}{m}}t)$
+
+
+```octave
+setdefaults
+g=9.81; % acceleration due to gravity
+m=linspace(50, 200,100); % possible values for mass 50 to 200 kg
+c_d=0.25; % drag coefficient
+t=4; % at time = 4 seconds
+v=36; % speed must be 36 m/s
+f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m
+df_m = @(m) 1/2*sqrt(g./m/c_d).*tanh(sqrt(g*c_d./m)*t)-g/2./m*sech(sqrt(g*c_d./m)*t).^2;
+```
+
+
+```octave
+newtraph(f_m,df_m,140,0.00001)
+```
+
+ ans = 142.74
+
+
+## Secant Methods
+
+Not always able to evaluate the derivative. Approximation of derivative:
+
+$f'(x_{i})=\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}$
+
+$x_{i+1}=x_{i}-\frac{f(x_{i})}{f'(x_{i})}$
+
+$x_{i+1}=x_{i}-\frac{f(x_{i})}{\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}}=
+ x_{i}-\frac{f(x_{i})(x_{i-1}-x_{i})}{f(x_{i-1})-f(x_{i})}$
+
+What values should $x_{i}$ and $x_{i-1}$ take?
+
+To reduce arbitrary selection of variables, use the
+
+## Modified Secant method
+
+Change the x evaluations to a perturbation $\delta$.
+
+$x_{i+1}=x_{i}-\frac{f(x_{i})(\delta x_{i})}{f(x_{i}+\delta x_{i})-f(x_{i})}$
+
+
+```octave
+mod_secant(f_m,1e-6,50,0.00001)
+```
+
+ ans = 142.74
+
+
+
+```octave
+Amt_numerical=mod_secant(@(A) car_payments(A,30000,0.05,5),1e-6,50,0.001)
+car_payments(Amt,30000,0.05,5)
+```
+
+ Amt_numerical = 563.79
+ ans = -160.42
+
+
+
+
+
+
+
+```octave
+Amt*12*5
+```
+
+ ans = 3.3968e+04
+
+
+Amortization calculation makes the same calculation for the monthly payment amount, A, paying off the principle amount, P, over n pay periods with monthly interest rate, r.
+
+
+```octave
+% Amortization calculation
+A = @(P,r,n) P*(r*(1+r)^n)./((1+r)^n-1);
+Amt=A(30000,0.05/12,5*12)
+```
+
+ Amt = 566.14
+
+
+## Matlab's function
+
+Matlab and Octave combine bracketing and open methods in the `fzero` function.
+
+
+```octave
+help fzero
+```
+
+ 'fzero' is a function from the file /usr/share/octave/4.0.0/m/optimization/fzero.m
+
+ -- Function File: fzero (FUN, X0)
+ -- Function File: fzero (FUN, X0, OPTIONS)
+ -- Function File: [X, FVAL, INFO, OUTPUT] = fzero (...)
+ Find a zero of a univariate function.
+
+ FUN is a function handle, inline function, or string containing the
+ name of the function to evaluate.
+
+ X0 should be a two-element vector specifying two points which
+ bracket a zero. In other words, there must be a change in sign of
+ the function between X0(1) and X0(2). More mathematically, the
+ following must hold
+
+ sign (FUN(X0(1))) * sign (FUN(X0(2))) <= 0
+
+ If X0 is a single scalar then several nearby and distant values are
+ probed in an attempt to obtain a valid bracketing. If this is not
+ successful, the function fails.
+
+ OPTIONS is a structure specifying additional options. Currently,
+ 'fzero' recognizes these options: "FunValCheck", "OutputFcn",
+ "TolX", "MaxIter", "MaxFunEvals". For a description of these
+ options, see *note optimset: XREFoptimset.
+
+ On exit, the function returns X, the approximate zero point and
+ FVAL, the function value thereof.
+
+ INFO is an exit flag that can have these values:
+
+ * 1 The algorithm converged to a solution.
+
+ * 0 Maximum number of iterations or function evaluations has
+ been reached.
+
+ * -1 The algorithm has been terminated from user output
+ function.
+
+ * -5 The algorithm may have converged to a singular point.
+
+ OUTPUT is a structure containing runtime information about the
+ 'fzero' algorithm. Fields in the structure are:
+
+ * iterations Number of iterations through loop.
+
+ * nfev Number of function evaluations.
+
+ * bracketx A two-element vector with the final bracketing of the
+ zero along the x-axis.
+
+ * brackety A two-element vector with the final bracketing of the
+ zero along the y-axis.
+
+ See also: optimset, fsolve.
+
+ Additional help for built-in functions and operators is
+ available in the online version of the manual. Use the command
+ 'doc ' to search the manual index.
+
+ Help and information about Octave is also available on the WWW
+ at http://www.octave.org and via the help@octave.org
+ mailing list.
+
+
+
+```octave
+fzero(@(A) car_payments(A,30000,0.05,5,0),500)
+```
+
+ ans = 563.79
+
+
+## Comparison of Solvers
+
+It's helpful to compare to the convergence of different routines to see how quickly you find a solution.
+
+Comparing the freefall example
+
+
+
+```octave
+N=20;
+iterations = linspace(1,400,N);
+ea_nr=zeros(1,N); % appr error Newton-Raphson
+ea_ms=zeros(1,N); % appr error Modified Secant
+ea_fp=zeros(1,N); % appr error false point method
+ea_bs=zeros(1,N); % appr error bisect method
+for i=1:length(iterations)
+ [root_nr,ea_nr(i),iter_nr]=newtraph(f_m,df_m,200,0,iterations(i));
+ [root_ms,ea_ms(i),iter_ms]=mod_secant(f_m,1e-6,300,0,iterations(i));
+ [root_fp,ea_fp(i),iter_fp]=falsepos(f_m,1,300,0,iterations(i));
+ [root_bs,ea_bs(i),iter_bs]=bisect(f_m,1,300,0,iterations(i));
+end
+
+semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))
+legend('newton-raphson','mod-secant','false point','bisection')
+```
+
+ warning: axis: omitting non-positive data in log plot
+ warning: called from
+ __line__ at line 120 column 16
+ line at line 56 column 8
+ __plt__>__plt2vv__ at line 500 column 10
+ __plt__>__plt2__ at line 246 column 14
+ __plt__ at line 133 column 15
+ semilogy at line 60 column 10
+ warning: axis: omitting non-positive data in log plot
+ warning: axis: omitting non-positive data in log plot
+ warning: axis: omitting non-positive data in log plot
+
+
+
+
+
+
+
+```octave
+ea_ms
+```
+
+ ea_ms =
+
+ Columns 1 through 7:
+
+ 43.43883 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+
+ Columns 8 through 14:
+
+ 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+
+ Columns 15 through 20:
+
+ 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+
+
+
+
+```octave
+N=20;
+f= @(x) x^10-1;
+df=@(x) 10*x^9;
+iterations = linspace(1,50,N);
+ea_nr=zeros(1,N); % appr error Newton-Raphson
+ea_ms=zeros(1,N); % appr error Modified Secant
+ea_fp=zeros(1,N); % appr error false point method
+ea_bs=zeros(1,N); % appr error bisect method
+for i=1:length(iterations)
+ [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i));
+ [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i));
+ [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i));
+ [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i));
+end
+
+semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))
+legend('newton-raphson','mod-secant','false point','bisection')
+```
+
+ warning: axis: omitting non-positive data in log plot
+ warning: called from
+ __line__ at line 120 column 16
+ line at line 56 column 8
+ __plt__>__plt2vv__ at line 500 column 10
+ __plt__>__plt2__ at line 246 column 14
+ __plt__ at line 133 column 15
+ semilogy at line 60 column 10
+ warning: axis: omitting non-positive data in log plot
+ warning: axis: omitting non-positive data in log plot
+ warning: axis: omitting non-positive data in log plot
+
+
+
+
+
+
+
+```octave
+ea_bs
+newtraph(f,df,0.5,0,12)
+```
+
+ ea_bs =
+
+ Columns 1 through 6:
+
+ 9.5357e+03 -4.7554e-01 -2.1114e-01 6.0163e-02 -2.4387e-03 6.1052e-04
+
+ Columns 7 through 12:
+
+ 2.2891e-04 -9.5367e-06 2.3842e-06 8.9407e-07 -2.2352e-07 9.3132e-09
+
+ Columns 13 through 18:
+
+ -2.3283e-09 -8.7311e-10 3.6380e-11 -9.0949e-12 -3.4106e-12 8.5265e-13
+
+ Columns 19 and 20:
+
+ -3.5527e-14 8.8818e-15
+
+ ans = 16.208
+
+
+
+```octave
+df(300)
+```
+
+ ans = 1.9683e+23
+
+
+
+```octave
+
+```
diff --git a/lecture_07/lecture_07.pdf b/lecture_07/lecture_07.pdf
new file mode 100644
index 0000000..4216be6
Binary files /dev/null and b/lecture_07/lecture_07.pdf differ
diff --git a/lecture_07/lecture_07_files/lecture_07_14_1.svg b/lecture_07/lecture_07_files/lecture_07_14_1.svg
new file mode 100644
index 0000000..9d54798
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_14_1.svg
@@ -0,0 +1,146 @@
+
\ No newline at end of file
diff --git a/lecture_07/lecture_07_files/lecture_07_22_1.svg b/lecture_07/lecture_07_files/lecture_07_22_1.svg
new file mode 100644
index 0000000..502c4c4
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_22_1.svg
@@ -0,0 +1,197 @@
+
\ No newline at end of file
diff --git a/lecture_07/lecture_07_files/lecture_07_24_1.svg b/lecture_07/lecture_07_files/lecture_07_24_1.svg
new file mode 100644
index 0000000..a8b614c
--- /dev/null
+++ b/lecture_07/lecture_07_files/lecture_07_24_1.svg
@@ -0,0 +1,182 @@
+
\ No newline at end of file
diff --git a/lecture_07/mod_secant.m b/lecture_07/mod_secant.m
new file mode 100644
index 0000000..a3caa10
--- /dev/null
+++ b/lecture_07/mod_secant.m
@@ -0,0 +1,31 @@
+function [root,ea,iter]=mod_secant(func,dx,xr,es,maxit,varargin)
+% newtraph: Modified secant root location zeroes
+% [root,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...):
+% uses modified secant method to find the root of func
+% input:
+% func = name of function
+% dx = perturbation fraction
+% xr = initial guess
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by function
+% output:
+% root = real root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+if nargin<4 || isempty(es),es=0.0001;end
+if nargin<5 || isempty(maxit),maxit=50;end
+iter = 0;
+while (1)
+ xrold = xr;
+ dfunc=(func(xr+dx)-func(xr))./dx;
+ xr = xr - func(xr)/dfunc;
+ iter = iter + 1;
+ if xr ~= 0
+ ea = abs((xr - xrold)/xr) * 100;
+ end
+ if ea <= es || iter >= maxit, break, end
+end
+root = xr;
+
diff --git a/lecture_07/newtraph.m b/lecture_07/newtraph.m
new file mode 100644
index 0000000..94ca667
--- /dev/null
+++ b/lecture_07/newtraph.m
@@ -0,0 +1,29 @@
+function [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,varargin)
+% newtraph: Newton-Raphson root location zeroes
+% [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,p1,p2,...):
+% uses Newton-Raphson method to find the root of func
+% input:
+% func = name of function
+% dfunc = name of derivative of function
+% xr = initial guess
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by function
+% output:
+% root = real root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+if nargin<4 || isempty(es),es=0.0001;end
+if nargin<5 || isempty(maxit),maxit=50;end
+iter = 0;
+while (1)
+ xrold = xr;
+ xr = xr - func(xr)/dfunc(xr);
+ iter = iter + 1;
+ if xr ~= 0
+ ea = abs((xr - xrold)/xr) * 100;
+ end
+ if ea <= es || iter >= maxit, break, end
+end
+root = xr;