diff --git a/HW3/README.md b/HW3/README.md new file mode 100644 index 0000000..494e0d5 --- /dev/null +++ b/HW3/README.md @@ -0,0 +1,74 @@ +# 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`, `incsearch.m`, `newtraph.m` and `mod_secant.m` + to solve for the angle needed to reach h=1.72 m, with an initial speed of 1.5 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 | | | | + |incsearch | | | | + |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: + + `![Plot of convergence for four numerical solvers.](convergence.png)` + + 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. + + 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/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": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t250\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t350\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t400\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "iterations = linspace(1,400,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f_m,df_m,200,0,iterations(i));\n", + " [root_ms,ea_ms(i),iter_ms]=mod_secant(f_m,1e-6,300,0,iterations(i));\n", + " [root_fp,ea_fp(i),iter_fp]=falsepos(f_m,1,300,0,iterations(i));\n", + " [root_bs,ea_bs(i),iter_bs]=bisect(f_m,1,300,0,iterations(i));\n", + "end\n", + " \n", + "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n", + "legend('newton-raphson','mod-secant','false point','bisection')" + ] + }, + { + "cell_type": "code", + "execution_count": 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": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "f= @(x) x^10-1;\n", + "df=@(x) 10*x^9;\n", + "iterations = linspace(1,50,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i));\n", + " [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i));\n", + " [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i));\n", + " [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i));\n", + "end\n", + " \n", + "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n", + "legend('newton-raphson','mod-secant','false point','bisection')" + ] + }, + { + "cell_type": "code", + "execution_count": 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/.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": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t25000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tprinciple amount left ($)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\ttime (years)\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "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": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t400000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t500000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t600000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t700000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t25\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tprinciple amount left ($)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\ttime (years)\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "Amt_numerical=mod_secant(@(A) car_payments(A,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": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t250\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t350\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t400\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "iterations = linspace(1,400,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f_m,df_m,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": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-14\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-12\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t102\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t104\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t30\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tnewton-raphson\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tnewton-raphson\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmod-secant\n", + "\n", + "\t\n", + "\t\tmod-secant\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tfalse point\n", + "\n", + "\t\n", + "\t\tfalse point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbisection\n", + "\n", + "\t\n", + "\t\tbisection\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=20;\n", + "f= @(x) x^10-1;\n", + "df=@(x) 10*x^9;\n", + "iterations = linspace(1,50,N);\n", + "ea_nr=zeros(1,N); % appr error Newton-Raphson\n", + "ea_ms=zeros(1,N); % appr error Modified Secant\n", + "ea_fp=zeros(1,N); % appr error false point method\n", + "ea_bs=zeros(1,N); % appr error bisect method\n", + "for i=1:length(iterations)\n", + " [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i));\n", + " [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i));\n", + " [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i));\n", + " [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i));\n", + "end\n", + " \n", + "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n", + "legend('newton-raphson','mod-secant','false point','bisection')" + ] + }, + { + "cell_type": "code", + "execution_count": 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 + + + +![svg](lecture_07_files/lecture_07_15_1.svg) + + + +```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 + + + +![svg](lecture_07_files/lecture_07_16_1.svg) + + + +```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 + + + +![svg](lecture_07_files/lecture_07_24_1.svg) + + + +```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 + + + +![svg](lecture_07_files/lecture_07_26_1.svg) + + + +```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 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -5000 + + + + + 0 + + + + + 5000 + + + + + 10000 + + + + + 15000 + + + + + 20000 + + + + + 25000 + + + + + 30000 + + + + + 0 + + + + + 1 + + + + + 2 + + + + + 3 + + + + + 4 + + + + + 5 + + + + + + + + + principle amount left ($) + + + + + time (years) + + + + + gnuplot_plot_1a + + + + + + + + + + + + + + + \ 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 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 10000 + + + + + 15000 + + + + + 20000 + + + + + 25000 + + + + + 30000 + + + + + 0 + + + + + 1 + + + + + 2 + + + + + 3 + + + + + 4 + + + + + 5 + + + + + + + + + principle amount left ($) + + + + + time (years) + + + + + gnuplot_plot_1a + + + + + + + + + + + + + + + \ 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 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + 100000 + + + + + 200000 + + + + + 300000 + + + + + 400000 + + + + + 500000 + + + + + 600000 + + + + + 700000 + + + + + 0 + + + + + 5 + + + + + 10 + + + + + 15 + + + + + 20 + + + + + 25 + + + + + 30 + + + + + + + + + principle amount left ($) + + + + + time (years) + + + + + gnuplot_plot_1a + + + + + + + + + + + + + + + \ 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 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 10-14 + + + + + 10-12 + + + + + 10-10 + + + + + 10-8 + + + + + 10-6 + + + + + 10-4 + + + + + 10-2 + + + + + 100 + + + + + 102 + + + + + 104 + + + + + 0 + + + + + 50 + + + + + 100 + + + + + 150 + + + + + 200 + + + + + 250 + + + + + 300 + + + + + 350 + + + + + 400 + + + + + + + + + + + + + newton-raphson + + + + + newton-raphson + + + + + + mod-secant + + + mod-secant + + + + + + false point + + + false point + + + + + + bisection + + + bisection + + + + + + + + + + + + + + + \ 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 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 10-14 + + + + + 10-12 + + + + + 10-10 + + + + + 10-8 + + + + + 10-6 + + + + + 10-4 + + + + + 10-2 + + + + + 100 + + + + + 102 + + + + + 104 + + + + + 0 + + + + + 50 + + + + + 100 + + + + + 150 + + + + + 200 + + + + + 250 + + + + + 300 + + + + + 350 + + + + + 400 + + + + + + + + + + + + + newton-raphson + + + + + newton-raphson + + + + + + mod-secant + + + mod-secant + + + + + + false point + + + false point + + + + + + bisection + + + bisection + + + + + + + + + + + + + + + \ 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 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 10-14 + + + + + 10-12 + + + + + 10-10 + + + + + 10-8 + + + + + 10-6 + + + + + 10-4 + + + + + 10-2 + + + + + 100 + + + + + 102 + + + + + 104 + + + + + 0 + + + + + 10 + + + + + 20 + + + + + 30 + + + + + 40 + + + + + 50 + + + + + + + + + + + + + newton-raphson + + + + + newton-raphson + + + + + + mod-secant + + + mod-secant + + + + + + false point + + + false point + + + + + + bisection + + + bisection + + + + + + + + + + + + + + + \ 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