diff --git a/HW3/README.md b/HW3/README.md new file mode 100644 index 0000000..45214a4 --- /dev/null +++ b/HW3/README.md @@ -0,0 +1,73 @@ +# Homework #3 +## due 2/15/17 by 11:59pm + + +1. Create a new github repository called 'roots_and_optimization'. + + a. Add rcc02007 and pez16103 as collaborators. + + b. Clone the repository to your computer. + + c. Copy your `projectile.m` function into the 'roots_and_optimization' folder. + *Disable the plotting routine for the solvers* + + d. Use the four solvers `falsepos.m`, `bisect.m`, `newtraph.m` and `mod_secant.m` + to solve for the angle needed to reach h=1.72 m, with an initial speed of 15 m/s. + + e. The `newtraph.m` function needs a derivative, calculate the derivative of your + function with respect to theta, `dprojectile_dtheta.m`. This function should + output d(h)/d(theta). + + + f. In your `README.md` file, document the following under the heading `# + Homework #3`: + + i. Compare the number of iterations that each function needed to reach an + accuracy of 0.00001%. Include a table in your README.md with: + + ``` + | solver | initial guess(es) | ea | number of iterations| + | --- | --- | --- | --- | + |falsepos | | | | + |bisect | | | | + |newtraph | | | | + |mod_secant | | | | + ``` + + ii. Compare the convergence of the 4 methods. Plot the approximate error vs the + number of iterations that the solver has calculated. Save the plot as + `convergence.png` and display the plot in your `README.md` with: + + `![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 for f(x)=x*exp(-x^2). + + b. Add the results to a table in the `README.md` with: + + ``` + ### divergence of Newton-Raphson method + + | iteration | x_i | approx error | + | --- | --- | --- | + | 0 | 2 | n/a | + | 1 | | | + | 2 | | | + | 3 | | | + | 4 | | | + | 5 | | | + ``` + + c. Repeat steps a-b for an initial guess of 0.2. (But change the heading from + 'divergence' to 'convergence') + +3. Commit your changes to your repository. Sync your local repository with github. Then +copy and paste the "clone URL" into the following Google Form [Homework 3](https://goo.gl/forms/UJBGwp0fQcSxImkq2) diff --git a/README.md b/README.md index 8059993..73ca14a 100644 --- a/README.md +++ b/README.md @@ -74,8 +74,8 @@ general, I will not post homework solutions. |3|1/31||Consistent Coding habits| | |2/2|5|Root Finding| |4|2/7|6|Root Finding con’d| -| |2/9|7|Optimization| -|5|2/14||Intro to Linear Algebra| +| |2/9|7| **Snow Day**| +|5|2/14|| Optimization | | |2/16|8|Linear Algebra| |6|2/21|9|Linear systems: Gauss elimination| | |2/23|10|Linear Systems: LU factorization| diff --git a/lecture_07/.lecture_07.md.swp b/lecture_07/.lecture_07.md.swp new file mode 100644 index 0000000..0998c86 Binary files /dev/null and b/lecture_07/.lecture_07.md.swp differ diff --git a/lecture_07/.newtraph.m.swp b/lecture_07/.newtraph.m.swp new file mode 100644 index 0000000..9759dd2 Binary files /dev/null and b/lecture_07/.newtraph.m.swp differ diff --git a/lecture_07/lecture_07.ipynb b/lecture_07/lecture_07.ipynb index d7ff495..08d29c4 100644 --- a/lecture_07/lecture_07.ipynb +++ b/lecture_07/lecture_07.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 13, "metadata": { "collapsed": true }, @@ -40,12 +40,12 @@ "\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$" + "Use Newton-Raphson to find solution when $e^{-x}=x$" ] }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 2, "metadata": { "collapsed": false }, @@ -54,7 +54,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "error_approx = 1\r\n" + "x_r = 0.50000\n", + "error_approx = 1\n" ] } ], @@ -63,13 +64,14 @@ "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_r=x_i;\n" + "x_i=x_r;\n" ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 3, "metadata": { "collapsed": false }, @@ -78,8 +80,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "x_r = 0.50000\n", - "error_approx = 1\n" + "x_r = 0.56631\n", + "error_approx = 0.11709\n" ] } ], @@ -91,7 +93,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 4, "metadata": { "collapsed": false }, @@ -100,8 +102,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "x_r = 0.56631\n", - "error_approx = 0.11709\n" + "x_r = 0.56714\n", + "error_approx = 0.0014673\n" ] } ], @@ -113,7 +115,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 5, "metadata": { "collapsed": false }, @@ -123,7 +125,7 @@ "output_type": "stream", "text": [ "x_r = 0.56714\n", - "error_approx = 0.0014673\n" + "error_approx = 2.2106e-07\n" ] } ], @@ -149,7 +151,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "metadata": { "collapsed": true }, @@ -176,12 +178,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "ans = 142.74\r\n" + "root = 142.74\n", + "ea = 8.0930e-06\n", + "iter = 48\n" ] } ], "source": [ - "newtraph(f_m,df_m,140,0.00001)" + "[root,ea,iter]=newtraph(f_m,df_m,140,0.00001)" ] }, { @@ -217,7 +221,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 11, "metadata": { "collapsed": false }, @@ -226,17 +230,24 @@ "name": "stdout", "output_type": "stream", "text": [ - "ans = 142.74\r\n" + "root = 142.74\n", + "ea = 3.0615e-07\n", + "iter = 7\n" ] } ], "source": [ - "mod_secant(f_m,1e-6,50,0.00001)" + "[root,ea,iter]=mod_secant(f_m,1,50,0.00001)" ] }, + { + "cell_type": "raw", + "metadata": {}, + "source": [] + }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 15, "metadata": { "collapsed": false }, @@ -245,23 +256,159 @@ "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" + "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": [ - "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)" + "car_payments(400,30000,0.05,5,1)" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 17, "metadata": { "collapsed": false }, @@ -270,12 +417,195 @@ "name": "stdout", "output_type": "stream", "text": [ - "error: 'Amt' undefined near line 1 column 1\r\n" + "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*12*5" + "Amt_numerical*12*30" ] }, { @@ -287,7 +617,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 19, "metadata": { "collapsed": false }, @@ -398,7 +728,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 20, "metadata": { "collapsed": false }, @@ -428,7 +758,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 23, "metadata": { "collapsed": false }, @@ -611,7 +941,7 @@ "\t\n", "\n", "\n", - "\t\n", + "\t\n", "\t\n", "\tmod-secant\n", "\n", @@ -667,19 +997,20 @@ "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_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", + "\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": 16, + "execution_count": 22, "metadata": { "collapsed": false }, @@ -688,34 +1019,30 @@ "name": "stdout", "output_type": "stream", "text": [ - "ea_ms =\n", - "\n", - " Columns 1 through 6:\n", + "ea_nr =\n", "\n", - " 2.3382e+03 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14\n", + " Columns 1 through 8:\n", "\n", - " Columns 7 through 12:\n", + " 6.36591 0.06436 0.00052 0.00000 0.00000 0.00000 0.00000 0.00000\n", "\n", - " 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14\n", + " Columns 9 through 16:\n", "\n", - " Columns 13 through 18:\n", + " 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n", "\n", - " 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14\n", + " Columns 17 through 20:\n", "\n", - " Columns 19 and 20:\n", - "\n", - " 1.9171e-14 1.9171e-14\n", + " 0.00000 0.00000 0.00000 0.00000\n", "\n" ] } ], "source": [ - "ea_ms" + "ea_nr" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 26, "metadata": { "collapsed": false }, @@ -953,7 +1280,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 27, "metadata": { "collapsed": false }, @@ -962,30 +1289,26 @@ "name": "stdout", "output_type": "stream", "text": [ - "ea_bs =\n", + "ea_nr =\n", "\n", - " Columns 1 through 6:\n", + " Columns 1 through 7:\n", "\n", - " 9.5357e+03 -4.7554e-01 -2.1114e-01 6.0163e-02 -2.4387e-03 6.1052e-04\n", + " 99.03195 11.11111 11.11111 11.11111 11.11111 11.11111 11.11111\n", "\n", - " Columns 7 through 12:\n", + " Columns 8 through 14:\n", "\n", - " 2.2891e-04 -9.5367e-06 2.3842e-06 8.9407e-07 -2.2352e-07 9.3132e-09\n", + " 11.11111 11.11111 11.11111 11.11109 11.11052 11.10624 10.99684\n", "\n", - " Columns 13 through 18:\n", + " Columns 15 through 20:\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", + " 8.76956 2.12993 0.00000 0.00000 0.00000 0.00000\n", "\n", "ans = 16.208\n" ] } ], "source": [ - "ea_bs\n", + "ea_nr\n", "newtraph(f,df,0.5,0,12)" ] }, @@ -1008,6 +1331,99 @@ "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, diff --git a/lecture_07/lecture_07.md b/lecture_07/lecture_07.md index ffdc48d..3aacefb 100644 --- a/lecture_07/lecture_07.md +++ b/lecture_07/lecture_07.md @@ -18,7 +18,7 @@ $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$ +Use Newton-Raphson to find solution when $e^{-x}=x$ ```octave @@ -26,12 +26,14 @@ 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_r=x_i; +x_i=x_r; ``` - error_approx = 1 + x_r = 0.50000 + error_approx = 1 @@ -41,8 +43,8 @@ error_approx = abs((x_r-x_i)/x_r) x_i=x_r; ``` - x_r = 0.50000 - error_approx = 1 + x_r = 0.56631 + error_approx = 0.11709 @@ -52,8 +54,8 @@ error_approx = abs((x_r-x_i)/x_r) x_i=x_r; ``` - x_r = 0.56631 - error_approx = 0.11709 + x_r = 0.56714 + error_approx = 0.0014673 @@ -64,7 +66,7 @@ x_i=x_r; ``` x_r = 0.56714 - error_approx = 0.0014673 + 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. @@ -90,10 +92,12 @@ 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 ```octave -newtraph(f_m,df_m,140,0.00001) +[root,ea,iter]=newtraph(f_m,df_m,140,0.00001) ``` - ans = 142.74 + root = 142.74 + ea = 8.0930e-06 + iter = 48 ## Secant Methods @@ -119,32 +123,46 @@ $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) +[root,ea,iter]=mod_secant(f_m,1,50,0.00001) ``` - ans = 142.74 + root = 142.74 + ea = 3.0615e-07 + iter = 7 ```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) +car_payments(400,30000,0.05,5,1) ``` - error: 'plot_bool' undefined near line 12 column 6 - error: called from - car_payments at line 12 column 3 - mod_secant at line 22 column 8 - error: 'Amt' undefined near line 1 column 14 - error: evaluating argument list element number 1 + ans = 1.1185e+04 + + + +![svg](lecture_07_files/lecture_07_15_1.svg) ```octave -Amt*12*5 +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) ``` - error: 'Amt' undefined near line 1 column 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. @@ -257,12 +275,13 @@ 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_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') ``` @@ -281,31 +300,27 @@ legend('newton-raphson','mod-secant','false point','bisection') -![svg](lecture_07_files/lecture_07_22_1.svg) +![svg](lecture_07_files/lecture_07_24_1.svg) ```octave -ea_ms +ea_nr ``` - ea_ms = + ea_nr = - Columns 1 through 6: + Columns 1 through 8: - 2.3382e+03 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 + 6.36591 0.06436 0.00052 0.00000 0.00000 0.00000 0.00000 0.00000 - Columns 7 through 12: + Columns 9 through 16: - 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 + 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 - Columns 13 through 18: + Columns 17 through 20: - 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 1.9171e-14 - - Columns 19 and 20: - - 1.9171e-14 1.9171e-14 + 0.00000 0.00000 0.00000 0.00000 @@ -344,32 +359,28 @@ legend('newton-raphson','mod-secant','false point','bisection') -![svg](lecture_07_files/lecture_07_24_1.svg) +![svg](lecture_07_files/lecture_07_26_1.svg) ```octave -ea_bs +ea_nr newtraph(f,df,0.5,0,12) ``` - ea_bs = - - Columns 1 through 6: + ea_nr = - 9.5357e+03 -4.7554e-01 -2.1114e-01 6.0163e-02 -2.4387e-03 6.1052e-04 + Columns 1 through 7: - Columns 7 through 12: + 99.03195 11.11111 11.11111 11.11111 11.11111 11.11111 11.11111 - 2.2891e-04 -9.5367e-06 2.3842e-06 8.9407e-07 -2.2352e-07 9.3132e-09 + Columns 8 through 14: - Columns 13 through 18: + 11.11111 11.11111 11.11111 11.11109 11.11052 11.10624 10.99684 - -2.3283e-09 -8.7311e-10 3.6380e-11 -9.0949e-12 -3.4106e-12 8.5265e-13 + Columns 15 through 20: - Columns 19 and 20: - - -3.5527e-14 8.8818e-15 + 8.76956 2.12993 0.00000 0.00000 0.00000 0.00000 ans = 16.208 @@ -383,6 +394,55 @@ df(300) +```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 index 79b130c..70d8fc2 100644 Binary files a/lecture_07/lecture_07.pdf and b/lecture_07/lecture_07.pdf differ 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_24_1.svg b/lecture_07/lecture_07_files/lecture_07_24_1.svg index a8b614c..b511b10 100644 --- a/lecture_07/lecture_07_files/lecture_07_24_1.svg +++ b/lecture_07/lecture_07_files/lecture_07_24_1.svg @@ -43,47 +43,47 @@ - + 10-14 - + 10-12 - + 10-10 - + 10-8 - + 10-6 - + 10-4 - + 10-2 - + 100 - + 102 @@ -98,28 +98,43 @@ - - 10 + + 50 + + + + + 100 - - 20 + + 150 - - 30 + + 200 - - 40 + + 250 + + + + + 300 + + + + + 350 - 50 + 400 @@ -141,7 +156,7 @@ - + mod-secant @@ -150,7 +165,7 @@ - + false point @@ -159,7 +174,7 @@ - + bisection @@ -168,7 +183,7 @@ - + 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_08/.ipynb_checkpoints/lecture_08-checkpoint.ipynb b/lecture_08/.ipynb_checkpoints/lecture_08-checkpoint.ipynb new file mode 100644 index 0000000..2fd6442 --- /dev/null +++ b/lecture_08/.ipynb_checkpoints/lecture_08-checkpoint.ipynb @@ -0,0 +1,6 @@ +{ + "cells": [], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lecture_08/Auchain_model.gif b/lecture_08/Auchain_model.gif new file mode 100644 index 0000000..2f8c78c Binary files /dev/null and b/lecture_08/Auchain_model.gif differ diff --git a/lecture_08/Auchain_model.png b/lecture_08/Auchain_model.png new file mode 100644 index 0000000..dce9cd7 Binary files /dev/null and b/lecture_08/Auchain_model.png differ diff --git a/lecture_08/au_chain.jpg b/lecture_08/au_chain.jpg new file mode 100644 index 0000000..492fee8 Binary files /dev/null and b/lecture_08/au_chain.jpg differ diff --git a/lecture_08/goldenratio.png b/lecture_08/goldenratio.png new file mode 100644 index 0000000..9493c8c Binary files /dev/null and b/lecture_08/goldenratio.png differ diff --git a/lecture_08/goldmin.m b/lecture_08/goldmin.m new file mode 100644 index 0000000..fb4ce0b --- /dev/null +++ b/lecture_08/goldmin.m @@ -0,0 +1,36 @@ +function [x,fx,ea,iter]=goldmin(f,xl,xu,es,maxit,varargin) +% goldmin: minimization golden section search +% [x,fx,ea,iter]=goldmin(f,xl,xu,es,maxit,p1,p2,...): +% uses golden section search to find the minimum of f +% input: +% f = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by f +% output: +% x = location of minimum +% fx = minimum function value +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +phi=(1+sqrt(5))/2; +iter=0; +while(1) + d = (phi-1)*(xu - xl); + x1 = xl + d; + x2 = xu - d; + if f(x1,varargin{:}) < f(x2,varargin{:}) + xopt = x1; + xl = x2; + else + xopt = x2; + xu = x1; + end + iter=iter+1; + if xopt~=0, ea = (2 - phi) * abs((xu - xl) / xopt) * 100;end + if ea <= es | iter >= maxit,break,end +end +x=xopt;fx=f(xopt,varargin{:}); diff --git a/lecture_08/lecture_08.ipynb b/lecture_08/lecture_08.ipynb new file mode 100644 index 0000000..32b973c --- /dev/null +++ b/lecture_08/lecture_08.ipynb @@ -0,0 +1,2762 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "function h = parabola(x)\n", + " y=1;\n", + " h=sum(x.^2+y.^2-1);\n", + "end\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 0\r\n" + ] + } + ], + "source": [ + "fzero(@(x) parabola(x),1)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Optimization\n", + "\n", + "Many problems involve finding a minimum or maximum based on given constraints. Engineers and scientists typically use energy balance equations to find the conditions of minimum energy, but this value may never go to 0 and the actual value of energy may not be of interest. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Lennard-Jones potential is commonly used to model interatomic bonding. \n", + "\n", + "$E_{LJ}(x)=4\\epsilon \\left(\\left(\\frac{\\sigma}{x}\\right)^{12}-\\left(\\frac{\\sigma}{x}\\right)^{6}\\right)$\n", + "\n", + "Considering a 1-D gold chain, we can calculate the bond length, $x_{b}$, with no force applied to the chain and even for tension, F. This will allow us to calculate the nonlinear spring constant of a 1-D gold chain. \n", + "\n", + "![TEM image of Gold chain](au_chain.jpg)\n", + "\n", + "Computational Tools to Study and Predict the Long-Term Stability of Nanowires.\n", + "By Martin E. Zoloff Michoff, Patricio Vélez, Sergio A. Dassie and Ezequiel P. M. Leiva \n", + "\n", + "![Model of Gold chain, from molecular dynamics simulation](Auchain_model.png)\n", + "\n", + "[Single atom gold chain mechanics](http://www.uam.es/personal_pdi/ciencias/agrait/)\n", + "\n", + "### First, let's find the minimum energy $\\min(E_{LJ}(x))$\n", + "\n", + "## Brute force" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.08\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "setdefaults\n", + "epsilon = 0.039; % kcal/mol\n", + "sigma = 2.934; % Angstrom\n", + "x=linspace(2.8,6,200); % bond length in Angstrom\n", + "\n", + "Ex = lennard_jones(x,sigma,epsilon);\n", + "\n", + "[Emin,imin]=min(Ex);\n", + "\n", + "plot(x,Ex,x(imin),Emin,'o')" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 3.2824\n", + "ans = 3.2985\n", + "ans = 3.3146\n" + ] + } + ], + "source": [ + "x(imin-1)\n", + "x(imin)\n", + "x(imin+1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Golden Search Algorithm\n", + "\n", + "We can't just look for a sign change for the problem (unless we can take a derivative) so we need a new approach to determine whether we have a maximum between the two bounds.\n", + "\n", + "Rather than using the midpoint of initial bounds, here the problem is more difficult. We need to compare the values of 4 function evaluations. The golden search uses the golden ratio to determine two interior points. \n", + "\n", + "![golden ratio](goldenratio.png)\n", + "\n", + "Start with bounds of 2.5 and 6 Angstrom. " + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "current_min = -0.019959\r\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.08\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_5a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_6a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "% define Au atomic potential\n", + "epsilon = 0.039; % kcal/mol\n", + "sigma = 2.934; % Angstrom\n", + "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n", + "\n", + "% calculate golden ratio\n", + "phi = 1/2+sqrt(5)/2;\n", + "% set initial limits\n", + "x_l=2.8; \n", + "x_u=6; \n", + "\n", + "% Iteration #1\n", + "d=(phi-1)*(x_u-x_l);\n", + "\n", + "x1=x_l+d; % define point 1\n", + "x2=x_u-d; % define point 2\n", + "\n", + "\n", + "% evaluate Au_x(x1) and Au_x(x2)\n", + "\n", + "f1=Au_x(x1);\n", + "f2=Au_x(x2);\n", + "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", + "hold on;\n", + "\n", + "if f2\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.08\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_5a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_6a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "% Iteration #2\n", + "d=(phi-1)*(x_u-x_l);\n", + "\n", + "x1=x_l+d; % define point 1\n", + "x2=x_u-d; % define point 2\n", + "\n", + "% evaluate Au_x(x1) and Au_x(x2)\n", + "\n", + "f1=Au_x(x1);\n", + "f2=Au_x(x2);\n", + "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", + "hold on;\n", + "\n", + "if f2\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.08\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_5a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_6a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "% Iteration #3\n", + "d=(phi-1)*(x_u-x_l);\n", + "\n", + "x1=x_l+d; % define point 1\n", + "x2=x_u-d; % define point 2\n", + "\n", + "% evaluate Au_x(x1) and Au_x(x2)\n", + "\n", + "f1=Au_x(x1);\n", + "f2=Au_x(x2);\n", + "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", + "hold on;\n", + "\n", + "if f2\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.08\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_5a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_6a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "% Iteration #3\n", + "d=(phi-1)*(x_u-x_l);\n", + "\n", + "x1=x_l+d; % define point 1\n", + "x2=x_u-d; % define point 2\n", + "\n", + "% evaluate Au_x(x1) and Au_x(x2)\n", + "\n", + "f1=Au_x(x1);\n", + "f2=Au_x(x2);\n", + "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", + "hold on;\n", + "\n", + "if f2\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.06\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.08\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_5a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "% define Au atomic potential\n", + "epsilon = 0.039; % kcal/mol\n", + "sigma = 2.934; % Angstrom\n", + "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n", + "\n", + "% set initial limits\n", + "x_l=2.8; \n", + "x_u=3.5; \n", + "\n", + "% Iteration #1\n", + "x1=x_l;\n", + "x2=mean([x_l,x_u]);\n", + "x3=x_u;\n", + "\n", + "% evaluate Au_x(x1), Au_x(x2) and Au_x(x3)\n", + " \n", + "f1=Au_x(x1);\n", + "f2=Au_x(x2);\n", + "f3=Au_x(x3);\n", + "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n", + "x_fit = linspace(x1,x3,20);\n", + "y_fit = polyval(p,x_fit);\n", + "\n", + "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n", + "hold on\n", + "if f2x2\n", + " plot(x4,f4,'*',[x1,x2],[f1,f2])\n", + " x1=x2;\n", + " f1=f2;\n", + " else\n", + " plot(x4,f4,'*',[x3,x2],[f3,f2])\n", + " x3=x2;\n", + " f3=f2;\n", + " end\n", + " x2=x4; f2=f4;\n", + "else\n", + " error('no minimum in bracket')\n", + "end\n", + "hold off" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.08\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_5a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n", + "x_fit = linspace(x1,x3,20);\n", + "y_fit = polyval(p,x_fit);\n", + "\n", + "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n", + "hold on\n", + "if f2x2\n", + " plot(x4,f4,'*',[x1,x2],[f1,f2])\n", + " x1=x2;\n", + " f1=f2;\n", + " else\n", + " plot(x4,f4,'*',[x3,x2],[f3,f2])\n", + " x3=x2;\n", + " f3=f2;\n", + " end\n", + " x2=x4; f2=f4;\n", + "else\n", + " error('no minimum in bracket')\n", + "end\n", + "hold off\n" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.08\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_5a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n", + "x_fit = linspace(x1,x3,20);\n", + "y_fit = polyval(p,x_fit);\n", + "\n", + "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n", + "hold on\n", + "if f2x2\n", + " plot(x4,f4,'*',[x1,x2],[f1,f2])\n", + " x1=x2;\n", + " f1=f2;\n", + " else\n", + " plot(x4,f4,'*',[x3,x2],[f3,f2])\n", + " x3=x2;\n", + " f3=f2;\n", + " end\n", + " x2=x4; f2=f4;\n", + "else\n", + " error('no minimum in bracket')\n", + "end\n", + "hold off\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Parabolic interpolation does not converge in many scenarios even though it it a bracketing method. Instead, functions like `fminbnd` in Matlab and Octave use a combination of the two (Golden Ratio and Parabolic)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Using the solutions to minimization for the nonlinear spring constant\n", + "\n", + "Now, we have two routines to find minimums of a univariate function (Golden Ratio and Parabolic). Let's use these to solve for the minimum energy given a range of applied forces to the single atom gold chain\n", + "\n", + "$E_{total}(\\Delta x) = E_{LJ}(x_{min}+\\Delta x) - F \\cdot \\Delta x$\n", + "\n", + "$1 aJ = 10^{-18} J = 1~nN* 1~nm = 10^{-9}N * 10^{-9} m$" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: called from\n", + " goldmin at line 17 column 1\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "xmin = 0.32933\n", + "Emin = -2.7096e-04\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.0004\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.0002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0004\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0006\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.7\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tLennard Jones Potential (aJ/atom)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tbond length (nm)\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "epsilon = 0.039; % kcal/mol\n", + "epsilon = epsilon*6.9477e-21; % J/atom\n", + "epsilon = epsilon*1e18; % aJ/J\n", + "% final units for epsilon are aJ\n", + "\n", + "sigma = 2.934; % Angstrom\n", + "sigma = sigma*0.10; % nm/Angstrom\n", + "x=linspace(2.8,6,200)*0.10; % bond length in um\n", + "\n", + "Ex = lennard_jones(x,sigma,epsilon);\n", + "\n", + "%[Emin,imin]=min(Ex);\n", + "\n", + "[xmin,Emin] = goldmin(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)\n", + "\n", + "plot(x,Ex,xmin,Emin,'o')\n", + "ylabel('Lennard Jones Potential (aJ/atom)')\n", + "xlabel('bond length (nm)')\n", + "\n", + "Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;\n" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\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\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.01\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.03\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.035\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tForce (nN)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tdx (nm)\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": [ + "N=50;\n", + "warning('off')\n", + "dx = zeros(1,N); % [in nm]\n", + "F_applied=linspace(0,0.0022,N); % [in nN]\n", + "for i=1:N\n", + " optmin=goldmin(@(dx) Etotal(dx,F_applied(i)),-0.001,0.035);\n", + " dx(i)=optmin;\n", + "end\n", + "\n", + "plot(dx,F_applied)\n", + "xlabel('dx (nm)')\n", + "ylabel('Force (nN)')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For this function, it is possible to take a derivative and compare the analytical result:\n", + "\n", + "dE/dx = 0 = d(E_LJ)/dx - F\n", + "\n", + "F= d(E_LJ)/dx" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "F =\n", + "\n", + "@(dx) 4 * epsilon * 6 * (sigma ^ 6 ./ (xmin + dx) .^ 7 - 2 * sigma ^ 12 ./ (xmin + dx) .^ 13)\n", + "\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.01\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.03\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.05\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "dx_full=linspace(0,0.06,50);\n", + "F= @(dx) 4*epsilon*6*(sigma^6./(xmin+dx).^7-2*sigma^12./(xmin+dx).^13)\n", + "plot(dx_full,F(dx_full),dx,F_applied)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "## Curve-fitting\n", + "Another example is minimizing error in your approximation of a function. If you have data (now we have Force-displacement data) we can fit this to a function, such as:\n", + "\n", + "$F(x) = K_{1}\\Delta x + \\frac{1}{2} K_{2}(\\Delta x)^{2}$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "function SSE = sse_of_parabola(K,xdata,ydata)\n", + " % calculate the sum of squares error for a parabola given a function, func, and xdata and ydata\n", + " % output is SSE=sum of squares error\n", + " K1=K(1);\n", + " K2=K(2);\n", + " y_function = K1*xdata+1/2*K2*xdata.^2;\n", + " SSE = sum((ydata-y_function).^2);\n", + "end\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "The nonlinear spring constants are K1=0.16 nN/nm and K2=-5.98 nN/nm^2\n", + "The mininum sum of squares error = 7.35e-08\n" + ] + } + ], + "source": [ + "[K,SSE_min]=fminsearch(@(K) sse_of_parabola(K,dx,F_applied),[1,1]);\n", + "fprintf('\\nThe nonlinear spring constants are K1=%1.2f nN/nm and K2=%1.2f nN/nm^2\\n',K)\n", + "fprintf('The mininum sum of squares error = %1.2e',SSE_min)" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\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\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.01\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.03\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.035\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\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", + "\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", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lecture_08/lecture_08.md b/lecture_08/lecture_08.md new file mode 100644 index 0000000..79edc0a --- /dev/null +++ b/lecture_08/lecture_08.md @@ -0,0 +1,626 @@ + + +```octave +%plot --format svg +``` + +# Optimization + +Many problems involve finding a minimum or maximum based on given constraints. Engineers and scientists typically use energy balance equations to find the conditions of minimum energy, but this value may never go to 0 and the actual value of energy may not be of interest. + +The Lennard-Jones potential is commonly used to model interatomic bonding. + +$E_{LJ}(x)=4\epsilon \left(\left(\frac{\sigma}{x}\right)^{12}-\left(\frac{\sigma}{x}\right)^{6}\right)$ + +Considering a 1-D gold chain, we can calculate the bond length, $x_{b}$, with no force applied to the chain and even for tension, F. This will allow us to calculate the nonlinear spring constant of a 1-D gold chain. + +![TEM image of Gold chain](au_chain.jpg) + +Computational Tools to Study and Predict the Long-Term Stability of Nanowires. +By Martin E. Zoloff Michoff, Patricio Vélez, Sergio A. Dassie and Ezequiel P. M. Leiva + +![Model of Gold chain, from molecular dynamics simulation](Auchain_model.png) + +[Single atom gold chain mechanics](http://www.uam.es/personal_pdi/ciencias/agrait/) + +### First, let's find the minimum energy $\min(E_{LJ}(x))$ + +## Brute force + + +```octave +setdefaults +epsilon = 0.039; % kcal/mol +sigma = 2.934; % Angstrom +x=linspace(2.8,6,200); % bond length in Angstrom + +Ex = lennard_jones(x,sigma,epsilon); + +[Emin,imin]=min(Ex); + +plot(x,Ex,x(imin),Emin,'o') +``` + + +![svg](lecture_08_files/lecture_08_3_0.svg) + + + +```octave +x(imin-1) +x(imin) +x(imin+1) +``` + + ans = 3.2824 + ans = 3.2985 + ans = 3.3146 + + +## Golden Search Algorithm + +We can't just look for a sign change for the problem (unless we can take a derivative) so we need a new approach to determine whether we have a maximum between the two bounds. + +Rather than using the midpoint of initial bounds, here the problem is more difficult. We need to compare the values of 4 function evaluations. The golden search uses the golden ratio to determine two interior points. + +![golden ratio](goldenratio.png) + +Start with bounds of 2.5 and 6 Angstrom. + + +```octave +% define Au atomic potential +epsilon = 0.039; % kcal/mol +sigma = 2.934; % Angstrom +Au_x= @(x) lennard_jones(x,sigma,epsilon); + +% calculate golden ratio +phi = 1/2+sqrt(5)/2; +% set initial limits +x_l=2.8; +x_u=6; + +% Iteration #1 +d=(phi-1)*(x_u-x_l); + +x1=x_l+d; % define point 1 +x2=x_u-d; % define point 2 + + +% evaluate Au_x(x1) and Au_x(x2) + +f1=Au_x(x1); +f2=Au_x(x2); +plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go') +hold on; + +if f2x2 + plot(x4,f4,'*',[x1,x2],[f1,f2]) + x1=x2; + f1=f2; + else + plot(x4,f4,'*',[x3,x2],[f3,f2]) + x3=x2; + f3=f2; + end + x2=x4; f2=f4; +else + error('no minimum in bracket') +end +hold off +``` + + +![svg](lecture_08_files/lecture_08_11_0.svg) + + + +```octave +p = polyfit([x1,x2,x3],[f1,f2,f3],2); +x_fit = linspace(x1,x3,20); +y_fit = polyval(p,x_fit); + +plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o') +hold on +if f2x2 + plot(x4,f4,'*',[x1,x2],[f1,f2]) + x1=x2; + f1=f2; + else + plot(x4,f4,'*',[x3,x2],[f3,f2]) + x3=x2; + f3=f2; + end + x2=x4; f2=f4; +else + error('no minimum in bracket') +end +hold off + +``` + + +![svg](lecture_08_files/lecture_08_12_0.svg) + + + +```octave +p = polyfit([x1,x2,x3],[f1,f2,f3],2); +x_fit = linspace(x1,x3,20); +y_fit = polyval(p,x_fit); + +plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o') +hold on +if f2x2 + plot(x4,f4,'*',[x1,x2],[f1,f2]) + x1=x2; + f1=f2; + else + plot(x4,f4,'*',[x3,x2],[f3,f2]) + x3=x2; + f3=f2; + end + x2=x4; f2=f4; +else + error('no minimum in bracket') +end +hold off + +``` + + +![svg](lecture_08_files/lecture_08_13_0.svg) + + +Parabolic interpolation does not converge in many scenarios even though it it a bracketing method. Instead, functions like `fminbnd` in Matlab and Octave use a combination of the two (Golden Ratio and Parabolic) + +## Using the solutions to minimization for the nonlinear spring constant + +Now, we have two routines to find minimums of a univariate function (Golden Ratio and Parabolic). Let's use these to solve for the minimum energy given a range of applied forces to the single atom gold chain + +$E_{total}(\Delta x) = E_{LJ}(x_{min}+\Delta x) - F \cdot \Delta x$ + + +```octave +epsilon = 0.039; % kcal/mol +epsilon = epsilon*6.9477e-21; % J/atom +epsilon = epsilon*1e18; % aJ/J +% final units for epsilon are aJ + +sigma = 2.934; % Angstrom +sigma = sigma*0.10; % nm/Angstrom +x=linspace(2.8,6,200)*0.10; % bond length in um + +Ex = lennard_jones(x,sigma,epsilon); + +%[Emin,imin]=min(Ex); + +[xmin,Emin] = fminbnd(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6) + +plot(x,Ex,xmin,Emin,'o') +ylabel('Lennard Jones Potential (aJ/atom)') +xlabel('bond length (nm)') + +Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx; + +``` + + xmin = 0.32933 + Emin = -2.7096e-04 + + + +![svg](lecture_08_files/lecture_08_16_1.svg) + + + +```octave +N=50; +dx = zeros(1,N); % [in nm] +F_applied=linspace(0,0.0022,N); % [in nN] +for i=1:N + optmin=goldmin(@(dx) Etotal(dx,F_applied(i)),-0.001,0.035); + dx(i)=optmin; +end + +plot(dx,F_applied) +xlabel('dx (nm)') +ylabel('Force (nN)') +``` + + warning: Matlab-style short-circuit operation performed for operator | + warning: called from + goldmin at line 17 column 1 + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + + + +![svg](lecture_08_files/lecture_08_17_1.svg) + + +For this function, it is possible to take a derivative and compare the analytical result: + + +```octave +dx_full=linspace(0,0.06,50); +F= @(dx) 4*epsilon*6*(sigma^6./(xmin+dx).^7-2*sigma^12./(xmin+dx).^13) +plot(dx_full,F(dx_full),dx,F_applied) +``` + + F = + + @(dx) 4 * epsilon * 6 * (sigma ^ 6 ./ (xmin + dx) .^ 7 - 2 * sigma ^ 12 ./ (xmin + dx) .^ 13) + + + + +![svg](lecture_08_files/lecture_08_19_1.svg) + + +## Curve-fitting +Another example is minimizing error in your approximation of a function. If you have data (now we have Force-displacement data) we can fit this to a function, such as: + +$F(x) = K_{1}\Delta x + \frac{1}{2} K_{2}(\Delta x)^{2}$ + + + +```octave +function SSE = sse_of_parabola(K,xdata,ydata) + % calculate the sum of squares error for a parabola given a function, func, and xdata and ydata + % output is SSE=sum of squares error + K1=K(1); + K2=K(2); + y_function = K1*xdata+1/2*K2*xdata.^2; + SSE = sum((ydata-y_function).^2); +end + +``` + + +```octave +[K,SSE_min]=fminsearch(@(K) sse_of_parabola(K,dx,F_applied),[1,1]); +fprintf('\nThe nonlinear spring constants are K1=%1.2f nN/nm and K2=%1.2f nN/nm^2\n',K) +fprintf('The mininum sum of squares error = %1.2e',SSE_min) +``` + + + The nonlinear spring constants are K1=0.16 nN/nm and K2=-5.98 nN/nm^2 + The mininum sum of squares error = 7.35e-08 + + + +```octave +plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2) +``` + + +![svg](lecture_08_files/lecture_08_23_0.svg) + diff --git a/lecture_08/lecture_08.pdf b/lecture_08/lecture_08.pdf new file mode 100644 index 0000000..c251694 Binary files /dev/null and b/lecture_08/lecture_08.pdf differ diff --git a/lecture_08/lecture_08_files/lecture_08_11_0.svg b/lecture_08/lecture_08_files/lecture_08_11_0.svg new file mode 100644 index 0000000..c61f790 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_11_0.svg @@ -0,0 +1,168 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.06 + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + gnuplot_plot_3a + + + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_12_0.svg b/lecture_08/lecture_08_files/lecture_08_12_0.svg new file mode 100644 index 0000000..b641032 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_12_0.svg @@ -0,0 +1,163 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + gnuplot_plot_3a + + + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_13_0.svg b/lecture_08/lecture_08_files/lecture_08_13_0.svg new file mode 100644 index 0000000..3e467e5 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_13_0.svg @@ -0,0 +1,163 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + gnuplot_plot_3a + + + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_16_1.svg b/lecture_08/lecture_08_files/lecture_08_16_1.svg new file mode 100644 index 0000000..ff27c29 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_16_1.svg @@ -0,0 +1,142 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.0004 + + + + + -0.0002 + + + + + 0 + + + + + 0.0002 + + + + + 0.0004 + + + + + 0.0006 + + + + + 0.2 + + + + + 0.3 + + + + + 0.4 + + + + + 0.5 + + + + + 0.6 + + + + + 0.7 + + + + + + + + + Lennard Jones Potential (aJ/atom) + + + + + bond length (nm) + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_17_1.svg b/lecture_08/lecture_08_files/lecture_08_17_1.svg new file mode 100644 index 0000000..1f3103b --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_17_1.svg @@ -0,0 +1,136 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + 0.0005 + + + + + 0.001 + + + + + 0.0015 + + + + + 0.002 + + + + + 0.0025 + + + + + -0.01 + + + + + 0 + + + + + 0.01 + + + + + 0.02 + + + + + 0.03 + + + + + 0.04 + + + + + + + + + Force (nN) + + + + + dx (nm) + + + + + gnuplot_plot_1a + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_19_1.svg b/lecture_08/lecture_08_files/lecture_08_19_1.svg new file mode 100644 index 0000000..2a891df --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_19_1.svg @@ -0,0 +1,135 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + 0.0005 + + + + + 0.001 + + + + + 0.0015 + + + + + 0.002 + + + + + 0.0025 + + + + + 0 + + + + + 0.01 + + + + + 0.02 + + + + + 0.03 + + + + + 0.04 + + + + + 0.05 + + + + + 0.06 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_23_0.svg b/lecture_08/lecture_08_files/lecture_08_23_0.svg new file mode 100644 index 0000000..8fa3bb1 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_23_0.svg @@ -0,0 +1,196 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + 0.0005 + + + + + 0.001 + + + + + 0.0015 + + + + + 0.002 + + + + + 0.0025 + + + + + 0 + + + + + 0.005 + + + + + 0.01 + + + + + 0.015 + + + + + 0.02 + + + + + 0.025 + + + + + 0.03 + + + + + 0.035 + + + + + 0.04 + + + + + + + + + gnuplot_plot_1a + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + gnuplot_plot_2a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_3_0.svg b/lecture_08/lecture_08_files/lecture_08_3_0.svg new file mode 100644 index 0000000..2c38150 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_3_0.svg @@ -0,0 +1,147 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_6_1.svg b/lecture_08/lecture_08_files/lecture_08_6_1.svg new file mode 100644 index 0000000..ecf7a18 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_6_1.svg @@ -0,0 +1,169 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + gnuplot_plot_3a + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + gnuplot_plot_6a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_7_1.svg b/lecture_08/lecture_08_files/lecture_08_7_1.svg new file mode 100644 index 0000000..e63011c --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_7_1.svg @@ -0,0 +1,169 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + gnuplot_plot_3a + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + gnuplot_plot_6a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_8_1.svg b/lecture_08/lecture_08_files/lecture_08_8_1.svg new file mode 100644 index 0000000..873b04d --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_8_1.svg @@ -0,0 +1,169 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + gnuplot_plot_3a + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + gnuplot_plot_6a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_9_1.svg b/lecture_08/lecture_08_files/lecture_08_9_1.svg new file mode 100644 index 0000000..28fd681 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_9_1.svg @@ -0,0 +1,169 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + gnuplot_plot_3a + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + gnuplot_plot_6a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lennard_jones.m b/lecture_08/lennard_jones.m new file mode 100644 index 0000000..d18a6ad --- /dev/null +++ b/lecture_08/lennard_jones.m @@ -0,0 +1,4 @@ +function E_LJ =lennard_jones(x,sigma,epsilon) + E_LJ = 4*epsilon*((sigma./x).^12-(sigma./x).^6); +end + diff --git a/lecture_08/octave-workspace b/lecture_08/octave-workspace new file mode 100644 index 0000000..8c437bb Binary files /dev/null and b/lecture_08/octave-workspace differ