From e16dc49abca343af1ab61c81ab9034c7c781485b Mon Sep 17 00:00:00 2001 From: Elizabeth Hunt Date: Fri, 13 Oct 2023 17:20:17 -0600 Subject: [PATCH] fix subst --- inc/lizfcm.h | 4 +- lib/lizfcm.a | Bin 14128 -> 15472 bytes notes/Oct-13.org | 8 ++++ src/matrix.c | 30 +++++++++++++- test/main.c | 100 +++++++++++++++++++++++++++++------------------ 5 files changed, 102 insertions(+), 40 deletions(-) create mode 100644 notes/Oct-13.org diff --git a/inc/lizfcm.h b/inc/lizfcm.h index 99acbf2..40046d2 100644 --- a/inc/lizfcm.h +++ b/inc/lizfcm.h @@ -30,7 +30,9 @@ extern Line *least_squares_lin_reg(Array_double *x, Array_double *y); extern void put_identity_diagonal(Matrix_double *m); extern Matrix_double *copy_matrix(Matrix_double *m); extern void free_matrix(Matrix_double *m); -extern Matrix_double **put_lu_decomp(Matrix_double *m); +extern Matrix_double **lu_decomp(Matrix_double *m); +extern Array_double *bsubst(Matrix_double *u, Array_double *b); +extern Array_double *fsubst(Matrix_double *l, Array_double *b); extern void format_matrix_into(Matrix_double *m, char *s); diff --git a/lib/lizfcm.a b/lib/lizfcm.a index 78ca79eefc769195e8b9beb2aafcacaeaa03fd2a..7ea59ed7ff0f7e88ae73cb986b2e9efbface855c 100644 GIT binary patch delta 2659 zcmcImU2GIp6h6arySpu=P+E2^ZFlHjie}61cIi)=sjHx=x+b;ZfoRG=7f@oM4Q)dU zm~8?g2@g;&(LO+0jA%l zcfb3cbM85J&Yd|kDfVB5Zvk?gA<;h;{bk|VniI#b(JnjS zCQ(!iJ`czGETXy;7yMKLydwO%@K=TJ2j`6?gm5ySKf&N4=+o*6p{aWOejf zd-o+e4|b_rrGE85yI6viR2JS6splGrglik>A}C?F_UZ?#&bVF9G?LZmlIh{x*Ytp` zUAe|vx7J40!@M1NRXH3VRmH|ws&Z6BZz^ccb||Y+FCy&y1@*RGlb5n^=~#>$^_AXK z+`QVMQx1@^hDgV}c^#2k9ju5}VNwaj)OwT{Q8_`2IKxCd+!+3=m-0t46T0^O+%bJ@ z?h6_mIT4t-qz7h3PgpY^Z)|oJ#oWWmT=eI3 zHqF`rvNI(sqwfqCR_TuKbS|%l%^g!uSGtn+hMBW(|A85? z&1(Gz9xeE?EysIil@^<$)6WD*qdB=4XK%=PZ+_mlVl1p(DHd@o|DI&YLNHDu7y2cF zO<|*R$F_~l>)PN=-MS0>A5LY>@?&|oprXzD&xbO<MP!@T~oKV_`F3P-2IRbZFbY^ur9U7?MsHx z{8-Q-X)RHP!xsOec&?uy3Kjtii8}lawygzjgx&*0Ylvq!Xig&@>nDK*@FuVbcpZq= z6u;)6IpFA%!9GmTYvMf)n%jY#DCS_B0o;TD8NgED=eQX3(BL=_uZ`d*K;({(I%tjn z(c$7395jc3SQvBQ%j`HRTdsa9@imlz<%9WveP~t8cL*AS50RMlQ$V~P zgFgWM$_&(ndFvstQs7>k8GFS9R|`(!_*lOa=vS8mHOi=|Pu_ub4Ms{_%n_*>)nnO* z!F-qS39-K}oZBhej|v|b`$^%KgugF*Nce=qlk~C#_(>ejio-qOHDb?^SI;T2=L>&J zxGB6x_%`8h3g0fAUk-~Gm?XNw3RF1tZhH<1bD?_!C(DAV&2QU`=p*5&X9V}TZBIQx zPYF*wU#Em0N0nH9$Kd$3AU}PBk>afUF^#vtxYQm^2kzF7KfiGuly zX#xAw?EHCy--0}=e%|y_B-M>9%i7<2C}COEsjnUBYaNRg`6PO>baTPNrw}c%tFF?! eSq;!sclR8wu2)TEts1GbWgFzf=StbB^?w5gGALpI delta 1452 zcmZvc-%FEG7{{O2Ip<8BOy|tpd~=d3neK=Cp)wPqh!8Oml4{MEiqO<$3RcDvBKiTf zIyYTdh|o4}uh z?(qhKK2`O*J*ujD-MFY8w>!|*hFHk`sMD6*B%-(IDF-fScM6DVn3LeDndluT^AT{+ zMD&?`j`QExTS0k$isQdI{)btrXzN&*4GQ&U7e0XrH6=vNoakU4XFlZxUom?*J_?r` zI>SCB{K^%PR$9a_#afn_nU2P1qQfJx$?=3JDtK1w@U(hEvIyKff%c#mH421Q)(h?% z^xCk$q*Kqs$EI_-`fT+7MQro>y`r>oz#KASD`*+)M7EGo^y|7TI|fZ<7=@vVeMIfz zxz1gsYK|Bp&qdm(S!zqs61Z$BW1FOLy5QXYFuQBY@;ul;h1~_?7 zoyPDsHY9)OhsMYzcmNA5gLPmACzRY^u7Xzh0$2|Yg8RWf&@Lhlt5Mn<4rUiO9M*Ol ztO|z1x@J*d=>traxs3-Vf?mdN)36~8}wrl9pb8l2R`xqV^(E!mp2c*i@k1lVR z!zpE3LYIr?S{j4h-U!+4#fo^*VNW*ZO$+6DbDlrQb9>&h$Mf>RynOXu7e5b$#ZXg| mZO?xW#iP;a=*&bc8g0pbf(uQb#4YD3@xf`~KbS4&?aIID3c1t( diff --git a/notes/Oct-13.org b/notes/Oct-13.org new file mode 100644 index 0000000..853a6d6 --- /dev/null +++ b/notes/Oct-13.org @@ -0,0 +1,8 @@ +* Root Finding +Finx x \in R such that f(x) = 0 + +If g(x) is a function and we want to find x such that g(x) is +extremal then we find x \in R such that g'(x) = 0 + + +Hanging Cable Problem y(x) = c_1 cosh(\frac{x- c_2}{c_2}) diff --git a/src/matrix.c b/src/matrix.c index 953d50f..438a468 100644 --- a/src/matrix.c +++ b/src/matrix.c @@ -1,7 +1,7 @@ #include "lizfcm.h" #include -#include #include +#include void put_identity_diagonal(Matrix_double *m) { assert(m->rows == m->cols); @@ -19,7 +19,7 @@ Matrix_double *copy_matrix(Matrix_double *m) { return copy; } -Matrix_double **put_lu_decomp(Matrix_double *m) { +Matrix_double **lu_decomp(Matrix_double *m) { assert(m->cols == m->rows); Matrix_double *u = copy_matrix(m); @@ -63,6 +63,32 @@ Matrix_double **put_lu_decomp(Matrix_double *m) { return u_l; } +Array_double *bsubst(Matrix_double *u, Array_double *b) { + assert(u->rows == b->size && u->cols == u->rows); + + Array_double *x = copy_vector(b); + for (int64_t row = b->size - 1; row >= 0; row--) { + for (size_t col = b->size - 1; col > row; col--) + x->data[row] -= x->data[col] * u->data[row]->data[col]; + x->data[row] /= u->data[row]->data[row]; + } + return x; +} + +Array_double *fsubst(Matrix_double *l, Array_double *b) { + assert(l->rows == b->size && l->cols == l->rows); + + Array_double *x = copy_vector(b); + + for (size_t row = 0; row < b->size; row++) { + for (size_t col = 0; col < row; col++) + x->data[row] -= x->data[col] * l->data[row]->data[col]; + x->data[row] /= l->data[row]->data[row]; + } + + return x; +} + void free_matrix(Matrix_double *m) { for (size_t y = 0; y < m->rows; ++y) free_vector(m->data[y]); diff --git a/test/main.c b/test/main.c index 75e60ee..2364b3d 100644 --- a/test/main.c +++ b/test/main.c @@ -10,19 +10,18 @@ int main() { printf("dmaceps(): %.10e\n", dmaceps()); printf("========\n"); + printf("Norm, Distance\n"); Array_double *v = InitArray(double, {3, 1, -4, 1, 5, -9, 3}); - char v_s[256]; - strcpy(v_s, ""); - format_vector_into(v, v_s); + char s[2048]; + strcpy(s, ""); + format_vector_into(v, s); + printf("v: %s", s); Array_double *w = InitArray(double, {-2, 7, 1, -8, -2, 8, 5}); - char w_s[256]; - strcpy(w_s, ""); - format_vector_into(w, w_s); + strcpy(s, ""); + format_vector_into(w, s); + printf("w: %s", s); - printf("Norm, Distance\n"); - printf("v: %s", v_s); - printf("w: %s", w_s); printf("l1_norm(v): %f\n", l1_norm(v)); printf("l2_norm(v): %f\n", l2_norm(v)); printf("linf_norm(v): %f\n", linf_norm(v)); @@ -38,52 +37,79 @@ int main() { printf("approx f'(1) w/ fw.d.: %f\n", forward_derivative_at(&f, 1, h)); printf("approx f'(1) w/ bw.d.: %f\n", backward_derivative_at(&f, 1, h)); printf("========\n"); + printf("Least Squares\n"); v = InitArray(double, {1, 2, 3, 4, 5}); - strcpy(v_s, ""); - format_vector_into(v, v_s); + strcpy(s, ""); + format_vector_into(v, s); + printf("v: %s", s); w = InitArray(double, {2, 3, 4, 5, 6}); - strcpy(w_s, ""); - format_vector_into(w, w_s); + strcpy(s, ""); + format_vector_into(w, s); + printf("w: %s", s); + Line *line = least_squares_lin_reg(v, w); - printf("Least Squares\n"); - printf("v: %s", v_s); - printf("w: %s", w_s); + printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a); v = InitArray(double, {1, 2, 3, 4, 5, 6, 7}); - strcpy(v_s, ""); - format_vector_into(v, v_s); + strcpy(s, ""); + format_vector_into(v, s); + printf("v: %s", s); w = InitArray(double, {0.5, 3, 2, 3.5, 5, 6, 7.5}); - strcpy(w_s, ""); - format_vector_into(w, w_s); - printf("v: %s", v_s); - printf("w: %s", w_s); + strcpy(s, ""); + format_vector_into(w, s); + printf("w: %s", s); + line = least_squares_lin_reg(v, w); printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a); printf("========\n"); printf("LU Decomp\n"); - char m_s[2048]; - Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0); - for (int i = 0; i < 8; i++) { - for (int j = 0; j < 8; j++) - m->data[i]->data[j] = (i + 1.0) + j * 3 + (rand() % 12); + Matrix_double *m = InitMatrixWithSize(double, 10, 10, 0.0); + for (int i = 0; i < 10; i++) { + for (int j = 0; j < 10; j++) + m->data[i]->data[j] = (100 - rand() % 200); } - format_matrix_into(m, m_s); - printf("m = %s", m_s); + strcpy(s, ""); + format_matrix_into(m, s); + printf("m = %s", s); - Matrix_double **u_l = put_lu_decomp(m); + Array_double *b = InitArrayWithSize(double, 10, 100.0); + Matrix_double **u_l = lu_decomp(m); Matrix_double *u = u_l[0]; Matrix_double *l = u_l[1]; - strcpy(m_s, ""); - format_matrix_into(u, m_s); - printf("u = %s", m_s); - strcpy(m_s, ""); - format_matrix_into(l, m_s); - printf("l = %s", m_s); + strcpy(s, ""); + format_matrix_into(u, s); + printf("u = %s", s); + strcpy(s, ""); + format_matrix_into(l, s); + printf("l = %s", s); + strcpy(s, ""); + format_vector_into(b, s); + printf("b = %s", s); printf("========\n"); - printf("Back Substitution\n"); + printf("Backward -> Forward Substitution\n"); + + Array_double *b_fsub = fsubst(l, b); + strcpy(s, ""); + format_vector_into(b_fsub, s); + printf("x: %s\n", s); + + Array_double *x_bsub = bsubst(u, b_fsub); + strcpy(s, ""); + format_vector_into(x_bsub, s); + printf("x: %s", s); + + free_vector(b_fsub); + + printf("Verifications (each should be approximately 100)\n"); + for (size_t row = 0; row < m->rows; row++) { + double curr = 0; + for (size_t col = 0; col < m->cols; col++) + curr += m->data[row]->data[col] * x_bsub->data[col]; + printf("Substitution for row %zu = %f\n", row, curr); + } return 0; }