From ccb1229b61cba9340f29eb4690b2e1d00f2d6f45 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Nils=20Forss=C3=A9n?= Date: Wed, 12 Oct 2022 16:05:17 +0200 Subject: [PATCH] push --- LorenzAttractor.m | 63 +++++++++++++++ Project/N_20.fig | Bin 0 -> 22253 bytes Project/N_20.svg | 180 ++++++++++++++++++++++++++++++++++++++++++ Project/N_5.fig | Bin 0 -> 22276 bytes Project/N_5.svg | 180 ++++++++++++++++++++++++++++++++++++++++++ Project/lufact.m | 49 ++++++++++++ Project/script.m | 59 ++++++++++++++ Project/thomas.asv | 21 +++++ Project/thomas.m | 21 +++++ RungeKutta4.m | 44 +++++++++++ SimpsonsRule.m | 72 +++++++++++++++++ adaptativeSimpson.m | 67 ++++++++++++++++ forwardEuler.m | 37 +++++++++ lab3_k.m | 30 +++++++ plot_composite_quad.m | 24 ++++++ 15 files changed, 847 insertions(+) create mode 100644 LorenzAttractor.m create mode 100644 Project/N_20.fig create mode 100644 Project/N_20.svg create mode 100644 Project/N_5.fig create mode 100644 Project/N_5.svg create mode 100644 Project/lufact.m create mode 100644 Project/script.m create mode 100644 Project/thomas.asv create mode 100644 Project/thomas.m create mode 100644 RungeKutta4.m create mode 100644 SimpsonsRule.m create mode 100644 adaptativeSimpson.m create mode 100644 forwardEuler.m create mode 100644 lab3_k.m create mode 100644 plot_composite_quad.m diff --git a/LorenzAttractor.m b/LorenzAttractor.m new file mode 100644 index 0000000..bc7ee53 --- /dev/null +++ b/LorenzAttractor.m @@ -0,0 +1,63 @@ +function LorenzAttractor(T) +%% +% Nonlinear ODE example from nonlinear dynamics to solve the +% Lorenz system using a Runge-Kutta 4 time integrator +% +% d(y1)/dt = 10*(y2-y1) +% d(y2)/dt = 28*y1 - y2 - y1*y3 +% d(y3)/dt = y1*y2 - 8/3*By3 +% +% INPUT: +% +% T - final time +% +% OUTPUT: +% +% NONE - this makes a movie +% +%% +% adhoc tactic to set the time step + num_steps = 2500; + dt = T/num_steps; +%% +% set the initial condition and allocate memory for the solution vector + y = zeros(3,num_steps); + ini = 5*ones(3,1); +%% +% save the initial condition and enter the Runge-Kutta loop to solve the ODE + y(:,1) = ini; + for i = 2:num_steps + k1 = RHS(y(:,i-1)); + k2 = RHS(y(:,i-1) + 0.5*dt.*k1); + k3 = RHS(y(:,i-1) + 0.5*dt.*k2); + k4 = RHS(y(:,i-1) + dt.*k3); + y(:,i) = y(:,i-1) + (dt/6).*(k1 + 2.*k2 + 2.*k3 + k4); + end +%% +% visualize the solution progessively to make a movie + clf; + hold on; + xlim([-20 20]); + ylim([-30 30]); + zlim([0 50]); + grid on; + for n = 1:(length(y)/7) + plot3(y(1,1:7*n),y(2,1:7*n),y(3,1:7*n),'b'); + drawnow + view(-37.5-n, 24) + pause(0.02) + end +end + +%%% +% auxiliary function to get the right hand side +%%% + function dy = RHS(y) +%% +% evaluate the right hand side of the nonlinear ODE for the Lorenz +% attractor. The input y is a vector with three components. + dy = zeros(3,1); + dy(1) = + dy(2) = + dy(3) = + end \ No newline at end of file diff --git a/Project/N_20.fig b/Project/N_20.fig new file mode 100644 index 0000000000000000000000000000000000000000..5eb85ad2dab4ce92c5055c278bde8ff17ada78c3 GIT binary patch literal 22253 zcma%i0&lP~6?!-QAtyuEpIc?(Qx{7g%8N#r3|wd;fu(GdYtp z`H*~?NuGHkFQWEKM3jPynUz9bM2*S9+RmJjLdnk9&BD>eo{vIFOjAbTCnqC?n2Wiw zo4FZn!7oLqNy7N8yf}NPd-jgK6X|Lc2;(FivLf^rv38&1WHC8{Xff(=F1na z+!}2JFPs%3R17MK-(^8R)G?K+I7I+qG|<6HRX6T(6a=mBhy#fNQL0gDNWaYkh zX`Ei&kgx7bf@>rLrvh5Z-fhA4T8^p#z0?j@N!LM9K}kmKTXiuW^%z2`D-R-+!IZ}B zeT-Y*hQgGo042Zn1JBqead7uZ39y{UUR`o23g)-A1jNK(VEp_K%Mo71pE{0#-$_91(^*cGLuE!4ae%$SDy}NPqcWE|e@VCEe7HIX@-q)@;JT7g& z>9q#>-yC8SVp6@LzOwArDwi7CJIe`i?c9qlbdT~me7=|8;2mO<(1|_N*OZJFysSAv zZXm7?wmw(x<7n;mV$5a_vyZOt+wAj=SCMqF5#|(2wS_(xM}(Of>mn^?{oo?pp{ebe zablfsRhK}Eb7_KwR{=#1s)zu$8PiTe>G7@42S1K8G=Hk57kCb$pXp2%oDBKPhlL|? zt6f8HLQa|c8GQEQi~e}2cRtiFr5bqMFGqjLB*{@|E>5ym9aG^-0-wh3>gwTQw}G0T zIi>@p7q?<~x#|{z?ibi734+f?+wX?DNi2PBYazP`?`I7Xx~#2fAjRHFJdzLljkJGu zD5HXZpL@0?@Se_8;RWk;Fd9FLa_BnYKHjV4!y@Y8pjg33zr8>al=l9=Bq-_jpKGU=&UMcLuajPETqqCn z=`cF;OD7zM;jAdkmOec?fyq;3Ls#`5%T19?Pi4Qoub<^CNUG@&egehU6#0MZTSnAs z5zWWFF6PAlo&Og%w`+$FnSJ=$jmu1$988AlBjD+0$Jf0RAe_mGYF+8pNUkd4a8eM` z5xk3lGc6deCcx{l_c+n98gO=&SuFzyz*_8^;o)D%kBTIdYcsK>(#SAy$RbQe%sx&M zzYFBE727b4_8PHqNb&s!#d(nn8|SY)`vZ#q@b-~<&vWTsrGN!HDWgYgnicAvR-lXJ z^d6Q|56k~$iW13g&MV&$An4FT@8SXRGSeg1Y`SRX;ndGi-|CYmIk|^vr0`UX9Vv~B z^!{As#G^VDKOUiIc!;0jco?PV7|TIG=*A5_1xv-?jupyXr`wF#j5d%fi1653-#)1( zHSg9*y1q$@>2Ozja30#t-Sb47V%p|(DgvOgaB1;LA_H^!E4;T7ik1=!h_|p4s5U=F zToo-zMYzwr2;PNq8cw^e(1Qm@+%V)_7DtdD#CW8%Ovtf94u=E8ACdI=1D33~J6ATJ zPjLFKl5Vo&uIrV%p7B{cr8nXTXM+&3PKRxVIjy#IG>P*`;j>8dbQD3}-Zxh* zlb?z_t)}1#D7kA5=Yy9L^7{R!HyjQW3l>B57hE0lh+lmy@5Fjmg1hCfVcT|VJsmfN z9@mjX2sagWmp-t(@1U4vKmvTU6k6?FR9(wfagM*((r)}fP%XaV@)iL~NbtRLg&2Iy znB`*45Ec&^_Kj#J2`^@!KO`mUj(WU?z9K))^5DS@Z6`TSYIHm+H-hzYDU?h;m;lgR z$O=b%Z1+}mv9GpFHB>f`f|rWs;^++TGr{@VALkUz4B4MlB>XRsSuq*ZBq0E>1O}-C zn+QB-L-s`YW9wjlfG}>JXaj4sGvnb-o^OBS&y|pU*FifsMQhcvnF_v>eIG`-rW`vm z@qi#8^Vd$8c7aTBVJH(UL}8aj2%!V z5i{fOSNrz_M&B@uZW;cG{>3nAd-b9)MXUpj4P}Xx9``$14LSR=Q{^EM7x+V@M#((& zP}kY3G*UX8efaaEBvXbE6dqzZHxe7A>eT(Gaoi|;t=yz;RUs^Y60*40=(Sz=2*Xu_ zo=#YK<2faA-B&?EiWD{Z`|TxKKGChwxHA50C|nm>cju3>2mxAis`Z7=qqP2Ll=C=l zkNn6`=k0+2H({Rv2ra_6y?{FjPRU(CbqNhr1{&r$`9rORu#W9EMd-U}(209S@B$3c z^Kd?!Q;F_I9jz>Cf+LLBk9obhocxq8M?Q0en()Au;OZgTz9mp@FvvW!oN$^J-a!Es ztBeyv^ozN0P-R0$a1la4ZyK)o-n%Z={MVDV6U!|${3}dj!laGPQ%o3Ig^4gO3~{st zOy6U=eU-c$QG{^E;#ub^*IY13BA;XrJ*J8!nfrU+fxji-jDmm-oK)lS6q4}54!|iA zqqf&u;5GP)LBh#?onJ0)1C#2oJ;J9-eN7yP-$=kA7eI%%CP3{n>KriQEI2Ht<2-?X z6Q{}S7pQ7BVCs)$jAyrSWJUF--Qhug(2aNIn@l0=X+hi?3ir=-4*+4w4`oFl07k<_O;eq%7+Ss1mZKg?H-Xz!pO%lJ<#43quhAeOZ;-)xbJ ztS{5uoQ=``4{gj%wcFUH$>_Vu_{K&|tcXROYDn6}_x|)kY&c0IX|x%Z{ca}Xd}d=( zRp?|WRJjEiw54hsoE=rUjlNO>1t^sm+UN$kjlj}EYwm_7=7!>rZi8KcPZkDXFfOrR zS_Uo-T4ZfptcDVSemj-r3%EoVpa)+(WjWaZd4&|0#(9@~rP7AON^^cociu@Y5?-9h zV*mBch#u&t!DE0sN(wh5{e`Zup`fav2RG@fuhf9wa~3KxuzZlGoHe2IX@C$Pc~(*; zL}_2Fzcj?ErexB)wwQ@f7H;5$6Zf4t0k858yJ1BfF+CVB%~{_3pNckZRiRUJWk+I^ z=tw>lp%&3tBD}<8)GQHL9rbV)lL$?kkXzd{uau44<~%U|S!M`Fh?gVYy)+`@h_%34 z%h@YkDTZhbn`l#p;6($bbcGppv5VzmArpa<(>3FF6XHtIelEf_w!m;>Inri+vFpNAf6H2i1J(f%v>ViN7p zvOW-?DqM&PmuMGqECyRVpht380g--govAl+84 zbHq883%`+u3I|b!gVEqb^^jsTgr^M33VS%gdg0E$MCQe_FpjBm3hU< zq=h<5b|7r{MF}f2L!2Qp*em=#aj>KAfHq6gUk!tFMu{*|{Y}(7EauXcTBU-*vvtLb zjOfqiwU6#l2FITwt|Sr_&V)~wghv;ptmcMX@A}+(H(p2ys(9CseGQ5|iPX~Kr&?(g)hV=NF$8Jm)%m}@^V(8i0>uy+ z1pkjeXw_j%VxVqziUq5(jDcqpCjn#gPi<`}?JQ0Gj+P{>@~?r?S*Nv5Zwz|0)n3WS zL}_4FO+n4WQz@+31nWZfQe$3-8ZS_c7c58&rYF`(a%#VM2tmD^L=Jvpe&u!ELSx^B z;zvy2;;-noJ~q;G6P_*8TTfhkaz%iIz;xM;Ly=CS@JS+SIL^yi8m!(_{>UfGEr;|* zpwan^kt}N0gE2&?NQy$u%=KaPr6%)WLhlm>UWWP^%$#6aP^r!@^fCr~CWZh0S;~W7 zlYy6}$z$i))${SBz}RujBf4YMS(16~JXy;1H7L0dSJnbPEp8rR)~4}+m#u-9u>pt^g-Sa3 zo}@d`SV{=%&7-dR;#`OySd#asdVWnqX%)GJ5V9kB(fw&`*<)H?EP$yH|>rI8>#(8KTgt+b3VB_PgrkTe{nN#p$XpeTb#oe0| zUbSmC`q!EHF*KQZb4*MEwY1CALde=P*soGuwv{H`-H{*;NIyqqj9B*H=%D=D8Zy^Z z;~Eb;Oy=7}E=-bY0NTUe`e6+07hAZ_|8NhG=&ew0uz2elRL1!3w?a2mMwd z&vCTgdg#Xo*7U=!gqkc3g5$6WjBdRD#b}X+z>ikDSAnylID{w5hO(w-u0J^kX3d+JvZ1YyeCDNtcn>*46DSp!(97QSDMxbrJls$ewqdqRB=JA|J z(PE)tC#Pjsq;(R=?uSTM==G_N4MRTc^-A;PzQ0kSk;`)Pm@quf<@(`I`C$!x#O?cM z5{zc3l9ZQw70_lsW))}kDc%N*Ij`=H5#T4@_CkzTQ6jwppT1nxUzR{(xEunna;xBk zrgn>sJaa=*mNP^&NTrg#Ys8^K)1m!A)%i4^n41@%wKKoFpo$Vrf-zOl2i;B`6UebnBG-tR&HwKhxchW!;g|xW5ChvWC81&wPSe z0=;c#WRGXNUuOl(ul#49p1tdCOhutBU*3Td;16SwfxtFFP9vM^n%mW@UTaa_SH&Ne zpGZgNLPfdnukh1oXy{L4@h@=VuyFwb6s%AAi7yc|$S9vF=r3>x5#b2pVPWA9s3`G> zh)AetI!e+=ZnDhOwQpOUqxxuqn-WM#E{wBQO9!o`3B(!(}B6fS9 z;gbgM{q_q|9QtVj%%9xi**;^BS2}&B{Cro}*M)o6U5}5o=*XzZ`Gz~4p5+V1`J7dziOZDKU}Z>s1G9wCO_>idtl!_`!u_a&br zUD5fk@E@_S^xMvcTsKu${%2eOMNK*q`z|dqf|cI29uKvTly~~wMnk3>W+Maf^}xx5 z+qawDwjoHb0n6MOoYxNbtI#tLy6TcJU#TSL`2AJ${K*qiB|6e~{Autx3aMJ?X=h(y+ zP*tqUTqrbqE0JQf;M~C%eE-%?C^C@}IlOBrVhtLjG7OjusWh~7=j_yB4bUu4B^zQS z_*kMHxI)lrsBt;(d~nopGA{Q^g3FqWv#j^C7>Kqct8zZ!M?pSZLesjI+THy(D~&=U zu=1MGZ-$ca-g~QJ&<0nMzwt28Yk~rOAMJ0T#uUYKcO6f?9fEsMhm`Rj48bY$&t_nF zYRpA*yCrD)PHxO|0B$tA_r7&fo^;2#iPw&@bfsSqW9C;{t#nbtu6ysBG-QV436pP#% zD(;^^x2`$em0hcdW^um%UqMa-TnB<6rR{lOG8#Vg3U4Nwx zPB^*=H)oO*5#y7ruuMT;QDhJ`m9}EkWiu_5rp_-Ne0Sdb$3&?!bn7gj9u444JOp(& zyP@vfq~wkjdwlwX=FX4VwD}zg+n_Zo=TJWpCMyp+QaSXGg!B6M#6BUkKv+X2lqWv4 z=Z&ebaOVoH-mQTr?gJ%DZ?l87ci?(Kqn#t8pnbQ;PHse_bQfQz+ZaAuMIT$q<{b@N zVLf>Qe7LU*Co_QEY53~CM7g_RXD1`3gtdC}-hK&9oU#SDCFRqD9Zr`M5rZvNWPkF2 zOj|APOfoA9YcDaF*_^a`rS^?C@Ec`MjxrP0;$}^Pr4ze}*!v_xxi@eE=kU~=mc zJ#g^F*{_|wcD~Mkd~7Ewetp5TuAlh&iMB!%a4$Rc$3F8|#w#gDpLp#Z<=vrUv*)`x^_iWUV=Hl*t!g#j-K)H7pnEJ0 zNc{a>X z-hSaUTenz*wSf+PWkRI!Nwm}5hy+xE#1cN0C;#;RYA9YdCZ?foJOT?KZ|;*nd4L7{ zRMAOB-2Vr95zL2e^!6#@Ldx0wcPqbSI?+);k`uDyR;Z@rCAT=QoxUMuF^zlDgxLF% zx}oF7N6LkyG`s(s@T0GeG;g3I1t|o93Ht(3sR?7@yjd5VVP?w4*W=urVFE9N7Gta^;2+z4HIvB>lTU52GmvJ3;ZR|pETMwhCgna8-{1yiwbo$YQ}ls%4NaP$Qn<6^mtY? z5@(fm8GHBU#IFH7azkZG+P{f})}TT@OVL`WC|oG^FZNbCHw7l@Kk7W6soMX}N zmS&V4^BoJN+CaQ60!KZ&$Q!$)ZOV2c>nsQywah)4#$}xX{d`R2xE<1_MEN ztE0W}V`D)-T$E4YBOa@WG5lxy?;jp|6f+2Oci~@Di=THK{Y9WtuYB{n*$K)~^~eDG zJ$2gf2C3HhCaeVn9c-S>ov*qzJdt6h zLuuFQl5;q{)G%tnG}BS=70(GUh*ohJoD|-1o#vMFf&R`uSD3D=MJMq=Dy|?~yN`z@ zXd#%*XcIPLF7C_8r$r|fwAR4n*=LiWy~4U+bRW4tKZTE0MRGZ0p!){AFynu79HkUH zw^oVI64ccb)J1AZvoL?sr{7FQUOC3ndlgqBSq$LWD`rxZxR>|>%A`3G$3hqXSVHeK zI8^ok)u1Dki)apL|LV(90ll%>xYNGo!%oOnuhO{omKwFWh*3x(gEZ;Iu~4dH#5kzl zIQ3>zhE3qpj(&28U%PsIT3B>98GS{9M%mHX@i_5ZS<+rBHJ!0!-(TW&03~_)PsF8Q zf`mEKebAV`leUfb+Y}`Xrkn5e!>G`i2U)e4(un_V*J|R!_O3P)H7mBfKfkDdoIdcG z5iPZ;CW1TDuBmR?=-bF-8G82p32nZaC^UB?nCOt&&zrQWn4LFtp9w?hgO{mu9Fow( zj&Hx}*tNQoKLyAgd;4$YqZcpi`OP-_sZkF85I$E`iNa;w#S5Cw(HeKx9DmY<&G-dw zkD*Xh73e~gvy`T330bTy!2P~CT0ov@I`V;eE36TakOtYF>@3Wj)1T%mCy;kA$gfH;HAdD>GX(+35>&H%DTWH=~xl&iOhPfi?P$B@NE~{+)tblVv0u56ehyAsed% zhdL53&1E~CT`{96fevmWh|k=Lsr^(BUiM1HoMe~%YZ)xv5Qq4>IY}RT(njllO{ZN& zYXc;atpmgO#wd5!v^cFvL>UF!T{rSvDxs-)3iMR=d2c)ual8Jlxen+Q@gdb5^U!1?QI za=LX516May`?KiKnJGuO%5q1+{Jh@aJUUTD&4W$Sg1>j)D4?)N#5}8#2pBc$qYaaG z#T*abZ2$ZvGeyj~`b-4rg{MT5*yFmR{IX({;^Gq%aYC^{XqEJXueM;cQ%)?8v>Z0Y zDl?#;?QwpH=+-=2W3^*eO=$zaa=M)!4;!;0cZ~i_?lG{Y7z-*5SuC5hxam?=JF#kd zm%g`#iM%3YvaxCMJ$}%yl9i0E@D5Prv|A(`#UQePe6r4f+mhser+p{z1)h`vsl zxTtt22a8TksUzFQE)u@-=iBbe^*!gQ6)JTkPCR6H;Q8`PE3X{Mv@hD~656p8ZIlNB zzcxM}Z(?qZP%e+N8Zs9?m=*t#lo_M;!rd95AII&$nVF*xGI8W4`)+(?(Q|Gg=oWft z@;&brCC|MjfsH!TZ5}vx(jc?W?pDVDgB5S_M6AaG|A_rFp&;HpHQswAQckYI>Ad5+ z`7vTR#Bg6wdxT?ozT#ix^VzdUUI5kp`5$;2)ia*Va*{Z7v{iqFRdoe`%;5F}5RR6f zp&IVs^$j#&k3a4osZYVPs#o|>7M$a}fNC8t_LMmJB*zU2A%3NwK;8ZifrvoVsNdH9 z9|G~$DxXZJ0-`4tiPQ5f@vLUeJ&GDLcwVHbDI(ZI?C7xHsK>P7gY}x8_urRzk+41i z5WAKlp_4Sl|-+H;H!1H&vV!$kX=qmz@xGe|j@y3?VLjUor`H+oMFA)#( zke@oj#1%>2jWoi*wI$z9A2qh@n9fFB+&r5>x-u~-nZ5JhJ(L*so7{Fy z-*!#lKf6d`h*_7x(v0g@e4nH}PiMo6gQvd)>##2~1AJ@2&f;Iu+-M7(_he*z z+c__nxGtNy7=Pbk8~`H~sKdCDHY({sm6QIZldnksGFg%WnOzFncd5VN(yGBhcTjW@%E`(KDUxD%E*cY1H4H4&p%9HY{z# za`Zw6C$OqB2%;PKfvksQWv5gNOU*EcpZnfenH+Kp^TO>$W+ z%)h}ugKzlAENvzcvV}R5@x~y(wemSeT|UMdJ`$kfWw49?10TL+gGd*4+!fzUI|59- z5->s~-3ql|Io1k(u=$T?D}vdUH{&zd{gk* zy8zzVC8qZYG`fLI*RHx6HpT}ku8;Wd3FdE+9^RMy^UOb-`_=>5-g;gyAgbiK%DvI1 zFBI=($$AE?yuOChUgG6zmNhf;jd`UEyyZ>2S~!!LLq0&t8L%d3KwajjSy9N7_j6%O zu4exn(;=@0CfI?24@*h+#`8b6!RO~i2D%cCIzfL145iWnEUE+&iS0a(B!)unvZC;L zgWC0OQ0U=Eos-9T-qb7L=>O+9YbQWMy?uKs!)6o2m&;iKEr0Y#O8Bxc}X z9b{;6iTLaljhT76%{^FG;CXmI)dOQtCUWp6(sM7R>2aid#}+n>j^?^FKuUEbrHVyq zW#|_Abfj%tD~p)nO0kt%PR%7-*Kx?(W8l7Rzq>L!#AyV!;G$R#JI;7gzr!n-nW-U( z?s5IOwb2KiB4)wX@}8{8hqQ)_eIm{M5UD0N2WG4`IE9}hp*A4!%rGn^w~G{%ARZ=l zeMsf)V(R6h4EXO7&W)nwPN&5mN+BMz?GvY(y~qFsGIL$yUbG=|QLG-n1U^RKouwoq z)Yh&IrZ3eNS?ju;d#^Nq$fGk@)#|}0$)?gk#cf)L2o(I0`#AD0kQ=+=*qnF=@d(Al z9!^Ny9Et<}5(IcSU-}6M=0}z5ivupG4y^Gp6}xyqLisHPz8z$j+8xw~rvd-cOWwU$ z)tNrJU?s~I$htCTGHi?~F3R4)dbeZhidl*HCV>Vg&m_V%BUs#L`?Io93PRSf=rqAS za2s3~2szBVVL9Ibm8p-&(|?9g(-aLmh2;7i9}CXkK7(p%C1a?AXR|Hs1Q9Gl`!p@= zk8aP2e(nk&+sewKU?r^dFAKUw1cc0>%NfxM2KkhRi^*q{_bhz^hM6enBX(Im9za$Z zb?5fPjqlueeEq!Xa0l4Q;D9xy#Q zcA?Nb%C@p(qp29MP%~tl84*+$zqvi&hw}6(c5CzMhx}GL&M!W_H{<8IZ6SSIf8bV= z&tg=T!RShU>G)y5oWjtBA)iu`sm`gL9ly`3!YinsJX1GkFkCTbqEQf4vhIS`g6S+U zN{_>C9U(4JBmH9_6!Zi5pRhUf;-JYVj@d=^-j(zxR`sH*(hO$+pFuWgP zT=om*6FZFim(a)$bI;wwK4#Qu!T{<&FiRo4t-LI5DQJWBAo!#XGC_&Bqnt=WqtS-- z<0esG1f# zi_Kx~xZJ>G>n_Cu6Dx}oU(D25ug-*I=wD;{vsIEMI0G_w2Vr8w9$FChE7HrFahp8j z-SN#keUxaE7sBTFG4_2%d7S!8o4&PXNv;#?iYMRq!}DOiBA;h@#OE{7SvmsgeYAdF zG_=9*GnK`i0+Q29W1U83QF3wlEDDpFdZ3MvNU-;%x7J7d3zU5p=9W7Fl(p?-+pi5V z;pG^(!AW*#H5AU)PM@>CKju$P6$ZoQjYVAFrLA)Gi*$*^L!^eB3sTg&vn~QX{nvkB z-Jbc;UZH%WgsbiKmyw1}KI;Qy#|Y8j12H`3_AXKX2^7+zY5xu6M?=yhsi#E~*yXuv zxf=6_0Y5?-fS>S3_3{J%;!-qYGC_xDKIZ;CEn*&*+(An9! z8>wm8WLQXT;tVVdWNefTOcD%qEPQk<6spLw;V2*iNz)(+S0Uo$SHeO+Bt!@^cJ?QA zFu@l}H1dT1WhKY*f0=PwrzyHpm<9L#!&2w#&QU1@Tf$$(d2@$gv5HrGF61@3Ze`t zJtx%o#S;=!<+)flot*<9ZVn)yMe?)vI=0SnM<^ooB0y;<_3xKE>+Db1oPv7G=M}N{ zfWnV1bSCo0EV1{OYUwtjfI!$g93HmjnZ;)QFdmmKQqMW78h#NEFY0k<SYJgxp8Xg>A*6aJc={;BHUau>}Z`NWBRhgw+b}A#R?wyh-y?ObI$KeJm+z`&0&CN(?S+@smCtOL zN@_YN>|JONMKzF2iMOJ1(L0SbVu2VO+rO2$QHw3${RHOs;DC058{8rfE5msviyxBMecfx}wwu>=eva1GP2E+6{1I`d@uhvYync>?%q&`ESTJg5l0JIm^Gs=LQ?>5b)Kym&Mk?kn3V( zoeG_8OJ`Gmz0LwB*zQQP6$rDZxOA&VCOWz)#4MG6(t<51sn{6Rmx+Q;*sZche#CI3 zsIxH_xn(;cPUp~Mki@h1beMUwlA3pHCq;20g-ZMqW0O~BNS@i3;9ofuFH7S58$<_l zpDkOw%y-{c8oGft7`JX19y7sLrjztBT(Yd} zit%A4s#*YxaEkX6=gD=3?}VJTu&v{@bOvQF37^oyBR#@;nnyv{3O-1({5 zOiq4Z!0H*+L)5$(SqzaS=dZJwo1pw1?$m*R2QSDh)!ei$WpJ-G7SgNbY}QH77u0dX z0Ch~aW&&B*otGJ5mWR7ZS+G_!=_r&N(L2jfsn|g5XboGHTXvdTZrKf7rMyChbnHB~ zyyu{K`pKd~);2!JkO;jP5825`4L>H)m40za_B!0eCgWPJYqqPN&iOJf^?%Q)Ca~^j zxaJI=ymT-rj-A?9PcR0{C=aeUlP;yeZ~k<3&-Pl;&XTtei^bhBO^IE&IXmI6p52W< z(zLJgIVy9?Dsw)<@Vh?}!nLoMPTKmH+5SmQvfGZy0~ozlOne7vM!efpwwCq)6RPm; zKK&YZN>_WoWzckjrn_xshRTbt3y6a~I6l4$vsk>r2tRUsqzR*QkO>8h5gI~en*jor z71(-_mam_JZY+awfMSKNWPN)ZhIga82k}A!UU49S{Et_GcXTfxx%t~#aK8B)cX023 z*EFbI;G-+vGt>$P-xiZ9`JAH{6UE@T`&i?9tZ9 z`?jHA`{YqaThQipOKOJ}=&?w&1$ai~Y+FgWD(3QgcF-}6q#bi@vkdwm(;3Rx9bVm= zgN(0wCJ5LBc%M1^b^!h%IUI=K^LPA{-w$r$ z6Y35F0_3`TRV-Zk6Z;Ghdt#k$mj)RM?^3=w7GAOK_A1_fBIZ1<*iEmn)hnk$K5@9& zpJ{F#MGjlLA`V=i{|z4M%0#&8$=6^%f0zRJ33eu~QNLbM@y|Zm#+zPnjW#-D@vMFp z960$OM>T*}j_ieu^{&3?eo{<0U*Kbmwi-$k*!M-TzBuz6-vKfH)6%eoM1Dk5R9r-0 zb|OYEidD$9zc4Lg-`5JsC+9-6eeXiKFkRvL@BiR^NR6i-yAb>EPh{ETPV@*mLq-f&A^4xrK_1&;!x9>|&_O+qJj_C)@Em^4l(bu6TI}6^S9Upw-^Au0) zQ9*C?98Z=BN7Ii>?1wjsQP)q3B~RzCjXV#&&PFpgo9_4V{u?XI$LlL#gJbT8<5$Gp zNf&K1!LPAzv+7iF_di7mz!yJ^Yuc3uwH(I~O?yab<2=vV2E9N5sjRE}rgr)n1^ik8 z$PD9ic?JQivE2_>79Ll|9_aP`mpw4Imv-B8b^;C1x>5k$kTMy~e^J%OA=|h1dnvPpO(` z0rIqH6jd1&S8xwc581u5Exh?gZaINRl^)>$0 z;bPbgtVd@%eq$aW;?N^_c&N-*YYhL`svsK&4WllqhY- z?;|}X)MPtmQt|DVKfQyVXZgYK<8 zB~fkF$sZ}LHHP2s`9*~u+%bDy5Ayw~v@yPx#gfy}n>%Sy*X-Nol(!9>b8qtq0^2b~ zZGA1n8eD5v;jM_uKsQ7$bcq6QRgfINc+gAXe2^bO(C9t%b}OT|_&#=C~NhCOYWmO>{{fIsZ+*NtL@2fR)aRh!f2wt0m?aqnc zGxzxtXwkhs>TdX9Fj3zlevCh4ytfuW1p$9yo2{pyf-qzffHV<*KjX@p1?17N_YFo* z*T;$e$60)`6!Nk5HKPAPktYN$iBinVns>l zYjsh87K@iE6FC2gUuSA6?B^VyMyi-SLVMM>PEcf{Vw_AZzcpV ze`i^Vo^U00rt48?_TqnrXin5~{1*>4a5rrHO_q|qPvyUCIH@3Va&Ztj3iiImiK#bs zr!)57ylEn43i_1n1_1x;dC`!cMF2{+V#ms`N!JptAceAk62?(i8o? zMu%~O^iP9sd$NA4Y{zi8P_3I;*Sq}< z_dop&MUHPgOt=p<&yW^g1qi7=N#4e@LVEqfQg-a3;%PW1y%Lw^jJoTJ-v2hp=iGDv zcf~k3keu$*vvC)7b(qzljpGE>!$z&i6f!tle?lKm+_BYIu4JhO5_;?P(e;f;Yi}(N zbf{o|stk0v$G=hUKN!Cr-t3Vzm|nbUUt8yh%LwRm?&I6&KGtg$IxrH>MX&Gu8T+TG z_C12SUBp_N_IagNqtv*u3{_udXkJnX{jY7A#AOf>if#KMmDan&T$+E26KTU~ za`v2N6Nv!-pTN@VoeUc0kGh!uz?_(&Ztz8+m;y0E$LPkEc=Xo%yc z?RhxYx_zJN8ImZO>)mrb))UZF`Okpv^zTN?W5En!C~&`)qk*>?i zM}@E`Y6_$iku$GF=wm34-pD9Rm*c|YH}B|+d7Mc`4{QBJ{*wK*k6ceLWo<{Tv^LTW%_{_4nUpHuGXDnZ#qW{? z(&abl#{XkPk5NJeJNFD?^=bn3>rdpmUaP$UK;9iMJBLP(;AVTkVgRAP>1P|y`jffV zj13`z_3NpY^l3g$Y|Ln=B%nE3V*emQ3IT7_C}Aa~Ka19swLZtUW0C_t<^5AaAVFlz z+qXK8BXQrm12LQPYM6y3qzJtedj`x+#Mnuwv0!Ay(WT$dyB7b`!nCAc%h~wUd)kKA za{u5RG`~NL#wP8Lw_wii+dda1Zh2}1OVSp<;eIksD1CH5kU+Chxg4rmI8pQ)G1P|{ z#sgY^)Xj!?kpd@=37tEK^&Ml;q5)ax%YRf&!&VvO;a~m(V#ungEGF44WcIyl*wb6j z+RR4hC-WxkV_d0oDF)r-u54brtZDdwSk4bi?t0+;U)*L$o@Oo6PU!KeLW*r!1jU%g zL@KkXsJEWt!36*J{V|)c5<=ZeoM^j9iPUV6vLUv&MZCmaY{y@L)=y@4iLYv}by6!^ zWGI!WPh&;tIE7wzrpBayj_)f5gz+u3LoF#Etx9Qw8tMCds5(bUKh1V_#5qYPlJuZa z9bN-wl7-{TLHN{^dyyBuf^pIjr8w(2cfUzw!&7&PIq!@(b?Zm`Huckmk7=t{!OEu< zv|OoERsi3V(AE+Tuid}Ul>(`HBD9r7W`nIO>AP$1so|w3AL@sD)fKduOvL8j=H7<7 zIEfh=^3@(rnQc!9P-DD(!CNq-=Y00obX+Yw^ndBZ@-D!; zu^N9m^2nkRG&Rqsy1j~yD#r#bD#bV@6Qv?<{<0pH_{#IX=RGQlSzXW^K_jY3v<)7j zbStdHYI+jIRpnP!Q&sHJ5!snAds$RI^Te;y|yYX2zURPVCbbsTg z@8dpSTPs0RAznD^@3q0?<#@@%jyst~{nwanVNBZ@xu|vO3C)(SS zf%Ans+k_?fBfuW=U1rcaJyiVjSTmLr3xmIQe~|bxUA=X4fG5>U_jgM&?_5OkzBk;3 zZz-&cwHud;3Ez$PZ#}O3RA~swqrMT^wm)fTljX*#D!)mrCEZ)yw8RFl?XBCW$4@`m zTJ7I>kM|tYX{8)p?6ot((N4_>-fC+gSJe}mxt*oPc$*V*8+Oy+Y8QN*-wS8N7?>|U z&YY#eKQ!r?s}jXze6;MC*%IkItYdTybLE}oBzW^bG^Lf9vR6;OH4+l=RW7bymb1_I zdaFywwN1ODRyE(kw!FS%bY z2&bfKB!236Yas8v>kKZn)c>fA?+(V=FYwIVq9NYDPki?i*Mz$+P&SXihk#ur_WJH&8Z!@ z+4mD6@=1i}7%Iru^3n}8s8wwvJUd?ObX~|Y=babAM-A*1l4|~jA39cho{pCvyMViy zyIy&$;>Uz0@%8dW#*rzR=vKk%CbQkN1g}T?uAbRypqazE6~ml4C9-vf>4>F$PK+2T z4luk%=G-$Y<$Vn_-x|7XK0)}l{yv2CHOu>YTl3r`nx3y^-?;OHV(+4;iz7~q&V@FrZm7GE!a8xNI#E9R)_}@(pM}ocQr!T z;ktT(@tX5eh!xQ97d{gK*tlT|yl)@U-2bPNs|t$(V78Qm0staM38 zcP&eI3evGONQZPTT}vpnbS~0MF0m}X|Gh8w^*)`M^E3}L=V9i2^TbHL>{^GZ&r45p zx3`4zeKrR|(vDEr*cd`Ytn((^hfQu;Qc+<=7c>GGsBk;&R25FsicJp;joXqfw<;HR zWTdWA%A;2tjmm&UYo*UgjG(TvtG^yqO_;LFP0uLfeu_{o$JyhUBXN~ukw|9R@Kp~9 zT96^|o^QkZu*M?8<)0b`@2t9b{hOaaP|;nXZm&@!|24ir>!)%~8rA|E{kf;ub_%A$nIU)vU_!T zZNEpmSoIYV1KqM*HdhSxP&C$DvNj+3Cd}l&a`MZ?H5FMOlx1EHg)Cw2>pCHu^VuvF zt1+If;ft-!dxOK_1yACDLb&8DvLm{RoF>}{!=c72y{HAeF9mxIB?i&@9IZ5J zw1s4*pQFyEr8OS@$i*?sYq&RQrc*rIqCyE!74_rkE5hFH>(ylhrH#UPvnD$7Inp`*M8DjOX|$ z;#^7M=vSWHqgkKXe-jR{bHMEFd_tCva*VP0(R}#wPSv?ED-0)=>Kf*@#NWw(B2R< zF0awYW`Z61NM~P)&QxZElKfh=8genmfx6s;#J7C0=5z1W=1*DYPm#=Wz0o^QQcj&9 z`PJjd!jNUllQHr<3;&d@FPH$!WKMPCNQ_-jf;7aM6>L9X%TE>empN=GYT3NeKFS#cp~ald~1~ zG-k-|NOB&&?#48WI1Xi#0;~DCVLP9l4}C3oZ>k67`aO5J^ITxG)=GM0E{QNo z;6x=dX->?J z$|P@c5R3|FF~YJ+sZB9id|wfgr`BXtPh^jgjt|dYp(D&CTcQ>FG&Q)`bAZ~yfRqf^ z1Pm~)j-SHTTW5_(c8`}hiKTu+`B>?XbM@0MhMh=`InIO|5{;|4H-Smg4@$i%-KdAN z9sMolu@;^LlonmsYOa6kCI5VPk5%@w(tLiUy#{{12%yWm+Zx+hicqc?xHdW60I58P zKN;@I3f;$FN>tEfAicg?2C53T`dRHn!Vasp&tv%{F zHK{#iv8;2JA&g+Dds7 zTCB8z3-9-AAAP=~83@si{aG^1UPf1xw1xPnOK!>)_FrTlqYCzFej6o}?Y;OJsJu97 z2hUPP!)f*vo)NKAt#rk#%!ol>6j#V-FK6n z;vlE+h0iMJwaLs?TIv1r67{^-mA~4KAWA>a+i&ys?izT0w(_O1Pd=wBy}k3BNppZ> zi9TrUTEcW&-(&%#-X8$Dl-Q8Jm+WZ2VHspU=xX;dHyTbK;Lg$~xK11M2*`MdJHH|z zbYXXkq79IEt)&E0_*9Q!k3}s2>3g{nm|a`jzGnK{e~&V6d!e{Y_VgKA?>D5h zBX9`i@|J9B|Meknq*AT0BOqPqGMdqe}L5+2yxTfrVY;Q z?7Wjk@i+bto;L=Xp*=|dz*FQ6=+yL|x#HiJ`is9RojjE(L-;w!`0OY&NXz}Fw5Fj> zN?PK1a#ndH6J}^wiXr>;4$p#I@V6?6q-sS;nnh(Gfa%>g)tRApzviXCJ>VI{4B^vM zcY0N8FKutMZu_(~PPp}Z`MMux@{J;nkAN51>%a>rlD=`h2k`RD+*L~4i~X`i?~_!{ zT}=f`{?P{?8abUGmMxn_Y2Mm2>)qJorg0ln4{Dngc=`)&bIui^2BFJJ?HelQ60q5A zaLtb{{3(3igkE>bP39a7EMAfp@3d+E{ujS!swLIR#)^D%9B~pVafnCliG*Ij@|=ws z$PSZV9FpC2x`1CDecV!pDf2X>7}e>p(Bx@TpS%_hPOU`{xH?%)zJ zcbL=c2bH|^?=bx;C#b`6aEkeH=v6MQB0i(Tsgt~}O3XJpIJtY{H%Y1x$b$WV+L?k?gb=N*~l5fgWMpx0kW1dRq zL?Wfz{bs8LDk!tEV!Z-gYbYgoMM2wYBG~gG!x|oc2(8jv=+q>ky=eh$TF2akRBmF< zUP#yTCZY!hCoUriGU_zazL57_)%Wpy6f4f}eN2I2;U2J!%@HWqXVf7WKshQ~dN(As zP}CYQmYnJkZwjQB})Hxy2*Q3i|H#JY)C- zw0Oe9LMOyE>drlp81p90KH9TzE577;TK9lBmQ>8f%c!p;X%Ik!+1u=|#@E73j^4$L za!`7(QU|V^M&>8F^ritP_-^ZUukY>+LDPaum!&GVvTS>aBo>mfV%wF-Vie4ipw)V^oL)K5cu!@ z8DRJZyP5Ni*r#g3S2Ux8wNaDM(Y#;id4J`V+n{IZBf zw;ZUPcuTvEan>s*Hn`Yi-kPv4aN9|Eg3*|}YS205y+;!$5>Nxbp97wb02x14J9I08 zI#9sVjQf~>)XJ@iGPHX=Z+7**%k|*3O`YjCTCr$$qUW>H?Z zU*?c|#62aW5xuLr8-X}8_&ZYM1c%|EWql9-tJDOLqUAN}6<9RV9|BOfDl-$mHrdm%=}%c2X8Hz~ zEcEJZX8ls@mkUPu+^@R_UwAM(vV)aloam_UlVlCO+9TWBY~9K95c?V|9my~B>$ zeZ+L$?|CJ5I__Ny+ve~R81qFvt$r_UVXW;%5jaeWR%{9!&TL-?yv^yq*tm&u5g^FX zP%jV&daxx4)SwChb`x&vFdN)_sh4jFrMVU=m(4?B#;w)+aaPMa#&PiOllquEW?6z{ z7ULo~BR05ajX&UQ;w6PpU*enCL~vSzu<+$q8_ZjyN}Z{2O>s^ggX!LhFJqY&oitJ? z?cILaB%GRxx&6Z5LptG%6SQX&YAG&7dx=0fxJ2y<#oWeQfGjiCso+8hxAB&s?TmFC zxKQw|ngu8&W1UqvcJ&eaf5lqO)<~DUJJU)#DcQM z8;5uTodT)ofj@qTo3^8;?qFO7B}Uck3Hgk zNtq5)!`$199z?Jec()|s&S3b_G} z72aULs9;S=<1)(KT(Fgjc{dT+S*#GAEc)j`B+3>!l`eAtneQf(U zw$#&bT}!Be&)j`QF!MH0iVgh1T&M?hIzzQ7jLvZx_VyWS@f|{{{OnOT@qbcldg9mDip;GzNr;`$34tUM=T~hy4LCGdz3`fu$`i0 zYHA9jb#4MOz7$`-zDKOg9xayD|DC9O(Ks`WFycBG5*xX^t@yhQx+#%A1!-|U=3028 z6Jdp8QdjhU+KIRP4Z4@sFXGBkd0jCRJvMY}fQDVEWFb1buYJxuJfgy5VH)Z@BL5(= zUrOvBJge_Js^RNUwiX&aD$Q6WZ#1;>ik)zZhH>+E04JZ058k5U``3kVa~25{GhO0x zU49S@|KHHKgLA?u8Kx%C@zz>s@Q5OFm8H?p!t2vSO^1RO8_InCaJlf>FF^fsjXto} zsMc8xh?tou`N_%H?d8aEEeSoe2DA(hV%986oX_u$(sm~*Sc2uHGTQ|pmk4pslka6b z@yo|ybLHX^l>rW`!@{+mjP|AX0NGi|H~v)5?|Zbh=|G+B7f*AsS-00T1Wu5w+YJBU zkVj~=_3aS1fK#+<_~4GQ>#0p3T-kM}FP_tLc>#=O|c2%gH!I2~))l}^pn zb9XVGD?860(|}-=o3}hwJJu+#ZkCcsl=9wjSDD-X%~yH=v~zqzbl#cFO>|hXmL5U* zW0eqCb>Kl=hm0gdm7+$13ZFFOw?{TC)pD^nMRP5hTOSLUdlUbFozby3U5IF2?$Rx; z?Lw|OYahlAfxAtnzQ_6u^yB%5ndb=c(rrU&vro`R+JbOi~F0kcRhII z9*r7r9Q((i{5GmWQzrJfw&vic0E@D%Ks<7>EhmOB`TAEUa2de=Cs}zPvdE(&%=!1$!%>PQArBx7yB+~~Biw{h23=4{KVa^{@_yZG~G&ndx;@A>8g_rgD ze=RZh1zB{!1etdeBeOO#CaoFYsd=oeB4VjCH!T#b|0jiSxWIyo6Y13$m7WVFOC|xI zk-S;0smIfnUa4aoSYsZ3YD2l?27zVqqg<~0IY@>5=T()3#cU(#a-g2Z)RU4ZUF^vk zikZ;H`1p8;Lz4gWzh;Dk0079!6b`wTj?R5z&GXfDzN-;*1dLLeDd($ZA9bI1)d&|c z2UMBcO@&Kc7k^ZB@>*F{1?oCEGNBxtX&-FG^Mb?zu@i5jm&j&(2W(>y+ij%Tlm|+Z+hWn0e4k|Ep3KwjVAbRtR9vZDfdKDY4DfhaB8M(P+p*(vh`tl2z;w+@X^HP*pK%W6`xP0WORE^1?Q1DK z^~S2y<_FGRI>qZU;w*$dmW4U)gx7mN(MZf4k4fnQNjzNC=ICG$O&$}-1f=C;p>D1F z?P(v~&b3H0DtaXl3s#FU@ad^TTMe%y#!>Ah<6TAlF#3p^#K)YseUg6B(f6yA_QW8u zHYfJ^KLq)_w-z-xE`0@( p+VNShR5c~Bi^_OarQ9CWRXfP$)39f{C(#f-cWDd^!_ZON{{Z+RZVUhb literal 0 HcmV?d00001 diff --git a/Project/N_20.svg b/Project/N_20.svg new file mode 100644 index 0000000..64b7199 --- /dev/null +++ b/Project/N_20.svg @@ -0,0 +1,180 @@ + + +11.11.21.31.41.51.61.71.81.920.50.520.540.560.580.60.620.640.660.680.7 diff --git a/Project/N_5.fig b/Project/N_5.fig new file mode 100644 index 0000000000000000000000000000000000000000..d622a0591e97e2b96c180b2238c1c6a27c94385d GIT binary patch literal 22276 zcma%?<8viU(C9a|ZEQBSZEkGa=80`P8{4*Rdt+~KVx5@xdGCGyfZJ75-8~0-MbzX(L`gXqnMvhE)EF$R?ab*(mF$e&EF4|zc}bPTG-VVx+388eT+EH#%*{w0 z9e7FA0Oq8Mrf#GxETk-)yzFee%xt8r%&e@W|DVD`{p0@$n2bEie-;7Nj~|{nm)dY{ zm@94eq~c%_lvv{CDv+446BLozGUBMP(4LW`;U&L<5<@AdNkbwd$wQ*4$>HH)X~`$} zBjj=uEBQi06NBV-)K|o?$crOCd=yERg2vx(>ic(fZfg7e>sK(<^#K>#Y?r-+5e+YIhOYR6aeOk+F$(3+G&o(*zT^%H_2N-TUrjsFIwSv)FNG)@9`P{q?Q zEU2bqAhM#P->1IC(%b1_Z;J_=a~v`ZN-7%$89pANZ}^c+*;!<<$L#y`H1aml==H?L zBj$b7kHy1jiD9jJ0}m73lv4^ie4bTm1Uk#qAQlm>*wAr4fp2te{wOAe#1N&Lr_Dlm zR-VJ*zs+OKM<79D^M0nY@_i@S_fcyy%k|g5M5^Pc-Zj;*lgC>H4S#8jAbnjXe^byUs*WfkTuviW+0fePhx61NM(JTlr>y;I*ocizk1ybH+gN-=ltNfMM9u5*p*Q1_kKRBl1^zA5~2}kw&2n4DT{CnCCV9;+jZuC4)i{V@s%p@glHU1cd zu>51t9yG&#y#L(J8WKN)VI3LPeu3+f%u=X411;NzsFdU>-l&O+)oOhVXP*g1XtKw#m_na30C9OTXy3g1-)P5wezzZEmRC|Cy8aO>$Z z+Kny-GV!tN+lN)_@h(JnH)X}~E^kNI?^$Qxe*LlQA)tduj𝔫rBT@rbi68F*sO~ zGqYHTILrWp7u$2AKziZnF&l$3JGG92@K+cvb;r=4)RQ}YR84v`eX4A9K)4Q|n2ndm zuBqYQBd3}IX9am7ya}dtuiiyM)o<0-BGi`!&M-<3PBZcqp@PX&fKYVEGt+ndcoKe%pdxav+zl|^Qy*Tvt(D4d^#?Oo;Cbccy_+cYh0Tm z+?Z2IrM>sB!La z2OmG+7GJ0|U191^Ab$0^(k=4##xFbd!YH&Of_+U2C#?SHvH>{#HpV$e; z#`pW%+sW{RdOmZ_3#R#?>YolmoW8GQJPLjGkfbu6_<1&v%3d0wh&zQ+x&JipV!zaM zBf?lXh#2ErIKH+=NkI98+i;vf{3Zru1%PzkACNNs$$KVd$JdDuh1iGKIQrtav^T!8 zB(5wru28hua+6SLdT&-a@?Xw$etPjbNAlvA<_IP3ab%Z|71urpH+9B zcvi#@JT=AB==^q16Y`k#ItP3(#2p?ClG-cM=Q4(XWPW~}-!0&3Kc-wcAHMd-wemZO z{v~Z->~Dx^JV@1jc$VVV-clkVMw5OV)3q(22c!bsVAel%S?mO%X>RT4r#zq0-=$hB z-Qm(#tn4>>9~s|6tV+0Z=?C3A0`(N3e-3*cStUx-4&zWnDY2320y6`f!DA#RMXMk) zyStJPO{Gw=B8jSIl1=!K>7Q{*_l55b`{88?^>ED`k|nH}GEClTQF*1hDg zg9&PN^dfQj=4U+*E%^9jFNoudVOZfm13-a8tNP=DjlcIccx^D8&QK@6tZ(q4OWHM1 z*trFvrHSphn41g__eWBo>)f5~0(zGAKR9pse`VeHem0{Q_&$4qri^kpot#hUZl0&u zZgb>7av`rH1Fdk&VRT61Xu^1?Mzs76Mz_0r^}spS+eHAJ6||e?2b+C(^m}9j;)IR< zcVbZSxcMM{h631PlIOh)no91SC`@U+1dxA)XC5#oMf8(WHuP_%lwbc|&t7%{FRWw@ zFGV#GcX4r$9E!v#wG07{+Z+xABupHV@46v*D;RBiok{)~`U}!njCgDo+eM*NE7W1Pu@%VgTtXVHjr?Q(LOb#O-$M;}wc(k3Dt8=kSSLIgMn zzZ$hf%sU6+(LRFh1GMn!?QrFT2g2KN0RCt1T|Ay4f;?JZVRA2o@mn&O2Dw6yPTMl5 zapt+=J~9zkE{$>(#T8+l>wpBc-<_oRSMKmd6I-)}?(ku<=At>GdtI*t#qTH1y(fhl zxup!E5B#g@ua|qiGsy=ZY$Rl?SL3-Sp7~>JsRK}ntBM&01CLri$h<9dH}k>leJ1hT z-aLoX7~|PyhQns&!)9(INbxyvob#{>Ing9pl68NRY0UWQBCTzZbiMG#$6qOb?%s}T=%$Tr8Gb6q3z>0iq_JJ6i|GJqnP56Q9iaZSr~J?jiciJG zBIRLk558y%_UZ}-$U#X%vsZ~rb;>3DR`OP?IC{}Aq5dWPE01)j?&w85c*eoXQr}1_ zgh#;N``Xcc41#afpGL@CiI=c}jd4@N+HE5~SjI8l_W=whk}==iXluTF6bl9v4pF29 zh0enrK%~cooGQ)Bb2cCjP6@NtJoH+{d`*oj3N8&+`Ux5KSsQ4!{tkaKwN|5!D}*)1 z6XyJP9A&ph$;svq1uWJeq#kYL7;lB%8trCcu^MWPqS_+xQd;>pEhucNqp-Vqajme}V&B>k~slb2(5(zh)$6>ew%* z$B4WE<>R#PkUn1lLiV*>LpHN5h{q&QN;nh_glG@z1!fS4FDriGJTNnz;!PqkTPH zPa6rC573;sWe@r#gC)z1Bh_)`?-B{ewGH}}g1Iy=|L6qufWa?gg*=7@&70#JV!SD~ zdGr&+f}rv9ZOGz4hxZkNe`w->Lsu^ZbA;&orOPAgbuNSonH+vZYpG1MRI(G3E~ps&+) z_|pp4GcW)@QxDDI^?0o8l*$Py=n~SL_0tc6>EGMTb395tbMQA+dVUAd)h;y}e7a+c z8&1o5Drn7k_hsCbM7|G-@=L^dk{zDot-?UnXexy4+q_sO608A3j2@r)^K6grw(l{= z1#q5geRM)8(;p{g7I1%2!3b0g&QxZ>ct5tJ0Sz=({WeHF%q@e|+;dgw;%jR$z@N3!dl2ap&0?5q0L{Phav z2EW5Re&wOm6(0`6v~jgmXgY=rc$JmZ$To( zw~sHcLd3w9e8Y7(sF0~O`^LI}bbVl3r~0EvsxAk~>f^Vjx$DDb&BsOPj3a6qfXo_#o%YWjedfIRWWyQB*>Q_~FmymeBrP?c_RZKOozO3oByYI`q z9%5r#;Jd7?6=iBHX?AP1Yn^MZYx|wIUj7PZ9U@W{{Nz9qOvgitvIh|^$?@aK6n;qNQ z*Jj~mWMYl&D3)(9DRV?{visyHj`QB?!^S)6OeOCuE@NJbzHdmN%ssh{U0|zJ1oaIJ z-V{c{h=FJgIVmDv60C{vehK`W?ygQL6?iV zn0#L9O$K)p6Acrm{zHGi z=p`i%6(PRU*Kg_PB9_YWveCWq-;Au zuk&-br_un=YU?5;RU(8n@l0K-_2eFLhOHHi4)4~_a}>3^9q$#QdbVqEu+{d>Uy9tq>rw{1GdnA~`V#7~=5$M`8k!;piRuyv3M9blIwp#MT zPIdR8KS?X6&-jpc$0?40^55zzdn9%kjLa88@)DW|9^Hkj=IgZ{e678gz|PSq6YzO; zML-q;DP8sM=_--#wT%b&zkMILjA$_*(B$t0t1xiC0GndGc(7yrEtX-rqZweAAl`Z_ zbWnSr4=!lBpkzl?ztZU&eECKAZA1o2LoCzx8Ukk?rbeDykHGl%jxUa%(}cJ>K4&es z7>Qn+rAFXfSw7xLVcVX85V=}UoK(@UHAe9-aJ#m2$_8; z7jIr4EJs0A^NcSJGq0 zs~8&N2EYG(rX+ruj-A~pU=aJ0Z$#h7Y3aKcT3eLBnVrZ|K6o+O_FlnZ7&K>P1NCJg zGjp>e7eh)VqSh%B`uWfT!t`hmJ}A*%d1AuCqAD)dKx-vj=*j^yGn0-6zs&+D8+zeh zYBs7i*s>`|52rGg7!1dX7dLCNRI|rJ9QJTv;Gajha{|m3CWHt^3q1{nbD40)$!ktp zc2za9U^%1t; zl@q}TpP@IAB%GARU|fd5#@~+AOCzR-e+JJT^0s{z4l=rR1-u`E4GNJR0x}QwiY-#= z4ee8!;8eHaP+XQ9S0epdRx-}B_k#GZYI3M2&aKw$(6^4K_b}+b1g}ThZo;4SgFSxB zT_6WMp#gz785-QKuT@QscWM?#CGSfC0~?ee2)=}^QU9F~TNLZ}fVwu4gz4k>h{P?+ zI_lZt3QJx$BlEZn_wdIH+)HVO&+Y^4|Br{BcM(g_|+)_C*I@3ZxNWqlvL z=N*t$(SqLd<-aiU5mX9GwogBrtlX41C45rHaQLFd5)MBCuWlg|gaxCyidt$UHe}cj zD{J>3&hfUOii}ltYH2B~FbbP5H>q3wU6MED4rgH4)Ak)-k^4+A4ynyfcij z?G_Ij;8=Cme0~`xvxAh-spa1&4&#cD9@_}z$DC;Yro3o7q`N2I@k3?nN<%9KnUx+_ ztMKeCq1Nf)OeBrHrXeSd1lzDKB*O|pq~B9CX8eb9CK1s-wyy#*@apMtWrGKT+yA)X zKKa$ktpYN&UOVMp)Bn*N$$TB2Y|Dqq6Yp8c)~xAl=(CNQYM1qXo{XQvua?NcW}EN12cG#N ziof^uq%WvKfi40kqLF!>Wt^5G;|A{7S3;X^ttnUZASYswD>2ANFj7+No@f4EITyD0 z2b`h|nZ|&x?hjU`g2aou#Jgw09M;-BUKp<|9UfqVRhHQhxl$JrwP;k^Q#yeZDBkfV`?{>voB*D2r5`X9}q zK79FO4EswOFOO`9@$e-rXr|86T;E7y$7{#m`;J6o41#VofQ*Tbnh04352v<(y~{-6 zx`a!T|1}axJbiwv!EYR|(ZqDOK?Ji}n4ADuA9pbpE-FJ+ZkHiFV0P%O!EJ~&Mg$ci z{ez~UCBLAdA1e{jM{&^iHR}=nuzZNSoH^e6WiS^9USU27vXDG6P!?)YOE>L8M*_hm z6+ZYeiq<=yfb6%7d7ly)Ya{y{2~u<_+7-uY zH9mEY7{_G_>ztu!ZYc|i4PijstBja%LyMfLe&vCB=WIH@24N0kYSGwdXt)M?xc^^a z^H$Z6*lW6OjuSAXxL)%zo7kWD7v@N6eOv5xw)0=hzZmKlwScxW@A#qIm_cg`2|v_A z)Thp5@c#vyTNa2Eou-1}Uoav)wi9B>6&arszt{Al*iNKdpVvC)Ru(h~m9y?uX$gxL zxDZ!cQS5w2McY8C(IQ&AS4?#&nhR9122{UP<@pxE(@Hs}y%U7Jsu=UD2%mWgAM}X@ zp=t&fVJ$SK(^xU=@M;ENrR=0;zB7Z1!6w@VKoZ{a&;M?gYXbAwq6nVR#Oz3udZkA_dgiD>oPZ z+*mF~_PjLJ@8AT4gyaT~Mgw;ik}hGF#;m!|usL3xaV>N?-qXi%fYvWNJl?f5uqer8 zU!*#mkasw`@%bs* z2;a(h0mmi=PR5!XorN)8BI8Pf=Z{@yGt<0%sI$vrAh$f*)5Ic=zhJqKUC2ok5vk$ub?M#aS)X<>hNor5Pd0^U%D0ontI=NjB>aZ_Jb3INg;H@r{e^?*K!R-siqEV zSULijgygcb&D#&+L$kg~6VEBZ1R>Uk@i@ZsD_) zg+UN>5EZ1A`)a3qER{(^F%DQ9k>c;o?n>wUI_Laabr{jeW?&@nDOmH?%qzTaR~4izX=YLkSJg*wzSqZwG1~scP&O&JCpBfCz7w?^3$Fh zn-u6EKDTwodQ`{!#PFV1N@u!D*HM{*dQ&gTF&{RI&;=FYBvto+3!8>y{UX)&>dCn1 z!o@M)Xs12tp&MHMyMTRw1+S+Z77w^+cO`w((VlP}0$&9|#ZtX)qnsS@*o=U1HRqcI zC(L_YpdmH;EQKAmGa!Oj_L{xw(q>(!mkFQg`CAUW(6Z=Psmg^-Q)lGdilNz>O^8od zmo4&qGPR%01Gm$XptFYTWT7UbTWJ#Yr$Se#IvrB*>QkJ0;(ZI@#DC&N@sr5x^O%AI z1Nbj)JQp(9d5|HQbpLtoqG>i>#|QAQIr7h={aZ{$P95N-0dw)*W~WAD-#xbWR{Llge`&OOf; z-}EII902}5IsTQVtrCVNrU7_xMRAddeJv*6TqK|b-E^RB4h$^@r792*2nNWmNmgir zBWn|FO*rG;)WNEK$#E~(#Adj9DZs9ZN7Na?X%#&#v=T)pCP$a@RinAuz1-o~2<<3L zAlc5$o)z%vcb(^X&U2jKSqXSZkD($A`;^5Qi<_qj1~jEsH}fNi{G;-_^4vyydGznO zyV!wZ&LEAPlzX=D4k+_Mv%82GC)PR%Iu^@B>l7_-`zP+6BaaP!M6QInVJAte%m`Ti<{kKhoYp$GLmyBP)PH(&T3&Tujj z=tC%56`UX`3#y#7YI8}Dsu7d-X{Qepq=6?#n5KW0?TvjnH#}Fo)f8^_$Mi#jyU|;2 z{XVmQ@H(7jk5j@h*y99-Syn?>&v+o~&-!nr3I^Lb`VXRuUc?}@c5|R=Rl}DJu0mc- zhV+W?R%QxYpXT2o?Yqj5!3V6KU3G5acKgdQEP&^ZZ2>_9e=l+=?PP>d*oDZwy2xn= zy1UDM$4Nj-fykTd8{TY$Ku&K?}R! zZi&|XFBy(usg0zIw7A*ze1gO!yw`@-V`%o3BOy7B)`6nI6lA}}Rq*nR_O&I4;Fld4 zQ$wr%s82hnYOiPA#{TvJ7<<`qhINYpgozoGwZI9?``p;pk+yH{mIZ)Drk zDcNQs&%cN@qy_XC*!3UTzptTxIaZCo+KhZygg!6S$YGHFVMJsO0q8{G0WCvxKbE$< zCin6u8S{kt{16rBQK`@(Q{YF!;9_EIUmst*Rf3$8;A3OdjPy#t5C7o@@HYxx3_l*^ zfwu>q=8Y3MM~*{e&fR;k$Vp4a`XR2#%I_ywd5<8qzjSmF`1ryCeMO#*(4NNpxpekd zX*0TBfEu?=*bQqUD=(3i`xgZcIb5puE{}b1dl7Xe?{s7USGo|pA8H))drtWNaQp@G zJ2NXAaUHt?h$J}qvvIsju&&UB@V&5W5u*(u#A7|I1O9IMxH`~mf3leRJWE%~5=heP z{$0KFqOp7@_Fs>%^}uZZpd#e^Q%ala(yF0gt^F$G#rXOmHFeBbOXV^Y1g zk8$TjIQEtM+~jvCmRUEwp!A|nD>I#TDO{{>OI@oVQraA}s4}A>-W8rqBVEcMxG1pB z7mH)bZcz7&qte8yK}AIMEET)fXcb+nAvCij zo+3x*-(it-t4i&aC}!~=Rrf`?u;`*?ESonZm9S4&Z0+)K-QIlp>B=<11(6Z!*afRp zc)Uu}0)-RO01WR)O(BDxU^^>>xjpwSzeGQG?w$A0i?xvXKYny5$jqXDpFVoaTrd{0 zKOxv@EUdDh8eAorb^9fs3-afEiV272?r_Itua7PyUZH zJbi*(hQ5(f6(6Plr$sLQzvY>8Xli7%Z)9d<1hFVQGx7thMoI7|`U&{IShqb-2v-Qw zz)DkN)WE3_tLLJT;MkcXlkbw6zemX?(1Hc)@X(>J!Ac;!19pUv`sM)UNEB99u0{$f z78xcYn^*%20~s4-1Cw|I9Sd(A3xz7;EGSZtK;jHg!c~YM`5mvo7Y-K0h>`V85rq4L z44EYUf1(;@{GX_5w9Z`hq|itG8xJ92!|kMVhT$|dy8o1l)77@xYNLa_ zTYI$dJ#9Vb0ef42xq`xeKm|bDCd~^g6RIEd&1YOdpk6!eNhu323M@C*d$hP5(Hf|) zy#o~glXX0~N_d)E10MfsI{^;0l2lF8Om#r7A)#^dn7W0tqeU-w zzXlH zv<@zx%G#y~xopQ$8Af{ixHbcsX?4pY5NwKiu6QA`m|iLAZeO?3lY z3i>ZC0-}zWeRu*9fdz&rS)S@6b3Pi}t_(`%o6H|-g>XPVe5PlFaib2J<{Q!`=g2J{ zFlQFI&7d&Xl))mz#!m_g6J3ra8(elXKICkbZkm<~nVIJn9Pf#)Na7=hXCA#+9kT0> zCb=eehE?i>^gnXCCfeIXy;r8s;aqNozsPJFl6<|2)M32!lEbDOk+Z{6c&{l#igo{3 zi&cmE5=O|QPHCrANnlmkL5XAZ?o4&zq)@0KAzGQii_j=&C@EkFX&4-^nUVm zwMeo&QIjrk$PylU%k-o4hA9AhJ@^AhXf9< zjkZJ>t=%0WomwG0*PH>hUYa`E_9%>dIf6!zADW)nboaTTSo1fXQGx13F5>Qie!e>R zd!`x(iI%$FUic?-ReNdgVSFX_-*4ZB-O|;(q-uyAR|}ceKl`?7h`tWco8vF(U^~Mv znm>AT?#v6&4=n=~Siabma_?M)Eph7uH;Ma4UT~>RJ|qf84);zs0#eq#-0MGgV#i;5 z=9j+%%fChHKWWASU{1GUzlj$-@DlDm3de|Ah@GS?sIFaEM8EM)QZ5qk7ZhF7;Gxq5y6W8OIx6p+7_;mh~hL4UtZQpVHuxITIG zlVVtqv79r;_j^CbovkP0ARBjA);;HDjW0odjVEE&cE)u2+nTse;efYK>B!j!&6x0p zk@G7?Cf=0fSFCA&DD$R2jq$^w!_cWz5c8We(|@+$P)5i&?BSPy2<{Co#}LYoS#*C~ ztGrtb;VkA|lwG+Fs0ZvJ&je38=~O_>;hwXLcChKH9g7$u-yjqEmvBZI zO_BUKmOVhYCLr#Avk)t}j^Qc4XXXmXl$bq`i$XIMb^< z(%8;7T;rLQaX8aD1@tk!TGTPUT=LjXUU^*C)~PYTU1#=ac|-im_Ti{y%!KQ=fyJ|+ z#q&Dv<*pFv_1tR6%&OQ$17N%P*Mwp_V*>LQXy$2%*RFr8xi|gofKBd-dfc&?D&wP- zV*Dm{4Na0nPhNy%?mrDal)!2jW8!oj=PlUjJAA{E?)c)@rqRccgnT^>nwf+0J|sqh zJ0!kB%_H5xzN3O3*%7cEJ!|Cn_rM*qOcA3RHfkL`YMnJ&*(I5yP>WWMhh4tO`_Iig zj(fI=JJ;Ab3-nK|UoxkmuVhAui|S$iO7-EU^QrP-%v0d3KWJCWkj-W|!voHAh6!yf z!;sRi#*Z$vp4o7&b~t?>zi*%D$ev<+T|fr6|7uSJx4(5y26z6;1MKtvHac8t$S7io zpfX~@;*xnTg}@+x{lE}#QkGOC{zQSs&@3vTX1u(PE-2UspE@5q4}b5I*ye&n|A3)1 zU?y%jmj`rZJh;9MrjRk!r!c9rE$ct31N>0Xe~_X$!zKkzNYNVyAehj+>!SG#=xayX z9q>is^*G{ffseS{31TUaS)PziZPf4Qv{=SXv~v#mUQi}FhD@}x4Ec8E`GU~WE*+%^ z&HGW6N57E8Io|ou#>L|gKm2iu!asK+b#`6Bos;))$j;ZP=E~R0`S@qQj7RXE&MtpF z(TthLdqFJ@-q0A7pQdL;c#Gq_rE2nBNT@zb)^MD3T$RxgJsPpidR(Qkt)&=YRr&57KNUc1Jq-w?S4<9|ZyXeu| zId+aWWVr!?f+hpEKHG!*`31I>-ww9kH7uit$VH5W)Y%2)e;W`1emm{^?*jvx=P@pv zO?%^ez~EI)AuGRLu-f|eg2zYXnKy3~9YYeMp1#yF!Qj=uR0~;eh&rl=_xcYG3q>AcdAq+te3%A$H5(=AO+J~oAx33Q^8c%C6g|qFgn+&;zoCij>v;Xis`XCeD$Hm=6Xz-KKhCaIKA;FF z^SrAUL`3Iz+r8%W81bHq*p+9sQ@wZ7Z7396=UDW(_4IS~%}7_X#J%hM=hC;O$?D%xM3b+BQxL%mk@}gQi8)nLPE0pyP^H3z(SHJrPIj4RZ61V zFqI)oBH5uz#Hd+HE-6c{N|N(x*;1M?l(GU}_Jp{Oa-WggA`&b!U0uGKQ+D{;egCVg z!{|fZzNnW=MRfZZ^gez&ePLOg zgrH~kczA8%CTOu+e=8sLTsKC7Cg{&v>9Lz<;zo0|9;*6^q!g?Xmb%ljl>369y5$A( z4K^2pt#<_bF3%jnFLq z9hhUSZ~aB|ZjBE<*h<#Hv&TZ*+~Ix8+=YMh1>IGkNK$7Q^4bAnJTNgYi(qJ^)y4&# zUsa*Xi318g7IWcQPpu=cKGt6&KDJ4r`~kcMxM-!M`(w#B%gI;|KE4-A5E@WFx5an` z%4ss(Tiza}+k;2{X@)S-C<#7b2&pq4sR*ApWy70g5j(GI$HP+##$4n0)t`DWv3ILK zpO{10OL=&?s)*dV@|z#P!%r(Wn61$Pb;7ZM?<*HfF@G@$@#Mg;S~xe|{rsbYW1M@# zHSyK96Ns^X3pSR72Mee(=+EJ3IdvgIwUI6Wn<;lM=qFm0AsDC$KuUXmb+CQE42Jig zJx2fQ%*U6e*HL!>U1s{ToVK3&FD_Y$OHIZAaiMsi_oPGJH=iv|j$9*Mq(a>`ae}7f zJ~Utq;ejs)k*~Cz&PI87E4EVHccpmML^vdI*=FH#&DD5)ij;OvE1ZXb@=}9M4h*wR z8p;QsAXZjmuC#!KP!cJ9XT&WFQSW=-q(Znf=vF>-!JSZAgpfJ=3!8ub;m&==FD2Ji z?Nn>sF;6Gzm|}0P_-=8x4wHga1VRy?1}ui zaWQmzJ$ill+_4hPcYikx%DYCfKJ#ZRxv*f1QDkE)3HZhkiY?4gLJT|_A1qm z&fc3kH%kx@A~i4Wb9k|pT>fJQFV}St((fv!CE>TZ-&Jr{o~fTU=R2X)NMhwaX;Hl z+{pIHwlF*J2T?dL4K3cmm0DF($^+LITk3>>D>3th?cYNCm)A@_9$Jx5oS-+={2F6l z|M14ll=NB`-ymM+l#l7LOA(vejbR;R?g2l6SX-3S@qQWrA-%KnmXI0n|#2OR{o{a)#yOK6{E*9chPk(p`&cm@h1n{#1 zQVsJzu>-1wFxm3k#6ekpG%2pAvwlsW;Uv_-L)PY>jR(2%hJR5ahy0ov7f-1XqJs$Q z8L;gmEl6~j)hnj$FhMB`Sy;Fc(JD&+luhmGZ!{hijdo;0Q(*+q1XYxhnQI5j+gE4g z#|T^euZ|Fa%F5^rP+oqD{7u|kQ0!^HPt*`oLS%=Jbe=7_^i$^i<_8Wun20^C^(a>@ zYxJ8urLix}IPrmdOyuUK_V4K6n;L&@@)ieN`>i0iR3K6{=!p$ta-pxb0WJO7bC~O0f|av_Sb>#?|4@gDJ0K7wh5(Q^V35pRU%zbSp=M9XHPn{ za_9&`c-|9`^3#)JLiW(n?vwq9akcjH3CmD&0UF{AxW6{f; zTX8HVBNyblkd5Kc4}h({WeCtzI8l#zM>ro@JtAy`pksBY^_wMP@P%_h9lZ7h^v6=B zAi7^Ok}d?WIZ#R*vgep>w$Bsn)iUUBvfs4_x%e9*rR2X+DHB>=hk~uf8wbi6=7tYS zxGXRa2c_;efsyADV?5F}he~VGTY9wB0Z1NHF42Gas6PSC1|4Jvw(e{SJF%`@hl_E$ z!!l$$m&K#sOmF^V#5;+sCe9&x45TY@zaddp;%v+a_PX(ypi2nVI2u}h(zJrfwEroZ zlScOsc_fGQ8(&Njbs9}cWlMdv3UG5vrpkAgo^+J5nXbOoB>|H4sK_RiD0$jC)mXoAh8Yd%7E5WB8^=50Jb@jgQMA zU%2yh`IqWuiF~|u_)FqOM~YcCjEf;_P`4IKUg5a)f;dYdr3viQU`6gFT{S7l?eW0K z(fW>)+xw5RM||Z)Y3j`r;+?4O4AN#k<;M{5fpKQBr0XcW^B}{$bI5a<{;N5yqk?49IZ|V zx*HyizujK1SC{e{4-}hW+AIy^Je6sOTb}hTB9&4p4F@pz<2nQveMqVLb+gbLCob+a zcH-c+QmSx>d`G(tUV&Ee|AP};PKF|A8?ci z=RF)DxJEQeFV+wTG!`b;l;VjK{@d7vC0~84daA3qs9P^kO?v&ka@!Hh1H!ujd1bL?uzAn%Pf?|b1K zXl*0k4G~W!^Q=en+%s;-?r$TxQO-Ihz(0RsG12Uljwo5<9kAeLC{WRia_U*ZPNr^2 zwCH8HN_N(*YiWQpP9oLTqX85!wA@=G`yc&k6;zIe{N1XorHILU;T`9GC>S2P@6_-&LBJwzS7Bsxj-=tLJHdhfkN ziQa`E%ILj?sDlY12uAc^7$pqR!svbUG0gn^{`WrKm-}+g`p#Z!uXWaWIN#c9?+-1I za%a@f} zW{EZ1-_-H$T`-I!K(?UcwyR@n;jiNFUeHYkgl??gk1+Tz?qEU`lne)y2;`X9 z)lzxZ7SBQGW+etP_*&h+nyHc|swz>@d49AqsBbj_)lsML)fO9wriZ5+CM0lWI$ol% z@~y`q6Zs~cMFDKRg8J95Ot5rRNA7Cr z@vAn|y)H*ol;g3iJm%1n^KmF@ej1h%bXRZ{`|B}h0awj8>D2_Z#1>!M?32lPc@QQW z2N%$h`g40jxGePii;i3;K6HBN4h1B|QDkAYWU=GsI34>(LDiWU_C`d<2P9(i(w`#YDKiT__2vl!1W=iE0Y z**6?3GmBR>aIIKE677MuYKWEbwmzj;d`pRYr{Taz)fj1o?q`#qvccAS!IYpbSN7Lm z@f|mK-8FY!-c|#Y&9oj#$z_|Q{o6Ee7O*Xf$?@M);0}M zy>=THboJ7HY%-Jz67YpZQ!=yN?*s{=(xDS$dz^VPosCmF1Jjy6lM@qIJNjz+f`VL& z$(qV&(py|}$B*5~2slV()>@IBJJ}?qQWuMKzC4OpHSL?*1AAew$r2HUeeYYK5-JHt zc&k%RNsqLpUm`(9UX}f# zFM37i<>?K>i-Ye^i{r5l1ILz-@EOUJ;Q(tm4(ob!Ums)Ew zIt+aAMqIW@(>EOhlPUH00qGLXn}65&2_}^BhEiNR?FSX{wnXB|n2rfj@$l%o0Q3qA zboH|pjuVG*BK_x*|5OBbU#!JjlUo$ktF7kz5s@t*aJQ6}{quzZubF(}K3?;|WPYQ8 z;bYeao0k1rdx0895T2j3%9u~xAW$uI^>b0x=4A_;2FK*|Nr-!UDWHqax&a%A5VN2Y z-FZDZW7x12@>R7ueq3W~cE_629e#juOAZRYc;4L&$W+gknW%f*#deeeVKn&d0d}$2 zsag*(n4jt06jQc+G4JaA;%jV8G-20g#m{y53lRFDsTL{~7v;&)F7Ttp@-@e14(U1%5TA*@Oxw76?OSB0pGWWMev91s;;Y5Q3##HIm4d$gJ#V{3e&Z$vK9F6 zQ9>NDXl=QSy`j4R-#08K3&g&?eS7I+7aK2EzU;o4P$UExRWLPza(_)f$w+adhFl57 zjLaq^*|^@alV}@(tsXcYPV~|#RCG6O%78>)t{Hb4eEU%O>|@<%ez)_KR2xmAC=vbvJ~FYW zjmQFwpgdBiOU85OtZh=amfwvYE2h_W97tTd1`d_y+u|@5m-`%Sc52uuQRj6XB*t_t z1}R1JjQQ)|RY*lvG%3J8@kEGPIpYkLy?8QDpT^M;w#-S}!l=DZ7f;wCV~Ln5q?C!? z7NJZvzomQ zV_C8LJI#JN9#}{=<5h%H;4EjA! zNilaBFo29;8b*H*$YUZH@|Cv}L`wKnf2PB$W1-V37Yr7vI!$u%D06&VE(xv(>PcF^ z|`$;-<@CBTF~j* zu2Q?PXC2e&#E42|G$Yh29Wh5dOd4619*)g*4)sDJDd#nPFwft^`M`7UGV32N$0#-1 z#2WV%#F<|w2C&)uu3|LgJ|0||zfvEih-vL?QbS#2@_b#4k+Qvwz<(AgQWQ(;+B^F% z=_&lvo8-an(%OGR6|$U49HoIdSDVyw@+oj(=2g9d4yTj5aj8be{4R>p?EqnK83!LU zonohZsERneR?Eh7ggYJ_$Uti4=ZE?kY_kZ8LnS<2M#g@plFj?a7c9^BU;%pD!qymt z*8h#xWxoPF`f->m0oCE{^Z5CO!7{lk<%7Ke^ulEDLUZuW5bxph57O;w$s-1g|p}XNgnEhtIFP zh;#0M8pp4nFOX$P3Tx7PSvUf{^*c$u>b#+OKFX!;HaGe14_QgFFz|!y;V06jt|gMD z^*pvL7sa8~aKKC-DtJKz_-Dmj_nR{_7_SG?`aR^$`&HvS`HGU~s}Sg|6K$mCTQ>tC zKc|XqyY0_6*>{`@kLucXZkB^#HsK6P@O5G?EbAo)1LDriN{pt&M&F0iT?06h4c zV-Noe=5xpnxea4QC0`+lPjXF}4`Py*`t^H-lTRh{M(Y-n9N-^hO60L#7!SU@D^%D@WFeR)<;o|-x)mxD zGI5c$+l{gk%hw3xuZWa1U-&M}-B@6Bxp#Vob4G(1w|G>2w?QjChYANQH_1g>1=o$* z^Dr6*U;N0!w znrGQ1B9IISFRQfq(U%~fYqyEs)#K2}p5tzwK@f(PjC` zSoA9B;A-OZbi?48uFjSW*}I_Yn~??Sk6_)Xr7^X&tMU!dCuj@2{Icpdrj=!31XG8# z6fouEJ!A1WBD`08U3x_1)!;s%&{yKkYOOE8Ge^Sv76^o(+zNlIbxnKt(<7Sff#W5E zWZom8;=;|kGHs*Zd##4KdsrG&Pnond>23JaC;e`ED#Q)ca^}O<+eAqC7IHDG$$nVP zy{Pc8+H|3<>8G{Tt>Y@^C?dr_%zs3O{l7h@r3;fGHKenUHS<{j;|$F{~Y_i}dpURTM_k~-qpA)QDRwRF! zD|0y@)+(>q57xacjTH&4*4$973|OwPWC>WSVV6t!)%_Qycv!~~(X9fR<^ZArejD#s zWKz>a*Sn8L5E>Eep?oTF?WJ3ufBh{I39{*w_SRe^%T?8^q#I*(K>h1Jj<3_}mDg}m z`p*w(7-xpJd9ZMT&2CY(d^sN_8D+^%2Az~P$<*AMMc!@p4!LK>j)?b_&*562G)|%p zH0kp84&?Qfg{;(Y5BcnI1OQzgKozorKFBPTUyA)gvKXjRa)*amZu>U9ob$lF`|lcB z;Q9@Rd<2xKxa)ihd@^)>#d;NXvCY@@?vMOcF<`Tn%ePx|eaxJZC zBPgDZO~(hWDj!8qFw@4JSqd6$hdQVJITH*AKk0Dx5KQ6(b9QdS?oEWcXN`HUmcss3 z+PM5v2y%H%L|&_h2QGp|UtBKG{v4g+hKY2!k(2=Kc#z)VZQoM$p3ls6yt&tU)tOPR zY9E_=mRTrJ!<@UHZ$H^y*2LE0H0TwEW=9yda|2D>{r7+QWvFrcOz&F}jH&pzu`Mj# zvG_BVShBPc`LXz{FsELqW*ZyG=E~v;*O9`=P)j-{fgDz3IYDA4K^Lm6+u?FHUWxEj zT%_4?aJ4?`YRm;D@|Pb*{xrlCqek!=dBb;WQOWV-sFwfSACzJp;B!grS9GWPA*1Rn zqSkxG3eP#|VKs0q>>4MPx2}TPA_C#V1JWw(40s zU4=rX?&kf222PE3>IG}`*&vsLnJ&antT>kYeUWR|?Ev+yAHIZkul-*)2drpOPr+P( zcK>mEm_1;nD~|Eq$(RbhdhbV`l~;IR#A!Z*i(C@p^kbE4bpS+bpV)>_NMfQ7UV|AG z_`@FNNl>DU%w)J`$`yA0m*=(0l)QU6o*wS!uqv}4K~c8Y_JHmKn%3rB+k2dq=|gia zH0=oFY4iKU1M_F-oDqm)^ZU31b9QvIgLTfqSl@E{YDatN6)k4sX(&*hRpMs3HoVwZ zF0)9trL9#uYu}p(zWde-sN50HJSw~X_DC<@HhGNdGtBX;!?g_BkS1W z$q6lNTe+`Dw?_ziCPl29vbp1q^dn$(9f4ya>BGK0jr8`0{^Jxx6&G6r?n~k2X+U5o&A4e$1baGk`4{75|SlxErDU{(mW+ zXdmOx6#g&8=l>{BEGlNF5;NdibiHNI5%pHYU`K~F1+Q>?#c<4B^|>&pxp8^2jUwUY zixbqb@xN z!DL+CP!mC5)!=~kqBGR9C9OCaP=GXk1eGyL9|Dx_uA_7psv0e2sfwizKY zKP5Liy%7kGuiw&=$UJ*mg>)lhNi_r`Pf+3PZ~6<9dQwlY*VP?rjkg|7L)pVaR-L}O zFmq$4gW~(hWjhHVP!ACSZ--88&_u+OqTNY*@R(fP-O-+>+C8;BM-nN4%?d?WI`TDb ztf`184XnP|U2nlKab|U3S|mYa0626^d3sM^ZSG zuenPR?qIf(0y$?8f{Pt!6sFc?;H-`5Qob{*G2=S+CC*|lh* zH94rES_GyM(s{Xl*k5Q_`{~-UfP>XvAWwyD_#Tiju+drr= zw=J-a+B$xFkV{Lrm3EvE>SIHWN^Ph6|0_-)dINXYR`3Xn+zVgSp1(|N?lP(~F8icY z@NrujrGy8CYIz%n#E_brsO;so5%$Q?+;sQIkRRPN8UY2~<&Agmpl!&F{=5#;-gvBdUJC%quq5*xvVHR3D6R#aR|Hw8qM$nhy1ub zGDH^g_k_-Fw%=GIy4@N}dT~SQ$YjU)&)W|(b_#jQbYX4mH20?7@Dxz!pKJXotG8HI zgz$#EOs?d49olFR4Mt#+x=@7w33K+2_s;1Dc1CPYv38Gb@a!o6Gc)7i{+bU{IKVL^ zoQ9FmI5Ud@UCsp3p!Uilz;v8h8J!wQe3@_6aiL82XT9qzUR-DRWnsgg|JcRVluZ}? z;i$*h_1Td(cbY_tv!4;6lE{*JxYBMPX6obHv2#EfNp~MC$at&?O+*2X*yJRanEf(e z>Yr(by_^xzj5*Uv%Fri#bIp&_^!$$u9>LNs6G7UH6kgm76a7g%Hyxf^P|S5ON~I%F zD&Me{_mLm&$vo_l7~?|3LpG~-He>$bNSx5|JdCon*yVhwEFP}#feO4szMd-BQT`ot zkX&xQ-^Hy%!X@5ffBzTr<&e#{g_369?EdJ#H%*oD;sk%r3D@a1oVT*6^@V$EzUMLw zS~-0h$f7RDa{i6{?8PEYd*+Y^{g}-zqX30H7%6=2=W5zWxv!h9xW!_(|I}{GGc81~ zCRjxMyp$=6qI~(K-@EYp%PnV|g;6DGS@08|CTYIYXG&7$BYD>p%j$~K9*CaX&6z*t z;h@HcD2=uf2*p2sZoP!lN5{3JV~0>-NEqbt<_;)(;48M-RgjtK55OCtBNXcuc|tWBO2+k*!b{c>#o1ME>>tpET3 literal 0 HcmV?d00001 diff --git a/Project/N_5.svg b/Project/N_5.svg new file mode 100644 index 0000000..f5e0e04 --- /dev/null +++ b/Project/N_5.svg @@ -0,0 +1,180 @@ + + +11.11.21.31.41.51.61.71.81.920.50.520.540.560.580.60.620.640.660.680.7 diff --git a/Project/lufact.m b/Project/lufact.m new file mode 100644 index 0000000..5ca3886 --- /dev/null +++ b/Project/lufact.m @@ -0,0 +1,49 @@ +function [L,U,P] = lufact(A) +%% +% Computes the LU factorization of the matrix A. The factorization is done +% in-place and then the L and U matrices are extracted at the end for +% output. The factorization is PA = LU +% +% INPUT: +% A - n by n square, non-singular matrix +% +% OUTPUT: +% L - lower triangular matrix with ones along the main diagonal +% U - upper triangular matrix +% P - permutation matrix for pivoting +%% +% get the size of the system + n = length(A); +%% +% initialize the pivoting matrix + P = eye(n); +%% +% Vectorized in-place LU factorization (with row pivoting) that keeps +% track of the total permutations by scrambleing the matrix P + for j = 1:n-1 + A + [v i] = max(abs(A(j:n, j))); + A([j i+j-1],:) = A([i+j-1 j],:); + P([j i+j-1],:) = P([i+j-1 j],:); + A +%% +% Find the index of the largest pivot element in the current column + + +%% +% Swap the rows within the in-place array as well as the permutation matrix P + + + +%% +% Perform the in-place elimination and save the new column of L + i = j+1:n; % indices for the "active" matrix portion + A(i,j) = A(i,j)/A(j,j); + A(i,i) = A(i,i) - A(i,j)*A(j,i); + %A, return + end +%% +% Extract L and U from the in-place form + U = triu(A); + L = eye(n) + tril(A,-1); +end diff --git a/Project/script.m b/Project/script.m new file mode 100644 index 0000000..1790f1f --- /dev/null +++ b/Project/script.m @@ -0,0 +1,59 @@ +%% Section 1 + +N = 5; + +% create three random vectors +main = N*ones(N,1) + rand(N,1); +upper = rand(N-1,1); +lower = rand(N-1,1); +% use the ’diag’ command to create tridiagonal matrix +A = diag(main,0) + diag(upper,1) - diag(lower,-1); + +known_x = ones(N,1); % make known solution vector of 1s +manuf_b = A*known_x; % manufacture the right-hand-side + +A + +%% Section 2 + +N = 5; +x = linspace(1,2, N); +h = x(2) - x(1); + +% Build tridiag matrix A +main = ones(N-2,1); +for i=2:N-1 + main(i-1,1) = (2/(x(i)^2)) - (2/(h^2)); +end +lower = ones(N-3,1); +for i=2:N-2 + lower(i-1,1) = (1/(h^2)) - (2/(x(i+1)*h)); +end +upper = ones(N-3,1); +for i=2:N-2 + upper(i-1,1) = (1/(h^2)) + (2/(x(i)*h)); +end + +A = diag(main,0) + diag(upper,1) + diag(lower,-1); + +% Build constant matrix b +b = zeros(N-2, 1); +for i=2:N-1 + b(i-1) = (2*log(x(i)))/(x(i)^2); +end +b(1) = b(1) + -((1/(h^2)) - (2/(x(1)*h)))*0.5; +b(end) = b(end) + -((1/(h^2)) + (2/(x(end)*h)))*log(2); + +answer = thomas(A, b); +answer = [0.5;answer;log(2)]; + +plot(x, answer); +hold on; +real_values = zeros(N, 1); +for i=1:N + real_values(i) = (4/x(i)) - (2/(x(i)^2)) + log(x(i)) - (3/2); +end +plot(x, real_values); + + +norm(answer - real_values, inf) \ No newline at end of file diff --git a/Project/thomas.asv b/Project/thomas.asv new file mode 100644 index 0000000..f4c60ae --- /dev/null +++ b/Project/thomas.asv @@ -0,0 +1,21 @@ +function [x] = thomas(A,b) + % from wikipedia + n = length(b); + c_prime = zeros(n-1, 1); + d_prime = zeros(n, 1); + c_prime(1) = A(1,2) / A(1,1); + d_prime(1) = b(1) / A(1,1); + for i=2:n-1 + c_prime(i) = A(i,1 + i) / (A(i,i) - (A(i, i-1)*c_prime(i-1))); + d_prime(i) = (b(i) - (A(i, i-1)*d_prime(i-1))) / (A(i, i) - (A(i, i-1)*c_prime(i-1))); + end + d_prime(n) = (b(n) - (A(n, n-1)*d_prime(n-1))) / (A(n, n) - (A(n, n-1)*c_prime(n-1))); + + % backward sub + x = zeros(n, 1); + x(n) = d_prime(n); + for j=n-1:-1:1 + x(j) = d_prime(j) - (c_prime(j)*x(j+1)); + end + +end \ No newline at end of file diff --git a/Project/thomas.m b/Project/thomas.m new file mode 100644 index 0000000..f4c60ae --- /dev/null +++ b/Project/thomas.m @@ -0,0 +1,21 @@ +function [x] = thomas(A,b) + % from wikipedia + n = length(b); + c_prime = zeros(n-1, 1); + d_prime = zeros(n, 1); + c_prime(1) = A(1,2) / A(1,1); + d_prime(1) = b(1) / A(1,1); + for i=2:n-1 + c_prime(i) = A(i,1 + i) / (A(i,i) - (A(i, i-1)*c_prime(i-1))); + d_prime(i) = (b(i) - (A(i, i-1)*d_prime(i-1))) / (A(i, i) - (A(i, i-1)*c_prime(i-1))); + end + d_prime(n) = (b(n) - (A(n, n-1)*d_prime(n-1))) / (A(n, n) - (A(n, n-1)*c_prime(n-1))); + + % backward sub + x = zeros(n, 1); + x(n) = d_prime(n); + for j=n-1:-1:1 + x(j) = d_prime(j) - (c_prime(j)*x(j+1)); + end + +end \ No newline at end of file diff --git a/RungeKutta4.m b/RungeKutta4.m new file mode 100644 index 0000000..c7cecff --- /dev/null +++ b/RungeKutta4.m @@ -0,0 +1,44 @@ +function [t,y] = RungeKutta4(f,t0,T,y0,h) +%% +% Classical four-stage Runge-Kutta time integration technique +% to solve the initial value problem +% +% y' = f(t,y), y(t0) = y0 +% +% INPUT: +% +% f - right-hand-side function that can depend on +% the time t and the function y +% t0 - intial time +% T - final time +% y0 - initial solution value +% h - fixed time step size +% +% OUTPUT: +% +% t - vector of time values +% y - solution vector of approximate values +% +%% +% discretize the time variable and set the problem size + t = t0:h:T; + y = zeros(length(t),1); +%% +% save the intial values + y(1) = y0; + t(1) = t0; +%% +% use Runge-Kutta 4 stage with the given timestep h to solve the problem and +% save the information for later analysis + for j = 2:length(t) + k_1 = f(t(j-1), y(j-1)); + k_2 = f(t(j-1) + (h/2), y(j-1) + ((h/2)*k_1)); + k_3 = f(t(j-1) + (h/2), y(j-1) + ((h/2)*k_2)); + k_4 = f(t(j-1) + h, y(j-1) + ((h*k_3))); + y(j) = y(j-1) + ((h/6)*(k_1 + (2*k_2) + (2*k_3) + k_4)); + + + + + end +end diff --git a/SimpsonsRule.m b/SimpsonsRule.m new file mode 100644 index 0000000..cecf9ad --- /dev/null +++ b/SimpsonsRule.m @@ -0,0 +1,72 @@ +function [S,x] = SimpsonsRule(f,a,b,n) +%% +% Implementation of the (possibly) composite Simpsons's quadrature +% rule to approximate the definite integral +% +% b +% /\ +% | +% | f(x) dx +% | +% \/ +% a +% +% INPUT: +% +% f - anonymous function integrand +% a - lower integration bound +% b - upper integration bound +% n - number of intervals (must be an even integer) +% +% OUTPUT: +% +% In - approximation to the definite integral on the given interval [a,b] +% x - set of n+1 quadrature nodes +% +%% +% Simpson's rule needs at least two intervals + if n == 1 + error('Simpsons rule requires n >= 2'); + end +%% +% check if the number of intervals is even + if mod(n,2) ~= 0 + error('The number of intervals n must be even'); + end +%% +% Find the interval spacing h + h = (b-a)/n; + +%% +% build the set of uniformly spaced quadrature nodes + x = zeros(n+1,1); + x(1) = a; + for i=2:n+1 + x(i) = x(i-1) + h; + +%% +% Initialize the integral value to zero + S = f(x(1)) + f(x(end)); + for i=2:n + if mod(i, 2) == 0 + S = S + 4*f(x(i)); + else + S = S + 2*f(x(i)); + end + end + + S = (h/3) * S; + +%% +% Approximate the integral with (composite) Simpson's rule + + + + + + + + + + +end diff --git a/adaptativeSimpson.m b/adaptativeSimpson.m new file mode 100644 index 0000000..5b96887 --- /dev/null +++ b/adaptativeSimpson.m @@ -0,0 +1,67 @@ +function [In,t] = adaptativeSimpson(f,a,b,tol) +%% +% Adaptive Simpson's rule to approximate a definite integral. This uses +% a posterioi error analysis to recursively apply Simpson's rule where it +% is needed +% +% INPUT: +% +% f - integrand function +% a - lower bound of interval +% b - upper bound of the interval +% tol - error tolerance for adaptivity +% +% OUTPUT: +% +% In - adaptive approximation of the definite integral to +% within the error tolerance +% t - vector of the adapted quadrature nodes +% +%% +%% +% use recursive bisection with error estimation to compute the integral approximation + m = 0.5*(b+a); % find the current midpoint + [In,t] = do_integral(a,f(a),b,f(b),m,f(m),tol); +%% +% this is a recursive function within Matlab so it calls itself + function [In,t] = do_integral(a,fa,b,fb,m,fm,tol) + %% + % need the two midpoints of the sub-intervals and the function + % evalutions for the recursion + xL = 0.5*(a+m); + fxL = f(xL); + xR = 0.5*(b+m); + fxR = f(xR); + %% + % save the five nodes at the current level of the recursion + t = [a;xL;m;xR;b]; + %% + % get the h value for the whole interval + h = 0.5*(b-a); + %% + % compute the Simpson's rule on the whole interval + S_coarse = h/3*(fa + 4*fm + fb); + %% + % compute the Simpson's rule on the two subintervals + S_Left = h/6*(fa + 4*fxL + fm); + S_Right = h/6*(fm + 4*fxR + fb); + S_fine = S_Left + S_Right; + %% + % error estimate + E = S_coarse - S_fine; + %% + % check against the user set tolerance + if abs(E) < 10*tol + % exit when tolerence is met + In = S_fine; + else + %% + % bisect again if the tolernce is note met + [IL,tL] = do_integral(a,fa,m,fm,xL,fxL,tol); + [IR,tR] = do_integral(m,fm,b,fb,xR,fxR,tol); + In = IL + IR; + % append the node edges together (this avoids duplicates) + t = [tL;tR(2:end)]; + end + end +end diff --git a/forwardEuler.m b/forwardEuler.m new file mode 100644 index 0000000..fa91267 --- /dev/null +++ b/forwardEuler.m @@ -0,0 +1,37 @@ +function [t,y] = forwardEuler(f,t0,T,y0,h) +%% +% Forward Euler time integration technique to solve the +% initial value problem +% +% y' = f(t,y), y(t0) = y0 +% +% INPUT: +% +% f - right-hand-side function that can depend on +% the time t and the function y +% t0 - intial time +% T - final time +% y0 - initial solution value +% h - fixed time step size +% +% OUTPUT: +% +% t - vector of time values +% y - solution vector of approximate values +% +%% +% discretize the time variable and set the problem size + t = t0:h:T; + y = zeros(length(t),1); +%% +% save the intial values + y(1) = y0; + t(1) = t0; +%% +% use forward Euler with the given timestep h to solve the problem and +% save the information for later analysis + for j = 2:length(t) +% this is an explicit method so we always use values from the current solution + y(j) = y(j-1) + h*f(t(j-1),y(j-1)); + end +end diff --git a/lab3_k.m b/lab3_k.m new file mode 100644 index 0000000..372f02f --- /dev/null +++ b/lab3_k.m @@ -0,0 +1,30 @@ +%lab3 +N = 1000; +x = linspace(0,1,N+1); +h = x(2) - x(1); +f_prime = zeros(1, N+1); + +f_1 = @(x) (1); +f_2 = @(x) (x); +f_3 = @(x) (x.^2); +f_4 = @(x) (exp(sin(4*x))); + +for i=1:N-1 + f_prime(i+1) = (f_4(x(i+2)) - f_4(x(i))) / (2*h); +end + +f_prime(1) = (f_4(x(2)) - (f_4(x(1)))) / (h); +f_prime(end) = (f_4(x(end)) - f_4(x(end-1))) / h; + +%f_prime(1) = ((-3*f_4(x(1))) + (4*f_4(x(2))) - (f_4(x(3)))) / (2*h); +%f_prime(end) = ((f_4(x(end-2))) - (4*f_4(x(end-1))) + (3*f_4(x(end)))) / (2*h); + +f_prime_exact = zeros(1, N+1); +for i=1:N+1 + f_prime_exact(i) = 4*exp(sin(4*x(i))) * cos(4*x(i)); +end + +f_prime +f_prime_exact + +max(abs(f_prime_exact - f_prime)) \ No newline at end of file diff --git a/plot_composite_quad.m b/plot_composite_quad.m new file mode 100644 index 0000000..d854aaa --- /dev/null +++ b/plot_composite_quad.m @@ -0,0 +1,24 @@ +function plot_composite_quad(f,t) +%% +% Script that will plot a function and the interval boundaries for use +% with a compositie quadrature rule +% +% INPUT: +% +% f - integrand function +% t - set of nodes used for the composite quadrature rule +% +% OUTPUT: +% +% NONE - Produces plot to screen +% +%% + clf + p1 = fplot(f,[t(1) t(end)],'-k','LineWidth',2); + hold on + p2 = stem(t,f(t),'fill','color',[0.9100 0.4100 0.1700],'MarkerSize',9,'LineWidth',1.5); + set(gcf,'Position',[500, 60, 1250, 1250]) + set(gca,'FontSize',16); + legend([p1 p2],{'$(x+1)^2\cos\left(\frac{2x+1}{x-3.3}\right)$','Quadrature nodes'}... + ,'interpreter','latex','fontsize',24,'Location','northwest') +end \ No newline at end of file -- 2.30.2