From f541b45704f9748b2f361d40d4db90079734156f Mon Sep 17 00:00:00 2001 From: "Ryan C. Cooper" Date: Fri, 8 Sep 2017 09:41:12 -0400 Subject: [PATCH] added freefall.m to HW1 --- 01_Introduction/octave-workspace | Bin 849 -> 153 bytes 03_Intro to matlab-octave/octave-workspace | Bin 0 -> 7355 bytes 04_intro to m-files/octave-workspace | Bin 7319 -> 153 bytes ...06_roots and optimization-checkpoint.ipynb | 913 +++++++++++++++++ .../06_roots and optimization.ipynb | 934 ++++++++++++++++++ 06_roots and optimization/bisect.m | 37 + 06_roots and optimization/falsepos.m | 39 + 06_roots and optimization/incsearch.m | 37 + 06_roots and optimization/lecture_06.md | 373 +++++++ 06_roots and optimization/octave-workspace | Bin 0 -> 1332 bytes HW1/.README.md.swp | Bin 12288 -> 0 bytes HW1/README.md | 6 +- HW1/freefall.m | 21 + 13 files changed, 2357 insertions(+), 3 deletions(-) create mode 100644 03_Intro to matlab-octave/octave-workspace create mode 100644 06_roots and optimization/.ipynb_checkpoints/06_roots and optimization-checkpoint.ipynb create mode 100644 06_roots and optimization/06_roots and optimization.ipynb create mode 100644 06_roots and optimization/bisect.m create mode 100644 06_roots and optimization/falsepos.m create mode 100644 06_roots and optimization/incsearch.m create mode 100644 06_roots and optimization/lecture_06.md create mode 100644 06_roots and optimization/octave-workspace delete mode 100644 HW1/.README.md.swp create mode 100644 HW1/freefall.m diff --git a/01_Introduction/octave-workspace b/01_Introduction/octave-workspace index cded7e4bdb6b5ea7b33e75fa5c7ceaaf0dad4467..8c437bb6e55a5d1b6115661b3a20e86870909d32 100644 GIT binary patch delta 19 acmcb}Hj|OtKe;5aELGP~*N0)^)jj|}F$Sst literal 849 zcmeZIE=ep))iu=hVPIrnVDJTE28RD^Kq5D>q$soE-~a#r*?~M}AO_LMzyZP$fYK6B zS^-L{Kxqvq&A_0!>FCR&bB!F1mbUr2SUEUsO20QtJ;=qut8T@bh#WVE+tGJJj^wyI zY(Jm%blJj2$fJ{Jbi|;^c7Q!ihPvvfLbICBM@9ebe26$NlGXemxHd z16z+NSM59;K#on!D+ZE4zpwy_ZazT=BDPA*Z~EAYB_;K zabbKh(C2ySFbyCB(sD|&Nzw!gtMt6mf}H#kkfq5$3=P-f5vpVG*gqa!h#@fsyD^q|#)?WbYFsq-u8R$fB1<+=qsG{ifE^@= zhyi9OBF&irh9Y%3FqEN)Xrf4UfN0$JzV{~Q?C#k)XZX(i<-Pmv{qA?~do0|rcip~C zBS(*1VAI3K#%4PG+St78MgQzrx7T-@{~s@3zORU8LxP_+V+2~qwT@?fZj`t`PT)rZt@c^xwT@?fZjyLznn0@^ z)^V-#TF0|K_X%NPoI$_TgUp{BHu?ClAhOFn>9;@qPGJ{EZt~f*W7E#P)93}^ZNClc zeAfB)UU??7!+m9&q@cd1P&a+O;5~ zklbp~dL7D@lRZ-%N^K`7$XiE#`^>*h226I&RvdUVeDvDMW{_#)y1nbzxc|-Qne@_s z7=efULkKWX;5P!l75JUNlLF5Q{7K*ifmZ}x6L?+VO@Yw@WdaigCJRg#s1&FYm?yAM zpjx0_V41*5fi(i}3w$W>vA|}5?E*UmJ{M>c%Kc#{&|Y9afsO(P2y_xSM4*d6H-YW~ zKM*)d;5dOF37jI32%IT!uD}HX7YSS{&|Bb2folbB6u4EOkHFmm_X_kE7%1?Y&_BtU zksmf|Syeot9IQCdoB&de9-asi9!^gn=7p!ue|kERtk^i&T+(=(tWNdpx4tBaY^I+F zyOIdyX!7bOcS!${@mCjqkWAjAIcDfiAuQ9b>`x_M{^(?v6q!bDTvoZL22iN~ia zbJGarvc1o)bn*#h|M~@Y@%+lST#olP4^omh%8M;yrz?qM{OZ({FO}rN&esmh z%9SK|#tENDM%wU7N6xLot%3TC<&UEMh4j}>mSpe%GcwN%QW|E?vBMF{SA)&o37_>{w=}#B)tMl!S)q;O)o^ za_aof*KFU;grCO0nMEkQ7EV2sMcQBJ{XK0|BzOMA<8iqv_E+6v_(iLF5BX&iwC)%lqP`dUU5n71%6PgN%!T-_lVS|bZUG9vE;9}DK zukGn>hf2u1)KM$ziZP!|RZN&hXBQU{O6#ZlKQD%V54jYB7kMK}$S$hlx0aTWpJyet zof=z0sFych&ME=#-c3{E{i&-pr1AHtmcb|0WOh>93w=)w`E9Cmr(TlNKGBONqSsoDS>r zTc8Jz^r3pJZ+WmDe0lFp1N<7~p(oqZ_x`A~H6V{3xf@V~mIvvH#6e-zxzMv5ywSO^e9fPz4e5q^Rn4b`F&Q*SYD5Ey8<-Vi|+(ZFhJuOm;D z*TUhycgKb!P6vBl2Y+v`x(^Si8^>7(F1e@%8^(9P)1X&ocPESuex&vRC5B%o9G1zuU{;?|$7egwi9j*Qi*G z8*C%yoi0VbEp(5;I{eF{5syiMG0=^Iu~MF2!BX-(szv+C?pUm+k1iJc2^}c~A2Rz$ zN$i~3fS=Y%2`m0v?PQ3fcbtr{W$ww9Azo)a`8Y+lq=corTat`Wx%BM!NQ(LTy2-#R zuY4Ke_u+RkLgizm?@cLm{Ko7U=v2k^SllnRm2rKKmBQ|n9a7BWcvs4GZg33jJw8K5 zc2I(k9ufz>H!O+cdg&PVXOqKtcVT-hDF~aD zNvJdZ5|iP_ac!yaqkWjr_lt7q#XDX~=wyI>2IuF~4Dco(JQI2p<(&y0PWwKS=Y3=r z^rG--7S`K*OvUxAPz9b3i&7zet!XMWcOxcb^E_FVjkvWe&p|w{uHyL0FFDBD9-rqT z4t>vwrciYc&0nW&KGyMkRX*nbp5%kK4KMOJzrW+U*CFFN*A`Yp?(LI@Ubn1caKtJb-Qw^260`!i{rg{8pL6uP6NIl?bIMI)cZBC zBj{%hp}nq0#Ayx2Jz%Th>p83e&u`NAn0|GpXrK>0?`i&Qe$f}W6rdlGo+v_m@*)e# zSG3??>L@}yPd_Ulv<&b5%8%>xeGjgm7vEOHuY~G0t0#Sjt5G*H&6$^_^#0-o&oFi=qv1`dc?(fnI3uchQt7!aGGX7 z-nCU3c;7(nN8eFtM7>%vtQ6~g8fJtJ&Dv%JzqFo4@OTof3rx?y_cWm1{-d`6`cheE zfFCYA)T0hf>t%pFsg-)f*}k71`qI5fOK3YX@Ia9cy8ft(+rNFf9=e^7#p_CbqaJm* zvQ`gX^eEB)xhZ1(SWd#rZ4tX}K zT8qABL4}U5V}+jcY@eR%Gxaa3pnIVf6P|C*Gog<8*Hl8++up7uUsHkIaHRr$ZnCTrI_=$2MQ9sb+Z|Mm zdg$k9hRzwH_`J!*X5>w0r5Sb0VRa4S`B88U`rJl`TG;v6qZWQ2V_OTmR&>=s2TC5- zpwBICGh^K|8_dX;nfYei_ghs1|E`TSgSWCUGxYP*Yi5jp=3+HrW%NK{1$eftq>9&b z=W6gNX;d}xczJX+`uvgy)tEp3JooEo^f~KhnfW_$jS_v>UEd7Y>$@ofyzU*If&8C6 znd5|mndmp%-^xN=_MX7?+h%tb^0;`23VbLGRS{Yq|37uHX^?r2;^ef0qje7CjObzPnr*2yGr3FzM^wXh~ z&#W)~B}9vOom#9zTxEmwh}WVH9s22fw(Yas%|L_X7i-KUK>0Zrn+p*D@pT2XDFKxnO4ON$7 z{U1A*A&+NVEdwtdJ`$j@k} zO6b(&N-#0sRq0=AF&4XRe!S~Rw^r_58yzA|Y;4k7?44%BThwmQ_FH@s# zFZSZ~YV=4g>>M1gK|T?U20Ao0lJjHJ7dq(8OjkX4U-v=>9s6L44t)HkzmD(UuWJ$S zcT=^fWA{hvu)glAI;?w)m!A94qDOt$Szv(Q1KSM9=f?X6=|Bize$AU#iv31WRVn&G!-z8IXzyP*uK9a8>Yt&h9C>y8RyoGsHPHmW^qpnmeZi;- z#D8E2*V*^_R+19h<_CRP2|t`IHDMj9S50Io+d%o4pcB0ZnqZHb@{S#{?Rm=tJDpdV zu%5ys6Zh||3apnMuCaaA+&S^MuP;u(zTrbtBJ^{`dVank|1k;r>pwLGx;x*Qp99!S z(xJ~!ho(bk=ljaR?g6VF~{5ul1ANwphM1sTJSvd-}uk%{3RN(yK*ii}Fmu;LZ#UZYJbfAQt z5@y>IalNF_Y~|r%e>H_VCPwR`-&y({H8Qq ztzhRf^dv@teQ<j!jUox@fm-DEczrtwYd%)dt*`u=B1YKWF5Tee_V&d7s8E z2|M?3>(wP;=L)fNS`@5bSG)G&`@9TC=<@j0y}+X(D|%ty;63S8C5vnP)xO|m&CP%mM9bYh-S!p`YxyWN!RJWrWFP|5ZygFX5nUX$1S6?_l) z$%*%`FC3J}bC*HbFP`7o6Z^e}MNUebtN0H+8YyyWSlN=EeRa2lC&YW6<|qyl97YvBP+_ zzwLZbgT6vZ`O9VvXsHGN`fsd*PCc%#g?%r>>fyJ!8}6fD*|g<8^7Y#1_p$$(IH>_T zQEdAF`{R#4e}F#8Yege;`smb$ysj>Ig#BR6nqQ%3`NrR1w`100%%?uk1pe08HA7b$ z&NQRH82M{6>dgyz3+nmgfL7GYi2kj}pS*J|7=OpfR`4arz72K$!MAP5@2IDxJHD(1x*nzG@Bj05EAs#7jaJmZ-AS#e1C0yXu)kP5x*c{8OlU{HGv``6 z>|1iR9pkTWX@?#zTiSuXr9;sHUEMb68FYGX`!n?v!sU>%otb%IxEUv{E?=6QA^Zaz_+sN>OYoruGV?q|F{0q_-EBN8}`Bf9K}DNaKXk2cenrFyZ;8XOtyUh literal 0 HcmV?d00001 diff --git a/04_intro to m-files/octave-workspace b/04_intro to m-files/octave-workspace index b8411417e628af07b405b9df2b347680dd9b307a..8c437bb6e55a5d1b6115661b3a20e86870909d32 100644 GIT binary patch delta 19 acmbPkIg^pgKe;5aELGP~*Jttt&VB$t+yoxzwfF1h8>peBRC5GrA^&$rgMR*q)QoO7Pj`Qtp#?C070ZQXwF`~Kdw-gkXp zug0lkXC-)gdxk68C=`mZ(pRC-wwHcVW7RW~vz}|US_7#?zk1SF`qXM#OZHn{r_1pmttF++LiY=4 zEy?#WX)Q@l(CZBBU!5Sa`(sqiTQuoFZ2YB z`!`OILb1dB)5DR4XEh%@D6hZ1>Y^sO$EN)`4QDh>$L`rwTbBw+rpWz-&m&*R>-d>( z#%XTv4_h`oFs@7^J7|T+T-eo&`U`FJX`71(N(1$0Hn%a6H2CM7%VdXAS3v#Vd|`+1Pw*yiIe#tp563AeX7by`qEEk z#-vQ?)|WO#ET901De7te(IR?JoRT(kdP=&QDvZsNetvVQjF-BoNqPyn5q&Pj*AQ7g zBG%FB)4|qFB#HI6g@VFK~XI=8%;5ujho|oFkN8@A2!bh`Gxdj2WK|l z5sA|5s<)kC#MUX%yX02>GA|Xnv-Pg1N!j^b=XGv4&v0IWj`Lf!2ld=x#qhddm1Iq{ z+5_7dZdVji^;O;Q?F4G4C;Iih6rUU6g(t>TPY?8=-O)XQ{j?(x!zdJ^!lK_ z^EyeR53-qlU+kq^YklFFM(yZ@*%u|aKE1JX2=%KsW_}@cr}V)>Cu&z;EPlsJLh6ep z^^_Ms{5*vM`QiOj)XsjeC)Du%&@Z3{_s0TzN`yab=2N@DtW`IZ_q z0PFivKL%pc#$}R%fq2GvBRDmH;sFntUx4XWiVK{tp+o@1EHHzxER6a&2r&!k!a;a< zFXeeKb}@Z}(X9*hdoWgvqdsm}1>+{;Cm22EPz6KaT10sog0;6LwIxH4??Cw)ib{5=9g4flDQ`nD?`z5Z+M(FU z;v0fpDjJv&C^IRKA&6c>H$*6mk5fKFVQxw{L?|BhrwfK*%q$w;VYslF@-Pf}Gz^Wy zkfo$NhoLQtcNh}0si9%mdxY{3j;y-_n=519k^i6wKWmX!TU7u+W!czZHmFK|tp`?G?NUoLPMSZdy?!KZLykU(D-QE8|vk{X&lz7d#H;<3exS z8}hmZ^woldPk*Tt{0d$UaXSk;$m0#S<#ht5Bk!Nh{cXG=k9(NN^Tqmsb#i}u&Y#XB zxxJaU>Fj^}P}&k%hn6BG&PxE;c20hF*2Uu`w%hlk}9`ML%|3#XJU`T8fs?Zz*LOR}FeUj^9o;^!Ps0eQLtGHx5(8N+5kzCyIc z)sye$?eT2~S&uGg!d9Qwg=}Yg94=k&;+?JsS#JmI@ueE|A^WS(OT7VPJ38VcY4(&5 zLr|1pYa_fTp}XBMf^3&Da#%?+fy|33`~xRTCYW+NO;P+8N?0duZzs?~w9Lzl+iiw| z`fj>Yo6NZV=1?(1%()*HparAzIt%WPB`nzjOYWB?*h8M575C2y9oQAdn)_)D7M==g z?yn7=t=5rDvf+N)f_fFLw&niY;-195f1n-b!47mc*qYgMKJ3A6HdXeV7YESjR}?vL zejGugTl^PC&XXf(fy<3>;(R%Qy#{k|=DazB?xx5)&YV9b7|qd2&Z81`W{+kqJgnq= zDsiF6LSn#$^XdY2lgf4B{JJ2LZ3q|6vkFGkoyjWBw+hE-rTtFDc~@Z}+Yqjte^;>Z z<+<`WxPl!Tg&XMUo%P4J+0cVyZ;$#ZhbXriMHGhYB#RP4D@K)nOXh57@apPuk@)9L zxWDRl4LRr2XTQbnU&uEmIas={j`Z%4?09U=Gg9Hd19UuFxRU)xwbOmT!yjp55lT2)y>BX zg7IppXZ60KA;>iDR5`kOC^io|e(s1}C=#M3o$NSa7<|Hy9d+LlhE!LhgSHLfI9Tak zQayY)cFnZgHMi(hgk?QE|72Mt)Sm>NaS9p%a%A2~FU3d{p4)oNqiiJ3oI4;<{2Jam zqWQWiY!rU@e?HJB{G$ zNce2k=CGdQ5TmN!ptc^53yTM>TdXqyP1my5EP5~j=iEPfFXMVNoMU$UZAj%rgf1yu z)V^X8-WVb6gsCxDAGPbi$9H02I?6oYdRi==hbH&4tBZw|U-9nXHz!NC`p3Hx+r=TG z$;*4ria4a6o8Yy?As*L0p0{o3ws_Cdk?fCOynR`^QVsRTT#DN;qoC!%oGwhivj z5^-d}`VZ@Jr@;AIj_n?=B%H{aY?c0f65cfElj|@w8E&h!D>MdE@u{`OaP0P}n7LA| ze`CZnJkV7xOm3Wp85sxv`s2nFB%U9ujSNf09Pb;?zI>dD_6w7iP)`oAYVeqpuqyC`d`j>2g-^)_mi2!nYgo2#B6613 zl71{bDDV}_>TTf7vbhenuq>_vC3B!H&)V%^!?Lw57Sos4e(k}YC-n57W@%FoT9%DF z@GPv4I+lG6FrH=Ijy&6TL^R8?hNx!Q)d=i)E5aBh^p&TB35K&QYJyQLotwguWld9@ zW!bV5*0VI+3BinrPH<#d(F|oQ8=7G$%Yxvu!Suv%XTQ2YR;oscF3m5`>GuRSoyVw9V@5yD4<2h!U630?=1%eurla? zhqRaub3{2cXQv~^vGV1F0NOp@3E$J_!&^?+%F2#2CbA9T3@4^niKmQzC3Z0Vl4Scq ziFcSCO2o0axFCqx;{qkyZ!X*}7hIrsnJO3TWcI00$o7#6@39AA75rEnRS04Bs^G@< wvn%d1yIoPq;_Hgd%zjs_Wc%-hq*t}yEkbnQrT^2l)_)l&@b8}4|Cz?W10_~5RsaA1 diff --git a/06_roots and optimization/.ipynb_checkpoints/06_roots and optimization-checkpoint.ipynb b/06_roots and optimization/.ipynb_checkpoints/06_roots and optimization-checkpoint.ipynb new file mode 100644 index 0000000..bcfbe7b --- /dev/null +++ b/06_roots and optimization/.ipynb_checkpoints/06_roots and optimization-checkpoint.ipynb @@ -0,0 +1,913 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true, + "scrolled": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Roots and Optimization\n", + "## Bracketing ch. 5" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "Given a function, numerical or analytical, it's not always possible a given variable. \n", + "\n", + "### Freefall example:\n", + "If an object, with drag coefficient of 0.25 kg/m reaches a velocity of 36 m/s after 4 seconds of freefalling, what is its mass?\n", + "\n", + "$v(t)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)$\n", + "\n", + "Cannot solve for m \n", + "\n", + "Instead, solve the problem by creating a new function f(m) where\n", + "\n", + "$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$. \n", + "\n", + "When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "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-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\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\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t140\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t160\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t180\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\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": [ + "setdefaults\n", + "g=9.81; % acceleration due to gravity\n", + "m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n", + "c_d=0.25; % drag coefficient\n", + "t=4; % at time = 4 seconds\n", + "v=36; % speed must be 36 m/s\n", + "f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n", + "\n", + "plot(m,f_m(m),m,zeros(length(m),1))\n", + "axis([45 200 -5 1])" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = -0.056985\r\n" + ] + } + ], + "source": [ + "f_m(140)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0\n", + "\n", + "Better methods are the \n", + "1. Bracketing methods\n", + "2. Open methods\n", + "\n", + "Both need an initial guess. \n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Incremental method (Brute force)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function xb = incsearch(func,xmin,xmax,ns)\n", + "% incsearch: incremental search root locator\n", + "% xb = incsearch(func,xmin,xmax,ns):\n", + "% finds brackets of x that contain sign changes\n", + "% of a function on an interval\n", + "% input:\n", + "% func = name of function\n", + "% xmin, xmax = endpoints of interval\n", + "% ns = number of subintervals (default = 50)\n", + "% output:\n", + "% xb(k,1) is the lower bound of the kth sign change\n", + "% xb(k,2) is the upper bound of the kth sign change\n", + "% If no brackets found, xb = [].\n", + "if nargin < 3, error('at least 3 arguments required'), end\n", + "if nargin < 4, ns = 50; end %if ns blank set to 50\n", + "% Incremental search\n", + "x = linspace(xmin,xmax,ns);\n", + "f = func(x);\n", + "nb = 0; xb = []; %xb is null unless sign change detected\n", + "for k = 1:length(x)-1\n", + " if sign(f(k)) ~= sign(f(k+1)) %check for sign change\n", + " nb = nb + 1;\n", + " xb(nb,1) = x(k);\n", + " xb(nb,2) = x(k+1);\n", + " end\n", + "end\n", + "if isempty(xb) %display that no brackets were found\n", + " fprintf('no brackets found\\n')\n", + " fprintf('check interval or increase ns\\n')\n", + "else\n", + " fprintf('number of brackets: %i\\n',nb) %display number of brackets\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/me3255_S2017/course_git/lecture_06/incsearch.m\n", + "\n", + " incsearch: incremental search root locator\n", + " xb = incsearch(func,xmin,xmax,ns):\n", + " finds brackets of x that contain sign changes\n", + " of a function on an interval\n", + " input:\n", + " func = name of function\n", + " xmin, xmax = endpoints of interval\n", + " ns = number of subintervals (default = 50)\n", + " output:\n", + " xb(k,1) is the lower bound of the kth sign change\n", + " xb(k,2) is the upper bound of the kth sign change\n", + " If no brackets found, xb = [].\n", + "\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help incsearch" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "number of brackets: 1\n", + "ans =\n", + "\n", + " 142.42 143.94\n", + "\n" + ] + } + ], + "source": [ + "incsearch(f_m,50, 200,100)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Bisection method" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Divide interval in half until error is reduced to some level\n", + "\n", + "in previous example of freefall, choose x_l=50, x_u=200\n", + "\n", + "x_r = (50+200)/2 = 125\n", + "\n", + "f_m(125) = -0.408\n", + "\n", + "x_r= (125+200)/2 = 162.5\n", + "\n", + "f_m(162.5) = 0.3594\n", + "\n", + "x_r = (125+162.5)/2=143.75\n", + "\n", + "f_m(143.75)= 0.0206" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 0.020577\r\n" + ] + } + ], + "source": [ + "f_m(143.75)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Bisect" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Much better root locator, with 4 iterations, our function is already close to zero\n", + "\n", + "Automate this with a function:\n", + "`bisect.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)\n", + "% bisect: root location zeroes\n", + "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + "% uses bisection method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% xl, xu = lower and upper guesses\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by func\n", + "% output:\n", + "% root = real root\n", + "% fx = function value at root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "test = func(xl,varargin{:})*func(xu,varargin{:});\n", + "if test>0,error('no sign change'),end\n", + "if nargin<4|isempty(es), es=0.0001;end\n", + "if nargin<5|isempty(maxit), maxit=50;end\n", + "iter = 0; xr = xl; ea = 100;\n", + "while (1)\n", + " xrold = xr;\n", + " xr = (xl + xu)/2;\n", + " iter = iter + 1;\n", + " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n", + " test = func(xl,varargin{:})*func(xr,varargin{:});\n", + " if test < 0\n", + " xu = xr;\n", + " elseif test > 0\n", + " xl = xr;\n", + " else\n", + " ea = 0;\n", + " end\n", + " if ea <= es | iter >= maxit,break,end\n", + "end\n", + "root = xr; fx = func(xr, varargin{:});" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'bisect' is a function from the file /home/ryan/Documents/UConn/ME3255/course_git_F2017/06_roots and optimization/bisect.m\n", + "\n", + " bisect: root location zeroes\n", + " [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + " uses bisection method to find the root of func\n", + " input:\n", + " func = name of function\n", + " xl, xu = lower and upper guesses\n", + " es = desired relative error (default = 0.0001%)\n", + " maxit = maximum allowable iterations (default = 50)\n", + " p1,p2,... = additional parameters used by func\n", + " output:\n", + " root = real root\n", + " fx = function value at root\n", + " ea = approximate relative error (%)\n", + " iter = number of iterations\n", + "\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help bisect" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## False position (linear interpolation)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner\n", + "\n", + "$ x_{r} = x_{u} - \\frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "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-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\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\t \n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%xl=50; xu=200; \n", + "xu=xr;\n", + "xl=xl;\n", + "xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu));\n", + "plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## False Position" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Much better root locator, with 4 iterations, our function is already close to zero\n", + "\n", + "Automate this with a function:\n", + "`falsepos.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)\n", + "% falsepos: root location zeroes\n", + "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + "% uses false position method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% xl, xu = lower and upper guesses\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by func\n", + "% output:\n", + "% root = real root\n", + "% fx = function value at root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "test = func(xl,varargin{:})*func(xu,varargin{:});\n", + "if test>0,error('no sign change'),end\n", + "if nargin<4|isempty(es), es=0.0001;end\n", + "if nargin<5|isempty(maxit), maxit=50;end\n", + "iter = 0; xr = xl; ea = 100;\n", + "while (1)\n", + " xrold = xr;\n", + " % xr = (xl + xu)/2; % bisect method\n", + " xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method\n", + " iter = iter + 1;\n", + " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n", + " test = func(xl,varargin{:})*func(xr,varargin{:});\n", + " if test < 0\n", + " xu = xr;\n", + " elseif test > 0\n", + " xl = xr;\n", + " else\n", + " ea = 0;\n", + " end\n", + " if ea <= es | iter >= maxit,break,end\n", + "end\n", + "root = xr; fx = func(xr, varargin{:});" + ] + } + ], + "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/06_roots and optimization/06_roots and optimization.ipynb b/06_roots and optimization/06_roots and optimization.ipynb new file mode 100644 index 0000000..72dfcbf --- /dev/null +++ b/06_roots and optimization/06_roots and optimization.ipynb @@ -0,0 +1,934 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": true, + "scrolled": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Roots and Optimization\n", + "## Bracketing ch. 5" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## Not always possible to solve for a given variable. \n", + "\n", + "### Freefall example:\n", + "If an object, with drag coefficient of 0.25 kg/m reaches a velocity of 36 m/s after 4 seconds of freefalling, what is its mass?\n", + "\n", + "$v(t)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)$\n", + "\n", + "Cannot solve for m \n", + "\n", + "Instead, solve the problem by creating a new function f(m) where\n", + "\n", + "$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$. \n", + "\n", + "When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "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-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\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\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t140\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t160\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t180\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\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": [ + "setdefaults\n", + "g=9.81; % acceleration due to gravity\n", + "m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n", + "c_d=0.25; % drag coefficient\n", + "t=4; % at time = 4 seconds\n", + "v=36; % speed must be 36 m/s\n", + "f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n", + "\n", + "plot(m,f_m(m),m,zeros(length(m),1))\n", + "axis([45 200 -5 1])" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 0.0053580\r\n" + ] + } + ], + "source": [ + "f_m(143)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0\n", + "\n", + "Better methods are the \n", + "1. Bracketing methods\n", + "2. Open methods\n", + "\n", + "Both need an initial guess. \n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Incremental method (Brute force)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function xb = incsearch(func,xmin,xmax,ns)\n", + "% incsearch: incremental search root locator\n", + "% xb = incsearch(func,xmin,xmax,ns):\n", + "% finds brackets of x that contain sign changes\n", + "% of a function on an interval\n", + "% input:\n", + "% func = name of function\n", + "% xmin, xmax = endpoints of interval\n", + "% ns = number of subintervals (default = 50)\n", + "% output:\n", + "% xb(k,1) is the lower bound of the kth sign change\n", + "% xb(k,2) is the upper bound of the kth sign change\n", + "% If no brackets found, xb = [].\n", + "if nargin < 3, error('at least 3 arguments required'), end\n", + "if nargin < 4, ns = 50; end %if ns blank set to 50\n", + "% Incremental search\n", + "x = linspace(xmin,xmax,ns);\n", + "f = func(x);\n", + "nb = 0; xb = []; %xb is null unless sign change detected\n", + "for k = 1:length(x)-1\n", + " if sign(f(k)) ~= sign(f(k+1)) %check for sign change\n", + " nb = nb + 1;\n", + " xb(nb,1) = x(k);\n", + " xb(nb,2) = x(k+1);\n", + " end\n", + "end\n", + "if isempty(xb) %display that no brackets were found\n", + " fprintf('no brackets found\\n')\n", + " fprintf('check interval or increase ns\\n')\n", + "else\n", + " fprintf('number of brackets: %i\\n',nb) %display number of brackets\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/me3255_S2017/course_git/lecture_06/incsearch.m\n", + "\n", + " incsearch: incremental search root locator\n", + " xb = incsearch(func,xmin,xmax,ns):\n", + " finds brackets of x that contain sign changes\n", + " of a function on an interval\n", + " input:\n", + " func = name of function\n", + " xmin, xmax = endpoints of interval\n", + " ns = number of subintervals (default = 50)\n", + " output:\n", + " xb(k,1) is the lower bound of the kth sign change\n", + " xb(k,2) is the upper bound of the kth sign change\n", + " If no brackets found, xb = [].\n", + "\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help incsearch" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "number of brackets: 1\n", + "ans =\n", + "\n", + " 142.42 143.94\n", + "\n" + ] + } + ], + "source": [ + "incsearch(f_m,50, 200,100)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Bisection method" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Divide interval in half until error is reduced to some level\n", + "\n", + "in previous example of freefall, choose x_l=50, x_u=200\n", + "\n", + "x_r = (50+200)/2 = 125\n", + "\n", + "f_m(125) = -0.408\n", + "\n", + "x_r= (125+200)/2 = 162.5\n", + "\n", + "f_m(162.5) = 0.3594\n", + "\n", + "x_r = (125+162.5)/2=143.75\n", + "\n", + "f_m(143.75)= 0.0206" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_r = 125\n", + "interval left f(x_l)= -4.6,f(x_r)= -0.4\n", + "interval right f(x_r)= -0.4,f(x_u)= 0.9\n", + "ans = -0.40860\n" + ] + } + ], + "source": [ + "x_l=50; x_u=200;\n", + "x_r=(x_l+x_u)/2\n", + "fprintf('interval left f(x_l)= %1.1f,f(x_r)= %1.1f\\n',f_m(x_l),f_m(x_r))\n", + "fprintf('interval right f(x_r)= %1.1f,f(x_u)= %1.1f\\n',f_m(x_r),f_m(x_u))\n", + "f_m(x_r)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Bisect Function" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Much better root locator, with 4 iterations, our function is already close to zero\n", + "\n", + "Automate this with a function:\n", + "`bisect.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)\n", + "% bisect: root location zeroes\n", + "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + "% uses bisection method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% xl, xu = lower and upper guesses\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by func\n", + "% output:\n", + "% root = real root\n", + "% fx = function value at root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "test = func(xl,varargin{:})*func(xu,varargin{:});\n", + "if test>0,error('no sign change'),end\n", + "if nargin<4|isempty(es), es=0.0001;end\n", + "if nargin<5|isempty(maxit), maxit=50;end\n", + "iter = 0; xr = xl; ea = 100;\n", + "while (1)\n", + " xrold = xr;\n", + " xr = (xl + xu)/2;\n", + " iter = iter + 1;\n", + " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n", + " test = func(xl,varargin{:})*func(xr,varargin{:});\n", + " if test < 0\n", + " xu = xr;\n", + " elseif test > 0\n", + " xl = xr;\n", + " else\n", + " ea = 0;\n", + " end\n", + " if ea <= es | iter >= maxit,break,end\n", + "end\n", + "root = xr; fx = func(xr, varargin{:});" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'bisect' is a function from the file /home/ryan/Documents/UConn/ME3255/course_git_F2017/06_roots and optimization/bisect.m\n", + "\n", + " bisect: root location zeroes\n", + " [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + " uses bisection method to find the root of func\n", + " input:\n", + " func = name of function\n", + " xl, xu = lower and upper guesses\n", + " es = desired relative error (default = 0.0001%)\n", + " maxit = maximum allowable iterations (default = 50)\n", + " p1,p2,... = additional parameters used by func\n", + " output:\n", + " root = real root\n", + " fx = function value at root\n", + " ea = approximate relative error (%)\n", + " iter = number of iterations\n", + "\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help bisect" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## False position (linear interpolation)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner\n", + "\n", + "$ x_{r} = x_{u} - \\frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false, + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "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-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\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\t \n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%xl=50; xu=200; \n", + "xu=xr;\n", + "xl=xl;\n", + "xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu));\n", + "plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## False Position" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Much better root locator, with 4 iterations, our function is already close to zero\n", + "\n", + "Automate this with a function:\n", + "`falsepos.m`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "^C\r\n" + ] + } + ], + "source": [ + "function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)\n", + "% falsepos: root location zeroes\n", + "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n", + "% uses false position method to find the root of func\n", + "% input:\n", + "% func = name of function\n", + "% xl, xu = lower and upper guesses\n", + "% es = desired relative error (default = 0.0001%)\n", + "% maxit = maximum allowable iterations (default = 50)\n", + "% p1,p2,... = additional parameters used by func\n", + "% output:\n", + "% root = real root\n", + "% fx = function value at root\n", + "% ea = approximate relative error (%)\n", + "% iter = number of iterations\n", + "if nargin<3,error('at least 3 input arguments required'),end\n", + "test = func(xl,varargin{:})*func(xu,varargin{:});\n", + "if test>0,error('no sign change'),end\n", + "if nargin<4|isempty(es), es=0.0001;end\n", + "if nargin<5|isempty(maxit), maxit=50;end\n", + "iter = 0; xr = xl; ea = 100;\n", + "while (1)\n", + " xrold = xr;\n", + " % xr = (xl + xu)/2; % bisect method\n", + " xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method\n", + " iter = iter + 1;\n", + " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n", + " test = func(xl,varargin{:})*func(xr,varargin{:});\n", + " if test < 0\n", + " xu = xr;\n", + " elseif test > 0\n", + " xl = xr;\n", + " else\n", + " ea = 0;\n", + " end\n", + " if ea <= es | iter >= maxit,break,end\n", + "end\n", + "root = xr; fx = func(xr, varargin{:});" + ] + } + ], + "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/06_roots and optimization/bisect.m b/06_roots and optimization/bisect.m new file mode 100644 index 0000000..c09ffbf --- /dev/null +++ b/06_roots and optimization/bisect.m @@ -0,0 +1,37 @@ +function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + xr = (xl + xu)/2; + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); \ No newline at end of file diff --git a/06_roots and optimization/falsepos.m b/06_roots and optimization/falsepos.m new file mode 100644 index 0000000..0a3477c --- /dev/null +++ b/06_roots and optimization/falsepos.m @@ -0,0 +1,39 @@ +function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + xr = (xl + xu)/2; + % xr = (xl + xu)/2; % bisect method + xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); diff --git a/06_roots and optimization/incsearch.m b/06_roots and optimization/incsearch.m new file mode 100644 index 0000000..bd82554 --- /dev/null +++ b/06_roots and optimization/incsearch.m @@ -0,0 +1,37 @@ +function xb = incsearch(func,xmin,xmax,ns) +% incsearch: incremental search root locator +% xb = incsearch(func,xmin,xmax,ns): +% finds brackets of x that contain sign changes +% of a function on an interval +% input: +% func = name of function +% xmin, xmax = endpoints of interval +% ns = number of subintervals (default = 50) +% output: +% xb(k,1) is the lower bound of the kth sign change +% xb(k,2) is the upper bound of the kth sign change +% If no brackets found, xb = []. +if nargin < 3, error('at least 3 arguments required'), end +if nargin < 4, ns = 50; end %if ns blank set to 50 +% Incremental search +x = linspace(xmin,xmax,ns); +f = func(x); +nb = 0; xb = []; %xb is null unless sign change detected +%for k = 1:length(x)-1 +% if sign(f(k)) ~= sign(f(k+1)) %check for sign change +% nb = nb + 1; +% xb(nb,1) = x(k); +% xb(nb,2) = x(k+1); +% end +%end +sign_change = diff(sign(f)); +[~,i_change] = find(sign_change~=0); +nb=length(i_change); +xb=[x(i_change)',x(i_change+1)']; + +if isempty(xb) %display that no brackets were found + fprintf('no brackets found\n') + fprintf('check interval or increase ns\n') +else + fprintf('number of brackets: %i\n',nb) %display number of brackets +end diff --git a/06_roots and optimization/lecture_06.md b/06_roots and optimization/lecture_06.md new file mode 100644 index 0000000..b3d2266 --- /dev/null +++ b/06_roots and optimization/lecture_06.md @@ -0,0 +1,373 @@ + + +```octave +%plot --format svg +``` + +my_caller.m +```matlab +function [vx,vy] = my_caller(max_time) + N=100; + t=linspace(0,max_time,N); + [x,y]=my_function(max_time); + vx=diff(x)./diff(t); + vy=diff(y)./diff(t); +end +``` + +my_function.m +```matlab +function [x,y] = my_function(max_time) + N=100; + t=linspace(0,max_time,N); + x=t.^2; + y=2*t; +end +``` + +In order to use `my_caller.m` where does `my_function.m` need to be saved? +![responses](q1.png) + + +What cool personal projects are you working on? +While we delve deeper into Matlab functions, could we review some of the basic logic +operators it uses and command codes. + +I still dont know when these forms are technically due. + + -by the following lecture + +I'm having trouble interfacing Atom with GitHub. Is there a simple tutorial for this? + + -Mac? Seems there could be a bug that folks are working on + +What are the bear necessities of life? +please go over how to "submit" the homeworks because it is still confusing + +Do you prefer Matlab or Octave? + + -octave is my preference, but Matlab has some benefits + +Would you consider a country to be open-source? + + -?? + +Is there a way to download matlab for free? + + -not legally + +how do you add files to current folder in matlab? + + -you can do this either through a file browser or cli + +How should Homework 2 be submitted? By simply putting the function into the homework_1 +repository? + + -yes + +How can we tell that these forms are being received? + + -when you hit submit, the form says "form received" + +can you save scripted outputs from matlab/octave as an excel file? + + -yes, easy way is open a file with a .csv extension then fprintf and separate everything with commas, harder way is to use the `xlswrite` + + +Also, can you update your notes to show what happens when these things are run, as you do +in class?" + + -I always update the lecture notes after class so they should display what we did in class + +I have a little difficulty following along in class on my laptop when you have programs +pre-written. Maybe if you posted those codes on Github so I could copy them when you +switch to different desktops I would be able to follow along better. + +Kirk or Picard? + + -Kirk + +Who is our TA? + + -Peiyu Zhang peiyu.zhang@uconn.edu + +Can we download libraries of data like thermodynamic tables into matlab? + +-YES! [Matlab Steam Tables](http://bit.ly/2kZygu8) + +Will we use the Simulink addition to Matlab? I found it interesting and useful for +evaluating ODEs in Linear systems. + + -not in this class, everything in simulink has a matlab script/function that can be substituted, but many times its hidden by the gui. Here we want to look directly at our solvers + + +# Roots and Optimization +## Bracketing ch. 5 + +When you are given a function, numerical or analytical, it's not always possible to solve directly for a given variable. + +Even for the freefall example we first explored, + +$v(t)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)$ + +There is no way to solve for m in terms of the other variables. + +Instead, we can solve the problem by creating a new function f(m) where + +$f(m)=\sqrt{\frac{gm}{c_{d}}}\tanh(\sqrt{\frac{gc_{d}}{m}}t)-v(t)$. + +When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity) + + +```octave +setdefaults +g=9.81; % acceleration due to gravity +m=linspace(50, 200,100); % possible values for mass 50 to 200 kg +c_d=0.25; % drag coefficient +t=4; % at time = 4 seconds +v=36; % speed must be 36 m/s +f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m + +plot(m,f_m(m),m,zeros(length(m),1)) +axis([45 200 -5 1]) +``` + + +![svg](lecture_06_files/lecture_06_5_0.svg) + + + +```octave +f_m(145) +``` + + ans = 0.045626 + + +Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0 + +Better methods are the +1. Bracketing methods +2. Open methods + +Both need an initial guess. + + +## Incremental method (Brute force) + +You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. + +```matlab +function xb = incsearch(func,xmin,xmax,ns) +% incsearch: incremental search root locator +% xb = incsearch(func,xmin,xmax,ns): +% finds brackets of x that contain sign changes +% of a function on an interval +% input: +% func = name of function +% xmin, xmax = endpoints of interval +% ns = number of subintervals (default = 50) +% output: +% xb(k,1) is the lower bound of the kth sign change +% xb(k,2) is the upper bound of the kth sign change +% If no brackets found, xb = []. +if nargin < 3, error('at least 3 arguments required'), end +if nargin < 4, ns = 50; end %if ns blank set to 50 +% Incremental search +x = linspace(xmin,xmax,ns); +f = func(x); +nb = 0; xb = []; %xb is null unless sign change detected +for k = 1:length(x)-1 + if sign(f(k)) ~= sign(f(k+1)) %check for sign change + nb = nb + 1; + xb(nb,1) = x(k); + xb(nb,2) = x(k+1); + end +end +if isempty(xb) %display that no brackets were found + fprintf('no brackets found\n') + fprintf('check interval or increase ns\n') +else + fprintf('number of brackets: %i\n',nb) %display number of brackets +end +``` + + +```octave +help incsearch +``` + + 'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/me3255_S2017/lecture_06/incsearch.m + + incsearch: incremental search root locator + xb = incsearch(func,xmin,xmax,ns): + finds brackets of x that contain sign changes + of a function on an interval + input: + func = name of function + xmin, xmax = endpoints of interval + ns = number of subintervals (default = 50) + output: + xb(k,1) is the lower bound of the kth sign change + xb(k,2) is the upper bound of the kth sign change + If no brackets found, xb = []. + + + Additional help for built-in functions and operators is + available in the online version of the manual. Use the command + 'doc ' to search the manual index. + + Help and information about Octave is also available on the WWW + at http://www.octave.org and via the help@octave.org + mailing list. + + + +```octave +incsearch(f_m,50, 200,55) +``` + + no brackets found + check interval or increase ns + ans = [](1x0) + + +## Bisection method + +Divide interval in half until error is reduced to some level + +in previous example of freefall, choose x_l=50, x_u=200 + +x_r = (50+200)/2 = 125 + +f_m(125) = -0.408 + +x_r= (125+200)/2 = 162.5 + +f_m(162.5) = 0.3594 + +x_r = (125+162.5)/2=143.75 + +f_m(143.75)= 0.0206 + + +```octave +f_m(143.75) +``` + + ans = 0.020577 + + +Much better root locator, with 4 iterations, our function is already close to zero + +Automate this with a function: +`bisect.m` + +```matlab +function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + xr = (xl + xu)/2; + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); +``` + +## False position (linear interpolation) + +Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner + +$ x_{r} = x_{u} - \frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$ + + +```octave +xl=50; xu=200; +xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); + +plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0) +``` + + +![svg](lecture_06_files/lecture_06_18_0.svg) + + +Much better root locator, with 4 iterations, our function is already close to zero + +Automate this with a function: +`falsepos.m` + +```matlab +function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin) +% falsepos: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses false position method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +test = func(xl,varargin{:})*func(xu,varargin{:}); +if test>0,error('no sign change'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +iter = 0; xr = xl; ea = 100; +while (1) + xrold = xr; + % xr = (xl + xu)/2; % bisect method + xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method + iter = iter + 1; + if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end + test = func(xl,varargin{:})*func(xr,varargin{:}); + if test < 0 + xu = xr; + elseif test > 0 + xl = xr; + else + ea = 0; + end + if ea <= es | iter >= maxit,break,end +end +root = xr; fx = func(xr, varargin{:}); +``` + + +```octave + +``` diff --git a/06_roots and optimization/octave-workspace b/06_roots and optimization/octave-workspace new file mode 100644 index 0000000000000000000000000000000000000000..0918300ce10589a9faf7b968c7df71a640e71a70 GIT binary patch literal 1332 zcmbu9OK4L;6o#*eMi5QGMG+Sr6lxksiXdIpcu*?VK9lCrq)Bfk_ZdiTHMuEOx=?(8 zsJKvYAu6c2sN$k4>Bj0t6hteys3NGiDC)vxSB*1s4oI=;!V4L==bt%`@4u5F!{?J$ zWKX1DTd!#v_v&KRwB|;!)D2$YwZO^4m*&s4zk4ioLpP>5x^?s?_dR4tPgYu=} z8Gd=fV=nj1inUdQpqqOW-tp>0qaI^>_sK}8+QI6_YCa1&EX=B`lZhf7EE;A$_gohG zrHMF;b~3SFXdNuVCjZsaI#Iw8zts)5o&6SnbaNL-)&EUZ0j));lYdYTDE04|34{b$ z<$kR^wbE=hWl2depb~4n@HEgi{o_e5oC(I3Kb}4Uhgu8%LQfxTlHouL9Acrh`l4*)r;*o_rBFoaw9C68kCDJc` z$q}D1*ewg(8zWA6P-Wdb@ya8IaDliLpehIZULbzs$l>8QDC|KUIULnNjk}=ZJll)J zwTN?_D1yQ}Q^fh+7l|{+IlDRW<~Z*SPTV=py~2rq3Fp^K|jB*hz%?dh)4-YtdL+28?s@+4@mgVt?HScws9;#NI})|=kDpc zb?!OecOG{%HSFBJa+S8c4T0ASLVW4>@7#U*xq}l=ofE>${FJ}`xYtOUM19;MTz{Jl zwlCcoGi(V2wf#rHh%|Ok- zV;LyyX#MQ-qTlOwcwJ-VCHmqQZawyv`eMyM%|OjS%|OjS%|OjS%|OjS&A|T`12#J$ zzKxJRRfTn@dOxxBUj10V)eO`O)C|-N)C|-N)C|-N)C|-N)C|-N)C|-Nd;%Hp@#po^ z_~rjcpN8}J`~Up+|9_kk;;+D8fFA-s0Nw)50GmJqSOH!DP5~!@6TtUR3NZyHz>C0_ zfj>Mi#QVT^fY*WRfB-)Fyb$jJzXIL?ehOR%KKh&x9|9i$e*?Y?On?Fy0DHhQz%M^5 z#4muK06zv2AO`LLw}B6z6XMUnpMc*0zXsk0z5<*EP5~!@pMOS(H-TN?S>R{S;xlj& zxBz_cj1a#A-Ukw(1-uOW@o6F61%3|jr9dZ{PWjmxQYj_v{Peear#->z0+iSxS<64TN7G0J5ifC;haPPk$%cIz6=~o+v5{}d z`uAl{BV8DpR91Z66dSZ?zcp(ev}|j28cA&Cesng4vpSKn4wRie9MgXo&)k#U$$o1x zYfT=C$Nz7~GxT#AQUPVkDA$j)$Q|!iCDJ3Y;n9^8tgDZH^@% zH}K=2h?(|7s#wEg`G6$oUnu!B(#fvP%+VWowBt`C%~+hNHuuSD)=$TwG;BMD%gxN_ zVJwrqe$)Bshnyq4(ngY+o$g-wY}OtdZ3Uf@rqv>sEM>wz4!%a2 zQK?;BT(UM;o2QfWjJL=GqcUvh;&m*9ZfnkS->4jnX0b-L=F}V^Sx4}D@93hJlzzih z>D(seGN6&s33WJ(euux%X_%Ku1?jIldoJs$9fs0Qr9^Sf2@mItD`fX%)vkdyLo8x$ zWrh+OQS{Z0*KHC0Azn0kzLxpv`e6QP0KKoRyY(JjMgg=a(V-g6TzsxXC}*CbG&I8I zDdm|A+={6#izsgqTcBJW$n`-qD0hodYLKAiEX|XKjC@ci15sh9G!Hii*9J^Os3t1p zdY9DgRg7TwxNxNmp|t29&&`%|%4O$*+lQNFM&{^3*bNz2`9p)bLSBOCucW0UgxA0X zL8sU0Yz=7j1|9*VZIqzanKWfz8XQD!tT^_E8ch(5 z3dFKaeod58uoSD{J|R^&zOnR;lWS#1TC37=zQW(9w6eWuWM#fSmLp4e8nN)`dZySF#< z?!DC~TUQoD6$w()huZKK{ed-c@m$u4bmiOV(ra3eW69AakQ?nLm)?TMXW{YLWsh&& z_{Ka`qsEctY3%+E3<4IMs!iBiYD5wsX`p%ht&wiJ4FlZg#< zRwni3n>VlTP!MB4EqJvBueO%Gnt*LfLv1Z86&|rSE!@gNHn-75(1n5sSs?S~Vu)v8 z?98&USf%JkX+A|`bL2fJP^g%H=jaZXSEyYwN60rWh8P)Ab>UpPk~Uu#js*Ot*fVrN zpCajVr<^K_Nl1gnfGeHs3`PbGI$v$4=}@58t6h=%w^F{9eOVx46TWXQ`<}@(P%*nD zL&U7?f+jAxuj#4tp64}c&9!sZektvjq;jr0aU7)(e8G_oII{7OBSuXi>Ri)ELJ7{3 zq28zFl@#X*KSmcS2|k)NbAHPGo5#GbItm8Kjl4VXFgp)BVDB=qPi z&KTfj!T&z|?^piMbvqJpP7FK$@_MeusS_@Uw8n&LUXxQ-J2)X>R8}($C!9;@4xk!$ zVU=s!EJ`(*9qW3UF5Kw?!K6DrBL{eBnqAr*BW;SIS77#}o(zkHj<9s4y;hNi*|e%g zPHgrXy%k5bHQM0_Tsv@HxJRa`G|y)NU8EM(>$`15&*?|0+~=9eX81vCLPL3yV3r5I zt2~@6nH@erfYAlzA7OEVlO`Wis~(fjxNR5Z8R*Cz?%>W{&eExGanw9+x+3JL*)CHy zA1BqL<2(GwxCc%i!}7nn