diff --git a/06_roots-1/06_roots-1.ipynb b/06_roots-1/06_roots-1.ipynb
index e95edef..1e1c23e 100644
--- a/06_roots-1/06_roots-1.ipynb
+++ b/06_roots-1/06_roots-1.ipynb
@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": 2,
"metadata": {
"collapsed": true,
"scrolled": true,
@@ -53,7 +53,7 @@
},
{
"cell_type": "code",
- "execution_count": 2,
+ "execution_count": 1,
"metadata": {
"collapsed": false,
"scrolled": true,
@@ -64,155 +64,9 @@
"outputs": [
{
"data": {
- "image/svg+xml": [
- ""
- ],
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAJNmlDQ1BkZWZhdWx0X3JnYi5pY2MA\nAHiclZFnUJSHFobP933bCwvssnRYepMqZQHpvUmvogJL7yxLEbEhYgQiiog0RZCggAGjUiRWRLEQ\nFBSxoFkkCCgxGEVUUPLDOxPn3vHHfX49884755yZA0ARBQBARQFSUgV8Pxd7TkhoGAe+IZKXmW7n\n4+MJ3+X9KCAAAPdWfb/zXSjRMZk8AFgGgHxeOl8AgOQCgGaOIF0AgBwFAFZUUroAADkLACx+SGgY\nAHIDAFhxX30cAFhRX30eAFj8AD8HABQHQKLFfeNR3/h/9gIAKNvxBQmxMbkc/7RYQU4kP4aT6edi\nz3FzcOD48NNiE5Jjvjn4/yp/B0FMrgAAwCEtfRM/IS5ewPmfoUYGhobw7y/e+gICAAh78L//AwDf\n9NIaAbgLANi+f7OoaoDuXQBSj//NVI8CMAoBuu7wsvjZXzMcAAAeKMAAFkiDAqiAJuiCEZiBJdiC\nE7iDNwRAKGwAHsRDCvAhB/JhBxRBCeyDg1AD9dAELdAOp6EbzsMVuA634S6MwhMQwhS8gnl4D0sI\nghAROsJEpBFFRA3RQYwQLmKNOCGeiB8SikQgcUgqkoXkIzuREqQcqUEakBbkF+QccgW5iQwjj5AJ\nZBb5G/mEYigNZaHyqDqqj3JRO9QDDUDXo3FoBpqHFqJ70Sq0ET2JdqFX0NvoKCpEX6ELGGBUjI0p\nYboYF3PAvLEwLBbjY1uxYqwSa8TasV5sALuHCbE57COOgGPiODhdnCXOFReI4+EycFtxpbga3Alc\nF64fdw83gZvHfcHT8XJ4HbwF3g0fgo/D5+CL8JX4Znwn/hp+FD+Ff08gENgEDYIZwZUQSkgkbCaU\nEg4TOgiXCcOEScICkUiUJuoQrYjexEiigFhErCaeJF4ijhCniB9IVJIiyYjkTAojpZIKSJWkVtJF\n0ghpmrREFiWrkS3I3uRo8iZyGbmJ3Eu+Q54iL1HEKBoUK0oAJZGyg1JFaadco4xT3lKpVGWqOdWX\nmkDdTq2inqLeoE5QP9LEado0B1o4LYu2l3acdpn2iPaWTqer023pYXQBfS+9hX6V/oz+QYQpoifi\nJhItsk2kVqRLZETkNYPMUGPYMTYw8hiVjDOMO4w5UbKouqiDaKToVtFa0XOiY6ILYkwxQzFvsRSx\nUrFWsZtiM+JEcXVxJ/Fo8ULxY+JXxSeZGFOF6cDkMXcym5jXmFMsAkuD5cZKZJWwfmYNseYlxCWM\nJYIkciVqJS5ICNkYW53txk5ml7FPsx+wP0nKS9pJxkjukWyXHJFclJKVspWKkSqW6pAalfokzZF2\nkk6S3i/dLf1UBiejLeMrkyNzROaazJwsS9ZSlidbLHta9rEcKqct5ye3We6Y3KDcgryCvIt8uny1\n/FX5OQW2gq1CokKFwkWFWUWmorVigmKF4iXFlxwJjh0nmVPF6efMK8kpuSplKTUoDSktKWsoByoX\nKHcoP1WhqHBVYlUqVPpU5lUVVb1U81XbVB+rkdW4avFqh9QG1BbVNdSD1Xerd6vPaEhpuGnkabRp\njGvSNW00MzQbNe9rEbS4Wklah7XuaqPaJtrx2rXad3RQHVOdBJ3DOsOr8KvMV6Wualw1pkvTtdPN\n1m3TndBj63nqFeh1673WV9UP09+vP6D/xcDEINmgyeCJobihu2GBYa/h30baRjyjWqP7q+mrnVdv\nW92z+o2xjnGM8RHjhyZMEy+T3SZ9Jp9NzUz5pu2ms2aqZhFmdWZjXBbXh1vKvWGON7c332Z+3vyj\nhamFwOK0xV+WupZJlq2WM2s01sSsaVozaaVsFWnVYCW05lhHWB+1Ftoo2UTaNNo8t1WxjbZttp22\n07JLtDtp99rewJ5v32m/6GDhsMXhsiPm6OJY7DjkJO4U6FTj9MxZ2TnOuc153sXEZbPLZVe8q4fr\nftcxN3k3nluL27y7mfsW934Pmoe/R43Hc09tT75nrxfq5e51wGt8rdra1LXd3uDt5n3A+6mPhk+G\nz6++BF8f31rfF36Gfvl+A/5M/43+rf7vA+wDygKeBGoGZgX2BTGCwoNaghaDHYPLg4Uh+iFbQm6H\nyoQmhPaEEcOCwprDFtY5rTu4bircJLwo/MF6jfW5629ukNmQvOHCRsbGyI1nIvARwRGtEcuR3pGN\nkQtRblF1UfM8B94h3qto2+iK6NkYq5jymOlYq9jy2Jk4q7gDcbPxNvGV8XMJDgk1CW8SXRPrExeT\nvJOOJ60kByd3pJBSIlLOpYqnJqX2pymk5aYNp+ukF6ULMywyDmbM8z34zZlI5vrMHgFLkC4YzNLM\n2pU1kW2dXZv9ISco50yuWG5q7uAm7U17Nk3nOef9tBm3mbe5L18pf0f+xBa7LQ1bka1RW/u2qWwr\n3Da13WX7iR2UHUk7fiswKCgveLczeGdvoXzh9sLJXS672opEivhFY7std9f/gPsh4YehPav3VO/5\nUhxdfKvEoKSyZLmUV3rrR8Mfq35c2Ru7d6jMtOzIPsK+1H0P9tvsP1EuVp5XPnnA60BXBaeiuOLd\nwY0Hb1YaV9YfohzKOiSs8qzqqVat3le9XBNfM1prX9tRJ1e3p27xcPThkSO2R9rr5etL6j8dTTj6\nsMGloatRvbHyGOFY9rEXTUFNAz9xf2pplmkuaf58PPW48ITfif4Ws5aWVrnWsja0Latt9mT4ybs/\nO/7c067b3tDB7ig5BaeyTr38JeKXB6c9Tved4Z5pP6t2tq6T2VnchXRt6prvju8W9oT2DJ9zP9fX\na9nb+aver8fPK52vvSBxoewi5WLhxZVLeZcWLqdfnrsSd2Wyb2Pfk6shV+/3+/YPXfO4duO68/Wr\nA3YDl25Y3Th/0+LmuVvcW923TW93DZoMdv5m8lvnkOlQ1x2zOz13ze/2Dq8ZvjhiM3LlnuO96/fd\n7t8eXTs6/CDwwcOx8DHhw+iHM4+SH715nP146cn2cfx48VPRp5XP5J41/q71e4fQVHhhwnFi8Ln/\n8yeTvMlXf2T+sTxV+IL+onJacbplxmjm/Kzz7N2X615OvUp/tTRX9KfYn3WvNV+f/cv2r8H5kPmp\nN/w3K3+XvpV+e/yd8bu+BZ+FZ+9T3i8tFn+Q/nDiI/fjwKfgT9NLOcvE5arPWp97v3h8GV9JWVn5\nBy6ikLxSF1/9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAAHXRFWHRTb2Z0d2FyZQBHUEwgR2hvc3Rz\nY3JpcHQgOS4xOJQFEHMAAA6+SURBVHic7d1BbuLm48dh+tdI7UjdMOtKrSz9TkCP4Cs4RyBHwKtK\nlbqA6QnwEeAI8RHGJ6jiRbeVwqZqVlX+CyQLkTRNOib5DnmeRQUJenmZBj74tQ1f3d3dTQDgtf3f\na08AACYTQQIghCABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACKM\nH6Ttdnt5eTn6sACctzGD1Pf9xcVF3/cjjgnAG/FuxLGKothsNpPJpK7rEYcF4C2wDwmACGNuIf03\n33///ddff72//O7du/fv37/iZG5vb193Ap/D5F+Lyb8Wk38xf//9919//TVcvbi4+OWXX0a/l9cP\n0u+///7aUwDgGf74449TDGvJDoAIggRABEEC4Hm++eabUww75j6kvu+bpplMJl3X7Y/8LopiPp+P\neBcAvLrvvvvuFMOOfB7ScrkccUAA3g5LdgBEECQAnufm5uYUwwoSAM9zolN6BQmA5xEkAM6ZIAEQ\nQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEEC\nIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiC\nBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgAR\n3o0+Yl3X+wtVVc1ms9HHB+AsjRykpmmGDtV1LUgAPNHIS3Zd1w0Rms1m2+123PEBOFcjB2k6nQ6X\nq6rqum7c8QE4Vw5qAOB5+r4/xbCCBMDzFEVximEFCYAIIwdpt9sNl7fbraPsAHiikYM0m82GAxm6\nrquqatzxAThXI5+HNJ/P67reH+2tRgA83fif1LBcLkcfE4Cz56AGACIIEgARBAmACIIEQARBAiCC\nIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRA\nBEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJ\ngAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACII\nEgARxg/Sdru9vLwcfVgAztuYQer7/uLiou/7EccE4I14N+JYRVFsNpvJZFLX9YjDAvAW2IcEQARB\nAiDC85bs+r5vmub+z5fL5UjzASDd7e3tKYZ9XpCKotAegDfu/fv3pxjWkh0AEQQJgAhjHvY97GHq\num5/5HdRFPP5fMS7AOBcjXwekj1MAPw3luwAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQ\nAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCC\nIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRA\nBEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAivBt3uK7r\nmqaZTqeTyaQsy7Isxx0fgHM1ZpB2u91qtdpsNvurq9WqKIqiKEa8CwDO1ZhLdn3fz+fz4episWia\nZsTxAThjY24hzWazw6tt2+7X7gDgX428D2lwtHwHwNm4vb09xbDPC1Lf9w+uwi2Xy8Oru92uruvN\nZmMLCeD8vH///iTj3o3t+vq6LMvr6+sn3v4kjwqAk/n48ePo7bi7uxv/sO+Li4vNZuPgOgCeZcyj\n7Lqu26/UDUc3bLfbEccH4IyNuYVU1/VsNttut4cdqqpqxLsA4FyNGaSrq6sRRwPgTfFZdgBEECQA\nIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBI\nAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBB\nkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIg\ngiABEEGQAIggSABEECQAIrwbd7iu65qmmU6nk8mkKIr5fD7u+ACcq5GDtN1u1+v1/nLbtm3blmU5\n7l0AcJZGXrJbLpfD5bIs27Ydd3wAztUJ9yHVdW3zCIAnOkmQ6rr+8OFDURSCBMATfXV3d/f0W/d9\n3zTN/Z8frtQNVqvVdDr91+Mavvrqq6dPAIBX9/PPP//000/jj3t3SovF4l9vM/5DAuCUPn78eIpk\nOA8JgAhjBqnrur7vh6uHlwHgcWOeh1QURV3X+7Ni9/99cN8SANw3ZpCm0+lwViwAPIt9SABEECQA\nIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBI\nAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBB\nkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIg\ngiABEEGQAIggSABEECQAIggSABFOFaSLi4umaU40OADn5yRBquu6qqq+708xOABnafwg9X3f931V\nVaOPDMAZGz9IdV0vl8vRhwXgvI0cpLZti6IoimLcYQHIcXt7e4ph3z3r1n3fP3iowrBJtFqtNpvN\nCPMCINX79+9PMezzglQUxSPLcW3bTiaT1Wq1v9p1Xdu2ZVl+zvwAeCvuTmaxWDzlZq/9DwDA83z8\n+PEU1XBiLAARThKktm3rut4v2Z1ifADOz/P2IT1RWZZ2HQHwLJbsAIggSABEECQAIggSABEECYAI\nggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIA\nEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAk\nACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIgg\nSABEeDfucHVdH14ty7Isy3HvAoDXdXNzc4phRw7SZDJZLpejjwlAjhMFyZIdABEECYAI4y/ZXVxc\nFEWxv7xYLKbT6eh3AcD5+eru7u7pt+77vmma+z9/cL/RbrdbrVb/ukvp119/HS5/++23//vf/54+\nHwBewJ9//vnbb78NV6uq+uGHH0a/l+cF6bnqunaMAwBPYR8SABHGDFLbtm3bDle7rhtxcADO28hL\ndk3T9H2/vzydTheLxYiDA3DGTrsPCQCeyD4kACIIEgARBAmACIIEQARBAiDC+J9l92U5PE798EMl\nuq7bbrf3fx5i/7FM+8tFUczn88PfBk5+u922bbter+//6vHZrlar3W43edWzCB6Z/PAFYNPpdD6f\nH31y4/Dbqqpms9mp53nfbrfbnx344OT3tttt0zRXV1dHP8+f/OEnmR09C8Inf/iyc/8DP1998vs/\n+P2s7s/h8afk507+7g1bLBafPn36p1/tL3z69Gm9Xr/gpJ5kPp/f3NzsLy+Xy6urq8PfRk3++vq6\nqqrlcjmfzx+8wSOzXa/Xw0O7urp6+cfy+OQXi8UwvZubm+GB7K3X6+Gv6+hXL+Pq6qqqqvV6/U//\n8nd3dzc3N2VZ3p9e/uQ/ffq0WCyGZ8Gh8MlfXV0tl8v95Zubm6PbJEz+cErr9fr6+vrw6iNPyc+f\n/NtdsmuapizLBxu+3W6Hn89ms8CPnJhOp8O7qsVicfgBGWmTL4pis9n80+e+Pz7bvu+Hbxwuy3J4\nU/liHp/8dDodpnf/Bl3XHT60YSvwxZRludls7m+3HVqtVg9ud+ZPfv/BzQ/+NnzybdsO/+aHT+S9\nV5983/eHCxXz+fzwA7Uff0p+/uTfbpAO/2WPdF1XVdVwNfAbNPabzHtt2x4+kPzJH/qyZnvk6KX8\n6Ml5+Fiqqnr1dwb3dV232+0efBaET75t28M/myPhky/L8vAd5OFzeRIw+aOOtm07fJ3Qv/r8yb/p\nfUi73a5pmt1u91prtf/ZfD7/8ccf9y8l9/ch8fLquv6n9zexVqvVI/uWku03MvY7aYqiqKrqC3or\nU5blarXabrfT6bTv++TPV+u6brVa3d+/eDpvdwtpt9tdXFzMZrPlctm27YPf8xRrv7N0uVxWVfXy\nC1kcaZrmi3tb0DTNbDb7gl7Hj+xf0xeLRVVVw272L8W+Q8vlMmFR/Z90XVfX9Wazeck7fbtBatt2\ns9ns39UuFosv6GW9aZphk242mx0t8vLCVqtV3/dfVo0mk8l2u93tdnVd13W9f+l57Rk9w3Q67bpu\nvV7v98HM5/OX39fynzVNs1gs9utgsa88Q41e+C3L212ye2Qzf787blikTnvz1ff94apuURSHf9Ph\nkz/yZc32vsvLy8lk8uDC1+FjOTx2I8ThOsz9L9IMn3xZloczPHoKhE/+6Pl7JGHybduuVqv/UKPP\nn/zb3UIqy3LYsDh6HTzcHXd43EiI+Xw+nIQ0mUxWq9Xh2/PwyR95fLZFUQy7f5+1c/VlXF5eTqfT\noUZHb9IPV2OOjt3IFz752WzW9/0QoaOnQPjkD5+/99+Bvfrkm6Y52jY6/MN+/Cn5+ZN/018/sd1u\nh3++o0N7A88tPXQ489lsdvQ/Pmryw9mLQ2+Odrcknxj7yOTbtj06kKHruqPdv697huPwhZnD5Muy\nvH/kxX7Jbr8z9ejn+wuZk3/83PDwySefGPvhw4f768+HfxsnPTH2TQcJgBxvd8kOgCiCBEAEQQIg\ngiABEEGQAIggSABEECQAIggSABEECYAIggRAhP8HJUscDjw38SYAAAAASUVORK5CYII=\n",
"text/plain": [
- ""
+ ""
]
},
"metadata": {},
@@ -349,7 +203,7 @@
},
{
"cell_type": "code",
- "execution_count": 6,
+ "execution_count": 3,
"metadata": {
"collapsed": false,
"slideshow": {
@@ -361,39 +215,21 @@
"name": "stdout",
"output_type": "stream",
"text": [
- "'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/course_git_F2017/06_roots and optimization/incsearch.m\n",
- "\n",
- " incsearch: incremental search root locator\n",
- " xb = incsearch(func,xmin,xmax,ns):\n",
- " finds brackets of x that contain sign changes\n",
- " of a function on an interval\n",
- " input:\n",
- " func = name of function\n",
- " xmin, xmax = endpoints of interval\n",
- " ns = number of subintervals (default = 50)\n",
- " output:\n",
- " xb(k,1) is the lower bound of the kth sign change\n",
- " xb(k,2) is the upper bound of the kth sign change\n",
- " If no brackets found, xb = [].\n",
- "\n",
- "\n",
- "Additional help for built-in functions and operators is\n",
- "available in the online version of the manual. Use the command\n",
- "'doc ' to search the manual index.\n",
+ "number of brackets: 1\n",
+ "ans =\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"
+ " 141.84 144.90\n",
+ "\n"
]
}
],
"source": [
- "help incsearch"
+ "incsearch(f_m,50,200)"
]
},
{
"cell_type": "code",
- "execution_count": 8,
+ "execution_count": 4,
"metadata": {
"collapsed": false,
"scrolled": true,
@@ -627,7 +463,7 @@
},
{
"cell_type": "code",
- "execution_count": 21,
+ "execution_count": 9,
"metadata": {
"collapsed": false,
"scrolled": true,
@@ -764,13 +600,13 @@
"\tgnuplot_plot_3a\n",
"\n",
"\t \n",
- "\t\n",
+ "\t\n",
"\n",
"\t\n",
"\tgnuplot_plot_4a\n",
"\n",
- "\t\t \n",
- "\t\n",
+ "\t\t \n",
+ "\t\n",
"\n",
"\t\n",
"\n",
diff --git a/06_roots-1/octave-workspace b/06_roots-1/octave-workspace
index 05454da..3a964d4 100644
Binary files a/06_roots-1/octave-workspace and b/06_roots-1/octave-workspace differ
diff --git a/07_roots-2/octave-workspace b/07_roots-2/octave-workspace
index 70eb62b..8c437bb 100644
Binary files a/07_roots-2/octave-workspace and b/07_roots-2/octave-workspace differ
diff --git a/08_optimization/.ipynb_checkpoints/lecture_08-checkpoint.ipynb b/08_optimization/.ipynb_checkpoints/lecture_08-checkpoint.ipynb
new file mode 100644
index 0000000..fa0f10b
--- /dev/null
+++ b/08_optimization/.ipynb_checkpoints/lecture_08-checkpoint.ipynb
@@ -0,0 +1,2836 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Optimization"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Lennard-Jones potential"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "The Lennard-Jones potential is commonly used to model interatomic bonding. \n",
+ "\n",
+ "$E_{LJ}(x)=4\\epsilon \\left(\\left(\\frac{\\sigma}{x}\\right)^{12}-\\left(\\frac{\\sigma}{x}\\right)^{6}\\right)$\n",
+ "\n",
+ "Considering a 1-D gold chain, we can calculate the bond length, $x_{b}$, with no force applied to the chain and even for tension, F. This will allow us to calculate the nonlinear spring constant of a 1-D gold chain. \n",
+ "\n",
+ "\n",
+ "\n",
+ "Computational Tools to Study and Predict the Long-Term Stability of Nanowires.\n",
+ "By Martin E. Zoloff Michoff, Patricio Vélez, Sergio A. Dassie and Ezequiel P. M. Leiva "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "\n",
+ "\n",
+ "[Single atom gold chain mechanics](http://www.uam.es/personal_pdi/ciencias/agrait/)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "### First, let's find the minimum energy $\\min(E_{LJ}(x))$\n",
+ "\n",
+ "## Brute force"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "setdefaults\n",
+ "epsilon = 0.039; % kcal/mol\n",
+ "sigma = 2.934; % Angstrom\n",
+ "x=linspace(2.8,6,200); % bond length in Angstrom\n",
+ "\n",
+ "Ex = lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "[Emin,imin]=min(Ex);\n",
+ "\n",
+ "plot(x,Ex,x(imin),Emin,'o')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Golden Search Algorithm\n",
+ "\n",
+ "We can't just look for a sign change for the problem (unless we can take a derivative) so we need a new approach to determine whether we have a maximum between the two bounds.\n",
+ "\n",
+ "Rather than using the midpoint of initial bounds, here the problem is more difficult. We need to compare the values of 4 function evaluations. The golden search uses the golden ratio to determine two interior points. \n",
+ "\n",
+ "\n",
+ "\n",
+ "Start with bounds of 2.5 and 6 Angstrom. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "current_min = -0.019959\r\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% define Au atomic potential\n",
+ "epsilon = 0.039; % kcal/mol\n",
+ "sigma = 2.934; % Angstrom\n",
+ "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "% calculate golden ratio\n",
+ "phi = 1/2+sqrt(5)/2;\n",
+ "% set initial limits\n",
+ "x_l=2.8; \n",
+ "x_u=6; \n",
+ "\n",
+ "% Iteration #1\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\n",
+ "\n",
+ "Gnuplot\n",
+ "Produced by GNUPLOT 5.0 patchlevel 3 \n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\tgnuplot_plot_1a\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_2a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_3a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_4a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_6a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Iteration #2\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\n",
+ "\n",
+ "Gnuplot\n",
+ "Produced by GNUPLOT 5.0 patchlevel 3 \n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\tgnuplot_plot_1a\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_2a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_3a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_4a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_6a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Iteration #3\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\n",
+ "\n",
+ "Gnuplot\n",
+ "Produced by GNUPLOT 5.0 patchlevel 3 \n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\tgnuplot_plot_1a\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_2a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_3a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_4a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_6a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Iteration #3\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\n",
+ "\n",
+ "Gnuplot\n",
+ "Produced by GNUPLOT 5.0 patchlevel 3 \n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\tgnuplot_plot_1a\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_2a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_3a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_4a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% define Au atomic potential\n",
+ "epsilon = 0.039; % kcal/mol\n",
+ "sigma = 2.934; % Angstrom\n",
+ "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "% set initial limits\n",
+ "x_l=2.8; \n",
+ "x_u=3.5; \n",
+ "\n",
+ "% Iteration #1\n",
+ "x1=x_l;\n",
+ "x2=mean([x_l,x_u]);\n",
+ "x3=x_u;\n",
+ "\n",
+ "% evaluate Au_x(x1), Au_x(x2) and Au_x(x3)\n",
+ " \n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "f3=Au_x(x3);\n",
+ "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n",
+ "x_fit = linspace(x1,x3,20);\n",
+ "y_fit = polyval(p,x_fit);\n",
+ "\n",
+ "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n",
+ "hold on\n",
+ "if f2x2\n",
+ " plot(x4,f4,'*',[x1,x2],[f1,f2])\n",
+ " x1=x2;\n",
+ " f1=f2;\n",
+ " else\n",
+ " plot(x4,f4,'*',[x3,x2],[f3,f2])\n",
+ " x3=x2;\n",
+ " f3=f2;\n",
+ " end\n",
+ " x2=x4; f2=f4;\n",
+ "else\n",
+ " error('no minimum in bracket')\n",
+ "end\n",
+ "hold off"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n",
+ "x_fit = linspace(x1,x3,20);\n",
+ "y_fit = polyval(p,x_fit);\n",
+ "\n",
+ "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n",
+ "hold on\n",
+ "if f2x2\n",
+ " plot(x4,f4,'*',[x1,x2],[f1,f2])\n",
+ " x1=x2;\n",
+ " f1=f2;\n",
+ " else\n",
+ " plot(x4,f4,'*',[x3,x2],[f3,f2])\n",
+ " x3=x2;\n",
+ " f3=f2;\n",
+ " end\n",
+ " x2=x4; f2=f4;\n",
+ "else\n",
+ " error('no minimum in bracket')\n",
+ "end\n",
+ "hold off\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n",
+ "x_fit = linspace(x1,x3,20);\n",
+ "y_fit = polyval(p,x_fit);\n",
+ "\n",
+ "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n",
+ "hold on\n",
+ "if f2x2\n",
+ " plot(x4,f4,'*',[x1,x2],[f1,f2])\n",
+ " x1=x2;\n",
+ " f1=f2;\n",
+ " else\n",
+ " plot(x4,f4,'*',[x3,x2],[f3,f2])\n",
+ " x3=x2;\n",
+ " f3=f2;\n",
+ " end\n",
+ " x2=x4; f2=f4;\n",
+ "else\n",
+ " error('no minimum in bracket')\n",
+ "end\n",
+ "hold off\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "warning: Matlab-style short-circuit operation performed for operator |\n",
+ "warning: called from\n",
+ " goldmin at line 17 column 1\n",
+ "warning: Matlab-style short-circuit operation performed for operator |\n",
+ "warning: Matlab-style short-circuit operation performed for operator |\n",
+ "xmin = 0.32933\n",
+ "Emin = -2.7096e-04\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "epsilon = 0.039; % kcal/mol\n",
+ "epsilon = epsilon*6.9477e-21; % J/atom\n",
+ "epsilon = epsilon*1e18; % aJ/J\n",
+ "% final units for epsilon are aJ\n",
+ "\n",
+ "sigma = 2.934; % Angstrom\n",
+ "sigma = sigma*0.10; % nm/Angstrom\n",
+ "x=linspace(2.8,6,200)*0.10; % bond length in um\n",
+ "\n",
+ "Ex = lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "%[Emin,imin]=min(Ex);\n",
+ "\n",
+ "[xmin,Emin] = goldmin(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)\n",
+ "\n",
+ "plot(x,Ex,xmin,Emin,'o')\n",
+ "ylabel('Lennard Jones Potential (aJ/atom)')\n",
+ "xlabel('bond length (nm)')\n",
+ "\n",
+ "Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=50;\n",
+ "warning('off')\n",
+ "dx = zeros(1,N); % [in nm]\n",
+ "F_applied=linspace(0,0.0022,N); % [in nN]\n",
+ "for i=1:N\n",
+ " optmin=goldmin(@(dx) Etotal(dx,F_applied(i)),-0.001,0.035);\n",
+ " dx(i)=optmin;\n",
+ "end\n",
+ "\n",
+ "plot(dx,F_applied)\n",
+ "xlabel('dx (nm)')\n",
+ "ylabel('Force (nN)')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "F =\n",
+ "\n",
+ "@(dx) 4 * epsilon * 6 * (sigma ^ 6 ./ (xmin + dx) .^ 7 - 2 * sigma ^ 12 ./ (xmin + dx) .^ 13)\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "dx_full=linspace(0,0.06,50);\n",
+ "F= @(dx) 4*epsilon*6*(sigma^6./(xmin+dx).^7-2*sigma^12./(xmin+dx).^13)\n",
+ "plot(dx_full,F(dx_full),dx,F_applied)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2)"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Slideshow",
+ "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/08_optimization/Auchain_model.gif b/08_optimization/Auchain_model.gif
new file mode 100644
index 0000000..2f8c78c
Binary files /dev/null and b/08_optimization/Auchain_model.gif differ
diff --git a/08_optimization/Auchain_model.png b/08_optimization/Auchain_model.png
new file mode 100644
index 0000000..dce9cd7
Binary files /dev/null and b/08_optimization/Auchain_model.png differ
diff --git a/08_optimization/au_chain.jpg b/08_optimization/au_chain.jpg
new file mode 100644
index 0000000..492fee8
Binary files /dev/null and b/08_optimization/au_chain.jpg differ
diff --git a/08_optimization/goldenratio.png b/08_optimization/goldenratio.png
new file mode 100644
index 0000000..9493c8c
Binary files /dev/null and b/08_optimization/goldenratio.png differ
diff --git a/08_optimization/goldmin.m b/08_optimization/goldmin.m
new file mode 100644
index 0000000..fb4ce0b
--- /dev/null
+++ b/08_optimization/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/08_optimization/lecture_08.ipynb b/08_optimization/lecture_08.ipynb
new file mode 100644
index 0000000..fa0f10b
--- /dev/null
+++ b/08_optimization/lecture_08.ipynb
@@ -0,0 +1,2836 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Optimization"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Lennard-Jones potential"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "The Lennard-Jones potential is commonly used to model interatomic bonding. \n",
+ "\n",
+ "$E_{LJ}(x)=4\\epsilon \\left(\\left(\\frac{\\sigma}{x}\\right)^{12}-\\left(\\frac{\\sigma}{x}\\right)^{6}\\right)$\n",
+ "\n",
+ "Considering a 1-D gold chain, we can calculate the bond length, $x_{b}$, with no force applied to the chain and even for tension, F. This will allow us to calculate the nonlinear spring constant of a 1-D gold chain. \n",
+ "\n",
+ "\n",
+ "\n",
+ "Computational Tools to Study and Predict the Long-Term Stability of Nanowires.\n",
+ "By Martin E. Zoloff Michoff, Patricio Vélez, Sergio A. Dassie and Ezequiel P. M. Leiva "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "\n",
+ "\n",
+ "[Single atom gold chain mechanics](http://www.uam.es/personal_pdi/ciencias/agrait/)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "### First, let's find the minimum energy $\\min(E_{LJ}(x))$\n",
+ "\n",
+ "## Brute force"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "setdefaults\n",
+ "epsilon = 0.039; % kcal/mol\n",
+ "sigma = 2.934; % Angstrom\n",
+ "x=linspace(2.8,6,200); % bond length in Angstrom\n",
+ "\n",
+ "Ex = lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "[Emin,imin]=min(Ex);\n",
+ "\n",
+ "plot(x,Ex,x(imin),Emin,'o')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Golden Search Algorithm\n",
+ "\n",
+ "We can't just look for a sign change for the problem (unless we can take a derivative) so we need a new approach to determine whether we have a maximum between the two bounds.\n",
+ "\n",
+ "Rather than using the midpoint of initial bounds, here the problem is more difficult. We need to compare the values of 4 function evaluations. The golden search uses the golden ratio to determine two interior points. \n",
+ "\n",
+ "\n",
+ "\n",
+ "Start with bounds of 2.5 and 6 Angstrom. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "current_min = -0.019959\r\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% define Au atomic potential\n",
+ "epsilon = 0.039; % kcal/mol\n",
+ "sigma = 2.934; % Angstrom\n",
+ "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "% calculate golden ratio\n",
+ "phi = 1/2+sqrt(5)/2;\n",
+ "% set initial limits\n",
+ "x_l=2.8; \n",
+ "x_u=6; \n",
+ "\n",
+ "% Iteration #1\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\n",
+ "\n",
+ "Gnuplot\n",
+ "Produced by GNUPLOT 5.0 patchlevel 3 \n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\tgnuplot_plot_1a\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_2a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_3a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_4a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_6a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Iteration #2\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\n",
+ "\n",
+ "Gnuplot\n",
+ "Produced by GNUPLOT 5.0 patchlevel 3 \n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\tgnuplot_plot_1a\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_2a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_3a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_4a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_6a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Iteration #3\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\n",
+ "\n",
+ "Gnuplot\n",
+ "Produced by GNUPLOT 5.0 patchlevel 3 \n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\tgnuplot_plot_1a\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_2a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_3a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_4a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_6a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Iteration #3\n",
+ "d=(phi-1)*(x_u-x_l);\n",
+ "\n",
+ "x1=x_l+d; % define point 1\n",
+ "x2=x_u-d; % define point 2\n",
+ "\n",
+ "% evaluate Au_x(x1) and Au_x(x2)\n",
+ "\n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n",
+ "hold on;\n",
+ "\n",
+ "if f2\n",
+ "\n",
+ "Gnuplot\n",
+ "Produced by GNUPLOT 5.0 patchlevel 3 \n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t-0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.02\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.04\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.06\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t0.08\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t2.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t3.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t5.5\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t6\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\tgnuplot_plot_1a\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_2a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_3a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\t\n",
+ "\t\n",
+ "\tgnuplot_plot_4a\n",
+ "\n",
+ "\t \n",
+ "\t\n",
+ "\n",
+ "\t\n",
+ "\tgnuplot_plot_5a\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% define Au atomic potential\n",
+ "epsilon = 0.039; % kcal/mol\n",
+ "sigma = 2.934; % Angstrom\n",
+ "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "% set initial limits\n",
+ "x_l=2.8; \n",
+ "x_u=3.5; \n",
+ "\n",
+ "% Iteration #1\n",
+ "x1=x_l;\n",
+ "x2=mean([x_l,x_u]);\n",
+ "x3=x_u;\n",
+ "\n",
+ "% evaluate Au_x(x1), Au_x(x2) and Au_x(x3)\n",
+ " \n",
+ "f1=Au_x(x1);\n",
+ "f2=Au_x(x2);\n",
+ "f3=Au_x(x3);\n",
+ "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n",
+ "x_fit = linspace(x1,x3,20);\n",
+ "y_fit = polyval(p,x_fit);\n",
+ "\n",
+ "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n",
+ "hold on\n",
+ "if f2x2\n",
+ " plot(x4,f4,'*',[x1,x2],[f1,f2])\n",
+ " x1=x2;\n",
+ " f1=f2;\n",
+ " else\n",
+ " plot(x4,f4,'*',[x3,x2],[f3,f2])\n",
+ " x3=x2;\n",
+ " f3=f2;\n",
+ " end\n",
+ " x2=x4; f2=f4;\n",
+ "else\n",
+ " error('no minimum in bracket')\n",
+ "end\n",
+ "hold off"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n",
+ "x_fit = linspace(x1,x3,20);\n",
+ "y_fit = polyval(p,x_fit);\n",
+ "\n",
+ "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n",
+ "hold on\n",
+ "if f2x2\n",
+ " plot(x4,f4,'*',[x1,x2],[f1,f2])\n",
+ " x1=x2;\n",
+ " f1=f2;\n",
+ " else\n",
+ " plot(x4,f4,'*',[x3,x2],[f3,f2])\n",
+ " x3=x2;\n",
+ " f3=f2;\n",
+ " end\n",
+ " x2=x4; f2=f4;\n",
+ "else\n",
+ " error('no minimum in bracket')\n",
+ "end\n",
+ "hold off\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n",
+ "x_fit = linspace(x1,x3,20);\n",
+ "y_fit = polyval(p,x_fit);\n",
+ "\n",
+ "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n",
+ "hold on\n",
+ "if f2x2\n",
+ " plot(x4,f4,'*',[x1,x2],[f1,f2])\n",
+ " x1=x2;\n",
+ " f1=f2;\n",
+ " else\n",
+ " plot(x4,f4,'*',[x3,x2],[f3,f2])\n",
+ " x3=x2;\n",
+ " f3=f2;\n",
+ " end\n",
+ " x2=x4; f2=f4;\n",
+ "else\n",
+ " error('no minimum in bracket')\n",
+ "end\n",
+ "hold off\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "warning: Matlab-style short-circuit operation performed for operator |\n",
+ "warning: called from\n",
+ " goldmin at line 17 column 1\n",
+ "warning: Matlab-style short-circuit operation performed for operator |\n",
+ "warning: Matlab-style short-circuit operation performed for operator |\n",
+ "xmin = 0.32933\n",
+ "Emin = -2.7096e-04\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "epsilon = 0.039; % kcal/mol\n",
+ "epsilon = epsilon*6.9477e-21; % J/atom\n",
+ "epsilon = epsilon*1e18; % aJ/J\n",
+ "% final units for epsilon are aJ\n",
+ "\n",
+ "sigma = 2.934; % Angstrom\n",
+ "sigma = sigma*0.10; % nm/Angstrom\n",
+ "x=linspace(2.8,6,200)*0.10; % bond length in um\n",
+ "\n",
+ "Ex = lennard_jones(x,sigma,epsilon);\n",
+ "\n",
+ "%[Emin,imin]=min(Ex);\n",
+ "\n",
+ "[xmin,Emin] = goldmin(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)\n",
+ "\n",
+ "plot(x,Ex,xmin,Emin,'o')\n",
+ "ylabel('Lennard Jones Potential (aJ/atom)')\n",
+ "xlabel('bond length (nm)')\n",
+ "\n",
+ "Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=50;\n",
+ "warning('off')\n",
+ "dx = zeros(1,N); % [in nm]\n",
+ "F_applied=linspace(0,0.0022,N); % [in nN]\n",
+ "for i=1:N\n",
+ " optmin=goldmin(@(dx) Etotal(dx,F_applied(i)),-0.001,0.035);\n",
+ " dx(i)=optmin;\n",
+ "end\n",
+ "\n",
+ "plot(dx,F_applied)\n",
+ "xlabel('dx (nm)')\n",
+ "ylabel('Force (nN)')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "F =\n",
+ "\n",
+ "@(dx) 4 * epsilon * 6 * (sigma ^ 6 ./ (xmin + dx) .^ 7 - 2 * sigma ^ 12 ./ (xmin + dx) .^ 13)\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "dx_full=linspace(0,0.06,50);\n",
+ "F= @(dx) 4*epsilon*6*(sigma^6./(xmin+dx).^7-2*sigma^12./(xmin+dx).^13)\n",
+ "plot(dx_full,F(dx_full),dx,F_applied)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2)"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Slideshow",
+ "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/08_optimization/lennard_jones.m b/08_optimization/lennard_jones.m
new file mode 100644
index 0000000..d18a6ad
--- /dev/null
+++ b/08_optimization/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/08_optimization/octave-workspace b/08_optimization/octave-workspace
new file mode 100644
index 0000000..8c437bb
Binary files /dev/null and b/08_optimization/octave-workspace differ
diff --git a/HW2/README.md b/HW2/README.md
new file mode 100644
index 0000000..a3b064d
--- /dev/null
+++ b/HW2/README.md
@@ -0,0 +1,99 @@
+# Homework #2
+## due 10/6/17 by 11:59pm
+
+
+**1\.** Create a new github repository called '02_roots_and_optimization'.
+
+ a. Add rcc02007 and zhs15101 as collaborators.
+
+ b. submit the clone repository URL to:
+ [https://goo.gl/forms/svFKpfiCfLO9Zvfz1](https://goo.gl/forms/svFKpfiCfLO9Zvfz1)
+
+
+**2\.** You're installing a powerline in a residential neighborhood. The lowest point on the
+cable is 30 m above the ground, but 30 m away is a tree that is 35 m tall. Another
+engineer informs you that this is a catenary cable problem with the following solution
+
+ $y(x)=\frac{T}{w}\cosh\left(\frac{w}{T}x\right)+y_{0}-\frac{T}{w}$.
+
+ where y(x) is the height of the cable at a distance, x, from the lowest point, $y_{0}$,
+ T is the tension in the cable, and w is the weight per unit length of the cable. Your
+ supervisor wants to know which numerical solver to use when they have to install these
+ powerlines in similar places.
+
+ a. Use the three solvers `falsepos.m`, `bisect.m`, and `mod_secant.m`
+ to solve for the tension neededi, T, to reach y(30 m)=35 m, with w=10 N/m, and $y_{0}$=30 m.
+
+ b. 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 | | | |
+ |mod_secant | | | |
+ |bisect | | | |
+ ```
+
+
+ c. Add a figure to your README that plots the final shape of the powerline
+ ($y(x)~vs.~x$) from x=-10 to 50 m.
+
+**3\.** 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-1)*exp(-(x-1)^2). The root is at 1, 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=3 for f(x)=(x-1)*exp(-(x-1)^2).
+
+ b. Add the results to a table in the `README.md` with:
+
+ ```
+ ### divergence of Newton-Raphson method
+
+ | iteration | x_i | approx error |
+ | --- | --- | --- |
+ | 0 | 3 | n/a |
+ | 1 | | |
+ | 2 | | |
+ | 3 | | |
+ | 4 | | |
+ | 5 | | |
+ ```
+
+ c. Repeat steps a-b for an initial guess of 1.2. (But change the heading from
+ 'divergence' to 'convergence')
+
+
+
+**4\.** Determine the nonlinear spring constants of a single-atom gold chain. You can assume
+the gold atoms are aligned in a one dimensional network and the potential energy is
+described by the Lennard-Jones potential as such
+
+ $E_{LJ}(x)=4\epsilon
+ \left(\left(\frac{\sigma}{x}\right)^{12}-\left(\frac{\sigma}{x}\right)^{6}\right)$.
+
+ Where x is the distance between atoms in nm, $\epsilon$=2.71E-4 aJ, and $\sigma$=0.2934
+ nm. The energy term that must be minimized is
+
+ $E_{total}(\Delta x)=E_{LJ}(x_{0}+\Delta x)-F\Delta x$.
+
+ Where $x_{0}$ is the distance between atoms with no force applied and $\Delta x$ is the
+ amount each gold atom has moved under a given force, F.
+
+ a. Determine $x_{0}$ when F=0 nN using the golden ratio and parabolic methods. *Show
+ your script and output in your README and include your functions*
+
+ b. Solve for $\Delta x$ for F=0 to 0.0022 nN with 30 steps. *Use the golden ratio
+ solver or the matlab/octave `fminsearch`
+
+ c. create a sum of squares error function `sse_of_parabola.m` that calculates the sum of
+ squares error between a function $F(x)=K_{1}\Delta x+1/2K_{2}\Delta x^{2}$ and the
+ Forces used in part B for each $\Delta x$.
+
+ d. Use the `fminsearch` matlab/octave function to determine $K_{1}$ and $K_{2}$.
+
+ e. Plot the force vs calculated $\Delta x$ and the best-fit parabola using $K_{1}$ and
+ $K_{2}$ in part d.
+