From c07d13fe4f902600949cd48f2d542c509fbebf9f Mon Sep 17 00:00:00 2001
From: Lucas Christesen Ahler <201907296@post.au.dk>
Date: Wed, 23 Jun 2021 13:46:47 +0200
Subject: [PATCH] root-finding homework
---
homework/root-finding/Makefile | 70 +++++++
homework/root-finding/func.png | Bin 0 -> 6772 bytes
homework/root-finding/functions.c | 249 ++++++++++++++++++++++++
homework/root-finding/functions.h | 33 ++++
homework/root-finding/improved_wave.png | Bin 0 -> 7944 bytes
homework/root-finding/main.c | 75 +++++++
homework/root-finding/out.txt | 58 ++++++
homework/root-finding/simple_wave.png | Bin 0 -> 8631 bytes
8 files changed, 485 insertions(+)
create mode 100755 homework/root-finding/Makefile
create mode 100644 homework/root-finding/func.png
create mode 100755 homework/root-finding/functions.c
create mode 100755 homework/root-finding/functions.h
create mode 100644 homework/root-finding/improved_wave.png
create mode 100755 homework/root-finding/main.c
create mode 100644 homework/root-finding/out.txt
create mode 100644 homework/root-finding/simple_wave.png
diff --git a/homework/root-finding/Makefile b/homework/root-finding/Makefile
new file mode 100755
index 0000000..5eeaabc
--- /dev/null
+++ b/homework/root-finding/Makefile
@@ -0,0 +1,70 @@
+CFLAGS = -Wall -O3 -std=gnu11 $(shell gsl-config --cflags)
+LDLIBS = -lm $(shell gsl-config --libs)
+
+default: simple_wave.png improved_wave.png func.png
+
+simple_wave.png: paths exact.txt Makefile
+ echo '\
+ set terminal png;\
+ set output "$@";\
+ set key top right;\
+ set tics out;\
+ set yrange [0:0.4];\
+ set grid;\
+ set xlabel "r (Bohr radius)";\
+ set ylabel "f(r)";\
+ set title "Shooting method with simple boundary condition";\
+ plot\
+ "simplepath2.txt" using 1:2 with line title "rmax = 2",\
+ "simplepath4.txt" using 1:2 with line title "rmax = 4",\
+ "simplepath6.txt" using 1:2 with line title "rmax = 6",\
+ "simplepath8.txt" using 1:2 with line title "rmax = 8",\
+ "simplepath10.txt" using 1:2 with line title "rmax = 10",\
+ "exact.txt" using 1:2 with line title "exact solution",\
+ ' | gnuplot
+
+improved_wave.png: paths exact.txt Makefile
+ echo '\
+ set terminal png;\
+ set output "$@";\
+ set key top right;\
+ set tics out;\
+ set yrange [0:0.4];\
+ set grid;\
+ set xlabel "r (Bohr radius)";\
+ set ylabel "f(r)";\
+ set title "Shooting method with improved boundary condition";\
+ plot\
+ "improvedpath1.txt" using 1:2 with line title "rmax = 1",\
+ "improvedpath2.txt" using 1:2 with line title "rmax = 2",\
+ "improvedpath3.txt" using 1:2 with line title "rmax = 3",\
+ "improvedpath4.txt" using 1:2 with line title "rmax = 4",\
+ "exact.txt" using 1:2 with line title "exact solution",\
+ ' | gnuplot
+
+func.png: Fe.txt Makefile
+ echo '\
+ set terminal png;\
+ set output "$@";\
+ set key top right;\
+ set tics out;\
+ set xlabel "Energy";\
+ set ylabel "Fe(rmax)";\
+ set title "Zoom on Fe to show instability of ODE";\
+ set xrange [-0.35:-0.32];\
+ set yrange [-4.6:-4.2];\
+ plot\
+ "Fe.txt" using 1:2 with line title "rmax=8",\
+ ' | gnuplot
+
+paths exact.txt Fe.txt: main
+ ./$< > out.txt
+ $(RM) N
+
+main: main.c functions.c
+
+
+.PHONEY: clean
+
+clean:
+ $(RM) main *.o *.txt *.png
\ No newline at end of file
diff --git a/homework/root-finding/func.png b/homework/root-finding/func.png
new file mode 100644
index 0000000000000000000000000000000000000000..c68e21e16175177fd8c646b0638142dee73cabf2
GIT binary patch
literal 6772
zcmaJ`2{csw-@g+XGS-qcyYxgz%HEJA%OgwHvL(urtl2e_eaV`oWQ&l!M0Q5@Vx;n@
z?1NA-c1FW}@4Y?m|GekC?|IJ|GxvUfzu$NHeAj!$8eO@}bcE*!002z7I+`W`a0m|o
zsA+l{h;oQU@*n6%*643jEf54DN`8Jm{6hc$JV*dQU}c4%2K4m->Mk=NK@h;J5%8dl
zJw6@))bId^2XWx%qSL3LqiWM1pmRJtJm=1xGcYi?b?X)mhkNqm$?Mmz2?WBxz`*?c
zJS1XORyKJhx$?{RRkPJ;5L}rE1K)I1udINh03;1?guZdY@xoH^7RY@(H@7FG3IGg7
z6Mq7*iXcE729G@gaYn0nJg5ci;k}nKf>ysjY^`@|h+^E@8qdcIf<?ee*h(7kX(E3j
z*+|XIO&Wg;z#lt?SA#xyygdPL58`3{)xdn{V-Mzgf&>Bx*#i;=>vLDZ)%-Vu9bm^G
z;bjP@mJbTHwY9akw@*z?<>lqs+1X*SSOEb6X_oXt$Yusl9V;IIVD5tdQ0WrPya2$p
zsH>@N8j!u5!<bQC!J~PPkJ#%a_>NlgA&gY0C33`4Z#mckTvE{~DrY{$)>a@G!Fq{a
z=>Pxv`Y}y}q%D%(%=BukIJG@X^2Cs3-Duh%)pWuZ%-fg0)3mfyVCH1;<Z0#c%4_>E
z`}8ByH$-JV=AXn>8KrATfZv?QS~;oA_7Z$MXU3k=e8HO4kAk;BCJCo8x8IxnZ%h*!
zC*rHp4x(Rgu<7v3wc;r0-jku?c<h*42I$zcy*>4Exrf(#@fvRt*%>{N^KGAAl%e3W
z-=>Yacw9G=*}bCFgo7J-imIZ`HFusaO_u(ph;7DcH{N+8Bl3)U&2(f+857!BXj5+X
zSU{ujmi39Ln;|1#Y_VZwdgzw|`dbNs^cM`dYQNx&-cSxwp4p6iO_D4D;@TJ;fjrG6
zROuOs;UB}!P0#jtTX*()ce^M;BXMQ5Cy2jnr_#4~SQzenyZf~wu&_#Cxp`*Zq<#Ei
z+4UFiv<#nLbwUK?eymOo1qM4ugz?`Qlnkyubl>RIuPa{qVp4Y%xM$5xK%b>>RV+18
zqSjtm`RT_;H?2iI%VM53$aD;KADh>YbHH3a-_~N&l6gP+O5Q2on=kH_DgFMkd-~kl
zH+vlWV_Y}(mF43fv2M_wxxwePAYZTXGMnzM5BXFrn6JEGEQ;zxO|!~c5#^646AXH<
zgZyRRfkgCbsCicLE2;YBb+O31mzge8w&ynrm9&B^YgN+@LPI`oo(_%?uO;<3YJpTv
zwTaAx%7H=ULoW^vCyNWQoaF-&xtGp7y?cjdm;6i-??+9&iIqir-jAht;AH;kSHy1=
z+Kjd9T6`6kqETplH>bZDjnB!VIR&CWghYTr?f|ETXRi<!jSmp;8eO|$KskFSNNwYo
z>W8e2G&7b{>8AUtSM}cm-z=Cknl7GMVh$o4+*|f7UTa6s>C*FU)80_AHQ7jT2SC+Y
zrs-AMV35S+84zaLoy^TkpbQuPN;-(h5wbn(88pJT*nQUi6@xu4OES}Y!2Z1pmyzRF
z28~$-#sdKYj>Vnk3)<l#<I29x8}r9$3|~td(bg)iCg<{KD-uz?7E9ZE<FReZzp{!Y
z!MwvP9$#+z{_t2B`aI}=M!X>M@yUW{Y<MUxl|FcTbZp@<TJqqq=YFyXevjc(#SL7;
zn}FEHlAkQoSq{w8ROj-kn^|NE4&!X=^zG5!VH2~^Ooc0oQf+bPJS@{TLTcUbGfyW~
zlfEyoc_&2q{57Y0jh~pZ;H)<qDHrkllBd*@|0`9K$s)y!Bb#)`$D?DP7d@%3k<s^m
ztd^Q~_sDc@L2|OAY0kusqvd&#_sd1iXJxs?7ChLpwxaDs-u2MEl<D``!dFwA|B!}H
zDK1kTJy#egSxMW+pIPB{If9%vh4BX0-z45{{a)Cr<@}Gpea*$3jpvg`9GLTsKIL}5
zrcG5&VgIsOv4cHET&dRaI;AIAdQ*97G1S9UzRJoh_^Hb|Jh#bDP~2-0(_uosWd3AA
zp?=67IQM)-my7U*TA#GYPA7x?A#BC6%gf3wDpW^dAgRP)Qd#MecMIxUF|~&+ROeW_
z16jG}l3lHW7fxaX*LtkqeO)#v3N#+$6(olD^|6Raj8wSICY)dTJEwRnU0(z6AfL41
z#0YX5+u)4q%Mxz8eSC5PH=0>puP3*DaIxavAdt#dw3BW-()0TKQk-_?wl(%T=TwX8
zq0#TS<E#tnj&spVq=U>N)=n{>8onc*%!?*5dc<xcsjEY`1rI;m5R_fN$*Xq#Q59G7
zWgmmsTP^wfcGZ@@Wm+7zg8m$@p{Zh;Hhn}K*sBa`o_VOL(59<uxG#-61uL9n0F-3;
zHceaOqFT7Zx<AlNu9sDDtJM%jXqeA0eg8P&M=z!1g9*58=HA+f?t+Y`ykda4VHzCv
zTYD?3N#>t$M`7Cq*-GU*tLfk^qWmk&tT0fb{E4eBYT4)%by5eAiuO8jz>VBFKzHI_
z;{{;#bws(;*1427q(%O4c_Gg!D2ejo*Q8@G4E3ZlAt<m2hes?3hVa1VB)|NO-H3+$
z!hU>b3oU`L3tN}?cQs!IKW7Z_<BNWcJ9jo`DkFy1?}e|V`z-+5OWt77&|EknMmk;O
zzsRBMlzdx(qx}U@S+A=un0HJDZ%5i4n>xRi?VlbqHt{VUMP=Vw8*P^;P$WyTl{9%V
zzI2q;3Hlqo>kCYXk@hxYFPK#-tg~)y&eFch_RoxY)Fgf9sM&i_7OxAzQL98yl5MBm
z)^ou1hNxl*dg=lF3OO^zt!tH+aU_Im3g|JOkzx~ncY}olarLH>I7cb&u?@r(Y_<uV
zZsZU(Pb>%+6(enK>j}TYq{i)a9ex2Kpfq+SjR@RD%wF+yqSc0hYJ6ApFU(b>wPo@j
zpyI8SsORu7c-C{*cI6OoId1^{KGv%LO`n<j!6J_=STwL<xnOu1$KMIWL0*Xw&s@k8
zN5rT!dQ&<Oo~(vHg8vP=7=fLd+W;(FacI92OA`Z~=?e@o&l3P(=Xmo}k1U`O>DedM
zwlRnhafaa?ra!a@l&B%Q<b4@9ZWio;{+qN};A>=wdU9ZMiWXsKSS5}1rh(b5KT4kY
zge>zN_1irE=gMd2n*<tR+)JB;gi@c>gWGl|5vBbYe78oTM2SOMvzv736-U>#FB!I5
znkZm76WsJ$wW<TO=^7tcEVBrhD72sADGA1QyFMX{Xm4W}vh9u~{M@XfRE2qA1Qn#2
zr)!+YSC#$yalCBbLCN79&_eD?MwJP)cwv~7NuYpQvyR;kkscxr2M$W)(>x`m2WEN>
zl?uKEbTS3MpB8w1Joa~HMFvIquE)W8_!nYzTo1z|jep{TkKXG#%4)#A;6B_?aMi?A
zh4DU{ac7Axr`pw$rg;l5vd};O9p4$_Q&ci7Vk)D9&Ha~1k~ylUvT?yl=BQ5n*ZY!Q
zRTu9e_&da^{FX<cbI?8Z$FssC{n3Zv?0l{xw=a)^VKu=l2w(&Q<BXb<@@SF7onLfA
zoT{t1$AjRQItW_97b&W-pDYMaG>N!<B>qb}|4O>e9~Sa3JFu8W6AA`xq6+IpCfIr2
zQ7F-*luan1V72?wZxgXd5cH*$ye)wPSvyL_cU>7);~_>jcx&oUt-GAn-HwEf2sW4h
zfiWytASuca5i+22e=Uy&E&`4{K2Ulij2|w{qGur<Apz96PWDew5-GTB$KQIEwGbGz
z=OJ9Mk;v|qo~kNbiXjOvlcs2dmB^>qM*WtYSkbYZAV2F+phc{L#gm3U*%-!{U)=vD
zdzI<<bE=D)ww*?;@OdJ;SLohexO&EqAg%}sR$p=f%Kzv5@z>4uf{tfza-7Wut%R>+
z3|>3r)?Gw%>tKkzV}GCBD}iE0Wz0`uGe7iowj-WAq3JzoBYgLT+2A#=A=D^J``-QQ
z!c|v;w6p3+{4^3FS&Jy^%$QG8{jpuLeSh1cHUZ6^UX=2mRGpob!8hakb}!)9!XfN>
zEJ``C9k8Q%EY7~hB~vuu%gbknHB<medkC~OW0^vE&wk3Xg^T|Upzwjxc4C0C^mVh0
zB3B!WEwMCpV%YXP8LN_$S(cQ`T9EY2+xnkOJvv1eUusAW`P$iuyMCMjnwjb<yQaX?
z@G16pf(wM^ZKwG9jSKfkR@`&6GVwT;Rloa|!@~_#<%WNSX9Ks^$f$$c9eS}dy@j%U
z9-hZRe$Ng{nQ2v|d#O|M%~fS=XC0t4x_vT3JrR`$zJ~Pi-7BkM{OwLp!<EIEjpL#k
zUb=RowRl#S{zA2momUZ>4}4!u@(k>x=FOLF7<~D}N+%=g2@CeOHP0fxYLG9*)VzPH
zq=sAA5`Fk1&CSb0kJ~YCG{SIpxutBtoR>J#%HWi&YAXxgUc)jTk6etlV_P+OEFs&+
zPgbO!vtB~Q*Zf6WaziH23cW}UDR+GKt*5Tj2x1cL-%XC&aCdtHKr<eHwbK5H#_u5_
zhpP1*01cpdcL#u@u;pp58)!QdwlPTfB-<Ajw6}JFDy?%uern!IhJA~d1;iYDlF3<J
zK}C-Lp0B~wkeJG`<JM1Nwu`JX#P3+w9pNTC{p`_r;2U1Qkpoio>j$iRfvEYqa55gO
z`P)plbID>rYxV&2)}VXU;=g&m3(H`|k?`A2eWeED5ZL+cJ&ja$#I6CMBCcn@T+m;G
zjytd$A{OLD{11xPHs(NFsHD1k>Ue^WX|jZ+#0qmYe<;)BiiML~RX;YR(_pw^MB0t}
zB70xgb#c2})FjlIh1|B&JoS6?xtM~$;JrBB=9?55j)bT88_$Hdgd7Y(X=K~L>M6xz
zA%lMwK%wS!<IJqTbPHNIISxrww|Y-w*(OTM*7r_|#y$K+>7(b{7ZuHTSV7VI(}e;U
zJiZ0Za>D5wfVTSqsODs%+r`pS6!RGJdWI-^W`6uHQ{J7=OV!l=rv)6V=g1|-zcJk_
zZ!k|+_m92-7;fZrV#7gsp_*mWQN60T;Hac~czBp0d~bD=W@#XWOrP(7bm<KqwdJ2C
zw-`AEI=PtGb91k+LuqJjdGzDEnima<#zAfCTjV&2-!vDBa+3qiwWC6$f}<qdF8KnD
z>P{wY0h{ukP6qVBgt^Zi<AKY^I1^D4#usEg8!yqVJ~re(;FkK#U>sy)Jnj7iGc7bV
zkUJB%aVct2Fdd~rDle?=7@|yfG@leXchMIrZrFvvUv?Zqk(lX=Q<fpMy^#Q2Zg`uZ
zi1fwmwlEa!Lh}ahI2~)2bmf6(Cm_7B&zy{k5DQ;*QK7LsT{OT_#h5sF`OV~R<U+yo
z!z}Om9+2---XL8-%~AGHAyk%WV)*{lRXU8-!?vi$L=OpHrt|!KjNiBk=P~vGA4k*a
zB6+`~x=t-AcSxKgftMqz{5<^i={~IYLnXJ%XVUE1`$km{Yi^H1q9atp1_N(LN{!qr
z@_FP!_aXH?G>~L$tQ@fk-P8Qpl%0o?Jbo32nGPN2Xd1}SuAEinCXn9iPchSzsY6F>
zYOL){22E^48WY;0pt;9l(RHAjWCtWXe93OB>sMGUlnS$YYu!JhB=!6TaBd)1sCZ-S
zh)qt}-Xqpt7a5wzGc$QY&;Z2IRF+d2%1s!u3vNe7j9ihSO%8Z!2w<jLwNi6b(U2O=
zRf79T&$$rX@;nNOs)!0@wPR3*ii(r}bjFly$@T*qq}hCa`KMF4av@uuN6a8Yl8Bg<
zVcP&aScnMy0?<gyV~z3v;I7!;7rO3*pxr+`lZW)sQhT%l!j5*Xas8$MTn<F1A2|(q
zz}@#DZAi^83KPtOP2dxF7LT;57)G=mGOolhLOFu$YjAAb>XAnFaoid@4EIgpx5J}I
z0~LN7Up<cK!ymZZ{1XK)k~%h*k%Iv;cx7abxn=hNDMnondU1amrvi0*nNv4AFRsVQ
zLPF|w**Q*~8LG9MRQ)wC>iV+yU)g!_J@e)gi6X0p?C|W9d9_j@ihOJ>$L@H-`%-hh
z>Q*H!-vtY;S=lm<-p+O(P5=95a#DrYyf9OW!b<IDc;u7SD;2J^WuGjDr<<>z<7WLm
zN@AC(WzN<Y>PP_xURk0<L+E%ZNAG5ixD7Da*t_G_?X@?=0>ea$caz?STbDMV04|_m
zvsOBWcMidhYd~L~VubCQ<>w2PMl8MV`LkaFA<^-LRwIz}u8!5aAta^0r?g=7yds0v
z9t9vng2ZPPN&^xqRYpDw;Rs5zhxi_lk)w<<K|8SQDEaLDW=2?c%gpeU4a{SP@smF<
ztl%yozB&w+jS&b~He*HBP}QwR5T;k&4$Jewswyzl=U@=+i5;o!C?o<E19mxCVHVlp
zGmP{w_bY^9c}ZkDsy&1e5qO2L8;OKaujigLVzok)vVbNMXb?#ZA$@V;fuf$!<qORS
zD|K}wg}t6XqfB)rNAIxwakd=)w3y(Vl?sV%#<2@6<(lKX5W<@YdU6XnVkBEj38F6^
zV=q_0QR<Zu6MO^S=M_8O(yTfD-nDzJJZU$=|KN6P!w%%TVY?9FhITvQ6(>l{PAXo9
z*sR~xT4B!DnTdz02y|U;&PVeYkBle5#xN^gWb%L>mKk8vED5_sC^$_-vO#?W5$%Qz
zbtM?LvxFQL_Cs(+psuKiBNL_yRi7C~W_QcO8{X0}EUzd?6zQ8@{7p+w<m8f7*t>`d
z+c+G7No-33rd#OL=sI`k0is)O1{5WK<@(*9bKxG}A|qGo+$xmNV|2+D{3?3esER;f
z_Y}sygLxHto)1HXq*l=Mt*RPMqUL_RpPl1V^TPvkJZe!a5qD_%zG5-5OoGy(>OxaM
z4NZh#Cm<W5GsmTNr42~21HPgzrC0#m7Zqd~1v&YjrF2nBovnv!e{?88wP6iOrz+!|
zirlXS^-^4FnkKZK#&TE4hLeg%dc^1j+xPII>Hw$KY?@T+=R4(cU6~!#7u8|~2wKAp
z#m9)dMGS(<YMKSg?UFp=dNzb&{Ye_{&2H#<FH>t)vs6>a^v6v?R-v>MTMpE-c~6Te
zx^<^y$EFe7hZQ#ROzTJmAr<xXgL|FQ&@{kH3@g{1zxH}t@Fz?Ajnh2h@AXVekC@+H
zR)=1YbXv<@@pH=jQQlCY+PBwv^nblx?dJ%!347UY>e`^Ao)wQC+qVdOQ^ID<%RDl<
z-%jvq?%?@FK5sD+dE9l(AoA4M!I@tDt7A3nm2KDROlmIkwSFjXeOslM_>X7)_6m6T
zel2fV$2(hNPD|%H{~}w-fgkp~fj*^euF)K#LSquY%I2Jz(5~bIyXngiXE`wY{puz6
z=2+(e<)@K9_N8w<%WiXh*}Pfu&BNnbI$!0C%<=@tNOhGr4KJ1-c^`C$G3T2>L}(|4
z|Fwhez5NtIx<*FMcClIiFc0&Ht=ES1YWai_Q&z-c?&+5wN+&ci5|6%j?L^l$U0uPN
z;B-7^bnHu7UCg{|qiE{<f5W>ynX=M3*dh{_Pm<anYG;WUmEt_4smfnm{>WEbo;ZNz
z+lRXcE=cc{PwYP4dOMcbfUV^52;gVntF4FBd^Iz+yx%+{z8TUi#B)q9b}1Ks$zr)j
ztD4jHM4X<_r4cvRuKZ7ABj5J_$R<mB7T937c|Sio)p_;p5y+|@Pgj|TOmDfqy0Pnf
zp}sve7rX4HU=)lxIMlmrOZYl1S9H{JvAeW@mtK{o6fN6o+whacd9Wuv=i^4I#g%X4
z%Pvtak1hWL?gPkepErg&B8kf<JMRnKOuRWZFMnCsKDb-XowIxfvJ1l1S^n>EZ<TH5
zLyR`H`v2d1_lj4HOC~FTR-dL>mYIBNIe(wXZWXh*=qVo{{nxvi^Sq`HfB!yZRa%UT
zDl78~5sr3c&(H$QB+9%`=c(L@6^l35E#2Xhi?0in_Y#~#Z`|eABVT@2?;<xw^Vl8r
zxxpj%>pe3w{|lc*yZS@e9lR~>HI=(%X_Lu?WVTy0oULDMcG5>LoD+1+>}FM=^PFT-
zn-VDjwmUuER$kIQd51<pVBxPuj6F!3ND!491DJW!%_dk06fPrj=3BPMw8f{p3C9Sp
z<R0)HdNCBHb=&W+qfSk&sk%~MUQgd;&Y!e*%|^=&8R;=A-z}P~NJ;kk<aCty{hAZ^
zrIq}q;?jY@{<;$e42_YC?@Y*_Cel<63mJMlEx7zAmqH$3l^Hx-uv2}l&_6-|U9Bsc
JRhR6;{tE{mt{eaW
literal 0
HcmV?d00001
diff --git a/homework/root-finding/functions.c b/homework/root-finding/functions.c
new file mode 100755
index 0000000..e836089
--- /dev/null
+++ b/homework/root-finding/functions.c
@@ -0,0 +1,249 @@
+#include "functions.h"
+
+static double E;
+void f(gsl_vector* x, gsl_vector* fx) {
+ double x0 =gsl_vector_get(x,0);
+ double x1 =gsl_vector_get(x,1);
+ gsl_vector_set(fx,0,400*pow(x0,3)-400*x0*x1+2*x0-2);
+ gsl_vector_set(fx,1,200*(x1-pow(x0,2)));
+}
+
+void radial(double r, gsl_vector* y, gsl_vector* dydr) {
+ double y0 = gsl_vector_get(y,0);
+ double y1 = gsl_vector_get(y,1);
+ gsl_vector_set(dydr,0,y1);
+ gsl_vector_set(dydr,1,2*(-1.0/r-E)*y0);
+}
+
+double Fe_simple(double energy, double rmax, char* save_path) {
+ double rmin = 1e-3;
+ gsl_vector* ya = gsl_vector_alloc(2);
+ gsl_vector* yb = gsl_vector_alloc(2);
+ gsl_vector_set(ya,0,rmin-rmin*rmin);
+ gsl_vector_set(ya,1,1.0-2*rmin);
+ E = energy;
+ //printf("%g\n",E);
+ int steps = driver(radial,rmin,ya,rmax,yb,0.1,1e-4,1e-4,save_path);
+ double result = gsl_vector_get(yb,0);
+ gsl_vector_free(ya);
+ gsl_vector_free(yb);
+ return result;
+}
+
+void print_vector(char s[], gsl_vector* v) {
+ printf("%s\n",s);
+ for(int i=0;i<v->size;i++) {
+ printf("%15g",gsl_vector_get(v,i));
+ printf("\n");
+ }
+}
+
+void print_matrix(gsl_matrix* M) {
+ for(int i = 0; i < M->size1; i++) {
+ for(int j = 0; j < M->size2; j++) {
+ printf("%10g\t", gsl_matrix_get(M, i, j));
+ }
+ printf("\n");
+ }
+}
+
+void GS_decomp(gsl_matrix* A, gsl_matrix* R) {
+ int n = A -> size1;
+ int m = A -> size2;
+ for(int i=0;i<m;i++) {
+ gsl_vector_view ai = gsl_matrix_column(A,i);
+ double Rii = gsl_blas_dnrm2(&(ai.vector));
+ gsl_matrix_set(R,i,i,Rii);
+ gsl_vector_scale(&(ai.vector),1/Rii);
+ for(int j=i+1;j<m;j++) {
+ gsl_vector_view aj = gsl_matrix_column(A,j);
+ double Rij;
+ gsl_blas_ddot(&(ai.vector),&(aj.vector),&Rij);
+ gsl_matrix_set(R,i,j,Rij);
+ gsl_vector* qiRij = gsl_vector_alloc(n);
+ gsl_vector_memcpy(qiRij,&(ai.vector));
+ gsl_vector_scale(qiRij,Rij);
+ gsl_vector_sub(&(aj.vector),qiRij);
+ gsl_vector_free(qiRij);
+ }
+ }
+}
+
+void GS_solve(gsl_matrix* Q, gsl_matrix* R, gsl_vector* b, gsl_vector* x) {
+ assert(Q->size1==Q->size2);
+ int m = Q -> size2;
+ gsl_blas_dgemv(CblasTrans,1,Q,b,0,x);
+ for(int i=m-1;i>=0;i--) {
+ double xi = gsl_vector_get(x,i);
+ for(int j=i+1;j<m;j++) {
+ xi += -gsl_matrix_get(R,i,j)*gsl_vector_get(x,j);
+ }
+ xi /= gsl_matrix_get(R,i,i);
+ gsl_vector_set(x,i,xi);
+ }
+}
+
+void rkstep12(void f(double t,gsl_vector* y,gsl_vector* dydt),double t,gsl_vector* yt,double h,gsl_vector* yh,gsl_vector* err) {
+ gsl_vector* k_zero = gsl_vector_alloc(yt->size);
+ f(t,yt,k_zero);
+ gsl_vector_memcpy(err,k_zero);
+ gsl_vector_scale(err,-1.0);
+ gsl_vector_scale(k_zero,h/2);
+ gsl_vector_add(k_zero,yt);
+ f(t+h/2,k_zero,yh);
+ gsl_vector_add(err,yh);
+ gsl_vector_scale(err,h);
+ gsl_vector_scale(yh,h);
+ gsl_vector_add(yh,yt);
+ gsl_vector_free(k_zero);
+}
+
+int driver(void f(double,gsl_vector*,gsl_vector*),double a,gsl_vector* ya,double b,gsl_vector* yb,\
+ double h,double acc,double eps, char* save_path) {
+ FILE* path_stream = fopen(save_path,"a");
+ if(save_path!="N") {
+ fprintf(path_stream,"%g ",a);
+ for(int i=0;i<ya->size;i++) {fprintf(path_stream,"%g ",gsl_vector_get(ya,i));} fprintf(path_stream,"\n");
+ }
+ gsl_vector_memcpy(yb,ya);
+ int n = 1;
+ double x = a;
+ double step_size = h;
+ while(x<b) {
+ gsl_vector* yx = gsl_vector_alloc(ya->size);
+ gsl_vector* err = gsl_vector_alloc(ya->size);
+ rkstep12(f,x,yb,step_size,yx,err);
+ double loc_tol = (eps*gsl_blas_dnrm2(yb)+acc)*sqrt(step_size/(b-a));
+ double loc_err = gsl_blas_dnrm2(err);
+ if(loc_tol>loc_err) {
+ n++;
+ x += step_size;
+ gsl_vector_memcpy(yb,yx);
+ if(save_path!="N") {
+ fprintf(path_stream,"%g ",x);
+ for(int i=0;i<ya->size;i++) {fprintf(path_stream,"%g ",gsl_vector_get(yb,i));} fprintf(path_stream,"\n");
+ }
+ }
+ step_size *= pow(loc_tol/loc_err,0.25)*0.95;
+ gsl_vector_free(yx);
+ gsl_vector_free(err);
+ }
+ fclose(path_stream);
+return n;
+}
+
+int search(int n, double* x, double z){
+ assert(x[0]<=z && z<=x[n-1]);
+ int i=0, j=n-1;
+ while(j-i>1){
+ int mid=(i+j)/2;
+ if(z>x[mid]) i=mid; else j=mid;
+ }
+return i;
+}
+
+void cubic_interpol(int n, double x[], double y[], cubic_coefs* a) {
+ gsl_matrix* A = gsl_matrix_alloc(n,n);
+ gsl_vector* b = gsl_vector_alloc(n);
+ gsl_vector* B = gsl_vector_alloc(n);
+ gsl_vector_set(B,0,(3*(y[1]-y[0])/(x[1]-x[0])));
+ gsl_vector_set(B,n-1,(3*(y[n-1]-y[n-2])/(x[n-1]-x[n-2])));
+ for(int i=0;i<n-2;i++) {
+ double B_i1 = 3*((y[i+1]-y[i])/(x[i+1]-x[i])+(y[i+2]-y[i+1])/(x[i+2]-x[i+1])*(x[i+1]-x[i])/(x[i+2]-x[i+1]));
+ gsl_vector_set(B,i+1,B_i1);
+ }
+ for(int i=0;i<n;i++) {
+ for(int j=0;j<n;j++) {
+ if(abs(i-j)>1) {
+ gsl_matrix_set(A,i,j,0);
+ }
+ if((i-j)==1) {
+ gsl_matrix_set(A,i,j,1);
+ }
+ if((i-j==-1)) {
+ if(i==0) {gsl_matrix_set(A,i,j,1);}
+ else {
+ double Qi = (x[i]-x[i-1])/(x[i+1]-x[i]);
+ gsl_matrix_set(A,i,j,Qi);
+ }
+ }
+ if((i-j)==0) {
+ if(i==0 || i==n-1) {gsl_matrix_set(A,i,j,2);}
+ else {
+ double Di = 2*((x[i]-x[i-1])/(x[i+1]-x[i]))+2;
+ gsl_matrix_set(A,i,j,Di);
+ }
+ }
+ }
+ }
+ gsl_linalg_HH_solve(A,B,b);
+ for(int i=0;i<n-1;i++) {
+ double pi = (y[i+1]-y[i])/(x[i+1]-x[i]);
+ double bi = gsl_vector_get(b,i);
+ double bi1 = gsl_vector_get(b,i+1);
+ (*a).b[i] = bi;
+ (*a).c[i] = (3*pi-2*bi-bi1)/(x[i+1]-x[i]);
+ (*a).d[i] = (bi+bi1-2*pi)/pow((x[i+1]-x[i]),2);
+ }
+ gsl_matrix_free(A);
+ gsl_vector_free(b);
+ gsl_vector_free(B);
+}
+
+double cubic_eval(cubic_coefs coefs, int n, double* x, double* y, double z) {
+ int i = search(n,x,z);
+return y[i]+(coefs.b)[i]*(z-x[i])+(coefs.c)[i]*pow((z-x[i]),2)+(coefs.d)[i]*pow((z-x[i]),3);
+}
+
+void newton_root(void f(gsl_vector* x, gsl_vector* fx), gsl_vector* x, double eps) {
+ int dim = x->size;
+ gsl_matrix* J = gsl_matrix_alloc(dim,dim);
+ gsl_matrix* JR = gsl_matrix_alloc(dim,dim);
+ gsl_vector* fx = gsl_vector_alloc(dim);
+ gsl_vector* fdx = gsl_vector_alloc(dim);
+ gsl_vector* Dx = gsl_vector_alloc(dim);
+ f(x,fx);
+ //printf("Dim:%i\n",dim);
+ //printf("x: %g\n",gsl_vector_get(x,0));
+ //printf("f(x): %g\n",gsl_vector_get(fx,0));
+ double stepsize = 1;
+ while(gsl_blas_dnrm2(fx) > eps && stepsize > DX) {
+ for(int j=0;j<dim;j++) {
+ double* xk = gsl_vector_ptr(x,j);
+ for(int i=0;i<dim;i++) {
+ f(x,Dx);
+ double fi = gsl_vector_get(Dx,i);
+ //printf("%.10g\n",*xk);
+ //printf("fi: %.10g\n",fi);
+ *xk += DX; f(x,Dx);
+ double dfi = gsl_vector_get(Dx,i);
+ //printf("%.10g\n",*xk);
+ //printf("dfi: %.10g\n",dfi);
+ *xk -=DX;
+ gsl_matrix_set(J,i,j,(dfi-fi)/DX);
+ }
+ }
+ //printf("Derivative: %g\n",gsl_matrix_get(J,0,0));
+ gsl_vector_scale(fx,-1.0);
+ GS_decomp(J,JR);
+ GS_solve(J,JR,fx,Dx);
+ gsl_vector_scale(fx,-1.0);
+ double l = 1;
+ gsl_blas_daxpy(l,Dx,x);
+ f(x,fdx);
+ while(gsl_blas_dnrm2(fdx) > (1-l)*gsl_blas_dnrm2(fx) && l > (1.0/64.0)) {
+ l /= 2;
+ gsl_blas_daxpy(-l,Dx,x);
+ f(x,fdx);
+ }
+ f(x,fx);
+ //printf("x: %g\n",gsl_vector_get(x,0));
+ //printf("f(x): %g\n",gsl_vector_get(fx,0));
+ stepsize = gsl_blas_dnrm2(Dx)*l;
+ }
+ gsl_matrix_free(J);
+ gsl_matrix_free(JR);
+ gsl_vector_free(fx);
+ gsl_vector_free(fdx);
+ gsl_vector_free(Dx);
+}
\ No newline at end of file
diff --git a/homework/root-finding/functions.h b/homework/root-finding/functions.h
new file mode 100755
index 0000000..8035849
--- /dev/null
+++ b/homework/root-finding/functions.h
@@ -0,0 +1,33 @@
+#ifndef FUNCTIONS_H
+#define FUNCTIONS_H
+
+#include<math.h>
+#include<time.h>
+#include<stdio.h>
+#include<float.h>
+#include<assert.h>
+#include<gsl/gsl_vector.h>
+#include<gsl/gsl_matrix.h>
+#include<gsl/gsl_blas.h>
+#include<gsl/gsl_linalg.h>
+#define RND (double)rand()/10000
+#define DX sqrt(DBL_EPSILON)
+typedef struct {double* b; double* c; double* d;} cubic_coefs;
+
+void f(gsl_vector* x, gsl_vector* fx);
+void radial(double r, gsl_vector* y, gsl_vector* dydr);
+double Fe_simple(double energy, double rmax, char* save_path);
+void aux(gsl_vector* x, gsl_vector* fx);
+void print_vector(char s[], gsl_vector* v);
+void print_matrix(gsl_matrix* M);
+void GS_decomp(gsl_matrix* A, gsl_matrix* R);
+void GS_solve(gsl_matrix* Q, gsl_matrix* R, gsl_vector* b, gsl_vector* x);
+void rkstep12(void f(double t,gsl_vector* y,gsl_vector* dydt),double t,gsl_vector* yt,double h,gsl_vector* yh,gsl_vector* err);
+int driver(void f(double,gsl_vector*,gsl_vector*),double a,gsl_vector* ya,double b,gsl_vector* yb,\
+ double h,double acc,double eps,char* save_path);
+int search(int n, double* x, double z);
+void cubic_interpol(int n, double x[], double y[], cubic_coefs* a);
+double cubic_eval(cubic_coefs coefs, int n, double* x, double* y, double z);
+void newton_root(void f(gsl_vector* x, gsl_vector* fx), gsl_vector* x, double eps);
+
+#endif
\ No newline at end of file
diff --git a/homework/root-finding/improved_wave.png b/homework/root-finding/improved_wave.png
new file mode 100644
index 0000000000000000000000000000000000000000..3c8fb2f4557dc16cc650cf0fb2dab8d7aec9abd3
GIT binary patch
literal 7944
zcmaiYc{Egi`2Wl`#@NTcjBUsgLb61bu|y(7gt8`Dgvh>K`;t&HmTaZ$BC?B-eJ9GE
zDP^azl)d?0ea`v*^ZlLkn{#L8+~+>;XL&u#`?>E}V?$kf8g3dG3`VbaN!tVl13(xI
zPDD{a6o66k3v{DlY;f%&nM{T#xw*L~e*_o|Bokm{*v19{2kY*JY25urCJ4Yl904TD
z*ny8>FdPUYgX9qM*c)+iXeo|346Wnl=011s+?6X=oSd9OLPFBh)2phg2n0fZfB(;)
zKOq$)S=pqGq_T+xk{O9eCU4A!k!xDZH#W%QFi0EB5&8`gd@LyS`1-&5$6Q?QkSQ?O
zK=j}$?8HPeY;c6knoJHcCV?RNBJ2S4SjPL47UDnD-mi<I{r!6;7Zf1Bfo+6sJco_W
z=FTP=<IL_$gRC%+l@-K67YN!BKsz#cVn2?Y3te{PTz4{oK!*Gw69#ItNn}!P%|I);
zb%0P2NXF%o1)7_ib8>Rj)zw>CT5N4?=jP`4`T5=4++f+B=I@sma$uw%nqq8yrVs?w
z?w2e*VK9dFlP^3)f&l}A@jlSg*0|=Ixt>MqUs}knBR}e{rpkb;S1Q_PBm&L)*Rr=r
zwzi9p%mm{q<LL3YQAp1J?*lg2J+EzIkv^|97bSls;P2wLw&>xhqrD6D<Xg<Ec`FI$
z-iKQQiEXiTeO09Xy~lQEgiGg-ObE>!>ir%=vQHyH4|3P3FS?(37xmy5qkdZYIvrg!
z7ANzzQ~lAw-It_n()+Frz%gAVzpd(+E0!(V8WU<iNQ!TU6sr+NX{Tz+;!bDAq>hR>
zY%gnWC)8xpZSn8p<E5_&4yXkT9@b~xl`gUSop^)_6n@&x*y+8i)z`^RA<D~-pZ&8U
z>D$5+<iqpyQq!t}vz*SQ^n0prgf&{;>nZ<Lua2&HQ7L5j!gXjNtwDH2atxkuBa|X&
z`pR{kn;JT=>Aznub(=hJSn4pl4#ZtiuBI0K{%R128_<5&7t($~@T-@@t(2uYqbcuC
zPtV>K{X5sA_`L%gt4UKgdn|C-WUZNnpXc(IS!T<wv%eiw4mQqf^>u0AS1^Y!np!dM
z+;kNhGUw6UFb!5&cr0ubVBzCwFw5-4z6I+lBRqZSc5ySwshIX+I$g2*(n2JS_^qcd
ze`O_sKcg5k<KnHQeO9BE-6ERz-;zODc#5!iv1WaTLQf2<d*k6|uJwfthk<Vv-%54o
zc%2R9)@u1FU(a!2q;(YE;H4L43<|D%ndLFn-}Fz}oHxhy;T-$9IxAO~cKk%w4;<Ls
zC$&Xh>!Vyd)%$AjVsZAIj!#~@9DH>yc)qldkKcJ58`4Sk*84(V)5x!HhIH7UmX9vZ
zZ2Y~{I8S$(ar0Nhu{iNjoi%Wk(q}toZz@lCeUI#=*fmxHLXMq&CEav0e^d5xUv1NQ
z)TKenQ|sL}gxl9`floHCUuLO96@hvB+Q;4l-*)^g6W+2Vl}_8D#4zd4hk)W$nTtfF
z=~dUa2o+b#g!tQeukiO8Uwj5R0IBsi)%pQfDVU|WH+Zv3(jPo~y3K9;H)P=D-UoxW
zKX2>C!0;8RVb_hsZP&Yu2AAmz;S?WrA8VhxEwx{?{&Adx3f?jW>S`%B)ZU%^+)OVP
zEG^0Y+09vsIF@PpqD+Z32U$Ps@#q1@z3Rit_|_9pKWqfkj1rp?b?hmd%D>^`c`3++
z!UN??&GbgydmaBweJAW9wjv9x<meY0rxxq=TZ!(~?h5*3o`2J|w&>G?#o%SB*zaP*
z8Cv{{vc{Cw>$Cip0j;YK@!Bt3x>@ch0#eiwWd?umHVAwy^RBm8dGGUq_-pJkZm`+E
z6-@#}8U<umN`=Wm{Kg88{{X!h2W?E4*{ic^aoQO{y%;AKrixU(v2{#+c+j)_ktn;V
zx1UEFreD6<mi}@Gx=^=n&$jv<hSt2=nx&~bCE&p_*_p}k?n8*IsO<)c*W_`0Z`H&c
zRRs7F>zU@CQ9r=#a6OTq2@yVOzZBm>r0L#P-m#e*3_Tx{^ks`xzHJnUTM+Ww0R=p>
z-+FO1f290>HI!UOR95!69J^Cyh6IfLgK{0}0lKu$aFy7McKdV0FijjuVBYtWmzbt*
zkvYM9SKm`5SyWBMP>rumZqUoWZwd}L;hRV^b23KDG77*?Rv22Q>PO__->46u2T2am
zFX{KKzQvv+EH)FlBW6eY)zAR0SM^A<Vc9)nq~fVMmB0;$)TUxm{<I7Qi+eF9YbE^&
zf6I^N<;RBEd!%=ijvnU1T>d?*{(0}=eNq%C&gA`)x?p7!R=CDHN0Y5|M+y}~5#dDi
zk-2~!xFgM+6q!9-!T@hI9Uf-!bNa_{Q#OWA7%>)#c8Rsha&3|N`o1_bLnvoSH2yX#
zuo4R53qLP=(~h4h>HU3luXBy)d3PKuP~mt}*|pK#+9=~{xrWrHUp5*DgT?FDZKS?V
zGfzn_dERertq>G2x6TrED*2f%Dj1&!_7!`)sLpA(75PJ5d*mE0B8ry59RQT425bmk
zW=i#^=#`9OxX}or0I{Js^QQPyBQritnV9gmNJoP*@d9&+>>SQ|#@1%YEmZ4}2io1z
z%y1TPEYM%V;m)8dELPBHR(%oS96bfX{68KyN&tC@5k6a_g)^2+e;z0@2}?}Y>k&e<
zu%S2wj_)H*y{8WnZ2;ippi?0&t+Et`OH&D6p*gWqHn?wJkH{^?^>F{OwHj-S+cX^R
zl^f27ouT8Ab_Hj#+t-|@HU;0Zo+18nnYkJ9mZsSV2a$VK#z0!KT*DPO*Ufnqk=NKe
z6oX6Fhlja}9*L`INiP{c^lH=tt|cuwRqC8u&FD1+k%)(Xk8_FT^ck6~lEj7I&8T(y
zz#Q=PLr#{Df|6S%rQmh^x2S$z*ISyt5f`HzsOxSB8V2TT7wA_@AsVVG20N>7$&iqc
z?$Y&`Q+EsnPZ-{&j_&kEJ?Ogq*Q<@*rE2r`1UAN!I|l6{D-vPw+n#zC?xWvcBocXK
z+~cMB(Q&?nB2ZHE0)?%pomxASd#!RC?x(KKAbnX5{}dJzvg>n*j_eb_gN?%9xh;!W
zyvRXEQ(cM}5uX<i)hkfnQFXR2^{e@=1%$JVaE{ID)%JD11jZ=OhFfjPexO9PZm0+m
zabu<`+ReYbqTly*YMgn%WW|YsCvRRux#E9uFdCT_W~Ycn@Q>yUg*Z)9Q#{$VBp6*O
zU1W1EtFso3C~PjTcxW6HkiXh_<92JEla}tO5Qo{QNl7Nfuh*CrD$DR2lfw~Tst#J;
z3x%OG<L@ZGQq*t1I<mX;I`C2OLcatigNF1TMK7U`y3=`-VvchNVJm3hInH$R0~!L>
z6ecT{XQv*}xmq-TqPgV9C^T{U_r7yFi~HbKig|Ysa`wf!g~{PBcR4ir3%T|hPGaV`
ze_2AHgSxN`CLK2{hH5dwjeBY4XFHz1%tFfX%3~P7(Z9W=c?siR^!gLx(J-|w4gW(v
z7R|3}GS<K5ywHt4e}b?RURg;4*R6!eF2g@w%fOAq39Eul<??f)0(aNCA|!3^Mncrt
zDbG8U!SizeG@dF)qmjY$HHO`gVYpsql4ROf%6E+d@o1MC!%Mb<*B|0Jo#y-VrUfIw
zdXwj1%iv2^MGDV1VsKTtk9nc#5UJoV7vi114^vq3)o0)$(?l*x`o0k=E&dr|=YID$
zOugi3J={<N;sC2;5)JBTti1v*q6~F30DoEpoPuy^;+XAJt`DkX02dO)VKlN@&WPBF
zYQ9b=Ip~v%Q3^@%x#EYNxN8uTn->}#(PV%dkLGd}Gse|I$;taw!Y*kqRzT%D1RvT=
z7_kE!Cow@!Jffz@X*ick<WhAH+)m3N8sn+4TN=x5@PEK?!nITl>@K{Wh=AI}ljF1`
z*J8%-!GWvMAiDrcKeID}&VX=Q!X6)ZBbx<%jE<d7L-*X62FKY_ub@ZQsj(9-2A8&N
z;N`1lN+48}pmpZCY-ubblw(rd&RI`#Kj8@RF=ujvN6xRh%%JXYwfE!wWT>zcSrfnt
z$CIi-or#AZbZBQSVDQ;kNn+2dhhHX982J&J{^%v37YwkF8&tsj88xw?ps2!Df})Q*
zC!K`$W~6HSK3IjK*9O~w?12iu<R^_ZeGRCV_~3MKZAW3-^^9q~C9%FDL4_}@#d$)`
zN@9a0aJ|}yc&tidN%*@n9iIhyRH-c`Q_;Q~kGnA}&F#}(g?zvGxeP}O&to}Oh<6=1
zKkYIz_1#!~&D+uvC2}vVOS$0|_XakPF*Ok*>6O64D9GPtg|B!9BzPft6;0i~w8YK*
z*5HrwBkH_|R`<EIr_K9`dRYU^4+nnU&eruQZBw-pY+>W{SJrDkjm>mz#!k=_DqxZa
zLVkXjcxX@ic?rU!Pv^P2r5oD3I;+BsZ+0YkGqqUt6gmA?NwLu6PV}adWF@qQ#Z4jP
z&FEv0byG%@$Lt=P@f~@(@c0j_OAPWex9$cF0Sq43oT|Q-4Aa)_m-yQq1OZ>?qP%Gw
zUWiu|qYd)A(xUwyb*aqUswY?MNTuBwjZKovIDK)H8sgR=B+_Th;05end#T-HbGLpL
z%rpCrZa2+<+k5uLHzLSQQpUV<xc#P#@NPgvC0)lyvQ4$~MdkMnNr7+oJQb?9_D<z;
zpY^!|O+1{L%9Jh*MI0(mHY}n8b~3wZjE+KFxEZ!i2a5B|{thdex=bUcFufawB~mA7
zxD){;G_B%09-k!hE-ZZG>!Vj45#v#Oe$I3F*C810!r6<KRGglp3M>GTPd`<oB-n;u
zRP1I*Lp#rP@duQ6Z8*^R90BC#Gr+ULC$3QnB1GqZ0)|LUg5wWqyBTk8tI$}#Mo>oP
z>^|)s!5Ce%DZtJf79`|5)z)HzyKOgu)F8loX_hD*+%H#h$OL2P%LvcMu(`IJd(U~b
zn8GOBv+46zdVUEW1;r_#03^W_XM_$Pg~KuL0aK`!o3D18;v(#PJUE6+J+{tJ42h$G
zhU7+&ky01&?HUe8i7>dx3moWTr74)gW<x<!e^i=?-vaPgjh&k*Vl=59WV{E5Cqw^H
z<`FTPd~byijTG?I9kk5FEoB59;NVX4-_~TvCE@JFDaue~S-5Z&PFD4mkzvv!%rZZm
zcT~D}k6F+zW8J9Iz5#b7P3(t@fM&&or-&J8CehWVx4fcRRlQMnDLF3@)XRTl>l&pr
zU54-Kw_mv*Nykz0Ab^Gw^&ZF(#e;dCqZjrQOW{gpt5j*u^LK^pi1Y1Cq<~1N>PO&t
zPSMg@%i`5Z^x6!KHp+XE3m$D#P4W5p+J}zIdPcNM7L=MepKMes{e0m1f$M6=S=xg8
z9mkQjLK_RZf3!@>sdwMzoWI|(u8iHv?%%ti&^+4TQZ%F)Muixr2~4x7wO7z@?}@Es
zb18VEp%DLI?_KYVnSX(K?&7ymk#bqYFMJG3`Ojl#!A08~_A}@DbpUu)F%Gq4xK{pM
zJM~tT$B@cwgS~=DyVhdk2qCtVU+2;%0WAsA)fKdiBdbRt0Kwm5iQ6`GAC|576^+CB
zc(BWu$!Md4(WtbchqijW05U<d1}S)gzx?upLo9kY3yJ`J_FxDds!&xH5f~Rw`QDcm
z>&ses4h?lO=&lgfMtkQoO$4|xhPCr`msr3G5UFKIT#!yp*2=q3HbTONKB@coqSE-T
zC}h9jHfD6~P|<7#L~_nDBh@sTfvG(XVVd;0I1Q9H6_8z^_!hwCUhqP3QneGQK#;>R
zDRQMqi(qkeDSermDQq;k1a#8igi@>w?D_Rktg>7Cv!jF^-^E1xm3I)I@71N)&kiSi
zWuVSJ_|GnF_j58*Xi}7lves`XW6$O#X|Tnh*AkM!jMIjyN{(%BzsFi3cac!|N-0mN
zx=rJE+aCqYW?1h|c8if-vT?HL!>Kt|&UQvb`6oAzluvr)Jc&m{gr0rJlJsZwD8Ft0
z<Be(7Dfz~YDtl48A7YPExuAkrL-0s_U1oXzdZhGZK@*~>=)|)x4{(hP-97(=pGF9w
zbla_P$5w96ng6)@5$)GMzN=ubICppd%zva(O_KFLbB`Zz_BOiS{L==O$T#98@8DOS
z)xyu<0}nQD%C4v#-wqCzwb&A$OkAw4wJCl)#7T20V?M)Lx4A0mI>Hqlk;imca(?}v
zUPIiDmY}_;K)r8)Zj00S0|iOSz+A3J!gf>UV-%wNXX~>--e33)gi`rX&C>W<L!VAC
zV8exiLmjcT)Xi_s%Q~cH=|-;0500?So8~-gZraJ&9#hE`M(wNu0#>0{yLbEubXugy
zP-sN62;t!QYl*v};L^s6ko-gP25L%2ys$_6y*dd@bsFe|(4X@bH+SxeIyecL!m9?Z
zw^~FI-?|;frcL_#zEhar5j8t0<}X?3j9aMSQo5a!Fm0_-I%8I#u*#u_8FsV)(`b-U
z2W$NW4W+mVZWXut8qE0Dhu(nuIz*(#c^UL@Mw91tZZFi{9ad%Y6T@)0!=30g4V=#k
znq==WWUGh6Q4KRv>+%9>ypheaP=YWZ9$ju3!Xh4t>cVOFKFX7M@R|t#-niio68VY&
z&Z7OY>_0(WP`5IGx|IPu%pEO|doI%PoF`NPdaPC@C#_K+H{hV-q*8eSkymN)P6kl`
zyF0%>aEuO-WgTv9xpwuAEL2atsHK#@EZ<3J833AxGop&d=jl=!l-N`#4Ul0K(oj1)
z!z7w}55JVq$K~(Easw|eK({0O>|}UGetYHb*W7Tk2<zK_214s0;>ap1)DSSqv)t-q
zHg@6w8d&skT+TfW+wWw5eAGTRvF(VmL)aUANBLR^y0RKibRA4rf2U!os5@Hd$ODSg
z2M%9+mbKo{DtMyF54pfc-x}BE-)eimBPiRpy#L$>)JLw;$v^n&?{v^h))aL*CC}>=
z)AhX$5n4(eIXQ~iSTBrIqxQ_vv!MG$13IAqd-|e5jJJ@#wV<JF!=?JYvPzO5Mk;c7
zl3~*g>BqGcwDWpdHM#4rxa*3CL5jUk%|{M|0lvuue&)0aA~j=`4g^ER<feO|66bl9
z)Cl3e9123MJhfQ?-uIy(N?=Vh7!c`?6s2@l<@G|P>9Kmqp|R&xQ&4-UTT8K@v<UHH
z%G-BAP~-H4-P46yguqwEVJ|N2XExEEyedA}TI4SnSo*;WJ&6}nB^xEb0Ldvf;Lm=S
zWU6}Vy?rmeY*W@s#S@VwpkYsyczD6}3xg{^${>T9aI7^jKnI4*fby{ZBDS##hcYy}
z4NRGZqGeRQA|aqX1?+h2a<PzM_N<jMCsKY>-?0}!7z|}T31Pb^rESk5%TiG+Te9Nl
z0F9FAI@q=<4M<29z@5#2FfdBN))oB<{73=MGnw9xC}4nZKhYEsjX-B+-=maM5<+pZ
ztiO6#Fj*+$Xy$tv8+$T$){VN1Sy;(uo1S-|Ewq)}ki@k7jwCL8#YGT(GbOow@bu%o
zg4?onN@q}<i6w0J<t7W|tv>tPz+<3sG{f3?valw%SoN$Ml*nYUM_x{-)=e4e))&>^
zZ+?l@Z!yoDYMoWc8N1)`;a*wY6Y#Z?6t~aoKj@jUAXh%^7W%v!avO3OGY*r_qKgJc
zn30By^MjRF-7v-;-rjdH_FlE{A@(O1dUR%i>?}SA28^b=DS2NB2Iud1$afsqa~bDM
z?v2${hK7dGZEm<1G|%k#(pWIZe7qn)oF@>@c%dxJSq~IWpe$#u*k3bZ1@afKU@f+?
z&wZc(^1ody@V~+fuyaN#&}&2xcbbd{PfnsyPTM;0Qew-*gPG)|*llO$1SpG2%hLos
z<5M*_h|?q;B6uW%puNfn-uLr2y#^{R!zcT<CM1YAQVq&EzbdfV%WhKCHX4lFUWgJv
z^oQ+c`||>L-bhD}#LE?H$<ZME6P4cIV(E|9SI`wSU;PfasP_he2~uzKrGI>1ei)|a
z{>K-szhOd+588r!VRCDvLCYxGd|@utB2?xL`*`S6_Nn<nhY3NXALgmA){+y2iLSAA
z`}9}One2*A`d(kEUr>E`{$Kgy6?sL0{7F6_cpG2x-%!-|BHQLR1=ctEp|HiDp00oI
z-Um@DuYNtFN9{$v47;5ZC<*TpnH$F#-tFtLRQSsLNYYosR3-szWNY}M_)+c4{qa7R
zZ#mg3<-K_T^Zu)G-cP0{hTYyeJaIEb4k4i+!C`30IUeU}E&4X8`uj}C>oBfT>g!`x
z4-X8x1ukZaw@ahL$_wjHbFM1%d!N`oeb293r-NfM$Db`8Y6ex>Gh0y}o0^t|J;*6A
zHB;im)W70}UP}01Kb4#J!!{%^(7t^RLSjFVPp43!J0h<WM+blWN{w<AU9l$Sm~vDI
zMS~Ge@j66dzxKo{UZ6QX12~J2eCE&EuY)U>u<~a0xct|WNAmCo7CNVjeHIed|BwML
zoDDRm812901lSFi8NOicCa-oY))$09^BZEvV8e`FlYqod<d|rkMU+Mh@<LImjV4lZ
zP%(vJP@FRRUFMpo2>}FEey$xIT%w~v&_G)X_$!ErHga8n_N&}0{?hU{!jF<M*CNmr
zo9yVDPrc-#z+c4Vl1MKUbY2o+qZvU&0&1y87HlV)_UvV8y3wTSz{=qY1~^k#*5(&9
zwk{O>x~wW00ScgVST(qSFln{KIE!cy)ynh<6BmP?!4V3V)pPXO`k38oBW4sZAwvrk
zb?$kuScYioZfHtbD}d5DXKVueYkv6XcT@HwHQjsLLe%)w)WnIdzoU{SFWDd|MHz$l
zh7?nK_}PtX$T16Ps8%*l8ISu5dBQC>TJKXK0)<RIDkJ*O(4;&+x{`QnS<9XF>!2s_
zOysg?`!3xlK0Z-0XQS{Bz-j*TZPH3|t)r#9l)HkBHL_gI?}lMZ-1hW(NgwgEDC7S<
z0$oxH&~D{iNU`Fv*^Nm%&Fgl0Qz^J@CYj|o&iPmH+RMe6eUgiR?%4j;p=~2NGl5+p
z%71V4(Sxq`#QR$|Y7PrFh$n%D<JP$eETM-ZREBy5*(~0d&oCPLFIKKLm0P0?b!lIm
zMFj*Qu64{Qsk6VnoER{sLH`?FQNTy#n%^~AG-o)Y{<nCzXsDXDbwB*TNoGZusAAl)
zX{utJKD;RCJtkD`wtSirnX@#CUVQcfv2#S{=c#5j!I>j+Gw*G40KJ&~poRfG+fGHy
zW22k^jG=Xa`f>hu&%8w!0lsl5y2Wt$C!5R_S9NzM!k18W1$+@|d;ra%wnXnfsdQbN
zKkPtkV*i@mSUDA2)DFP#Vcd)}#)v=4zvHT4xz(?mK~?ywsE3Tvf#NvmB!A&6CC$YZ
zv?L99=-sY7?u~w`^*#{evBCfK-8ymmn~&PIlaybMDx$TP<&!T>cT=(Yt9y)#iAIY)
zq7RIh)2eWJz{c|N0Xc6^O_#B4`HuaqBj$k5>IQ5tx4WF}0K@cKTrzWv<qaNusxsv$
z6`Kn!e{E3^gqV4u-17|$St}XY2er2XUEAFKB*u@I+J8wk+4kLS@NQPU6w(9*nkmUZ
zV|brpI3OhJhZo1wL+cBE>#+eILzr)VfoxR&H~jxg5jdwT%h2Sy?@Fxwrh@eNNIPJW
zi2bd?p391qUr)?vigVMwX^EP+j||vx>S<kF9*_0$Htn)N6^6gz>!h!YN_H*XK+5++
za}0fUR&@EMVbe<VK>l_~wh_ziKbgCxD88;^6Rv<8Ql{xNbYH0pMrGRGAx)vv0^*)*
z;+v$$OiEmBmCHNDInjHcYLb&{HYk5wi}0wgRoiz|6S%d<gVi-!QrYZ0Pez%E$S0cw
z<`kyhh!J<slm>rMUTblz4BkmgU>qBg)eg>7c(fJ~li`*PxMb6XyXuZgm#1ZlxMr3X
zOw*o!pBcZIA)M0#IQD=4wC9j<Y#L$l=E)f88crkwe4c}vkhz*(oRk$T#h-LmpCi<j
za&12!`TT|$yhmyz*I#7EfS{vvJ{O5dl{CNdfKZRjoC1$-60VMDVT}DHGoBRxEbt$_
zLlH*HUXpJk;~dX4S!*58>OmVs->uE7OEuz&jW+~(_w(ANJ8Vx2U7Ib{kmk&+vXGyp
zx|T&8soz)2_4i~o7LyQ8Uf1yBPkI6{R4yx4!@J6N-1*U*OfFOJf^u%*5J(;EgTYcG
q^(-N=@NDgg)7eq*{|zMU>7KgdM5cX#-ZSJ0QSYLm_FGN6u>S!oOJ$n?
literal 0
HcmV?d00001
diff --git a/homework/root-finding/main.c b/homework/root-finding/main.c
new file mode 100755
index 0000000..456dd1e
--- /dev/null
+++ b/homework/root-finding/main.c
@@ -0,0 +1,75 @@
+#include "functions.h"
+#define nr 2*64
+
+static double xs[nr]; static double ys[nr];
+static double my_cubic_b[nr-1];
+static double my_cubic_c[nr-1];
+static double my_cubic_d[nr-1];
+static cubic_coefs coefs = {.b=my_cubic_b,.c=my_cubic_c,.d=my_cubic_d};
+
+void aux(gsl_vector* x, gsl_vector* fx) {
+ double energy = gsl_vector_get(x,0);
+ double M = cubic_eval(coefs,nr,xs,ys,energy);
+ gsl_vector_set(fx,0,M);
+}
+
+int main() {
+ printf("Finding root of Rosenbrock's valley function:\n");
+ gsl_vector* x0 = gsl_vector_alloc(2);
+ gsl_vector_set(x0,0,-3);
+ gsl_vector_set(x0,1,3);
+ print_vector("Starting guess:",x0);
+ newton_root(f,x0,1e-9);
+ print_vector("Found root:",x0);
+ gsl_vector_free(x0);
+
+ printf("\n\nFinding bound states of hydrogen with simple boundary:\n");
+ //Making spline of Fe(rmax)(e) since its instability makes newtons method inable to calculate proper derivates (See func.png)
+ for(int rmax=2;rmax<=10;rmax+=2) {
+ int index = 0;
+ for(double e=-2;e<0;e+=2.0/(nr)) {xs[index]=e;ys[index]=Fe_simple(e,rmax,"N");index+=1;}
+ cubic_interpol(nr,xs,ys,&coefs);
+ gsl_vector* x00 = gsl_vector_alloc(1);
+ double guess = -1.0;
+ gsl_vector_set(x00,0,guess);
+ printf("Rmax = %i\n",rmax);
+ print_vector("Starting guess:",x00);
+ newton_root(aux,x00,1e-9);
+ print_vector("Found root:",x00);
+ double energy = gsl_vector_get(x00,0);
+ char filename[30];
+ sprintf(filename,"simplepath%i.txt",rmax);
+ Fe_simple(energy,rmax,filename);
+ gsl_vector_free(x00);
+ }
+ printf("\n\nFinding bound states of hydrogen with improved boundary:\n");
+ for(int rmax=1;rmax<=4;rmax++) {
+ int index = 0;
+ double rd = (double)rmax;
+ for(double e=-2;e<0;e+=2.0/(nr)) {xs[index]=e;ys[index]=Fe_simple(e,rmax,"N")-rd*exp(-rd*sqrt(-2*e));index+=1;}
+ cubic_interpol(nr,xs,ys,&coefs);
+ gsl_vector* x00 = gsl_vector_alloc(1);
+ double guess = -1.0;
+ gsl_vector_set(x00,0,guess);
+ printf("Rmax = %i\n",rmax);
+ print_vector("Starting guess:",x00);
+ newton_root(aux,x00,1e-9);
+ print_vector("Found root:",x00);
+ double energy = gsl_vector_get(x00,0);
+ char filename[30];
+ sprintf(filename,"improvedpath%i.txt",rmax);
+ Fe_simple(energy,rmax,filename);
+ gsl_vector_free(x00);
+ }
+ FILE* exact_stream = fopen("exact.txt","w");
+ for(double r=1e-3;r<8;r+=1.0/64) {
+ fprintf(exact_stream,"%g %g\n",r,r*exp(-r));
+ }
+ fclose(exact_stream);
+ FILE* plot_stream = fopen("Fe.txt","w");
+ for(double r=-0.35;r<-0.32;r+=1.0e-4) {
+ fprintf(plot_stream,"%g %g\n",r,Fe_simple(r,8,"N"));
+ }
+ fclose(plot_stream);
+return 0;
+}
\ No newline at end of file
diff --git a/homework/root-finding/out.txt b/homework/root-finding/out.txt
new file mode 100644
index 0000000..fb3a8a4
--- /dev/null
+++ b/homework/root-finding/out.txt
@@ -0,0 +1,58 @@
+Finding root of Rosenbrock's valley function:
+Starting guess:
+ -3
+ 3
+Found root:
+ 1
+ 1
+
+
+Finding bound states of hydrogen with simple boundary:
+Rmax = 2
+Starting guess:
+ -1
+Found root:
+ -0.127929
+Rmax = 4
+Starting guess:
+ -1
+Found root:
+ -0.483549
+Rmax = 6
+Starting guess:
+ -1
+Found root:
+ -0.499289
+Rmax = 8
+Starting guess:
+ -1
+Found root:
+ -0.499977
+Rmax = 10
+Starting guess:
+ -1
+Found root:
+ -0.5
+
+
+Finding bound states of hydrogen with improved boundary:
+Rmax = 1
+Starting guess:
+ -1
+Found root:
+ -0.500014
+Rmax = 2
+Starting guess:
+ -1
+Found root:
+ -0.500873
+Rmax = 3
+Starting guess:
+ -1
+Found root:
+ -0.500605
+Rmax = 4
+Starting guess:
+ -1
+Found root:
+ -0.500018
diff --git a/homework/root-finding/simple_wave.png b/homework/root-finding/simple_wave.png
new file mode 100644
index 0000000000000000000000000000000000000000..d032b1aa7191d94db1af02ef66ff14a135fdb74e
GIT binary patch
literal 8631
zcmai4XH-*7w<aNhgc=AX280A@3PPx&D1=_5Ns%rQ>CFI&Anj13D=5-iRHRE2X(9xv
z0@9J9LKHy+!B9mIxxx4S?w|Xub+b;Ab>{5ZvuB=X@7c4HcoQQXRwgtP9UUF3p01`T
z9UT~;qXSW3^t2W*yW~6Cjhcypxz^#~A+3c(A{~8^>F9t%GTkBF+8P;8*V{{{=KkrB
zjG+VYWZ)1-021lw@BrN*a2S3#{`Ax-T30+}nAQ)CM$5^`85$b8xVVIehi7GFJ%9e3
zOePNu3@j`x(1xJO%BHNPJeyjiT2Lv6hie}r59`~j*47Ru>1d<TInh4Dg%X9N5-tDU
zCmuUSq?tlT_cnHLh3?42L%P9{L$1`ra1$y39BR?+0e;`oZ&MeOUNv|&#xQShe;@&v
z!>4p>k!ug>#y*lhrkLO@Jf#6HI)IA{z|$@OAdmsVA#h|r{*Xkw2!|x%A(?zgo6jNn
zZ9_iwkV>k5+kV*omRucjh$kIlT3cHS3JO$JRomLy92^|x=H>(i1?lqrFV}rrETE%)
zQK$2(ZPvK=l(qzLqOP?+9UWWe(HoQ|!Nx~NhoRHcR5K6EUCo2vetHU_t?-_xYR$-2
znELV%Dj)Tn1l)=I^ELR}tyJ-hP9&2_$zz(H|MzA5v^40A#CoLMzKr3ehdbBb8`nO3
zv%BM28?t}$-7N>f*|fZC2L@?Go%{>t)BBe6qS>`BhF!lZ0dENv^w5I^Jppjv#%i(E
z+<_4&_Ehkp2U*j;KU{fTCj>(}lMfcib(}W>iKw-qaqVe<)&y5eS>SYEzt|s#a@FNO
zTK4Hr6XbH^GY+tysef(_l{bH>NImQ|&I)^EB>QmMZ<fun%EuX5@#lUxxZF&**8Ere
zZ~Ck@LrcZxucX;j)*sJWv9ph|C$8x9*5_M>S7<8lsGgd&D6?^XdXZDr@j}(Fo$rtO
zmRa{Yik<*SZ|m*0)9be0nRRL*rSAPMw?|%l))xnLVbjh<%dIC0K)7Pz^j}mz%ZaL#
zceXDVdS~6V*0bNr+DNe6K(0IbCrG_?G5B*=7xeY|SbVwUY>_FL?fS`Y`ju@muNH(U
zCI>Z~r!#tb3o0tJcIDj0=hP1q@2fhV;eLK0xTIE+pyYL)sd41R(?+L?rT6N-dFO7n
zFHbXbyKXsWhZgfez4ndAUzoA#_wNTbOPSfvNB#u&Aw7c9w9HB>KYxGthptrfDKxk3
zZrEm-Ir!dnYp(e-&80i3z;-j|9%EGu`i<p*?N3#qoP=o_Qq~Dd{EGb0+;5pl8hNcr
z79Y=mDsrzFfbYfjB3EP}ynF9!4@IPmCvO@k9`~*}ZKl)Yu@t6Pg<enobKN(gGKZzP
zBVw~*gx8M2Q8D=9fyluR^L;&X>ynm;roZj>oT<aIt?A`uv_Y5GP(EBqpuR+fsFif5
z$06^}MB?1iuHy;JbLL+Cx2EKpk~H9^#9*n!F9B^3GspEV+(%*3P7)kc+#o3mnl|4*
z4szD1U3>y8$G+nTEct9cPzlUmMGUitzA~pX`xw~3L)0NZVpiZAU8@6L3re-Bwq<vm
zJYo1uWX2uj8X8u`Q*G99=QgUp_e2b3j&I@-G>K)BC%&?uxUj#*QgFo!Vz&)ya53K=
zEvhMT7BmU5x01}B>k}fOOZAKE!8`UB23@XO8-$XrS0}2|CpX#D9N^Db3m{y19E#%x
z$IT48yapAwKxZLacjC*dVuR<?=8_`2Zo2b;`z{g9{8NF?Ge0#<gGTEInjgB}Gb1iX
zW^`EH^i9x&%~WUS*pR2i8*P#z`?tj#)6-77g72jv>J>L^KVk@G=^AgxlN(Xcsv&{$
z78FziZ;U<*xqsiSd1ECAKW$;dBrf(QQwCC=YJO=$SR>1i^Vu<|*TPAUJ#C)q^I8*?
zQWIXeSjgnqc?vcCjz?3BUi+uh%G1fYv6NCy1A*6p3WKggq9!q8JwfER(;8{Kn*f5!
zG=*qhb6k4aEFe&9MdW>Kv3%-l&J~V0C;i~%9sWj^v&yf3cu78niSH#RwVG3PAJ?_A
zu6+Xkw;)1=S$$Nkn;olu=>eYAHLLgV9>&KW94Duhtn$0`OM1XyQkv3|dGCbHZ+df8
z%ZIKD%k&KtJ+Aum$#qGKt>d1s-y3<zK>a|-sdF;lS4qkde(kFN-q&!MNIi)Grqi>G
zwZwR^nd^hX1o*6x0mxd0O}*^TDvPJXiNoW7Q+()xMbxBm$BnzdWRY{3OxZg3qXF?d
zi>59Q)vqy6^{%ElEjfoSy{cr-e5=s=bU+fyf4%F)fwH|ArZ+;fGDK*CE+e%!O$Qbx
zfuR0T-{h#CQ%QZm$P6nP$zd2nXZEHU7orInSBr0q%IaNYX29;+m~%EcE4HY#<<6Wg
zIMxt<+{g32&lwo#hx8U3$eFWAE4x{y&9491zLUh}FgS9vB#ARo{^n_ky1+bx@9wIp
zkJw>@Pg0wtke9Z8EjZn){^S!nHup*}Z^dR=YzEo10c)SE<t7eDpXuKfZXFv7pP{}L
zSesKpd4P|GNLSOX;Ue})i}{b?Hmwb#p%>VB5L44IrG<<Y{j=)tVxOtwlB-FQ_azE+
z6|O{;LcP#w4o1#30tU9xl%k5v=sMAOyP#-_C#}n3hS^xiU74$cLQyi@e`nvs%O63T
zyC3Fnck+U;-vUSE922Ca%>)*x8gEW-XC!2)a^C~!*Zb0U0@B6}3k-i`4cC3u;q3zc
z6AHQH2@b%{eNvPFg{u~XYk{EGqeEif=R$s`?&*DC(S+3_>qHBye-QbMv}xiAf_I)o
zAcl>mvJv}i4WFEZDLPG-y3C2;aj%PEwn#^sqTCS3?~t&VjMrobRA5ZC0(R%PPsbZ=
z^3dGSRy>aSV;CPzP*(F)Pg4HyXO~c-rvjGpkPb^9nrlBsQ#_oVl6NGpMQ(cv^pV|O
z1gX$u0|FXZ#nOR9F5l1B^;w!@T$>H^;yOG+-o0($XD%rJoCCK(KJG1KZ`9?O?DM|u
z++pwARMa&dDQIzx?|O0+PQ5yhq<~E@fY=qrnb|Y_(h<?Qbai4%9$-^Lffw=oAl~U?
z%{KMX8`3?0CjG@ppFA^JL%=mh;u20z2FHKDY+aqVW!0cxQsiVPMe2~9W26-JOb%EF
zHC=~q<cMjKpZLPmiEMT!2tDJ2l?Gw6*|UCa==vC+@N;PJBH04li25PzlV4}zFZWp}
z#9CEFb+PJ@%c0=aX8!iNa8ZGTQDtKnCc3nys=CJP+J5&|$&ebe+djkXh-suJ+9h4c
z$oLeTFQ+M8V4}Kb`eP|7n!+Ea$Z_E!hh<`oj>tlJtn1+Kva14-;ybq^5GlEp3yaoC
zM8~-VGmU%oZQR9cZ;pMP`U$NIXKskDDx6m+moz>V{h+-A_VVNTEwHracD~2TY6Zp*
zF)R>q>noTeay9xH{8HZDy{)U6l8_G<#%7tLfq=+&(Y<dJ91jGrnA;J5AkW|&C;Dq8
zR~lXbIzzKhcu{y7dPQFhN2m?kpj$72E+b?d!b{A}FfPBxew+R97Y%I&T5928PlA9}
zY28Q8Slzu<`&#b`d44zHEBF9DfPDNr#JX0U5hQbL)sRNAo$lxfHXm!pBLv9`*o{9t
zPoBMzgUAa&LH-mtSm1mnld0%UP~O06ryH$>Qj`%%(rZZ!4GCds(%84|EEFu0-k1X+
zqh!X(JG`wZ9fNJ!K!F+zk3<TkNb|oCRL&UhLD(A&^7ZyGP$M=3G*SKb81^iz<hwEs
z(UHPz@^dZjr)$4De%o&j-=Sd*YY#B0fCo)H7G4-9%tl1yh$6r<@P6dzM$xq@oAk4u
zbYG<nFbw;Mbg{yc)5?+sm{5?|>b-3I7dBr}oN;e;7ZiIIzBCOHg1kRst>7C6h2%C`
ztawld_@QsleqU9LP+ie+#gL9Dd)sZJNCB%){>3}YEG$G@$uOfmLr--T)pG40CGcqA
z9u20E_qng!40~P27FU@Xvzd-uwo&wpcf;r35W6!G2pppya@Z3E9u;c_x(j{pP3(Hv
zo-aVL$YSb&#X?`7iCsDPk_w9&Y%pR9oSoX;MDv-5#3)fBM--O8rpYt>Kz{HUTF@*w
zwv_D3VUIsXJfk=Nq3miL89{ab&tq*+vLm%F0fJRxl;Wv2Mt?PDee|FTNTe;!=*ZIj
zc}M-In0^iFo;y=;W#iD?&jze(_BYS@o*~QiCgtMyNdK<VVy<*Eba7msgimL#yiZ~c
z^tX&+fDTf|Dp%7`tilW|RU{u5r^0$oAN}<*t3NN5D0Mrj9oEkGIz3YU1rKo~+k_)r
zXy&2(U@004=0OmXtDj=JQXad!FE}^k-0(SFZ{Dtq;}piRI*t4k^9oqu5P`#0lc+UW
zi|DBKY??3>V&^vw*=d+ZM4o}1h8{DPkD;gvvRK@H3_u}uh&rx9ac$57D2r8{U&g#c
zmf4KoxD&(h9(IvkHurhvC`srlX!Z2*HAY6xdM+m=*X&V@Fs+3(&z40SHj8NFa`M{Q
zc**@W;Szpj#D~vWjm#fd1WXm4;CsML-E4T0HLrv!9HNKv)5uO_J8zCZ%gYM#$Jm2@
z3M+Hb7|binO|3@8L0@O36az+o8H~rP3SbG>@(eD^O>A}qtg(+#bwR$$J_YZZhc?ym
zK~@(-=C5f~f&~;K0wpKz#R4OGmgEB@zauqR1825!p~zX_pqc$wFL%b!@uZdEFRu~5
z&Zo~P=-1yg#UBh?JYzR2Fqq?Z@SI!6y}tB#c^<{@%AXqrR8)+Q=3LF{sat#E<Y(<X
z`R<iwi%ZhvbNeTY;IfV__l{dZ5=>z_MTG~j??RN|;E1vp$Q{=H^EO3oQ{Qw5Qz%;t
zzUi&GVRok)J^angU8kLo&4GqR!Z1|QaM{YaW5LDF0v;cGc5w$7f>?$t&zK|obKv=p
zGO=PT>GdPzQl{rrB8xh1rNDJ9##x}FIoWsOy;8HAW<jOPt1MIbRGKUr@|(7|tZG0X
zGK6)xf#ncH{V$I)f|6ELZl^YPQ7<b^(K1dUitY8DZz?&|LWf0}%KhQ;na)HteD}mB
z!@2(XSAfU24@K;~4x*G1+Qy_mksP`2EHKv0p0#4=BQ@^#QQ`ZS`P~=98!}9l6PI4-
z-)llcHOTQ#cs6u2>tSFQM5M*iLIo-#h-F2k8-8Kr<TzvB>2u55=LurUF|v^Uj`U67
z2M4ur{ahpOh6l)4LwU00Vm-LIlknl*#eyED-Q`kplCC+qL#;XraFB_?{(!afZg_#d
zUBf|yMX!N4q6+YmAl@9M-kzwv(@PMkcn`ZyUrt>2suEsT7{t3@x)Og$fwc=#<mz<|
zt~(1^<cxeCgpz##^mF;vzD>=VwO!#~EJGWyPTr^KcO**Hw%C!K6B@&=y&cRI9W-a<
z1UZB%OYu3vmF5O4v>d0w8_h$kM03_y8C_p97KXY`4&uwa7c(&1smVI2(q>)A?*uGC
z>eI#oHZWDF(&7NWz3aV_axz{UKLO$VF`k|%Da}iKQE<ECm;Ba8ug&Yd++WjwtnM#m
z3jVC6`i|-3+WGb$jOT<cS9xx^#<-q-cbYWzlK)I9Dn{72!aeg>*m8(nbQe>Dov$!B
ztSeA>F|oFKpl?7NB5m3kKk4@;>y4}TKD10{f_TGaj$DPBr`S0Ky#zM7e4ks`bUb;X
zGlAi#MnO7P8>96$u0ke#a%{yA$xn+0qM#2!QRf-Ha!NS}MqY>)!2}^Q6D1MD$<7KN
z79ZKG<9vFa8fAl<^B3RILN0uh0O~;)9E<XYXTaLSG?6LI@Bu9z*pf~iq5?|O?iADl
zWHK4#v_L;7RI{V%WRx<9@rds|kI@1@2Xmpnx2wjkJSzNRx)=<(k;@?YK{+C>y>uf4
zB*UVY&B&|-VKD{$qgoPrLv_d+Uv^bjs+Q*J>d+g>w1DTnfj$yVrHL}I^r3PuHfq9l
zPLu%;J5ubNB7iVW*7mv%1fB^3X?Q>a!fHl+D3XL>Rt>+u=GtK~djHN21s>dmi=*t+
zBf)#_p99{ry&{_jRRf!VnkWrMPV7{&Cqy`|y{Q)rI)L#Gi!sax*xhc?$KL>t&9{T4
zzzk_XQ*^&;EMf0P?Fh6`7RPIyFM;t4vQBStuI$P23s|-2xTUPx(^Je5HP&f{@4j=}
za16%!R`i+7*Y7js+UB6E1dm;~Rm0eyqPz?1?|!pf8?1cg5m*7c6nKX+243Wp5=2M7
z%o`dI8~7<w@K>KLM@~e!gQZ}U@!Digq9i+aHT5`Pb@V98fXp4&uJ`#g>Fi_<>%}8Y
zc(UcE(}ac5O}?5hM*~(MKJW!@_~_C&+Jw!VXeqlf!_K(l7jW0Ys%Uw`=-Kva;0I&E
zuU(W%tEuip7b8E^sY}syGM!nSs>JeWrWj^|*uR`B<2vc1@LQv>JYoqd8dm^#eU+12
z%(V}mz#xfZ*B&8o+GoTwWRXksms`-W$bs4Xw@1@0)xH5Trh#yO;Srd}&~~G~tn6L-
z)@q2t=R1T>3S_+?&c*xfO8STjhn&R;N~!bpBaCCe<TDl`iF`O1%{@PDJi<06u1ivC
z(K67#xNMk!`ipR^Sn+ZWq1fvz-7_gSs3H7`i89{Cg&vyRDT@L8m`*u}rT6rXW*I{s
zPz1<yQYf*tG$ydKwnv`JLA0D81WNzKjlf^Y+CC*Jdi;tgCI(x&4g<YqK+|_m|I4p=
zS<uid>zgB*_gLMnUM~%2>1O?;8Tbr9#T2Uffp1qbz1xkY4AZFkmpOW5(%1#ad76{l
zhl2MwH1Nz~Ie2`(9E3#{0WKm(LF=GNW;DXS7E!MjyahvTUHMD1#uJbe48`OeXeUk>
zBp^vrAS5v`Y=QrAhg|MqCrS+fG9>(U2S=D2qgli^ss|{McF;Hz4@Z3v{S{4-qyO}e
zRt<&7f=h6&;r#=@e{3Kh$01W$q<ImkBezmBfln%4tRI+B?KE|NO@0sA^7i+?zpxFK
zGVSEndU9pu_$kBSQ=PO-LJz+M{vq<<==>3Bv&u<3Txe9eU+XW2YD(j{h2Tf04mj_`
z#o&OZat$I8nR&4>BnWoJ*|j7Xw9;X}A&tb<@zNqn6|YIRwIX*VpN5oyeDPlMdB_vq
zOqeg^Lt)$)m}x&d;Ed%5IBz6tG;qtW6!37}43M@JeF`vwPE^1LMC^nhyfKK1sV~g=
z3kmzbgBFh2hJ*wQ9LlxwfM*M0Psb(;`*odV>P>J*$ik-dB=+C=YB}7I`&QVY907ZM
zhLd)T$QSYnoUVBB!F)amCI;!Wql9Ya#O0fI>eVGG@H$xV*p#WAqFhQADS;=Au=YJ2
zNZLD^@`F~H3a?|W7y^T{iwH>!0>0ftX^NPtSe3$%z#S<s1FL}^&<WB%Chq&b1S<bU
z$Oqn@8LtOeGRIw^9Gthp#{x3f!5?o&G(}&GCa8)-?D&>8RyKUV|KSKN$8F(*DMenP
z9H&VSA->mJT>6U*yZznN91AF{Q62%l;^1Zq&TpY*O-}A7De*lVYA6Je?RG>)>O;V_
z>Tnd_0DPKe{<~cJ_v3&=wlUaOA!f+di5TEs1<QGBC<P9&V~P_qlU{5GH#4wQiQ@i-
zXB}CvfMM5{BhmS8L<=Y471=WakO0-Fe6tY{p^W0c0-wfco2V2=<bkd(Ih6iCR^Pg8
zh>4`$DHNvsas}x^I+ecY)0i7fo`SixVKe?lXAS1S=1e^3LN@ttgziKPN2LL<*OMR+
z_Y!or>p8>@wS>gw=8*qk?i;vUu)G4VHHHGnxaT&ZJHdO^JP5VjQ;y7TMrhZ{?z>>0
z9DD_PQ_SClJ*(XJ8cb0>VsRh)cV+L29Cf<Iw1kd6ve9El_MKOm8pl?zV|6G$SUg^m
zB;uOUgl=G3;*SQMWak-{BX~jN8UAsbZV^aw(`OY(Q%`Hh(7|-b>QDdhu=<B88*jMJ
zAXB*C7y8<A!rF^O5o${3MvAW*=*>$`=o;Y-bb7@8Xx5!P8u{lAh-M&du6v%e+EER=
z+eNf%8G~MnGM1DY@z>2U$dq#o^Y`qWFcg;5h$V(@5iH2Pa!Z&LPR79+*<SxckWs|2
z7Y34Ky%YsyA&X<0BV-}~X8uSVW}zJ%14-)m7dXZ|AT`1qET#SWukJW30j3QaVNHTT
z>>Onf?8V)mY2r{8e}ur#XK#+YX#ajAP|LR){CGIj?<5j@&Q}WZgfCF$8%#ItB<>!7
z!*LEFiSjrc!3t`=Ir6V9ye0wkha3L}5|auB9zo5kew>^sWffN7$9$+LMrIU0?z>v!
z?w7v|=~X9Nd1D|g1rOD!VF+!mqJ;=>{uyF1yt5;s^Sl+A3S8yI>GQ=QLR8)1jn^oX
zs(slz&YT&%<p!`N-A*Ht0F&8UD=^IMY7ah1Q<t_ayMkjZ$W<_Jj5_t8P?R)*w-EQ`
zqm{74ULacdqL|MjYQ)AND(5eer}FULn7s&L+9>{@4o8B6cxeQz{by<pFp94Y7DC{_
zG@tGeGg#IAm|cGVGE9Z(b=+%BO-|x{A)nucpACjl_lAt|ql(n7|Fk_@oglSja$%>0
zC>18PiU@3Gu2w+Ix6~P~(0vbKG14?S=+SHZgX2gH=4|3``7yDU%8P>o6Bt(4j81|x
zCVx~_k5*wh)?8#99ct(~kYPP~MVAwFmBMXz3X>{rGiKY3pxz+k;Ekuuj|J{@<J56y
z(-9S!&QjmYNq$|S;{GY2V-Gj|Ly~aR6%Cu8(W*2ux(Y<!Fzswrv%BvE3y5nE<%`4d
z8_=@EpoGRxsV<cZ&FChk5*AGzT&5dx<T;n`yqE&7^us&vJ%pzi<&U|EZ_Q0>r;r|J
zJj}WzNK6y*X_2;EpJr`#R(QZg^d!{`r!uIM_pv-Frx;#BuYP0@VqpY<ZuugKOE7h{
zZ~;lmO>CeSeMb?5Bl3rs1?UWe4f6r-ER5DiPXY3h@#^tA?)gT6NCcJ66|}F2xS*zL
z%8?S-0IEf_C=wfn+wsOwTF<Y?z&i|kijfFvSP8U6gn{=ON}E4`RwChaqDXVyIB^UV
z&-}ZH8O01rX3pcf5lSCT+1Ils--(dAv%_$1=%NOBKOT!Yc15N!a|Yjel@UXcX$22*
zb-I0tV`SFl#knvv^PEF34KoQr=8R{+gDf|sw&u1r++5J$S@tHNP(}#@8|omE-`@X8
z2eOB8Aehdr6F?{FyWL0*kNlyl^BlARgGzwSd=fet6>e};P^4}0m=+KERkl`M*phY$
zzz21{j$kvQ1?*8k28J?SvT47me8lTF%t|QK<t6QQEoB-_W3S+sQ6WOGi;4^$jQ4M{
zft&q|@ERwqmaKknuziYx3bUbX$jCcrZppBD-8jT|9R{!my)n1P<?7u843-wYqYxpA
z3=)z|JQ^2wTLW1@A0Od}I_;NUTUzfDKyT<{xKc>NcOw#I!^8Dv)M?v`5DX>8U}hI}
z^F-rm|2o{12F%UnzPgm30_WZK2d|XQ`(Ei_X}nC<VrrV2An~RQ;&q@RZF1~_Zl7P!
z5`f;WN^LICx5vYUB(to(j=8qvCke22QJ!hO$Cf92LB%6+##wC_x?Bf$?g>ETYxJ#=
zO=u2?a8mn&@U@7*x8je#Jdde}S#L5YYcW(UO~hd_H4!SIYI1FD%-jmHCWqEaKLfUC
zIeiKAA=9X4{nLbZvfhr(8m8~>L)3}z*}<FO+Rx2ux(@{<G7?(%VjG_8VtHSpt6CXo
zh4|GZ>-ox(?((Yk^6=f`l`Y&qmp!S=f&%}y++VHom{@Bt=I>DKj`EGL?>~Dpi(A@r
zMwTx$RJUr}+~uLf!X=LP9-olr&J%I$qp~9Q)117$o|j7*Z+=Tr4Y7UjM{jc3QZM-F
zjjWo=hsIi7<#&up=S!Cgwo6N)=0maNb9Bd-Z&{C7)FI!M_G<(_?Uub}e=@nuN^o4&
ztmsOp`umCm2Xn(fd+uewfV%B*@TlI$aT}gK35B#{U-g3s4yQwkOn$NdnzJwq=mhr(
z=3A+j>W<0GR1bcwHAz4AP5&3cp(2P=`{VPaPPYgva#2NJZAE<ej?F_oTVoOC^|Wlq
zv%k_dgdPC)7eCrd%w*Sjin!~~a<LY#(%0uIbkgT!g0&k>uV}*CKR@={$U5GXn_!2(
zY;w5$&B0A;rt}_f<TbmKs*Wt`U!p<oK@Z9tv<-i)j84{JxlL9KY(<L&qFJ1^n@sbq
zESU=co!!!AZSqP?T@@+JXGSIPi7kFqc6;H`E7N8Y%CZl+;_=B$OBU{&g<f_?@mv+C
zzZ^;vh<%r^p}&35XJy$(zh~=Pk(|Eu{NgycOfa7#LoZyS^7@W$)x%SAvuXeJ=bC>p
zO*zdyCeeAZ)#}z-BR9{y)>W&S_al`UhM}isFWG#XC^18m>w&N-bh<;GPjS@bn2|C>
zVY6sJ1CL}~-E*hnl*(`TPnO=iVuOFNXIE*SdKz8TSze!W{9KDiP@CDc{-F1<gjZZG
z=`{^sHY6=0+Lb@uY!fhTAn22g)<lNca<rM6(4x2KQS7S@FYE-`oqB2>yo`~8et7{~
zsIA-ibLRubeiwm#A+}we6+1ie__8SXE7Oo5Q;B4RA!&9BkcOH&c#W&ubP@=z9~#tu
zsFpKgA20AJrlWU0mwdIgOR7RQ8E4H}l-Mv$rR9DL{D|q2tCnl-EWY5qaJZFjYN}BV
z++6~{`yfOZiR(7M<27R#OYT6$!nW<m=|eX7plHsN(?-kVy@mB~D{)5UI)IQdQr60@
zR0o6#ro4xaC;03Eqt`~f%`d1bWNQl|zT<3Y^_-wFd4qf2#h0qjzy1JLe0;ZDS;p_G
z;YViiGtHb1TyxL(E|<SSCq+i&jg*B6n6O5^@b5TJ+H;a1nOxQv-H-AUs2%|a)UHha
z1Z|(bFhCC$>)p|@y?6GUCV`?^80@n8%K8~zrA|`)Q6QRuE?UJyi-soR=)VB}>#3Tb
YRV~5#-{pTp!jIbYw2U+>)d`XR1F&CatpET3
literal 0
HcmV?d00001
--
GitLab