From fda0cd082dcac7057795000a1917e348619088ce Mon Sep 17 00:00:00 2001
From: nguyed99 <nguyed99@zedat.fu-berlin.de>
Date: Tue, 21 Nov 2023 15:58:30 +0100
Subject: [PATCH] Add UB4
---
UB2/MyIntegrator | Bin 33512 -> 33512 bytes
UB2/QuadratureRule.cpp | 7 +-
UB2/QuadratureRule.h | 2 +-
UB2/main.cpp | 6 +-
UB2/polynomial_new.py | 65 ++++++++++++
UB2/quadrature_full.py | 60 +++++++++++
UB2/shell_script.sh | 2 +
UB2/test_quad.py | 38 +++++++
UB2/test_quad_composite.py | 46 +++++++++
UB2/test_quad_transformed.py | 42 ++++++++
UB3/__pycache__/run_rw2d.cpython-310.pyc | Bin 0 -> 672 bytes
UB3/random_walk.py | 58 +++++++++++
UB3/run_rw2d.py | 48 +++++++++
UB4/UB4.py | 123 +++++++++++++++++++++++
shell_script.sh | 2 -
15 files changed, 492 insertions(+), 7 deletions(-)
create mode 100644 UB2/polynomial_new.py
create mode 100644 UB2/quadrature_full.py
create mode 100644 UB2/shell_script.sh
create mode 100644 UB2/test_quad.py
create mode 100644 UB2/test_quad_composite.py
create mode 100644 UB2/test_quad_transformed.py
create mode 100644 UB3/__pycache__/run_rw2d.cpython-310.pyc
create mode 100644 UB3/random_walk.py
create mode 100644 UB3/run_rw2d.py
create mode 100644 UB4/UB4.py
delete mode 100644 shell_script.sh
diff --git a/UB2/MyIntegrator b/UB2/MyIntegrator
index 026b9ee91cf3c24edc6d2e0cd1b732abd64b7a49..50c1a570180b960467c9c23fed2122ee6a925bda 100755
GIT binary patch
delta 6701
zcmaFS%Jia@X#)o%m#j1c7|2Rb7Gl)nViy4O7}y0Sdojv$#zJ@+v6G7!^<g}Z*vYFH
z&k61@fJnmW$(>9YFf|&nlYcSQiM)Kg_*e3dn3l>N%SvOuyk}AV*SdCd4|6ILW5VW7
ztlOFEp%!pJjAB?I1YsDUi6=nC9ni!ppl;wufarsn1LIpj)iW|MFt9Q(FhG6Ba3LO|
zzMhSNfk70)29r0S5+HGAun2^BfFur*gNQOPyg(9%84i;BfFureJ4l$}2a>o514BJn
z5JU(-eE@PPKS%(I9iZYM*NZ{LKvV(LAa(`@29O?*7zkG&i8Fx&ptu1^oCPWdw!Z^O
z0%RCWZUT}xELcHuGmykNK>|>`07;w+Dh8rfAc=E>1fX~WNZgfyfdOPG4^$9D?EpzY
z#rZ%2P<#MMTo5V-qE5tbe#NcCBInP*;L&=Zgz5hUkLDvBhe1*d|4ps@85sVn+W1e_
z5;$A$$H0)rFYm&@@R?uk*MHSCG=3Z!KMalUhsJk9<J+O}&CvKxdT0V_XnZ*|z8D&x
z&u{V$K`k|328J(bdi?S&U_W~F+A8`nFnIK`N`WZLA0?t5-K_U~CvynNaQK2GOkep<
z))A6n^quS}B&X@ez`#(d9qQ3+8|BNuz~IsOw1mf_+qTV@f#F5wzyJSVr2PN?{}}6D
zzsVDX0vLZxeki1(wx5ZC;S;|AYrQOp7IbX@ajKx44iG0*mVsfivaqMX8g`JCtpWf4
z|M%!T{$lUs8eusmCfUhzguMj}QRSQ_e-~D$w-*4JDZpAP1EK|6U;O+3A7o-Kl>Y(3
zZw1jF-Mt`bkKWJ)9=%{Y()cw_q<Qe?9rWnD2C^?&fPo>+gTL;8M`!DgfB*k`^s=@{
zgM4|+qnCA)H;8ro#i~Cb#~t$MZvF93?*D&~{sSJJ*IzLI{r}&iTT}(C3~b!S$u%M`
zMX&z`sXzYW{qO((Uz`GQ8jt+gTqT;!$RzHz`My{HBhyir$%YbV88bHjmuO{@-SHKo
z;R1;6c72iJq0OSI!N9;!!un#)=7rLJOpLypKgi}VYCQM|cH>L|h|3B=+-}z$9^IiJ
zO(HN&J2$t;Co{5ze)#|Y<l)H$iV~Aq6;l{DPtH;dmDYX#|No00??6s>ec*BM0W&BZ
zEZ+bBKjGx$e@Z%=*_4hjGB!;<tz2Jk!^ptk(e3)eqc=dpqucd{N9Qq+zf~C-7+%YP
zGKELyVX%DX@fQZ)L4LRaGVJynkm0dEP9B~BQWWaZ`P8HNO@K#d?G2C4(hHyX1zk^k
z;urKikjDT10EqYE6Tg7#flvG~AaMcS)!d+{jXPKmlKbIt-1P$}<-U#wJE1%D12pdR
zA*oE@6Te{S4zL`=nIE9e<Yxl8UcjT<^#Ul$4uPB|57rauy5keSfbWY>{1FF1cD(Rt
zcKyHzHlp+Ri}fI1Te^PWmzs2-<v?i{$Z;<|@khG;0ILOaC+n&1<bV40|NjXd-L5l0
zCZGF0`I~BDL@8LT^M*%v=#3N)ZFr2#d<6=g<1aS8`u`skD3Fi?>s$hgbxa-UP#yKJ
z{{MdwHF>>SlW6)&P{6vr@aT?^@aT1Y0J67#vYUDgW76cs>gqf~UqF_C;~;qQHT3|t
zM=!t$;DNg0WM2(U)=41l<bO(1lc#DpFg9+!rs2xSIBT+umbP%i50KVF9=)Jg2B+ou
zo1?UhSSJ}ws1F5$tM*_j^s7B6HW(N_?1$1&K7#|JG_&;+CI*IR1_p+V`~Uw-Ff%Y%
zJox|Lfti89<H7&`5zGt>91s8hXJBSv`1J7qe-jo429Zbq|4(6IU=Vov|Nj#f28QWR
z|NjqQWnftL^#A_?RtAPGPyg5dpTNq%@bu~b{~K5t7!02M|9^p%fnoi#|NsB6F)%!T
z_W!>OI|GB@^Z);2*cliip8x+phn<0;==uNud)OHmBwqage}sd9!RW>R|4%p=7#v>w
z|Ifn7z~J}d|9=%u28LNL{{Q#jWMDY`;{X2)P6mc&FaH1M;9_9-^5Xx04K4--otO3h
z|9fyTFu1+^|G$8lfg$GQ|Nj%185r_j{{O#$nSr6<<^TU5m>C%Mz5M@QgoT0O-OK;~
zLs%FX_+I`0U&6w`VE5|(|06663^}j<|7T%kU@&_9|Gy3^14Ga2|Np<RGBC_}{r|rR
z8w10U*Z=>UurV;ad;R}^2pa=~!<+yAE7%wq!s_4r|38C`fnma%|NnQeF)%c}{r~?C
z8w10XxBvg!urn|$d;9-?4LbwFk+=W<KVWBI;ClD}KLZB?L(IGX{}nhG7>eHg|L?%T
zz##Vi|Nj&Y28QJK|NpmeFfhz`|Ns9I4hDup@Bjbr;ACL<^8Ww-6`TwVr62zPKf%es
zu<pbE|1UTh80rsw`2RnHi-F<ChyVXOxEL6cKmPx}f{TG+%E$lzPjE3XJpTCq{~Im_
zhOZz0|L5UmU`YP-|Gxn@14HSj|NjHH85r6>{r_LU&A_nZ)BpbyxEUDEefs}@12+T1
z|4;w_UjWH}{{R02Hv_|}&;S1$@GvlleEI)Bf`@_O+n4|U1$Y?h8HB$6|6jqwz)<z=
z|NkXC3=E6D{r`W2hk@bPxBvg2@Gvl_eE<KSgO`E9?EC-!8oUe)bH4xoAHd7N@Z$Ub
z{{_4Z3>H8B|DV9iz>xRj|NjlV3=AiJ{QrM}mx1BskN^Ka@G>wc{rvx5fRBM8_~-xs
z27C+*i+=w9AHv7Lu;u6f|0R4p3=B%Y{{NrB$H3t7YjT~jmL$Wy|NqrMsow?C09IjO
zU@!qyj*|}>s~di}^Z!3`<4y!rkTWnaytwoK{{fILD3)PhVCZ3BU^q8f(nOx=$kWOH
z?8GJqm~b#{gD?wBgvCm=85pFcdDuB7FfuSGFfcHHtm>FN(?pqZ(d4}*@{HRi-!+kE
zoHY3_n3Oe@XRWyV|No!Kj;1mKkC}MbIT#o~_B>!@U?`fLYbwvA@?dhGsT^a%<h7>i
zjQ=KIHC1M`pZwEQo=NY)<UUI&wlz!)3^ETUTbao*>Q0U{lV^;WTxq7vIAQW!GkL~_
z$$QP@89OH5H3Qperp$O?vZ6W26^`cej29*+nkzHSdj|Cg%j5;-ER3rsFEy8EJUscR
zxja(<Txya93**|!tQPW2O3x;%S;#SNnVf4WJ~_sMgK^#DR10OsJCi#t<e8YBO<rXo
z$7nM7q=h`A*5sEK%8aWf^IFO?UY{&zEzNp>m4V^SWM4~pMxV*Kmhy~llY1@Y8BHdy
zwUlR4x<C1xr4pmp<d2r}jM0+?t>hVvCL3DGGrCL;w326BI=RqFo^jvgiB|HA*CubY
zl4m?J`J$CP<DJPLt>hVlCJS22GiFRSw3cVwJUP%BBwJ{$%=8}Q0jbFwtXY^IKAS9O
zFURJ?!N8#X9PHx{)*MW-5T<|)2NVDE$x=3QObX8@TiM7li9MejWh2L={(N$kjT~eD
z<e4_|jQ1w*w2@~lntannU2)6v|NkF>5{4_EKpT@YFB{WMaO8nn0c$uI7$!|twAE%j
zH`&(~<lbCcdB$6ldu^2&MJE5X6Xyb@iWp7?2DVp|FWIUwxxAYE%T|uDV6vp0Jkx?#
zll3g6SkG`VFl?I~X(!M4VREINJmbg7GwtM=xL!}*WhciZ{d)2(J0(W9$%*#jT%cqI
z%BB{tC(GF@F*Z$hw3laWnVe{^EV-PK^)w5^A|}=+EDV>JSf{cwEMN})%fj$r^IZF4
z-pLxi8AbxyAYU*rNJ1zE4QQVV+HnLGOkf!)Fv-AR5IgyhZ#`rDW+y*!M#k958=1r=
z7x>R)Oqk3S@LVzxq6*e!2?cjw#6g<40`(bVCtC%&F`k&*7$`1z0@^QYhv;CCWMp7q
zVh~{PpS&?pTrv*Y8G`i?)<R`-CO-@mmz)cg^@i&8V}j^?F<CH3oH2c}W{^AM{K<tu
z;`JxRAxhjK6oU~9s6Qpa0P8uz^1U+)!~_8eh(efe3!&m#Q1N9DX@;9n@e5FKn0gLY
zh&ity;*;5f1wDU2*bFegfT~H5PuZYB?gKS*D^za~RNN3MF2n{gvw3o2u(;$rs4UF1
zJy6*-lP3m?Gd`KTGI%-Tgvo&+`iwg!H-<>Sy5LZ=7@Rmj0U*T?0qt(V#3MKu7?c<U
z7>+@Gq7RiQfQnl{#qFTtT~P5jIfww1V%P~4kJ5sO*TcHU9GnmfxS<mM5P60$sQ3p^
z&S79+FolXA;ACJBVh~_R0_7D328K+i_*+g;KVO<b0@`zi3Nw6xnzK$7A_6K!85kIh
z7{s_B7FU7t30MT8pjiY|LQ69+pc&*j`C+IyW5i^^Fmc9O5NWy=YQi#T7|25mU^oU9
z=aYp9SVAa<7a;M;8<~VtnYkf$*FaqZ^*e()R2<z+QQQzupMi$jA7~#x9wrVA3~<TA
zz`!sWBo4C4_ZG}VMTjF{!N>xt4H*O&-l;>xQy~<CCRDsl86p5nY@s|13_>6i>+7K^
zM4-NgI*DO64+DcBlQaWt5C)>3fnfzyeE>9^VSYUT6%T=i0W_W&-a*BUpe}<_4BVg!
zo<X1<BPl9FBpALy9i#`f*_M}qL6J$CK>!+&u;GkIsQNWf2Q@;~H$uhjp@9hv0frk;
zaY$0<XJCM(^-sJE41%CA@l@u6#OEETV_{(u4i%S%S_KRJem({UAtpfvSRo5_0mIbE
z4<p14H$ipjLgNgUe(plW%b*^G<;Bk+UG-}G5I1;0Q{6|Xp(Rjp6==zQ3@W|?D&C+8
z5rA68@Ej_B3mQ={aZ3S+`7Y4-f(=)NL&TjJUO_F0f|$q9E&@rk2cW|?;FcBx1498s
z1H&n3j!cKv3cZsrMv60D1e1?JWc?4Q%GpW~BVoZRDhP2(6g0S@5yqeg6_<ua9ZbEq
zAjBo;Wkv;5eGRmf0*y<7Y+40Y&oFs2lc4W$sG2&cU2z}<3=9kpq2lPN?l)8%J(whg
zpc#nKWb(%-@p>yEh#!_ggEJ6nZWL4;-Az40knDwCzRiHDUjWTN0Z?<cK*bB7VF9fg
z8LmRa>ls9##T88A2SkEF4Vt=OfhQ^qap*aydeG<*$PNdncp5Z;2!jL}7#Q-Q;`UIB
zVI|0HVFm_41{ns(&}KaYg8($%H$pX_7lv1%;y<7v1Z@~F{Dg{chg#49wKzZo;!rhc
zf`+BTMiEHOfKs4?hG6@^5;6=C(AEjm9EQnIgVsWmT^>|oA5=U98q}}~`VLeaJp+7&
znvWi3jG_z-qM%Ck4>XIxg4$v7#8`3998rir{y}SWSOk?o#dDy#V9vS-6-TdzUV%(w
z{0AaU<-{Onq89=|P;vBBSR}^4pa3e`82cuFj1{k63RPzX%_%S!-++qiLgNA^{sk(&
z0h&0np+P7k4spjJs5wDUac`*jEoj`q63sz~coc&cH1#cpYWM(=U|0y%04pPXB_I}{
z*K0*kaUW=$!kQU<lOM*3o6d#GhCq`VG=dpUL&c?`x?s-YmxN_z)ya<W;+FQ35L?lU
zv>2#3dKo(#D(<s+Vtf_HX4{N!Op~W%3Qp$EkYh?o*=&%#jeqltst3$Wev*?<)#-7;
z>Inhu$xQWbj5{WK*6VX3W?&lY%{dnsLqcZ3<b(C*oKv8KIa4MxHduh#*A4cJv6HL7
n<SsDztHGWzezH}gJ!8V;DloaL@f_p%$(>E+jQNufHmL&u8~P^B
delta 6686
zcmaFS%Jia@X#)o%*Gee{Fjy%yS%^`O>pef1$MBwivKON~r!R!3;XApAQ6I+h@SVJh
z@toiS1BfJ?p4`cl0aK&lJNXw=oye;@_e-BXQ0P8TJ0t4f|1$9!Ug!Utdze$17y~wc
zV%^SM54C^;VidyyAqc|&O*{cA?tms<0d)gM07M_m92nmMs-BU7fq|8Qfk6PGjNyVm
zgi+7Nz`!61VS~vVPzjJYGgt&dJU|i$$w5RJ7+xTW!wd(>eLxb2x*a6U@B>L)gn^+R
zEC?b5pgsV(6dDI0VFm}NILP&4U=avW05yo6fq?;}2O`SAP=O@Q1QTRnU}!)RXMu~?
zGca@@Nq`K4$xT2KhXpH0ZU&M#CrAK_7a)mqLB&AS3M6rEkN^~K0ExRYFff2D<$(%<
zs2v~)s5l=;0E!PFi3>u-K-3BU&9AtXSmf*(7(7}Jlra6j;L&`9<1k2y;lHVsJp;pk
zRU7-sS^{V5?HCx+_~l&~7(VmM{raz(hQ^OW<A<T~{m}SsXnZ>~z8M<dNe@jx4UI2{
z#ur25^Vv<_A*iKh%fRp@O^;u`1?)$UURy;w1_qB_Rw)o=`J+VCqnq`f?PLxi84g>J
zgy}2W$vQ$ZjJA_Kh2%8t7#J8zwL?9cZKG@%7#KV{pO)}=blbMsGBCWz{P+L=i<JNW
z{~u%BYd3j<PypkL$q$88<m#C~A@}~kCw@WKhEMzgtYxwc3?NQN8vpwPleL9C4NBNS
zR<#EF|Nq~k^Z1Ka4hDt^Ae918l?@=3U@aigPcoC22zv|sWJA%!HTl1=LMpoeNHBDP
zM=#X0Yaj!o1sE98JoxJlcyzY@`1k+6M=$FUX^``edGxY=vIen^zj*ZrWcwkHZm^Ox
z5B|IZ9-Y@;IRE|s-=kYJ1*{5e+Q-RzL|%%DgOwkDq5tRq{})^!PUDdmn|FyOGctvn
zZI%}gU}R!7om?n!mT|>qd&yQN*&km)<{y6{@C{_3>x&c*Z5CAx1_p)_))#j+Ka}=k
zV(i^)AeY0asqhKx#G4?$bRK`P@grEn4v+3okR}nBrk@`tduvE;z9;X=$U61I|Nkc^
zFVHlZoS^8*{`vj?|0fSmn7mL^V)6+^bE(|-|Np<Rcn>z_fycoI%%DK2c>n+ZgvswT
zq$kf;R+yZnBC<I{X%Qpip~>7T^-dg&3=AIKt}i@#10+1UU2k}F9s>o&M+OFl*K(j7
z;L&*)EZ=$jg~4}_$8Lbklzt1cBlZWxa*=viynXOE?)n3i;$Fvtb#;gS@aO~wrUDZv
zF$f5J;uj3P0hWW<_ycMq0~5$00v_G24?t0J2xN-@SWl$ujZgdnz8^mEM;rvDjSn8p
zu0I&TMsyy3u>j-&eoNOM{F4r}94PJj#4q6bVe%r?o&1vD{{Nrg(d{||WE9hn$ysWN
zCM&;!#5!+ybcfzZ@z90`G4E@T*N(qXef|GG!q;G(HzxmA6XU!63ZmfitN;IBoSUqy
z-XwbaB`6GBUwCv!NO<(RJ^)#Bc=9gw7{+On#WmD<+`fRcfP?G!WG{^Xrm&Zj=V<sb
zvAvx9MkAJS>t-)aS4PHjlP73t3vc-WQhUgw7Zd~F@VdPDoR$&mWC5O%Nd^<3Rlo;(
zP|z|keAv$@&20UIiGd*+RBGS<|6hWcfg#|*|Njom3=9bm{{N3)W?)cw`2RlxGXn$H
zqyPU+SQr>|9{vA6g@u9P#gqU4pRh16<UjrYKY*2iq3-Gb{{^fJ3_VZ(|DV9hz_9n}
z|Nk3U85lU8{r`V~m4TuC+5h_gf7lop_CNdoUxuB5;pMac|6|x07!02O|38PFfx+wf
z|Nnc~85ll1|Ns982Ll7wi~s+ha4;}Py!ij0g_D6n?ZyB9Dx3@qMKAvU_uyn;Sp4Gu
z{|rtBhJ7#o|L5RhU^w&Q|9=fG1_l;TfO0V~$i4jkzkr#6!Q|!t{}Y%Q7~Edg|Np;%
znSmkV<^TU5m>C$Rz5M@QgoT0O*vtR_Ls%FXp1u74zl4Q>LG0E4|3_FD7+hZc|Ifn8
zz`*tT|9>4;28NW^|NnnsWnd_I{r|rR8w10f*Z=>UurV+kd;R}^2pa=~#GC*BE7%wq
z^xpjcKZA{dA>+;e|2x<i7^2?(|9^*#fgz{<?f?Ha><kQbZ~y<VVP{~N^Y;J$2kZ<C
zkKX?O&%nXJVDj$&e+3Q(2CsMj|2uFnFuZ&B|9=Vx1B3PZ|NmP!7#IrP|Np;)gMnez
z`~UwtI2jnuy#N1y1t$Z8?}z{YPjE6Yw0-#h{{<%l!>kYg|7UP9FkJZX|9=M;1B3O)
z|NmETF)-BUeEk3a1Q!Ft?vMZfzu{tFIQ#Mce;#fI2J27%{~K^KF!+A@|3846fg%3W
z|NjNtpi=e!{|Vd-49h<K|G$Bof#Le6|Nk$5<Ujxa|ACu<q3QGg{{}n^3~xUF{~y7_
zz%c8}|NjC!3=G@8{QqCU!@wZ-?f?HJJPZtO-~RtU!oyI{koE2V|0g^Q3@5(*|Ifk8
zz;Nr^|Nk1i3=9t6|Njr*Wnk#}{{MdgF9XA!@BjZ#;ALPC`tkq&23`h+oFD)HU*Kh6
z=>75k{|8<MhGRee{}<q6VEF#y|9=BM1_rmE|Nn>ZF)&2@{Qtj%kAdOX&;S2t@G&sF
z`1$|;4rt9T!Z3M}u^gks<b%fIj0ux(8mk*}-2ML_xhW+ADqcaE<L>|e2S79w%P=r7
z^e`|mESqd-BF{Kya-fMkW8dUL6D_e`X$A&qX&!ct35*O33JeSkDhvz^36nRPC^J?~
zzG@=Rcy97n6M4q0$)ctp($rL*b;8~M|MyIeG?fup%*4aa!N34A{Q)BbgV*F%Q+cL8
z_a`qim1DG+eAHB(@!I6Krpk=slLgJ>8QCTqn#nVMnapV=&3cB3fq`{$p_x3R!Q_c%
z%8VJ4x0=Z_MohkHCeN5K`KuYoE>Ux3#u<|>%|T9$G?!;wF}c!QnXzp0TyuHG=E-}_
z<r!yBzH2Ves4@AkxjbX*WLXP&re9Ad+gZpl_DqhokY{Y0Tx+4sxMlJ}3wgy`Pyhc9
z0)>VvpFkUvGcOyHC*18etPBi1lOI~hGcrwPv{YtnnyhLm&$xcFrnU6s1WOjiBa>4t
z<r!5b$6AO_R<LH7yugx$N$&pSQ&v)w4_LA=@=QKzDbMul-sESNN{mXAIj!UwjVEhb
z$un|I_Oy~`l$o4qCC^wpxzkFXaoXgSR`QJNCZDvDXPR?w@;-AZwiI>-hAsCd$61I?
zR<Pz^)S9emEzjsM+0$B{v3qi+H7K?^t(6&%PhM*+&$x5)S!-oR=E)ze<r%+Do@*|x
z+QY%X@ch~T|AxqZ%K&=~94~7)7#My`erO@?35umNpiGHJ^&1WbhX2q0|L;IHjj0}N
z97sP8Cj&$J<c&7+j7^g-+Q>6{P5x-3uGsVZ|NnI;I+lWUfGT^C{brLbZM7MTCa2oU
zGgeJ*wUuYAo4nLknQ6z%$@9#`CO@#{U|RKZ@+(_4#TPIC|Mx*L$r)@C$nG941_qJI
zhIaCdE|UZ8<QbzT7uv})PMJK>PM&e<<c)UnOzU1v&a#$b^WkP-*#ByBm$lqv0ecpv
zSFa|ovzD4{V9&yI=hftM)^e;{xEUDKCP&)KGikgAyI_JnkK}wt)@~MtSxl@;SQt(*
zv0i0ic)%3AmxW=$=B@U{yps<wNle!8%`lpy1#%t(gCvAv(17-uq@jFR^Ir<eH}IW&
z$hV%+f3uUHI3uI)<c&;XlMDQ3G6qcM3V1FV2vG&=l!SsiDB=tTzLU8E^%;F9TLro?
zc1&&z6qoFP_Pp95Iv6Aw85o!t1Q_@yZwwTdl!JD6pd|&vTBxkf<cERclCe-(Z>U~B
zCWzh@lLdps8Pz9i2DvlFPc94+ukRFxC~=2S3`Q)V9+LzEtZxV_Yn)jiCR~83hxxV;
zD*g&8z6@&4O{n+;35WnpJqIhqoK+C<$?U;`o;x6H2AE$!bu!4O$Drci^ufTuuobFT
z2pSdO9KgW9AjAeS(|K}Yu()I#R2F919;j^2<cY!Jj7ui33|`JyF*z_qpK-zD#t;cu
zCmZ5S1_mb%Pyk3VL_j-KF!2Zu1_mJp0frt>amm2I08ai43=CZy3=B#P0t{_XAL>Ij
zEP{&vfI0`%u>eW!fr`t?fp`oI^-!AOHB{o6CPV_(F}C7_STG$b4s#c%It5u=1Ily^
z3=F0aZ493{Ar4i7Iy4h1F3$z(;Y%|}K>N&4;Uoq%E{FwIst_Si)&sFTVB(-kf`Ng-
z2qehBz|bfHsyCz=7*GvjsODl|5M&TyV1N!gz&g@hlOKkPGp+}bruU&{ra?nX9%|+f
zsQ5Xk<(5!!DQ*S^rO6wagndmxECvAv8>k9su4IUWid&(%vWFWI3}>L>xCh$bp9m5M
zn|K4JQxRerEC88!AmLG^4iQg<Pz>r&aW`d%04(_h^Dr<7flREgg{rs$jdQ5g3^REc
z7zCN58DIl15d91c%b@B7pl*lxX&+Qv1nO643^Tleitkf{m;iMk11B%o#q}78P!S@*
zunijAdJugK*1QZ1icHcBuyGYw85Ryzp96JJBUF7oRQx~G>rjf}8dMyT$oUx<V5uC`
z;t>ReiKikT*li3MP{+c;Bor!s8=3@Qq2J5Lz#zmV$N(#Lp)O#UJo#aSxZwt<uD4JJ
z!qUlYsJI!_qp)1~5u~eLg&*Pu4q1q4piUl$RRk4(0@Za4BnWCWLB)BXSq_>&7@k7K
z)1VOr6E_!tn9l-@FMg;wp%8H=hE>qYA_^+eDgsHI2cSbW(D-D?gJ@s~gJ#8aXg$$A
z`C_Cv<9RUo5Jc90gQ|>Hf*1)4R$)PiQ=UPMgcbq}I#BW3XzD!$Aud5LEy|$kZJ@;y
zXe<h3(+aS9hRK_m1bvS})o?+hHV&kKfq~&ZR2)6k{e+662a~uEGy^djP5u}qUT+};
z@k1IkIH7eWLj+VD-A!FWknDw4vN249s@H+$A6R*{2`X*?O~%m5lHoE$yq@6-)YCAD
zZx9KFXK49aSQz5aK4=<-whS2TpyF!K1R@L#m|Up%e`qAYN|2etppsJtI-FV0z#ssP
z_w`T>=!M>8sQ3<O2tk_)4Bw&R<xro2TXhT!41OXIhi-y~1S}oai$LlGl%gEe=k5VZ
z$S_2}!XD~_iBN-bp+3liitmAnKY%6$Se1PXDvq83K10n%kFx(F3=E>6TCD?`#bD{$
zZt}!fanCGKh(Gp0BLWsdMNsh{&<qQ6)_JHndUf*xWE$fi5NRqS1~C)85D0*Zqo=|G
zF$M+&P|?QNGx=kzc>Q9ix?j-31?J*wQ1Q3Wpn{2if{GVF6Gt{Q2&Kd!?r4H0K5*ZK
zfq}sjDxL=Q8!XZ6hlocpyn<?23^nK-M1nyVssUC;dP_hoK(E&dpyE7Ghr${dJ(C~C
ziJQ)b%A(grC!ylEq27Qwi&qksnUyCy#*16pN<wT!FVdo*;^<}UOqlrQiSbn&lZCS7
zHh)j&VxC;VX0bUpyPALVv#LkTOt&N^pQ_X2f>jc8v?epvyD=`9>{+kRiI{R}tT*TE
zFopz6$K-?c=A1Rqfux$rj13l`=5~WUqwnM@Fu4m%{%Wvi^q*|iXwMihxe83~YCOky
Pc5-KvIb-hRgH7rHHge{Y
diff --git a/UB2/QuadratureRule.cpp b/UB2/QuadratureRule.cpp
index 3063ce3..7591a8d 100644
--- a/UB2/QuadratureRule.cpp
+++ b/UB2/QuadratureRule.cpp
@@ -5,13 +5,16 @@ QuadratureRule::QuadratureRule(const std::vector<double>& w, const std::vector<d
: weights(w), points(x) {}
// the integrate method takes a callable function f and integrates it over the domain R
-double QuadratureRule::integrate(const std::function<double(double)>& f) {
- double integral = 0.0;
+double QuadratureRule::integrate(const std::function<double(double)>& f, double a, double b) {
+ double integral = 0;
+ double stepSize = (b-a)/2;
// the integration is a weighted sum of function values at the specified points
for (size_t i = 0; i < weights.size(); i++){
integral += weights[i] * f(points[i]);
}
+
+ integral *= stepSize;
return integral;
}
\ No newline at end of file
diff --git a/UB2/QuadratureRule.h b/UB2/QuadratureRule.h
index 134ae1a..f90d8b3 100644
--- a/UB2/QuadratureRule.h
+++ b/UB2/QuadratureRule.h
@@ -12,7 +12,7 @@ private:
public:
QuadratureRule(const std::vector<double>& w, const std::vector<double>& x);
- double integrate(const std::function<double(double)>& f);
+ double integrate(const std::function<double(double)>& f, double a, double b);
};
#endif // QUADRATURERULE_H
diff --git a/UB2/main.cpp b/UB2/main.cpp
index 9709cf5..1501401 100644
--- a/UB2/main.cpp
+++ b/UB2/main.cpp
@@ -4,12 +4,14 @@
int main() {
std::function<double(double)> func = [](double x) { return x * x; };
- std::vector<double> weights = {1/3, 4/3, 1/3};
+ std::vector<double> weights = {1.0/3.0, 4.0/3.0, 1.0/3.0};
std::vector<double> points = {-1,0,1};
QuadratureRule quadrature(weights, points);
+ double a = -1;
+ double b = 1;
- double result = quadrature.integrate(func);
+ double result = quadrature.integrate(func,a,b);
std::cout << "Approximate integral of f(x) = x^2 over [-1, 1] is: " << result << std::endl;
diff --git a/UB2/polynomial_new.py b/UB2/polynomial_new.py
new file mode 100644
index 0000000..3107878
--- /dev/null
+++ b/UB2/polynomial_new.py
@@ -0,0 +1,65 @@
+import numpy as np
+
+class polynomial:
+
+ def __init__(self, coeff):
+ self.coeff = coeff
+ self.deg = len(coeff) - 1
+ self._diff_coeff = np.zeros(self.deg + 1)
+ self._diff_coeff[:-1] = np.arange(1, self.deg + 1) * self.coeff[1:]
+
+ def __call__(self, x):
+ if isinstance(x, float):
+ x = np.array([x])
+ X = self._evaluate_powers(x)
+ return (self.coeff @ X)[0]
+ else:
+ X = self._evaluate_powers(x)
+ return self.coeff @ X
+
+ def diff_p(self, x):
+ X = self._evaluate_powers(x)
+
+ return self._diff_coeff @ X
+
+ def _evaluate_powers(self, x):
+ X = np.zeros((self.deg + 1, len(x)))
+ for ii in range(self.deg + 1):
+ X[ii, :] = x**ii
+
+ return X
+
+class polynomial_new(polynomial):
+
+ def __init__(self, coeff):
+ # This method calls the constructor of the parent class:
+ super(polynomial_new, self).__init__(coeff)
+ # Overwrite array of derivative coefficients, they are replaced by
+ # a matrix of coefficients for all (relevant) derivatives:
+ self._diff_coeff = self._eval_diff_coeff()
+
+
+ # Overwrite derivative method:
+ def diff_p(self, x, order=0):
+
+ if order > self.deg:
+ return np.zeros(x.shape)
+ else:
+ # Evaluate all powers of x at all data sites:
+ X = self._evaluate_powers(x)
+
+ return self._diff_coeff[order, :] @ X
+
+ # New function to generate derivative coefficients:
+ def _eval_diff_coeff(self):
+ # Compute all derivative coefficients:
+ diff_coeff = np.zeros((self.deg + 1, self.deg + 1))
+ diff_coeff[0, :] = self.coeff
+ for ii in range(1, self.deg + 1):
+ if ii == 1:
+ upper_ind = None
+ else:
+ upper_ind = -ii + 1
+ diff_coeff[ii, :-ii] = np.arange(1, self.deg + 2 - ii) * diff_coeff[ii-1, 1:upper_ind]
+
+ return diff_coeff
\ No newline at end of file
diff --git a/UB2/quadrature_full.py b/UB2/quadrature_full.py
new file mode 100644
index 0000000..dbe04a0
--- /dev/null
+++ b/UB2/quadrature_full.py
@@ -0,0 +1,60 @@
+import numpy as np
+
+class QuadratureRule:
+
+ def __init__(self, x, w):
+ self.x = x
+ self.w = w
+ self._n = self.x.shape[0]
+
+ def integrate(self, f):
+ # Evaluate f on all grid points:
+ fx = np.array([f(self.x[i]) for i in range(self._n)])
+
+ return np.dot(self.w, fx)
+
+
+
+class LinearTransformed(QuadratureRule):
+
+ def integrate(self, f, T, dT):
+ # Evaluate f on transformed grid points:
+ fx = np.array([f(T(self.x[i])) for i in range(self._n)])
+ # Transform weights:
+ wT = np.array([self.w[i] * dT for i in range(self._n)])
+
+ return np.dot(wT, fx)
+
+
+class CompositeQuadrature:
+
+ def __init__(self, a, b, n, x, w):
+ """ Composite quadrature on interval [a, b] using n sub-
+ intervals, and a quadrature rule on reference interval [-1, 1].
+
+ Parameters:
+ a, b: lower and upper bounds of integration domain
+ n: number of sub-intervals
+ x: grid points for reference quadrature rule
+ w: weights for reference quadrature rule
+ """
+ self.a = a
+ self.b = b
+ self.n = n
+ # Create sub-division of [a, b]:
+ self.h = (b - a) / n
+ self.y = np.array([a + i * self.h for i in range(self.n+1)])
+ # Mid-points of the sub-intervals:
+ self.c = 0.5 * (self.y[:-1] + self.y[1:])
+ # Create quadrature rule object:
+ self.QR = LinearTransformed(x, w)
+
+ def integrate(self, f):
+ I = 0.0
+ # Iterate over sub-intervals:
+ for i in range(self.n):
+ # Linear transformation:
+ T = lambda x: 0.5 * self.h * x + self.c[i]
+ I += self.QR.integrate(f, T, 0.5*self.h)
+
+ return I
\ No newline at end of file
diff --git a/UB2/shell_script.sh b/UB2/shell_script.sh
new file mode 100644
index 0000000..1c30853
--- /dev/null
+++ b/UB2/shell_script.sh
@@ -0,0 +1,2 @@
+g++ main.cpp QuadratureRule.cpp -o MyIntegrator
+./MyIntegrator
diff --git a/UB2/test_quad.py b/UB2/test_quad.py
new file mode 100644
index 0000000..280f47f
--- /dev/null
+++ b/UB2/test_quad.py
@@ -0,0 +1,38 @@
+import numpy as np
+import quadrature
+import polynomial_new as poly
+
+""" Settings: """
+# Create an instance of the quadrature rule:
+x = np.array([-1.0, 0.0, 1.0])
+w = np.array([1.0/3, 4.0/3, 1.0/3])
+S = quadrature.QuadratureRule(x, w)
+
+# Boundary values of interval:
+a = -1.0
+b = 1.0
+# Number of test:
+ntest = 10
+
+""" Function to integrate polynomial: """
+def int_poly(P, a, b):
+ # Coefficients of the primitive function:
+ pc = P.coeff / np.arange(1, P.deg+2, dtype=float)
+ # Evaluate powers at both boundary points:
+ bx = np.array([b**i for i in range(1, P.deg+2)])
+ ax = np.array([a**i for i in range(1, P.deg+2)])
+
+ return (pc @ bx - pc @ ax)
+
+""" Main Test: """
+for ii in range(ntest):
+ # Generate random third order polynomial:
+ #c = np.array([1.0, 0.0, 0.0, 1.0])
+ c = np.random.randn(4)
+ print("Test %d: coeffs: "%ii, c)
+ P = poly.polynomial(c)
+ # Compute integral:
+ I = int_poly(P, a, b)
+ # Approximation with Simpson rule:
+ ISimp = S.integrate(P)
+ print("Test %d: diff = "%ii, np.abs(I - ISimp))
\ No newline at end of file
diff --git a/UB2/test_quad_composite.py b/UB2/test_quad_composite.py
new file mode 100644
index 0000000..e56026e
--- /dev/null
+++ b/UB2/test_quad_composite.py
@@ -0,0 +1,46 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+import quadrature_full
+
+""" Settings: """
+# Function to integrate:
+f = lambda x: np.arctan(x)
+# Anti-derivative:
+F = lambda x: x * np.arctan(x) - 0.5 * np.log(1 + x**2)
+
+# Interval:
+a = -1.0
+b = 3.0
+
+# Array of grid sizes:
+n_array = np.array([5, 10, 50, 100, 500, 1000, 5000, 10000])
+# Reference quadrature:
+x = np.array([-1.0, 0.0, 1.0])
+w = np.array([1.0/3, 4.0/3, 1.0/3])
+
+""" Demonstration: """
+# Compute reference:
+Iref = F(b) - F(a)
+# Convert grid sizes to step sizes:
+h_array = (b - a) / n_array
+# Apply composite quadrature:
+E = np.zeros(h_array.shape[0])
+for i in range(h_array.shape[0]):
+ n = n_array[i]
+ QC = quadrature_full.CompositeQuadrature(a, b, n, x, w)
+ I = QC.integrate(f)
+ E[i] = np.abs(I - Iref)
+
+""" Figure: """
+plt.figure(figsize=(6, 4))
+plt.plot(h_array, E, "o--", label="Error Simpson")
+plt.plot(h_array, h_array**4, "ro--", label=r"$h^4$")
+plt.xscale("log")
+plt.yscale("log")
+plt.xlabel(r"$h$")
+plt.ylim([1e-16, 1e-2])
+plt.legend(loc=2)
+plt.show()
+
+print(I, Iref)
\ No newline at end of file
diff --git a/UB2/test_quad_transformed.py b/UB2/test_quad_transformed.py
new file mode 100644
index 0000000..c03aadb
--- /dev/null
+++ b/UB2/test_quad_transformed.py
@@ -0,0 +1,42 @@
+import numpy as np
+import quadrature
+import polynomial_new as poly
+
+""" Settings: """
+# Parameters of the quadrature rule on reference interval:
+x = np.array([-1.0, 0.0, 1.0])
+w = np.array([1.0/3, 4.0/3, 1.0/3])
+
+# Define linear transformation of [-1, 1] to [1, 5]:
+a = 1.0
+b = 5.0
+T = lambda x: 2 * x + 3
+dT = 2.0
+# Number of tests:
+ntest = 10
+
+""" Function to integrate polynomial: """
+def int_poly(P, a, b):
+ # Coefficients of the primitive function:
+ pc = P.coeff / np.arange(1, P.deg+2, dtype=float)
+ # Evaluate powers at both boundary points:
+ bx = np.array([b**i for i in range(1, P.deg+2)])
+ ax = np.array([a**i for i in range(1, P.deg+2)])
+
+ return (pc @ bx - pc @ ax)
+
+
+""" Main Test: """
+S = quadrature.LinearTransformed(x, w)
+
+for ii in range(ntest):
+ # Generate random third order polynomial:
+ #c = np.array([1.0, 0.0, 0.0, 1.0])
+ c = np.random.randn(4)
+ print("Test %d: coeffs: "%ii, c)
+ P = poly.polynomial(c)
+ # Compute integral:
+ I = int_poly(P, a, b)
+ # Approximation with Simpson rule:
+ ISimp = S.integrate(P, T, dT)
+ print("Test %d: diff = "%ii, np.abs(I - ISimp))
\ No newline at end of file
diff --git a/UB3/__pycache__/run_rw2d.cpython-310.pyc b/UB3/__pycache__/run_rw2d.cpython-310.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..eb602d2a62b12a374abbab7324c6114c71e11aed
GIT binary patch
literal 672
zcmd1j<>g{vU|@I?7MyyOnStRkh=Yt-7#J8F7#J9e{TLV+QW#Pga~Pr+!8B7Ya}*0B
z11CcYQxt0oa}-+&OB8!5;{uKp)>QU1riF}AoGENkTq*2P+$kJUJSm(}yqpXvTq)cw
z3{iY3jKK_=JTE~u`K8aA2?96m)2l_+xj(F0VxRefnStRY69WT7I!FZwzZ6066b?Xz
z6&#XTL6T4mVzV(YFgSzU8^yrDP{Xi*VIcz}LkZ&orW%l|8CscA7}J;}8Q?r75N{z9
zSPgRua}ARi!$L+dkEMpWh9#W|%40?FY8ZkUG+F$L7#SECUV?}!>45x#(wxMS%=|os
z#1e&)%-mFk5(QfYRXyW0O}1OC$vLTsMYos>a!PKo7bho{l%y8jVlBxm$w|G%QjwFH
zdyAzqCo@-*^%iqlX8JA0h+B*$w^%B2GxI=XV#P{^A^`>lhF^O68Tq-X`US<AMX8B7
zCPw<jC8a5q`pNmZ1-iw_nfjql#`;C2dGST%Mk#s)mAClfOAAsGOH$+0GSf?oQpG^Q
z!o$G8z`>}($OncjOhrNr3=Eq5w^;K^a|<eOG36E963k63Dagq$$;nL8E2spqZm~h=
zB2ehwVoIsF#gtNVi_z~EYgKAde(^1)3WHlrB?h+`y^Hu57#NE9K?KNPu)C26K?Vkf
oTO2mI`6;D2sdk{SE0zK|h=GTZiG_uUg^`7c4@`pCOiVnC0D%vmwg3PC
literal 0
HcmV?d00001
diff --git a/UB3/random_walk.py b/UB3/random_walk.py
new file mode 100644
index 0000000..3fd6fbb
--- /dev/null
+++ b/UB3/random_walk.py
@@ -0,0 +1,58 @@
+import numpy as np
+import random
+import matplotlib.pyplot as plt
+
+from run_rw2d import _update_figure
+
+class RandomWalk2D:
+ def __init__(self, dx: float, N: int, x0: np.ndarray, dt: float):
+ """
+ :param x0: vector of initial positions
+ :param N: population size
+ """
+ assert x0.shape == (N,2), "x0 must have the shape (N,2)"
+ self.positions = x0
+ self.time = [0]
+ self.dt = dt
+ self.N = N
+ self.dx = dx
+
+ def simulate(self, K: int = 1):
+ """
+ param K: number of discrete steps
+ """
+
+ choices = [(-1, 0), (1, 0), (0, 1), (0, -1)]
+
+ for _ in range(K):
+ self.time.append(self.time[-1]+ self.dt)
+ directions = np.array([random.choice(choices) for _ in range(self.N)])* self.dx
+ self.positions += directions
+
+
+### testing ###
+N=100
+K=20
+dt=4e-2
+dx=4e-1
+RW = RandomWalk2D(dx=dx, N=N, x0 = np.zeros((N,2)), dt=dt)
+
+### visualisation ###
+
+# Axes limits for figure:
+xmin = -8.0
+xmax = 8.0
+X = RW.positions
+fig = plt.figure(figsize=(6, 4))
+fig = _update_figure(fig, X, 0.0, xmin, xmax)
+for _ in range(1, K+1):
+ # Update positions:
+ RW.simulate()
+ X = RW.positions
+ print(f'{X=}')
+ t = RW.time[-1]
+ # Update figure:
+ fig = _update_figure(fig, X, t, xmin, xmax)
+ plt.pause(0.25)
+
+plt.show()
diff --git a/UB3/run_rw2d.py b/UB3/run_rw2d.py
new file mode 100644
index 0000000..24ac93b
--- /dev/null
+++ b/UB3/run_rw2d.py
@@ -0,0 +1,48 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# import random_walk_2d as rw
+
+""" Settings: """
+# Mesh sizes:
+dx = 0.4
+dt = 0.04
+# Population size:
+N = 1000
+# Initial conditions:
+x0 = np.zeros((2, N))
+t0 = 0.0
+# Number of discrete time steps:
+K = 20
+# Axes limits for figure:
+xmin = -8.0
+xmax = 8.0
+
+""" Function to update the figure: """
+def _update_figure(fig, X, t, xmin, xmax):
+ fig.clear()
+ plt.scatter(X[:, 0], X[:, 1])
+ plt.title("Population at time t = %.3f"%t)
+ plt.xlim([xmin, xmax])
+ plt.ylim([xmin, xmax])
+
+ return fig
+
+# """ Instantiate the class: """
+# RW = rw.RandomWalk2D(dx, dt, N, x0)
+
+# """ Simulation: """
+# X = RW.positions()
+# fig = plt.figure(figsize=(6, 4))
+# fig = _update_figure(fig, X, 0.0, xmin, xmax)
+
+# for ii in range(1, K+1):
+# # Update positions:
+# RW.simulate()
+# X = RW.positions()
+# t = RW.time()
+# # Update figure:
+# fig = _update_figure(fig, X, t, xmin, xmax)
+# plt.pause(0.25)
+
+# plt.show()
diff --git a/UB4/UB4.py b/UB4/UB4.py
new file mode 100644
index 0000000..402e74d
--- /dev/null
+++ b/UB4/UB4.py
@@ -0,0 +1,123 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Ex1: Numerical integrators
+
+def implicit_euler(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
+ """
+ A-stable
+ :param f: function to be integrated
+ :param y0: initial value
+ :param t: time interval
+ :param dt: time step
+ """
+ no_of_steps = int(t // dt)
+ y = np.zeros((no_of_steps, *y0.shape))
+ y[0] = y0
+
+ for i in range(1, no_of_steps):
+ y[i] = y[i-1] + f(y[i-1]) * dt
+ y[i] = y[i-1] + f(y[i]) * dt
+
+ return y
+
+def explicit_euler(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
+ """
+ :param f: function to be integrated
+ :param y0: initial value
+ :param t: time interval
+ :param dt: time step
+ """
+ no_of_steps = int(t // dt)
+ y = np.zeros((no_of_steps, *y0.shape))
+ y[0] = y0
+
+ for i in range(1, no_of_steps):
+ y[i] = y[i-1] + f(y[i-1]) * dt
+
+ return y
+
+def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
+ """
+ Not L-stable, doesn't decay properly - oscillation
+ :param f: function to be integrated
+ :param y0: initial value
+ :param t: time interval
+ :param dt: time step
+ """
+ no_of_steps = int(t // dt)
+ y = np.zeros((no_of_steps, *y0.shape))
+ y[0] = y0
+
+ for i in range(1, no_of_steps):
+ y[i] = y[i-1] + dt * f(y[i-1] + dt/2 * f(y[i-1]))
+ y[i] = y[i-1] + dt * f(1/2 * (y[i-1] + y[i]))
+
+ return y
+
+# Ex2: dy/dt = lambda * y, y(0) = 1, lambda = -1
+y0 = np.array([1])
+lamda = -1
+t = 100
+dt = 1e-2
+f = lambda y: lamda * y
+y = [implicit_euler(f, y0, t, dt), explicit_euler(f, y0, t, dt), implicit_midpoint(f, y0, t ,dt)]
+integrators = ['implicit euler', 'explicit euler', 'implicit midpoint']
+
+# for i in range(3):
+# plt.figure()
+# plt.plot(y[i])
+# plt.title(integrators[i])
+# plt.show()
+
+# Ex3:
+def f(y: np.ndarray, k: float = 1) -> np.ndarray:
+ assert y.shape[0] == 3, 'y has the wrong dimension. It should be 3'
+
+ A = np.array([[-1, 1, 0],
+ [1, -1-k, k],
+ [0, k, -k]])
+
+ return np.dot(A, y)
+
+dt = [1, 1e-1, 1e-2, 1e-3]
+t = [10,100]
+y0 = np.array([1,0,0])
+
+# short time
+fig = plt.figure()
+for dt_i in dt:
+ y = implicit_euler(f, y0, t[0], dt_i)
+ plt.plot(y[:,0], '--',label = 'State 1')
+ plt.plot(y[:,1], 'x', label = 'State 2')
+ plt.plot(y[:,2], '>', label = 'State 3')
+ plt.legend()
+ plt.title(f'{integrators[0]}, dt = {dt_i}')
+ plt.show()
+
+
+
+
+
+
+fig = plt.figure()
+for dt_i in dt:
+ y = explicit_euler(f, y0, t[0], dt_i)
+ plt.plot(y[:,0], '--', label = 'State 1')
+ plt.plot(y[:,1], 'x', label = 'State 2')
+ plt.plot(y[:,2], '>', label = 'State 3')
+ plt.legend()
+ plt.title(f'{integrators[1]}, dt = {dt_i}')
+ plt.show()
+
+fig = plt.figure()
+for dt_i in dt:
+ y = implicit_midpoint(f, y0, t[0], dt_i)
+ plt.plot(y[:,0], '--', label = 'State 1')
+ plt.plot(y[:,1], 'x', label = 'State 2')
+ plt.plot(y[:,2], '>', label = 'State 3')
+ plt.legend()
+ plt.title(f'{integrators[2]}, dt = {dt_i}')
+ plt.show()
+
+
diff --git a/shell_script.sh b/shell_script.sh
deleted file mode 100644
index 6d1b172..0000000
--- a/shell_script.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-g++ main.cpp QuadratureRule.cpp -o MyApp
-./MyIntegrator
--
GitLab